/*
* SUMMARY:      Groundwater.c - Route subsurface flow
* USAGE:        Part of DHSVM
*
* AUTHOR:       Mark Wigmosta, Scott Waichler
* ORG:          PNNL
* E-MAIL:       scott.waichler@pnl.gov
* ORIG-DATE:    Apr-96
* Last Change: Tue Dec 10 10:28:42 2002 by Scott Waichler <waichler@tuff.pnl.gov>
* DESCRIPTION:  Route subsurface flow
* DESCRIP-END.
* FUNCTIONS:    InitGroundwater(), RouteGroundwater(), StoreGroundwaterState(),
*               ReadGroundwaterState()
* 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 "trans.h"
#include "soilmoisture.h"
#include "slopeaspect.h"
#include "DHSVMChannel.h"
#include "groundwater.h"
#include "fileio.h"
#include "varid.h"
#include "sizeofnt.h"
#include <malloc.h>
#include<string.h>

#ifndef MIN_GRAD
#define MIN_GRAD .3             /* minimum slope for flow to channel */
#endif
/*****************************************************************************
RouteGroundwater()

  Sources: 
  Wigmosta, M. S., L. W. Vail, and D. P. Lettenmaier, A distributed 
  hydrology-vegetation model for complex terrain, Water Resour. Res.,
  30(6), 1665-1679, 1994.
  
    Quinn, P., K. Beven, P. Chevallier, and O. Planchon, The prediction of 
    hillslope flow paths for distributed hydrological modelling using 
    digital terrain models, Hydrological Processes, 5, 59-79, 1991.
    
      This routine follows Wigmosta et al. [1994] in calculating the subsurface
      flow.  The local gradient is based on the local hydraulic head, consisting 
      of the height of the pixel surface minus the depth of the water table 
      below the water surface.  This has the disadvantage that the local gradients
      have to be recalculated for each pixel for each timestep.  In Wigmosta et
      al. [1994] the local gradient is taken to be equal to the slope of the land
      surface, which is a reasonable assunption for mountainous areas.  For the 
      flat boreal forest landscape it is probably better to use the slope
      of the water surface.
      
        Set the gradient with pixels that are outside tha basin to zero.  This 
        ensures that there is no flux of water across the basin boundary.  In the 
        current implementation water can only leave the basin as surface flow.  
        This may not be entirely realistic, and should be analyzed further.  
        One consequence of this could be that the soil in the basin is more 
        saturated than it would be if subsurface flow out of the basin would
        be allowed.
        
          The surrounding grid cells are numbered in the following way
          
            |-----| DX
            
              0-----1-----2  ---
              |\    |    /|   |
              | \   |   / |   |
              |  \  |  /  |   | DY
              |   \ | /   |   |
              |    \|/    |   |
              7-----*-----3  ---
              |    /|\    |
              |   / | \   |
              |  /  |  \  |
              | /   |   \ |
              |/    |    \|
              6-----5-----4
              
                For the current implementation it is assumed that the resolution is the 
                same in both the X and the Y direction.  If this is not the case an error
                message is generated and the program exits.  The reason is that the 
                formulation for the flow width in the diagonal direction changes if the
                grid is not square.  The method for defining the flow widths in the case
                of square grids is taken from Quinn et al [1991]
                
                  WORK IN PROGRESS
*****************************************************************************/
void RouteGroundwater(int Dt, MAPSIZE *Map, TOPOPIX **TopoMap, SOILTABLE *SType, 
                      SOILPIX **SoilMap, CHANNEL *ChannelData, GWPIX **Groundwater,
                      float *GwToChannel)
{
  
  int x, y, k;      /* counter */
  float GwKsLat, GwOut, GwDeepLoss;
  
  for (y = 0; y < Map->NY; y++)    
    for (x = 0; x < Map->NX; x++)
      if (InBasin (TopoMap[y][x].Mask))
      {
        /* update groundwater store and set hydraulic conductivity */
        
        GwKsLat = SType[SoilMap[y][x].Soil - 1].groundwaterKsLat * Dt; /* Lateral flow */
        /* SoilMap[y][x].GwRecharge is the vertical flux (downward percolation) from
           the upper saturated zone.  It is computed in RouteSubSurface.c  */
        Groundwater[y][x].depth += SoilMap[y][x].GwRecharge;
        
        /* Initialize deep loss to the vertical flux rate (unit gradient).  */
        GwDeepLoss = SType[SoilMap[y][x].Soil - 1].baseLayerKs * Dt;
        
        if ( !Groundwater[y][x].toChannel ) {
          /* GwOut is total outflow, lateral + vertical (deep loss) */ 
          GwOut = ( (GwKsLat * Groundwater[y][x].FlowGrad * Groundwater[y][x].depth) / 
            (TopoMap[y][x].Area) ) + GwDeepLoss;
          GwOut = (GwOut > Groundwater[y][x].depth) ? Groundwater[y][x].depth : GwOut;
          Groundwater[y][x].depth -= GwOut;
          
          if( GwOut > GwDeepLoss ) {
            GwOut -= GwDeepLoss;
          } else { /* all outflow is assigned to downward flux */
            GwDeepLoss = GwOut;
            GwOut = 0.0;
          }
          
          /* Assign the water to appropriate surrounding pixels */
          GwOut /= 255.0f;
          for (k = 0; k < NDIRS; k++)  {
            int nx = Map->xneighbor[k] + x;
            int ny = Map->yneighbor[k] + y;
            if (valid_cell (Map, nx, ny)) {
              /* Put some of the groundwater outflow onto the neighbor cell */
              Groundwater[ny][nx].depth += GwOut * Groundwater[y][x].Dir[k];
            } else if (Groundwater[y][x].Dir[k] > 0) {
              /* The neighbor cell is not valid, but the flow direction points to it.
                 In some applications, the user may wish to use the option 1 below
                 to avoid buildup in the pour point or other border cell that may be
                 intentionally set to flow outside the basin.  For the Hanford Recharge
                 project, I use option 2 to keep the water as deep groundwater.  Note
                 either option assumes the flow direction grid is defined correctly for
                 the user's intentions.  */ 
              /*  Option 1:  */
              /*  SoilMap[y][x].GwRecharge += GwOut * Groundwater[y][x].Dir[k]; */
              /*  Option 2:  */
              Groundwater[y][x].depth += GwOut * Groundwater[y][x].Dir[k];
            }
          }
        }  /* end if block, cell does not have a channel  */
        else   /*  cell has channel  */
        {
          float gradient = Groundwater[y][x].FlowGrad;
          if(gradient == 0 ) gradient = 2.0 * Groundwater[y][x].depth;
          if (gradient <= 0.0) gradient = 0;
          
          GwOut = ( (GwKsLat * gradient * Groundwater[y][x].depth) / 
            (TopoMap[y][x].Area) ) + GwDeepLoss;
          
          GwOut = (GwOut > Groundwater[y][x].depth) ? Groundwater[y][x].depth : GwOut;
          Groundwater[y][x].depth -= GwOut;
          
          if( GwOut > GwDeepLoss ) {
            GwOut -= GwDeepLoss;
          }
          else { /* all outflow is assigned to downward flux */
            GwDeepLoss = GwOut;
            GwOut = 0.0;
          }
          
          /* contribute to channel segment lateral inflow */
          channel_grid_inc_inflow (ChannelData->stream_map, x, y,
            GwOut * TopoMap[y][x].Area);
          
          SoilMap[y][x].ChannelInt += GwOut;
          *GwToChannel += GwOut * TopoMap[y][x].Area;
          
        } /* end else block, cell has channel  */
        
        Groundwater[y][x].cumDeepLoss += GwDeepLoss * TopoMap[y][x].Area;
        Groundwater[y][x].deepLoss = GwDeepLoss * TopoMap[y][x].Area;
      }  /* end block, cell in basin  */
      
}  /* end RouteGroundwater() */


