/*
 * SUMMARY:      WaterTableDepth.c - Calculate the depth of the water table
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Bart Nijssen
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    Apr-1996
 * LAST-MOD: Tue Dec 17 08:36:19 2002 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Calculate the depth of the water table below the ground
 *               surface corrected for road and channel effects
 * DESCRIP-END.
 * FUNCTIONS:    WaterTableDepth()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: WaterTableDepth.c,v 1.3 2001/04/04 16:43:51 waichler Exp $";
#endif /* lint */

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "settings.h"
#include "soilmoisture.h"

/*****************************************************************************
  Function name: WaterTableDepth()

  Purpose      :   This function calculates the depth of the water table
                   below the ground surface, based on the amount of soil
                   moisture in the root zone layers. 

  Required     :
    int NRootLayers  - Number of soil layers
    float TotalDepth - Total depth of the soil profile (m)
    float *RootDepth - Depth of each of the soil layers (m)
    float *Porosity  - Porosity of each soil layer
    float *FCap      - Field capacity of each soil layer
    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 

  Returns      :
    float TableDepth - Depth of the water table below the soil surface

  Modifies     :
    float *Moist - Amount of soil moisture in each layer

  Comments     :   
    This function completely overhauled by SRW 12/05/02.  In the original version,
    a thin lower layer could lead to negative soil moisture in that layer.  There
    is really no good way to fix that under the old logic, because of the bottom-up
    approach for redistributing soil moisture.  And in the end, the TableDepth in the old logic
    was set in a very simple way using the moisture of the entire soil column anyway.
    So, in this new version of WaterTableDepth(), we start with the TotalStorage, and
    satisfy field capacity of all soil layers first, if possible.  Then the remaining
    water content is added to create saturated thickness, starting with the bottom layer, 
    until it is used up.  This approach is a great simplification over the
    too-complicated old method.  It assumes that vertical percolation is not delayed, and
    it does change the moisture regime in the upper layers and thereby affect ET.
    We start with the assumption that each soil layer has at least a field capacity
    moisture content.  If there is not enough soil moisture to allow that, then the
    water table height is set to zero, and no redistribution is done in that timestep.
    If total soil moisture exceeds field capacity of the total profile, then we proceed
    to compute water table height assuming that all soil above the water table is at
    field capacity.

    Notes for original version:
    Since no unsaturated flow is allowed to occur when the moisture content
    is below the field capacity, and because lateral saturated flow is the
    only mechanism by which water can disappear from the soil below the
    deepest root layer, the soil moisture content in the soil below the
    deepest root layer can never fall below field capacity.  The water
    immediately above the water table is assumed to be at field capacity. 
    Changes have been made to account for the potential loss of soil storage
    in a grid cell due to a road-cut or channel.
*****************************************************************************/
float WaterTableDepth(int NRootLayers, float TotalDepth, float *RootDepth, 
		      float *Porosity, float *FCap, float *Adjust, float *Moist)  {
  int i;			/* counter */

  /* these arrays INCLUDE the layer below the deepest root layer (m) */
  float SoilWater[5] = {0., 0., 0., 0., 0.}; /* water content for soil layer (m) */
  float SoilThickness[5] = {0., 0., 0., 0., 0.};  /* thickness of soil layer (m) */
  float SoilFCap[5] = {0., 0., 0., 0., 0.};  /* field capacity soil layer (m/m) */
  float SoilPorosity[5] = {0., 0., 0., 0., 0.};  /* porosity of soil layer (m/m) */
  float Capacity[5] = {0., 0., 0., 0., 0.}; /* water content at saturation for soil layer (m) */
  float DiffCapacity[5] = {0., 0., 0., 0., 0.}; /* water content between saturation and 
                                                   field capacity for soil layer (m) */

  float DeepLayerThickness; /* thickness of soil layer below deepest root layer (m) */ 
  float TableDepth;		/* depth of the water table below ground surface (m) */
  float TableHeight = 0.0;		/* height of the water table above soil bottom (m) */

  float TotalStorage = 0.0; /* actual water content in entire soil profile (m) */  
  float TotalRootStorage = 0.0; /* actual water content in root zone layers (m) */  
  float TotalCapacity = 0.0; /* water content at saturation for entire soil profile (m) */
  float TotalCapacityAtFC = 0.0; /* water content at field capacity for entire soil profile (m) */
  float RemainingStorage = 0.0; /* convenience water content variable for redistribution (m) */
  float test = 0.0f; /* convenience variable for debugging  */
  float MakeupStorage = 0.0; /* water content needed to raise negative soil moisture to FC (m) */  
  float SoilWaterFC;  /*  water content at field capacity (m)  */

  assert(NRootLayers <= 4);  /* Size for local arrays is hardwired above */

  /* Initialize the local counterparts to the global variables that were previously
     set only for the rooting zone layers.  */
  DeepLayerThickness = TotalDepth; /* init only */
  for (i = 0; i < NRootLayers; i++) {
    SoilThickness[i] = RootDepth[i];
    SoilWater[i] = Moist[i] * SoilThickness[i] * Adjust[i];
    TotalRootStorage += SoilWater[i];
    DeepLayerThickness -= RootDepth[i];
    SoilFCap[i] = FCap[i];
    SoilPorosity[i] = Porosity[i];
  }
  SoilThickness[NRootLayers] = DeepLayerThickness;
  SoilWater[NRootLayers] = Moist[NRootLayers] * SoilThickness[NRootLayers] * Adjust[NRootLayers];
  SoilFCap[NRootLayers] = FCap[NRootLayers-1];
  SoilPorosity[NRootLayers] = Porosity[NRootLayers-1];
  
  /* Compute total actual water storage and storage capacities (m) in soil profile.
     Remember, Moist has units of volumetric water content (m/m).  */  
  for (i = 0; i <= NRootLayers; i++) {
    TotalStorage += SoilWater[i];
    Capacity[i] = SoilPorosity[i] * SoilThickness[i] * Adjust[i];
    DiffCapacity[i] = (SoilPorosity[i] - SoilFCap[i]) * SoilThickness[i] * Adjust[i];
    TotalCapacity += Capacity[i];
    TotalCapacityAtFC += SoilFCap[i] * SoilThickness[i] * Adjust[i];
  }

  if(TotalStorage <= TotalCapacityAtFC) { 

    TableDepth = TotalDepth;

    /* Redistribute water if bottom layer has negative Moist.  In that case, shift water
       to the bottom layer to raise its moisture content to field capacity, and reduce
       soil moisture in the other soil layers such that all have the same volumetric 
       water content (m/m).   */
    if (Moist[NRootLayers] < 0.0) {
      MakeupStorage = (SoilFCap[NRootLayers] - Moist[NRootLayers]) *
                      SoilThickness[NRootLayers] * Adjust[NRootLayers];
      Moist[NRootLayers] = SoilFCap[NRootLayers];

      for (i = 0; i < NRootLayers; i++) {
        SoilWater[i] -= MakeupStorage * (SoilWater[i] / TotalRootStorage);
        assert(SoilWater[i] >= 0.0);  /* get rid of this eventually */
        Moist[i] = SoilWater[i] / (SoilThickness[i] * Adjust[i]);
      }
    }

    return TableDepth;

  } else { /* enough water to redistribute and have a positive 
              water table height (saturated thickness) */

    RemainingStorage = TotalStorage;

    /* Initialize moisture content of all the soil layers to field capacity; at least one
       will have more than that after redistribution is complete */
    for (i = 0; i <= NRootLayers; i++) {
      SoilWaterFC = SoilFCap[i] * SoilThickness[i] * Adjust[i];
      SoilWater[i] = SoilWaterFC;
      RemainingStorage -= SoilWaterFC;
    }

    /* Distribute the water that is in excess of total soil profile field capacity */
    for (i = NRootLayers; i >= 0; i--) {
      if(RemainingStorage >= DiffCapacity[i]) {  /* layer is fully saturated */
        SoilWater[i] += DiffCapacity[i];
        /* Moist[i] = SoilPorosity[i]; */
        RemainingStorage -= DiffCapacity[i];
        TableHeight += SoilThickness[i];
      } else if (RemainingStorage > 0.0) {  /* layer is only partially saturated */
        SoilWater[i] += RemainingStorage;
        TableHeight += RemainingStorage / (SoilPorosity[i] - SoilFCap[i]);
        RemainingStorage = 0.0;
      } /* else RemainingStorage is 0, no further changes to TableHeight. */
      Moist[i] = SoilWater[i] / (SoilThickness[i] * Adjust[i]);
    }

    assert(RemainingStorage >= 0.0);  /* get rid of this eventually */

    /* if RemainingStorage > 0, water is ponded on surface */
    TableHeight += RemainingStorage;

    TableDepth = TotalDepth - TableHeight;

    return TableDepth;

  } /*  end else block--there is a positive water table height (saturated thickneses)  */

}  /*  end of WaterTableDepth()  */
