/*
 * SUMMARY:      EvapoTranspiration.c - Calculate evapotranspiration
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Bart Nijssen
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    Apr-96
 * Last Change: Tue May 22 13:38:36 2001 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Calculate evapotranspiration
 * DESCRIP-END.
 * FUNCTIONS:    EvapoTranspiration()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: EvapoTranspiration.c,v 1.6 2001/05/22 20:39:43 waichler Exp $";
#endif /* lint */

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

/*****************************************************************************
  EvapoTranspiration()
*****************************************************************************/
void EvapoTranspiration(int Layer, int Dt, PIXMET *Met, float NetRad, 
			float Rp, VEGTABLE *VType, SOILTABLE *SType, 
			float MoistureFlux, SOILPIX *LocalSoil, float *Int, 
      EVAPPIX *LocalEvap, float *Adjust, float Ra, int HeatFluxOption)
{
  float *Rc;			/* canopy resistance associated with 
				   conditions in each soil layer (s/m) */
  float DryArea=1.0;		/* relative dry leaf area  */
  float DryEvapTime;		/* amount of time remaining during a timestep
				   after the interception storage is depleted 
				   (sec) */
  float F;                      /* Fractional coverage by vegetation layer */
  float SoilMoisture;		/* Amount of water in each soil layer (m) */
  float WetArea = 0.0;		/* relative leaf area wetted by interception 
				   storage */
  float WetEvapRate;		/* evaporation rate from wetted fraction per 
				   unit ground area (m/s) */
  float WetEvapTime=0.0;		/* amount of time needed to evaporate the 
				   amount of water in interception storage 
				   (sec) */
  float Temperature = 0;  /* if HeatFluxOption == TRUE use soil temp in CanopyResistance,
                             otherwise use Tair */
  int i;			/* counter */

  /* Convert the water amounts related to partial canopy cover to a pixel depth 
     as if the entire pixel is covered.  These depths will be converted 
     back later on. */ 
  F = VType->Fract[Layer];
  *Int /= F;
  NetRad /= F;
  MoistureFlux /= F;

  /* allocate memory for the canopy resistance array */
  
  if (!(Rc = (float *) calloc(VType->NSoilLayers, sizeof(float))))
    ReportError("EvapoTranspiration()", 1);
  
  /* Calculate the evaporation rate in m/s */

  LocalEvap->EPot[Layer] = (Met->Slope * NetRad + 
                            Met->AirDens * CP * Met->Vpd/Ra)/
    (WATER_DENSITY * Met->Lv * (Met->Slope + Met->Gamma));

  /* 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 from the potential evaporation.
     Another mechanism that could be used to account for this would be to 
     decrease the vapor pressure deficit while going down through the canopy
     (not implemented here) */

  LocalEvap->EPot[Layer] -= MoistureFlux/Dt;

  if (LocalEvap->EPot[Layer] < 0)
    LocalEvap->EPot[Layer] = 0;

  /* SRW 1-17-01.  More changes below.  */  
  if (*Int > 0.0) WetArea = pow((double) (*Int/VType->MaxInt[Layer]), (double) (2.0/3.0));
  if (WetArea >= 1.0) {  /* avoid potential roundoff problems */
    WetArea = 1.0; 
    DryArea = 0.0;
  }
  else DryArea = 1.0 - WetArea;
  
  /* calculate the amount of water that can evaporate from the interception 
     storage.  Given this evaporation rate, calculate the amount of time it 
     will take to evaporate the entire amount of intercepted water.  If this 
     time period is shorter than the length of the time interval, the 
     previously wetted leaves can transpire during the remaining part of the 
     time interval */

  /* WORK IN PROGRESS:  the amount of interception storage can be replenished 
     during (low intensity) rainstorms */

  WetEvapRate = WetArea * LocalEvap->EPot[Layer];
  if (WetEvapRate > 0.0) {
    WetEvapTime = *Int/WetEvapRate;

    /* SRW 1-17-01.  WetEvapTime should never exceed Dt!  Not having this check
       was causing transpiration to exceed potential ET in the equation
       LocalEvap->ESoil[Layer][i] *= WetEvapTime * (1 - WetArea) + DryEvapTime; */ 
    if (WetEvapTime > Dt)  WetEvapTime = Dt; 
  }
  else if (*Int > 0.0)
    WetEvapTime = Dt;
  else
    WetEvapTime = 0.0;

  if (WetEvapRate > 0.0) {
    if (WetEvapTime < Dt) {
      LocalEvap->EInt[Layer] = *Int;
      *Int = 0.0;
      DryEvapTime = Dt - WetEvapTime;
    }
    else {
      LocalEvap->EInt[Layer] = Dt * WetEvapRate;
      *Int -= LocalEvap->EInt[Layer];
      DryEvapTime = 0.0;
    }
  }
  else {
    LocalEvap->EInt[Layer] = 0.0;
    if (*Int > 0.0)
      DryEvapTime = 0.0;
    else
      DryEvapTime = Dt;
  }

  /* Correct the evaporation from interception and the interception storage for
     the fractional overstory coverage */
  LocalEvap->EInt[Layer] *= F;
  LocalEvap->ETot += LocalEvap->EInt[Layer];
  *Int *= F;

  /* calculate the canopy conductances associated with the conditions in 
     each of the soil layers */
  for (i = 0; i < VType->NSoilLayers; i++) {
    if( HeatFluxOption == TRUE )
      Temperature = LocalSoil->Temp[i];
    else
      Temperature = Met->Tair;

    Rc[i] = CanopyResistance(VType->LAI[Layer], VType->RsMin[Layer], 
			     VType->RsMax[Layer], VType->Rpc[Layer], 
			     VType->VpdThres[Layer], VType->MoistThres[Layer], 
			     SType->WP[i], Temperature, 
			     LocalSoil->Moist[i], Met->Vpd, Rp);
  }
  
  /* calculate the transpiration rate for the current vegetation layer,
     and adjust the soil moisture content in each of the soil layers */
  for (i = 0; i< VType->NSoilLayers; i++) {

    /* SRW 1-11-01.  Only do transpiration if soil moisture > wilting point. Was
       getting continued transpiration at sometimes high rates below WP.
       Suspect problem lies in using large resistances instead of zero
       conductances.  */
    if (LocalSoil->Moist[i] >= SType->WP[i]) {

      LocalEvap->ESoil[Layer][i] = (Met->Slope + Met->Gamma) / 
        (Met->Slope + Met->Gamma * (1 + Rc[i]/Ra)) * VType->RootFract[Layer][i] * 
	    LocalEvap->EPot[Layer] * Adjust[i];

    } else LocalEvap->ESoil[Layer][i] = 0.0;

    /* calculate the amounts of water transpirated during each timestep based 
       on the evaporation and transpiration rates.  While there is still water
       in interception storage only the area that is not covered by intercep-
       ted water will transpire.  When all of the interception storage has 
       disappeared all leaves will contribute to the transpiration */
    LocalEvap->ESoil[Layer][i] *= ((WetEvapTime * (1 - WetArea)) + DryEvapTime);

    SoilMoisture = LocalSoil->Moist[i] * VType->RootDepth[i] * Adjust[i];
    
    if (SoilMoisture < LocalEvap->ESoil[Layer][i]) 
      LocalEvap->ESoil[Layer][i] = SoilMoisture;

    /* correct the evaporation for the fractional overstory coverage and update
       the soil  moisture */
    LocalEvap->ESoil[Layer][i] *= F;
    SoilMoisture -= LocalEvap->ESoil[Layer][i];
    LocalSoil->Moist[i] = SoilMoisture/( VType->RootDepth[i] * Adjust[i]);
  }  /* end of soil layer loop for transpiration */
  
  for (i = 0, LocalEvap->EAct[Layer] = 0; i < VType->NSoilLayers; i++) {
    LocalEvap->ETot += LocalEvap->ESoil[Layer][i];
    LocalEvap->EAct[Layer] += LocalEvap->ESoil[Layer][i];
  }

  /* clean up */
  free(Rc);
}





