/*
 * SUMMARY:      SoilResistanceEvap.c - Calculate evaporation from the soil 
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Mark Wigmosta
 * ORG:          Pacific Northwest National Laboratory
 * DESCRIPTION:  Calculate evaporation from the soil 
 * DESCRIP-END.
 * FUNCTIONS:    SoilResistanceEvap()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: SoilResistanceEvap.c,v 1.1 2001/05/15 17:47:40 waichler Exp $";
#endif /* lint */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include "settings.h"
#include "DHSVMerror.h"
#include "massenergy.h"
#include "constants.h"
#include "extreme.h"

/*****************************************************************************
  SoilResistanceEvap()
*****************************************************************************/
float SoilResistanceEvap(int Dt, float Temp, float Slope, float Gamma, float Lv, 
		      float AirDens, float Vpd, float NetRad, float RaSoil, 
		      float MaxEvap, float FieldCap, float Rmin, float Cs, float RootDepth, 
          float *MoistContent, float Adjust, int Overstory)
{
  float EPot;			/* Potential evaporation from soil during timestep (m/s) */
  float SoilEvap;		/* Amount of evaporation directly from the soil (m) */
  float SoilMoisture;		/* Amount of water in surface soil layer (mm) */
  float Rsoil;          /* Soil resistance to evaporation (s/m)
  
  /* PET at the soil surface w/ corrections for overstory if required */
  EPot = (Slope * NetRad + AirDens * CP * Vpd/RaSoil)/
    (WATER_DENSITY * Lv * (Slope + Gamma));

  if(EPot < 0) EPot = 0.0;
  
  /* The potential evaporation rate accounts for the amount of moisture that
     the atmosphere can absorb.  If we do not account for the amount of
     evaporation from overlying evaporation, we can end up with the situation
     that all vegetation layers and the soil layer transpire/evaporate at the
     potential rate, resulting in an overprediction of the actual evaporation
     rate.  Thus we subtract the amount of evaporation that has already
     been calculated for overlying layers to specify MaxEvap (m/s).
 */

  EPot = ( MaxEvap < EPot ) ? MaxEvap : EPot;
  
  /* Calculate soil resistance */
  Rsoil = Rmin * exp( -Cs * *MoistContent/FieldCap );

  SoilEvap = EPot * (Slope + Gamma) / (  Slope + Gamma*(1 + Rsoil/RaSoil) );

  if( SoilEvap > EPot )SoilEvap = EPot;

  /* adjust for width and convert from flux to depth */
  SoilEvap *= Adjust * Dt;
  SoilMoisture = *MoistContent * RootDepth * Adjust;

  if (SoilEvap > SoilMoisture) {
    SoilEvap = SoilMoisture;
    *MoistContent = 0.0;
  }
  else {
    SoilMoisture -= SoilEvap;
    *MoistContent = SoilMoisture/( RootDepth * Adjust );
  }

  return SoilEvap;
}
		    
