/*
 * SUMMARY:      RouteSubSurface.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: Tue Dec 10 09:03:16 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 "trans.h"
#include "soilmoisture.h"
#include "slopeaspect.h"
#include "DHSVMChannel.h"
#include <assert.h>

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

  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].

  SoilMap[y][x].SatFlow, the amount of water in a cell resulting from outflow from the upgradient
  cell, is made available for recomputing the water table in the NEXT timestep,
  in UnsaturatedFlow.c  In other words, water moves just one cell per timestep.
  In this scheme, the order in which you process the cells does not matter.

  WORK IN PROGRESS
*****************************************************************************/
void RouteSubSurface(int Dt, MAPSIZE *Map, TOPOPIX **TopoMap, VEGTABLE *VType, 
                     VEGPIX ** VegMap, ROADSTRUCT **Network, SOILTABLE *SType, 
                     SOILPIX **SoilMap, CHANNEL *ChannelData, AGGREGATED *Total)
{

  int x;			/* counter */
  int y;			/* counter */
  float BankHeight;
  float Width;
  float *Adjust;
  float fract_used;
  float depth;
  float OutFlow;
  float water_out_road = 0.0f, water_out_stream = 0.0f;
  float Transmissivity;
  float AvailableWater, channelLoss = 0.0f;
  int k;
  int i, NegSoilMoisture;
  int NRootLayers = 3;

  /* reset the saturated subsurface flow to zero */ 
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin(TopoMap[y][x].Mask)) {
        SoilMap[y][x].SatFlow = 0;
        SoilMap[y][x].ChannelInt = 0;
        SoilMap[y][x].RoadInt =  0;
        SoilMap[y][x].GwRecharge = 0;
        SoilMap[y][x].SatThickness = SoilMap[y][x].Depth - SoilMap[y][x].TableDepth; 
      }
    }
  }

  /* next sweep through all the grid cells, calculate the amount of
     flow in each direction, and divide the flow over the surrounding
     pixels */
  for (y = 0; y < Map->NY; y++) {
    for (x = 0; x < Map->NX; x++) {
      if (InBasin (TopoMap[y][x].Mask)) {
        BankHeight = (Network[y][x].BankHeight > SoilMap[y][x].Depth) ?
          SoilMap[y][x].Depth : Network[y][x].BankHeight;
        Width = Network[y][x].Width;
        Adjust = Network[y][x].Adjust;
        water_out_road = 0.0;
        
        if (!channel_grid_has_channel(ChannelData->stream_map, x, y)) {

	      /* only bother calculating subsurface flow if water table is above bedrock */          
          if (SoilMap[y][x].TableDepth < SoilMap[y][x].Depth) {
            
            /* get the depth to top of appropriate saturated thickness for
               subsurface flow */
            depth = ((SoilMap[y][x].TableDepth > BankHeight) ? 
                      SoilMap[y][x].TableDepth : BankHeight);
            
            Transmissivity = 
              CalcTransmissivity (SoilMap[y][x].Depth, depth, 
                                  SType[SoilMap[y][x].Soil - 1].KsLat, 
                                  SType[SoilMap[y][x].Soil - 1].KsLatExp);
            
            SoilMap[y][x].GwRecharge = SType[SoilMap[y][x].Soil - 1].groundwaterKs * Dt;

            OutFlow = ( (Transmissivity*TopoMap[y][x].FlowGrad*Dt)/
	                    (TopoMap[y][x].Area) ) + SoilMap[y][x].GwRecharge;
            
            /* check whether enough water is available for redistribution */
            AvailableWater = 
                CalcAvailableWater(VType[VegMap[y][x].Veg - 1].NSoilLayers, 
                                 SoilMap[y][x].Depth, 
                                 VType[VegMap[y][x].Veg - 1].RootDepth, 
                                 SType[SoilMap[y][x].Soil-1].Porosity,
                                 SType[SoilMap[y][x].Soil - 1].FCap, 
                                 SoilMap[y][x].TableDepth, Adjust);

            OutFlow = (OutFlow > AvailableWater) ? AvailableWater : OutFlow;

            if( OutFlow > SoilMap[y][x].GwRecharge ) {
              OutFlow -= SoilMap[y][x].GwRecharge;
            }
            else {
              SoilMap[y][x].GwRecharge = OutFlow;
              OutFlow = 0.0f;
            }
          }
          else { 
            depth = 0.0f; 
            OutFlow = 0.0f;
            SoilMap[y][x].GwRecharge = 0.0f;
          }

          SoilMap[y][x].CumGwRecharge += SoilMap[y][x].GwRecharge;

          /* compute road interception if water table is above road cut */
          if (SoilMap[y][x].TableDepth < BankHeight &&
              channel_grid_has_channel(ChannelData->road_map, x, y)) {
            fract_used = ((float) Network[y][x].fraction/255.0f);
            Transmissivity = 
              CalcTransmissivity(BankHeight, SoilMap[y][x].TableDepth,
                                 SType[SoilMap[y][x].Soil - 1].KsLat,
                                 SType[SoilMap[y][x].Soil - 1].KsLatExp);
            water_out_road = (Transmissivity*fract_used*
			      TopoMap[y][x].FlowGrad*Dt)/(TopoMap[y][x].Area);
            
            AvailableWater =
                CalcAvailableWater(VType[VegMap[y][x].Veg - 1].NSoilLayers,
                                 BankHeight,
                                 VType[VegMap[y][x].Veg - 1].RootDepth,
                                 SType[SoilMap[y][x].Soil-1].Porosity,
                                 SType[SoilMap[y][x].Soil - 1].FCap,
                                 SoilMap[y][x].TableDepth, Adjust);

            water_out_road = (water_out_road > AvailableWater) ? AvailableWater : water_out_road;

            /* save intercepted subsurface flow in Runoff to prevent recycling
               in grid cells that contain a culvert 
            SoilMap[y][x].Runoff +=  water_out_road; */ 

            /* increase lateral inflow to road channel */ 
            SoilMap[y][x].RoadInt =  water_out_road;
            channel_grid_inc_inflow(ChannelData->road_map, x, y, 
                                    water_out_road * TopoMap[y][x].Area);
            
          } /* end if block, grid cell has a road channel */

          /* Subsurface Component - Decrease water change by outwater */
          SoilMap[y][x].SatFlow -= OutFlow + water_out_road + SoilMap[y][x].GwRecharge;

          if(OutFlow > 0.0) {	  
            /* Assign the water to appropriate surrounding pixels */
            fract_used = 0.0f;
            for (k = 0; k < NDIRS; k++) {
              fract_used += (float)TopoMap[y][x].Dir[k];
            }

            OutFlow /= fract_used;

            for (k = 0; k < NDIRS; k++) {
              int nx = Map->xneighbor[k] + x;
              int ny = Map->yneighbor[k] + y;
              if (valid_cell(Map, nx, ny)) {
                SoilMap[ny][nx].SatFlow += OutFlow * TopoMap[y][x].Dir[k];
              }
            }
          }

        } else {   /* cell has a stream channel */
          
          if (SoilMap[y][x].TableDepth < BankHeight) {
            
            float gradient = 4.0 * (BankHeight - SoilMap[y][x].TableDepth);
            if( gradient < 2*TopoMap[y][x].FlowGrad ) gradient = 2*TopoMap[y][x].FlowGrad;
            if (gradient < 0.0) gradient = 0.0;
            
            Transmissivity =
              CalcTransmissivity (BankHeight, SoilMap[y][x].TableDepth, 
                                  SType[SoilMap[y][x].Soil - 1].KsLat, 
                                  SType[SoilMap[y][x].Soil - 1].KsLatExp);

            water_out_stream = (Transmissivity*gradient*Dt) / TopoMap[y][x].Area;
          } else water_out_stream = 0.0;

          SoilMap[y][x].GwRecharge = SType[SoilMap[y][x].Soil - 1].groundwaterKs * Dt;

          /* check whether enough water is available for redistribution */
          if (SoilMap[y][x].TableDepth < SoilMap[y][x].Depth) {

            AvailableWater = CalcAvailableWater(VType[VegMap[y][x].Veg - 1].NSoilLayers, 
                                 SoilMap[y][x].Depth, 
                                 VType[VegMap[y][x].Veg - 1].RootDepth, 
                                 SType[SoilMap[y][x].Soil-1].Porosity,
                                 SType[SoilMap[y][x].Soil - 1].FCap, 
                                 SoilMap[y][x].TableDepth, Adjust);

            /* If AvailableWater is less than potential outflow, satisfy GwRecharge first. */
            if(AvailableWater < SoilMap[y][x].GwRecharge) {
              water_out_stream = 0.0;
              SoilMap[y][x].GwRecharge = AvailableWater;
            } else if(AvailableWater < (SoilMap[y][x].GwRecharge + water_out_stream)) {
                water_out_stream = AvailableWater - SoilMap[y][x].GwRecharge;
            }
          } else { /* not enough water for redistribution */
            water_out_stream = 0.0;
            SoilMap[y][x].GwRecharge = 0.0;
          }

          /* add water infiltrated from channel bed */
          channelLoss = channel_grid_get_recharge(ChannelData->stream_map, x, y) / 
                (float) channel_grid_get_numcells(ChannelData->stream_map, x, y);
          SoilMap[y][x].SatFlow += channelLoss / TopoMap[y][x].Area;

          /* remove water that is going to a channel or percolating to
             deep groundwater zone from the grid cell */
          SoilMap[y][x].SatFlow -= water_out_stream + SoilMap[y][x].GwRecharge;
            
          /* contribute to channel segment lateral inflow */
          channel_grid_inc_inflow(ChannelData->stream_map, x, y, 
                                  water_out_stream * TopoMap[y][x].Area);
            
          SoilMap[y][x].ChannelInt = water_out_stream;
          Total->channelLoss += channelLoss;
          SoilMap[y][x].CumChannelLoss += channelLoss / TopoMap[y][x].Area;
          SoilMap[y][x].CumGwRecharge += SoilMap[y][x].GwRecharge;

        }  /* else block--cell does have a stream channel */
      }  /* end if in basin */
    }  /* end loop thru cols (x) */
  }  /* end loop thru rows (y) */

} /* end RouteSubsurface() */
