/*
* SUMMARY:      Groundwater.c - Route subsurface flow
* USAGE:        Part of DHSVM
*
* AUTHOR:       Bill Perkins
* ORG:          PNNL
* E-MAIL:       nijssen@u.washington.edu
* ORIG-DATE:    Apr-96
* Last Change: Mon Dec  9 17:10:32 2002 by Scott Waichler <waichler@tuff.pnl.gov>
* DESCRIPTION:  Route subsurface flow
* DESCRIP-END.
* FUNCTIONS:    RouteSubSurface()
* COMMENTS:     
*/

#ifndef lint
static char vcid[] = "$Id: RouteSubSurface.c,v 1.2 2001/03/31 00:10:34 waichler Exp $";
#endif /* lint */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "constants.h"
#include "frozensoil.h"
#include "fileio.h"
#include "varid.h"
#include "sizeofnt.h"
#include "Calendar.h"
#include <malloc.h>
#include<string.h>

static float K_L, B, N, thresValue;
static int startJDay;



/*****************************************************************************
InitFrozenSoil()
*****************************************************************************/
void InitFrozenSoil(LISTPTR Input, OPTIONSTRUCT *Options, MAPSIZE *Map, 
                     FRZPIX ***FrozenSoil)
{
  enum KEYS {thres_value, frozen_soil_k_l, frozen_soil_b, frozen_soil_n, frozen_soil_jday};
  const char *Routine = "InitFrozenSoil";
  int i;      /* Counter */
  int x;      /* Counter */
  int y;     /* Counter */
  
  /* Process the [FROZENSOIL] section in the input file */
  STRINIENTRY StrEnv[] = {
    {"FROZENSOIL", "THRESHOLD VALUE", ""  , ""},
    {"FROZENSOIL", "SOIL K_L", ""  , ""},
    {"FROZENSOIL", "SOIL B", ""  , ""},
    {"FROZENSOIL", "SOIL N", ""  , ""},
    {"FROZENSOIL", "START JDAY", ""  , ""},
    {NULL       , NULL            , ""  , NULL}
  };
  
  /* Read the key-entry pairs from the input file */
  for (i = 0; StrEnv[i].SectionName; i++) {
    GetInitString(StrEnv[i].SectionName, StrEnv[i].KeyName, StrEnv[i].Default,
      StrEnv[i].VarStr, (unsigned long) BUFSIZE, Input);
    if (IsEmptyStr(StrEnv[i].VarStr))
      ReportError(StrEnv[i].KeyName, 51);
  }
  

  /* get cumulative flux that indicates frozen soil */ 
  if (!CopyFloat(&thresValue, StrEnv[thres_value].VarStr, 1))
    ReportError(StrEnv[thres_value].KeyName, 51);

  /* get ratio of soil thermal conductivity to representative depth */ 
  if (!CopyFloat(&K_L, StrEnv[frozen_soil_k_l].VarStr, 1))
    ReportError(StrEnv[frozen_soil_k_l].KeyName, 51);

  /* get proportionality constant B */ 
  if (!CopyFloat(&B, StrEnv[frozen_soil_b].VarStr, 1))
    ReportError(StrEnv[frozen_soil_b].KeyName, 51);

  /* get constant N */ 
  if (!CopyFloat(&N, StrEnv[frozen_soil_n].VarStr, 1))
    ReportError(StrEnv[frozen_soil_n].KeyName, 51);

  /* get constant julian day to reset cumulative flux to zero */ 
  if (!CopyInt(&startJDay, StrEnv[frozen_soil_jday].VarStr, 1))
    ReportError(StrEnv[frozen_soil_jday].KeyName, 51);

  /* build array of structures */
  if (!(*FrozenSoil = (FRZPIX **) calloc(Map->NY, sizeof(FRZPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map->NY; y++) {
    if (!((*FrozenSoil)[y] = (FRZPIX *) calloc(Map->NX, sizeof(FRZPIX))))
      ReportError((char *) Routine, 1);
  }
  
 
  for (y = 0; y < Map->NY; y++)
    for (x = 0; x < Map->NX; x++) {
      (*FrozenSoil)[y][x].airTemp = 0;
      (*FrozenSoil)[y][x].oldAirTemp = 0;
      (*FrozenSoil)[y][x].heatFlux  = 0;
    }
}

/*****************************************************************************
UpdateSoilPhase()
*****************************************************************************/
void UpdateSoilPhase(int NDaySteps, int JDay, MAPSIZE Map, SNOWPIX **SnowMap, 
                     FRZPIX **FrozenSoil){

  int x, y;
  float grd_flux, up, snow_depth;
  
  up = 2.5 * sin( JDay + 80 );


  for (y = 0; y < Map.NY; y++)
    for (x = 0; x < Map.NX; x++) {
      snow_depth = SnowMap[y][x].Swq;
      FrozenSoil[y][x].airTemp /= (float)NDaySteps;

      if( snow_depth > 0.1 ) FrozenSoil[y][x].airTemp = 0;
      if( FrozenSoil[y][x].heatFlux < 0 && FrozenSoil[y][x].oldAirTemp > 0 )
        FrozenSoil[y][x].oldAirTemp = 0;

      grd_flux = K_L*(FrozenSoil[y][x].airTemp - FrozenSoil[y][x].oldAirTemp/B)
                   * ( 1.0 - snow_depth/(snow_depth + N) );
     
      if( JDay == startJDay ) FrozenSoil[y][x].heatFlux = 0.0;
      
      FrozenSoil[y][x].heatFlux += grd_flux + up;

      FrozenSoil[y][x].oldAirTemp = FrozenSoil[y][x].airTemp;
      FrozenSoil[y][x].airTemp = 0;
    }

}

/*****************************************************************************
FrozenSoilReduction()
*****************************************************************************/
float FrozenSoilReduction(FRZPIX **FrozenSoil, int y, int x){

  if( FrozenSoil[y][x].heatFlux < thresValue )
    return 0.0;
  else
    return 1.0;

}

/*****************************************************************************
FrozenSoilReduction()
*****************************************************************************/
float FrozenSoilReductionWithMoisture(FRZPIX **FrozenSoil, int y, int x, 
                                      float UpperSoilMoisture, float UpperFCap, 
                                      float UpperPorosity) {
/*   This approach uses a combination of soil heat flux and soil moisture of the  */
/*   upper layer to determine the infiltration reduction factor.  Both of these */
/*   individual factors involve linear change between end points corresponding */
/*   to 0 and 1.0.  If the heat flux is less than some threshold value, then the */
/*   infiltration factor will be set as the product of the individual factors.  */
/*   This is obviously very ad hoc and crude, but it is a reasonable first cut */
/*   at addressing the interaction of the two soil properties in limiting  */
/*   permeability of the soil surface.  */
  float lower_heatflux_endpoint = -33.6; 
  float upper_heatflux_endpoint = -20.0;
  float lower_moisture_endpoint = 0.5 * (UpperFCap + UpperPorosity); 
  float upper_moisture_endpoint = UpperPorosity;
  float heatflux_threshold;
  float fmoist, fhf;

  heatflux_threshold = upper_heatflux_endpoint;

  if( FrozenSoil[y][x].heatFlux < heatflux_threshold) {
    if (UpperSoilMoisture <= lower_moisture_endpoint) fmoist = 1.0;
    else if (UpperSoilMoisture >= upper_moisture_endpoint) fmoist = 0.0;
    else fmoist = (UpperSoilMoisture - lower_moisture_endpoint) / 
                  (upper_moisture_endpoint - lower_moisture_endpoint);

    if (FrozenSoil[y][x].heatFlux <= lower_heatflux_endpoint) fhf = 0.0;
    else if (FrozenSoil[y][x].heatFlux >= upper_heatflux_endpoint) fhf = 1.0;
    else fhf = (FrozenSoil[y][x].heatFlux - lower_heatflux_endpoint) / 
                  (upper_heatflux_endpoint - lower_heatflux_endpoint);
    return fmoist*fhf;
  } else return 1.0;

}

/*****************************************************************************
FrozenSoilReductionFromSnow()
*****************************************************************************/
float FrozenSoilReductionFromSnow(float swq) {

  /* "all or nothing"  0 or 1 */
/*   if(swq > 0.010) */
/*     return 0.0; */
/*   else */
/*     return 1.0; */

  /* threshold more than which no water is infiltrated */
  float swq_threshold = 0.017;  
  if(swq == 0.0) {
    return 1.0;
  } else if (swq >= swq_threshold) {
    return 0.0;
  } else {
    return 1 - (swq/swq_threshold);
  }
}

