/*
 * SUMMARY:      UnsaturatedFlow.c - Calculate the unsaturated flow
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Bart Nijssen and Mark Wigmosta (*)
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    Apr-1996
 * Last Change: Fri Dec  6 15:39:03 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Calculate the unsaturated flow in the soil column (vertical
 *               flow) 
 * DESCRIP-END.
 * FUNCTIONS:    UnsaturatedFlow()
 * COMMENTS:     (*) Mark Wigmosta, Batelle Pacific Northwest Laboratories, 
 *               ms_wigmosta@pnl.gov
 */

#ifndef lint
static char vcid[] = "$Id: UnsaturatedFlow.cpp,v 1.6 2001/05/15 17:30:27 waichler Exp $";
#endif /* lint */


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "constants.h"
#include "settings.h"
/*  #include "soilmoisture.h" */

#include "infiltration.h"
#include "unsatflow.h"

/*****************************************************************************
  Function name: UnsaturatedFlow()

  Purpose      : Calculate the unsaturated flow in the soil column, and
                 adjust the moisture in each soil layer

  Required     :
	int SpatialExtent  - BASIN or POINT value
    int Dt             - Time step (seconds)
    float DX           - Grid cell width (m)
    float Infiltration - Amount of infiltration entering the top of the soil
                         column (m)
    float SatFlow      - Amount of saturated flow entering the soil column
                         from neighbouring pixels (m) 
    int NSoilLayers    - Number of soil layers
    float TotalDepth   - Total depth of the soil profile (m)
    float Width        - Width of channel or road surface (m)
    float *RootDepth   - Depth of each of the soil layers (m)
    float *Ks          - Vertical saturated hydraulic conductivity in each
                         soil layer (m/s)
    float *PoreDist    - Pore size distribution index for each soil layer
    float *Porosity    - Porosity of each soil layer
    float *FCap        - Field capacity of each soil layer
    float *Perc        - Amount of water percolating from each soil layer to
                         the layer below (m)
    float *PercWidth   - Width of the bottom of each soil layer as a fraction
                         of the gridcell width DX
    float *Adjust      - Correction for each layer for loss of soil storage
                         due to channel/road-cut.  Multiplied with RootDepth
                         to give the layer thickness for use in calculating
                         soil moisture 
    int CutBankZone    - Number of the soil layer containing the bottom of
                         the cut-bank 
    float BankHeight   - Distance from ground surface to channel bed or
                         bottom of road-cut (m) 

  Returns      : void

  Modifies     :
    float *TableDepth - Depth of the water table below the ground surface (m)
    float *RunOff     - Amount of runoff produced at the pixel (m)
    float *Moist      - Moisture content in each soil layer

  Comments     :
    Sources: 
    Bras, R. A., Hydrology, an introduction to hydrologic science, Addisson 
        Wesley, Inc., Reading, etc., 1990.
    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.

    This function is based on Wigmosta et al [1994], and assumes a unit 
    hydraulic gradient in the unsaturated zone.  This implies a 
    steady-state situation and uniform moisture distribution.  This is not
    unreasonable for situations where the groundwater is fairly deep 
    [Bras, 1990], but may be unrealistic for the boreal forest system, or 
    other similar situations where the groundwater table is relatively close
    to the surface. The current formulation does not allow for upward
    movement of water in the unsaturated zone, by capillary and matrix
    forces.  No unsaturated flux is assumed to occur if the water content
    drops below the field capacity.

    The unsaturated hydraulic conductivity is calculated using the
    Brooks-Corey equation see for example Wigmosta et al. [1994], eq.41

    The calculated amount of drainage is averaged with the amount calculated
    for the previous timestep, see eq. 42 [Wigmosta et al., 1994].
  
    The residual moisture content is assumed to be zero (easy to adapt if 
    needed).

    If the amount of soil moisture in a layer after drainage still 
    exceeds the porosity for that layer, this additional amount of water is
    added to the drainage to the next layer.
  
    CHANGES:
  
    Changes have been made to account for the potental loss of soil storage
    in a grid cell due to a road-cut or channel.  Correction coefficents are
    calculated in AdjustStorage() and CutBankGeometry()
*****************************************************************************/
extern "C" {
float WaterTableDepth(int NRootLayers, float TotalDepth, float *RootDepth, 
		      float *Porosity, float *FCap, float *Adjust, 
                      float *Moist);

void UnsaturatedFlow(int SpatialExtent, int Dt, float DX, float Infiltration, float SatFlow,
                     int NSoilLayers, float TotalDepth, float Width, 
                     float *RootDepth, float *Ks,  float *PoreDist, 
                     float *Porosity,  float *FCap, float *Perc, 
                     float *PercWidth, float *Adjust, int CutBankZone, 
                     float BankHeight, float *TableDepth, float *Runoff, 
                     float *Moist, int infiltrationRate, 
                     float infiltCoeff, float infiltExp, float *QuickVol)
{
  float DeepDrainage;           /* amount of drainage from the lowest root 
                                   zone to the layer below it (m) */
  float DeepLayerDepth;         /* depth of the layer below the deepest
                                   root layer */
  float SubSurfaceRunoff = 0.0;	/* amount of water leaving saturated bottom layer(m) in point model mode */
  float Drainage;               /* amount of water drained from each soil 
                                   layer during the current timestep */
  float Exponent;               /* Brooks-Corey exponent */
  float FieldCapacity;          /* amount of water in soil at field capacity
                                   (m) */
  float MaxSoilWater;           /* maximum allowable amount of soil moiture
                                   in each layer (m) */
  float SoilWater;              /* amount of water in each soil layer (m) */
  int i, start_loop = 0;                        /* counter */

  float quick_flow_vol = 0.0, influx_vol, before = 0, old_water;
  
  *QuickVol = 0;
  
  before = Moist[0];
  if( infiltrationRate == VARIABLE ) { /* if infiltration rate variable w/ soil moisture */

    MaxSoilWater = RootDepth[0] * Porosity[0] * Adjust[0];
    SoilWater = RootDepth[0] * Moist[0] * Adjust[0];
    FieldCapacity = RootDepth[0] * FCap[0] * Adjust[0];
    old_water = SoilWater;

    influx_vol = Infiltration;

    if (0 != CutBankZone) { /* doesn't contain bottom of cutbank -- what if above it? */
      influx_vol *= 1 - Width/DX;  
       SoilWater = CalcInfiltration( SoilWater, infiltCoeff, infiltExp, MaxSoilWater, 
                        FieldCapacity, 0.5*Ks[0], (float)Dt, influx_vol, &Drainage, 
                        &quick_flow_vol  );
       Moist[0] = SoilWater / ( RootDepth[0] * Adjust[0] );
    }
    else {  /* contains bottom of cutbank */ 

      if (*TableDepth <= BankHeight) {  /* water table above road cut */
        *Runoff += Infiltration * Width / DX;
        influx_vol *= 1 - Width/DX;
      }
 
       SoilWater = CalcInfiltration( SoilWater, infiltCoeff, infiltExp, MaxSoilWater, 
                        FieldCapacity, 0.5*Ks[0], (float)Dt, influx_vol, &Drainage, 
                        &quick_flow_vol  );

      Moist[0] = SoilWater / (RootDepth[0] * Adjust[0]);
    }

    *Runoff += quick_flow_vol;
    *QuickVol = quick_flow_vol;

    // check mass balance
    float error = SoilWater - old_water - influx_vol + quick_flow_vol + Drainage;


    Perc[0] = Drainage;
    if( NSoilLayers > 1 ) Moist[1] += Perc[0] / (RootDepth[1] * Adjust[1] );

    start_loop++;   // skip upper layer in loop below
    
  }


  for (i = start_loop; i < NSoilLayers; i++) {
    
    /* Update soil moisture in top layer and check for Runoff */
    
    if (i != CutBankZone) {
      if (i == 0) 
        Moist[i] += Infiltration / RootDepth[i];
    }
    else {
      if (*TableDepth <= BankHeight)
        *Runoff += Infiltration * Width / DX;
      else
        if (i != 0)Moist[i] += Infiltration * Width / (DX * RootDepth[i] * 
                                                       Adjust[i]);
      
      /* if required add throughfall/snowmelt to upper zone */
      if (i == 0) {
        if (*Runoff > 0.0) 
          Moist[i] += Infiltration * (DX - Width) / 
            (DX * RootDepth[i] * Adjust[i]);
        else
          Moist[i] += Infiltration / (RootDepth[i] * Adjust[i]);
      }
    }
    
    /* No movement if soil moisture is below field capacity */
    if (Moist[i] > FCap[i]) {
      Exponent = 2.0/PoreDist[i] + 3.0;
      
      if (Moist[i] > Porosity[i])
        /* this can happen because the moisture content can exceed the 
           porosity the way the algorithm is implemented */
        Drainage = Ks[i];
      else
        Drainage = 
          Ks[i] * pow((double) (Moist[i]/Porosity[i]), (double) Exponent);
      
      Drainage *= Dt; 
      Perc[i] = 0.5 * (Perc[i] + Drainage) * PercWidth[i];
      
      MaxSoilWater = RootDepth[i] * Porosity[i] * Adjust[i];
      SoilWater = RootDepth[i] * Moist[i] * Adjust[i];
      FieldCapacity = RootDepth[i] * FCap[i] * Adjust[i];
      
      /* No unsaturated flow if the moisture content drops below field 
         capacity */
      
      if ((SoilWater - Perc[i]) < FieldCapacity) 
        Perc[i] = SoilWater - FieldCapacity;
      
      /* WORK IN PROGRESS */
      /* If the moisture content is greater than the porosity add the 
         additional soil moisture to the percolation */
      
      SoilWater -= Perc[i];
      if (SoilWater > MaxSoilWater)
        Perc[i] += SoilWater - MaxSoilWater;
      
      /* Adjust the moisture content in the current layer, and the layer 
         immediately below it */
      
      Moist[i] -= Perc[i]/(RootDepth[i] * Adjust[i]);
      if (i < (NSoilLayers-1)) 
        Moist[i+1] += Perc[i]/(RootDepth[i+1] * Adjust[i+1]);
    }
    else
      Perc[i] = 0.0;
    
    /* convert back to straight 1-d flux */
    Perc[i] /= PercWidth[i];
  }
  
  DeepLayerDepth = TotalDepth;
  for (i = 0; i < NSoilLayers; i++) 
    DeepLayerDepth -=  RootDepth[i];
  
  if ((NSoilLayers == CutBankZone) || (BankHeight == TotalDepth)) {
    if (*TableDepth <= BankHeight)
      *Runoff += Infiltration*Width / DX;
    else
      *Runoff += 0.0;
    
    /* if required add throughfall/snowmelt through channel/road bed */
    if (*Runoff == 0.0)
      Moist[i] += Infiltration * Width / 
        (DX * DeepLayerDepth * Adjust[NSoilLayers]);
  }
  
  DeepDrainage = (Perc[NSoilLayers-1] * PercWidth[NSoilLayers-1]) + SatFlow;
  Moist[NSoilLayers] += DeepDrainage/(DeepLayerDepth * Adjust[NSoilLayers]);

  /* SRW 1-3-01.  Add logic to deal with point mode.  Drainage below root
  zone becomes SubSurfaceRunoff to the extent that it exceeds porosity of
  the bottom soil layer (which is assumed equal to porosity of lowest root
  zone layer).  SubSurfaceRunoff is added to surface Runoff
  to simplify accounting (in point mode, "routing" of surface runoff
  is just wholesale removal from the grid cell). 
  Modification here and in Aggregate.c, MassEnergyBalance.c, soilmoisture.h    */
  if (SpatialExtent==POINT) {
	  if (Moist[NSoilLayers] > Porosity[NSoilLayers-1]) {
		  SubSurfaceRunoff = (Moist[NSoilLayers] - Porosity[NSoilLayers-1])
		                     * DeepLayerDepth * Adjust[NSoilLayers];
		  Moist[NSoilLayers] = Porosity[NSoilLayers-1];
		  *Runoff += SubSurfaceRunoff;
    }
  } 
  

  /* Calculate the depth of the water table based on the soil moisture 
     profile and adjust the soil moisture profile, to assure that the soil 
     moisture is never more than the maximum allowed soil moisture amount,
     i.e. the porosity.  A negative water table depth means that the water is 
     ponding on the surface.  This amount of water becomes surface Runoff */
  
  *TableDepth = WaterTableDepth(NSoilLayers, TotalDepth, RootDepth, Porosity,
                                FCap, Adjust, Moist);
  
  if (*TableDepth < 0.0) {
    *Runoff += -(*TableDepth);
    *TableDepth = 0.0;
  }
}


}