/*****************************************************************************
InitGroundwater()
*****************************************************************************/
void InitGroundwater(LISTPTR Input, OPTIONSTRUCT *Options, MAPSIZE *Map, 
                     TOPOPIX **RealTopoMap, GWPIX ***Groundwater)
{
  FILE *GwToChannelFile;
  const char *Routine = "InitGroundwater";
  char VarName[BUFSIZE+1];  /* Variable name */
  char FlowDirFileName[BUFSIZE+1];  /* File with flow direction */
  char data_line[120];
  int i;      /* Counter */
  int x;      /* Counter */
  int y, n;     /* Counter */
  int NumberType, row, col;   /* Number type of data set */
  unsigned char *Mask = NULL;   /* Basin mask */
  float *Elev;                  /* Surface elevation */
  float initialDepth;
  TOPOPIX **TopoMap = NULL;
  
  /* Process the [GROUNDWATER] section in the input file */
  STRINIENTRY StrEnv[] = {
    {"GROUNDWATER", "GWATER DEM FILE"        , ""  , ""},
    {"GROUNDWATER", "GWATER FLOW DIRECTION FILE", ""  , ""},
    {"GROUNDWATER", "GWATER TO CHANNEL FILE", ""  , ""},
    {"GROUNDWATER", "INITIAL GROUNDWATER DEPTH", ""  , ""},
    {NULL       , NULL            , ""  , NULL}
  };
  
  /* don't look for flow direction file name if we are not reading 
  flow directions 
  if( Options->ReadFlowDir != TRUE )StrEnv[1].SectionName = 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);
  }
  
  /* Assign the attributes to the correct map pixel */
  if (!(*Groundwater = (GWPIX **) calloc(Map->NY, sizeof(GWPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map->NY; y++) {
    if (!((*Groundwater)[y] = (GWPIX *) calloc(Map->NX, sizeof(GWPIX))))
      ReportError((char *) Routine, 1);
  }
  
  /* get initial groundwater depth */
  if (!CopyFloat(&initialDepth, StrEnv[initial_groundwater_depth].VarStr, 1))
    ReportError(StrEnv[initial_groundwater_depth].KeyName, 51);
  
  /* In elevation band mode set groundwater flow directions equal to surface */
  if( Options->Extent == ELEVBAND ){  
    for (y = 0; y < Map->NY; y++)
      for (x = 0; x < Map->NX; x++) {
        (*Groundwater)[y][x].Dem  = RealTopoMap[y][x].Dem;
        (*Groundwater)[y][x].depth  = initialDepth;;
        (*Groundwater)[y][x].TotalDir = RealTopoMap[y][x].TotalDir;
        (*Groundwater)[y][x].FlowGrad = RealTopoMap[y][x].FlowGrad;
        (*Groundwater)[y][x].toChannel = 1;
        (*Groundwater)[y][x].deepLoss = 0.0;
        (*Groundwater)[y][x].cumDeepLoss = 0.0;
        for (n = 0; n < NDIRS; n++){
          (*Groundwater)[y][x].Dir[n] = RealTopoMap[y][x].Dir[n];
        }
      }
  }
  else {  /* In basin mode allow groundwater flow direction to differ from surface */
    
    /* assign name to flow direction file if required */
    if( Options->ReadFlowDir == TRUE ) {
      if (IsEmptyStr(StrEnv[gwater_flowdir_file].VarStr))
        ReportError(StrEnv[gwater_flowdir_file].KeyName, 51);
      else
        strcpy(FlowDirFileName, StrEnv[gwater_flowdir_file].VarStr);
    }
    else
      strcpy(FlowDirFileName, "none");
    
    /* Read the elevation data from the DEM dataset */
    GetVarName(001, 0, VarName);
    GetVarNumberType(001, &NumberType);
    if (!(Elev = (float *) calloc(Map->NX * Map->NY, SizeOfNumberType(NumberType)))) 
      ReportError((char *) Routine, 1);
    //  Read2DMatrix(StrEnv[0].VarStr, Elev, NumberType, Map->NY, Map->NX, 0,
    //         VarName);
    Read2DMatrix(StrEnv[gwater_demfile].VarStr, Elev, NumberType, Map->NY, Map->NX, 0,
      VarName);
    
    for (y = 0, i = 0; y < Map->NY; y++) 
      for (x = 0; x < Map->NX; x++, i++){ 
        (*Groundwater)[y][x].Dem  = Elev[i];
        (*Groundwater)[y][x].depth  = initialDepth;
        (*Groundwater)[y][x].toChannel = 0.0;
        (*Groundwater)[y][x].deepLoss = 0.0;
        (*Groundwater)[y][x].cumDeepLoss = 0.0;
      }
      
      free(Elev);
      
      
      /* build fake topomap with groundwater info  */
      if (!(TopoMap = (TOPOPIX **) calloc(Map->NY, sizeof(TOPOPIX *))))
        ReportError((char *) Routine, 1);
      for (y = 0; y < Map->NY; y++) {
        if (!((TopoMap)[y] = (TOPOPIX *) calloc(Map->NX, sizeof(TOPOPIX))))
          ReportError((char *) Routine, 1);
      }
      
      for (y = 0; y < Map->NY; y++)
        for (x = 0; x < Map->NX; x++) {
          TopoMap[y][x].Area = RealTopoMap[y][x].Area;
          TopoMap[y][x].Dem = (*Groundwater)[y][x].Dem;
          TopoMap[y][x].TotalDir = 0;
          TopoMap[y][x].Travel = 0;
          TopoMap[y][x].Grad = 0;
          TopoMap[y][x].Slope = 0;
          TopoMap[y][x].Aspect = 0;
          TopoMap[y][x].FlowGrad = 0;
          TopoMap[y][x].Mask = RealTopoMap[y][x].Mask;
          for (n = 0; n < NDIRS; n++)TopoMap[y][x].Dir[n] = 0;
        }
        
        /* Calculate slope, aspect, magnitude of subsurface flow gradient, and 
        fraction of flow flowing in each direction based on the land surface 
        slope.  The slope is also used in the inline radiation calculations. */
        if (Options->FlowGradient == TOPOGRAPHY ) 
          ElevationSlopeAspect(Map, TopoMap, Options->ReadFlowDir, FlowDirFileName);
        
        /* copy back to Groundwater */
        for (y = 0; y < Map->NY; y++){
          for (x = 0; x < Map->NX; x++) {
            (*Groundwater)[y][x].TotalDir = TopoMap[y][x].TotalDir;
            (*Groundwater)[y][x].FlowGrad = TopoMap[y][x].FlowGrad;
            for (n = 0; n < NDIRS; n++){
              (*Groundwater)[y][x].Dir[n] = TopoMap[y][x].Dir[n];
            }
          }
        }
        
        /*  Cleanup TopoMap  */
        for (y = 0; y < Map->NY; y++) {
          if( TopoMap[y] != (TOPOPIX *)NULL )
            free(  TopoMap[y] );
        }
        if (TopoMap != (TOPOPIX **)NULL)
          free(TopoMap);
        
        
        /* read groundwater flow to surface channel map file name */
        if (IsEmptyStr(StrEnv[gwater_to_channel_file].VarStr))
          ReportError(StrEnv[gwater_to_channel_file].KeyName, 51);
        else
          strcpy(FlowDirFileName, StrEnv[gwater_to_channel_file].VarStr);
        
        /* open groundwater flow to surface channel map file */
        OpenFile(&GwToChannelFile, StrEnv[gwater_to_channel_file].VarStr, "r", FALSE);
        
        /* loop through header info */
        for( i = 0; i < 10; i++ ){
          if( fgets( data_line, 120, GwToChannelFile ) == NULL ){
            printf("Error Groundwater to channel map File: exit program\n");
            exit(1);
          }
        }
        
        while( fgets( data_line, 120, GwToChannelFile ) != NULL  ){
          if( sscanf(data_line, "%d %d",  &col, &row ) != -1 ) {
            (*Groundwater)[row][col].toChannel = 1;
          }
          else {
            break;  /* done reading file */
          }
        }
        
        
  }
  
}

/*****************************************************************************
StoreGroundwaterState()

  Store the current groundwater state.
  
    *****************************************************************************/
    void StoreGroundwaterState(char *Path, DATE *Now, MAPSIZE *Map,  SOILPIX **SoilMap, 
      TOPOPIX **TopoMap, GWPIX **Groundwater)
    {
      char OutFileName[BUFSIZ+1] = "";
      char Str[BUFSIZ+1] = "";
      FILE *OutFile = NULL;
      float *Array;
      const char *Routine = "StoreGroundwaterState";
      int y, x;
      
      printf("storing groundwater state \n");
      /* Create storage file */
      sprintf(Str, "%02d.%02d.%04d.%02d.%02d.%02d", Now->Month, Now->Day,
        Now->Year, Now->Hour, Now->Min, Now->Sec);
      sprintf(OutFileName, "%sGroundwater.State.%s", Path, Str);
      OpenFile(&OutFile, OutFileName, "wb", TRUE);
      
      /* Store data */
      
      if (!(Array = (float *) calloc(Map->NY * Map->NX, sizeof(float))))
        ReportError((char *) Routine, 1);
      
      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask) )
            ((float *)Array)[y*Map->NX + x] = SoilMap[y][x].GwRecharge;
          else
            ((float *)Array)[y*Map->NX + x] = NA;
        }
      }
      fwrite (Array, sizeof (float), Map->NY * Map->NX, OutFile);  /* recharge  */
      
      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask) )
            ((float *)Array)[y*Map->NX + x] = Groundwater[y][x].depth;
          else
            ((float *)Array)[y*Map->NX + x] = NA;
        }
      }
      fwrite (Array, sizeof (float), Map->NY * Map->NX, OutFile);  /* recharge  */
      
      free(Array);
      
      /* Close file */
      fclose(OutFile);
    }
    
    /*****************************************************************************
    ReadGroundwaterState()
    
      Read the current groundwater state.
      
    *****************************************************************************/
    void ReadGroundwaterState(char *Path, DATE *Now, MAPSIZE *Map,  SOILPIX **SoilMap, 
      TOPOPIX **TopoMap, GWPIX **Groundwater)
    {
      char InFileName[BUFSIZ+1] = "";
      char Str[BUFSIZ+1] = "";
      FILE *InFile = NULL;
      float *Array;
      const char *Routine = "ReadGroundwaterState";
      int y, x;
      
      /* Re-create the storage file name and open it */
      sprintf(Str, "%02d.%02d.%04d.%02d.%02d.%02d", Now->Month, Now->Day, 
        Now->Year, Now->Hour, Now->Min, Now->Sec);
      sprintf(InFileName, "%sGroundwater.State.%s", Path, Str);
      OpenFile(&InFile, InFileName, "rb", TRUE);
      
      if (!(Array = (float *) calloc(Map->NY * Map->NX, sizeof(float))))
        ReportError((char *) Routine, 1);
      
      /* read groundwater recharge */
      fread( Array, sizeof(float), Map->NY * Map->NX, InFile );
      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask))
            SoilMap[y][x].GwRecharge = ((float *)Array)[y*Map->NX + x];
          else
            SoilMap[y][x].GwRecharge = NA;
        }
      }
      
      /* read groundwater depth */
      fread( Array, sizeof(float), Map->NY * Map->NX, InFile );
      for (y = 0; y < Map->NY; y++) {
        for (x = 0; x < Map->NX; x++) {
          if (InBasin(TopoMap[y][x].Mask))
            Groundwater[y][x].depth = ((float *)Array)[y*Map->NX + x];
          else
            Groundwater[y][x].depth = NA;
        }
      }
      
      free(Array);
      
      /* Close file */
      fclose(InFile);
      
      
 
}
    
    
    
    
    
    
    
    
