/*
 * SUMMARY:      CanopyResistance.c - Calculate canopy resistance
 * 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-MOD: Tue May 22 12:47:39 2001 by Scott Waichler <waichler@tuff.pnl.gov>
 * DESCRIPTION:  Calculate canopy resistance
 * DESCRIP-END.
 * FUNCTIONS:    CanopyResistance()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: CanopyResistance.c,v 1.4 2001/06/01 19:29:02 markw Exp $";
#endif /* lint */

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

/*****************************************************************************
  CanopyResistance()
*****************************************************************************/
float CanopyResistance(float LAI, float RsMin, float RsMax, float Rpc, 
		       float VpdThres, float MoistThres, float WP, 
		       float TSoil, float SoilMoisture, float Vpd, float Rp)
{
  float MoistFactor;		/* multiplier for resistance due to soil
				   moisture feed-back */
  float Resistance;		/* Canopy resistance (s/m) */
  float RpFactor;		/* multiplier for resistance due to light
				   level feed-back */
  float TFactor;		/* multiplier for resistance due to soil 
				   temperaure feed-back */
  float VpdFactor;		/* multiplier for resistance due to vapor
				   pressure deficit feed-back */

  if (TSoil <= 0) {
    Resistance = DHSVM_HUGE;
    return Resistance;
  }

  /* for OBS */
  TFactor = 1.0/(0.176 + 0.0770 * TSoil - 0.0018 * TSoil*TSoil);

  /* for OJP */
/*   TFactor = 1.0/(.0705 * TSoil - 0.0013 * pow(TSoil, (double) 2.0)); */

  if (TFactor <= 0) {
    Resistance = DHSVM_HUGE;
    return Resistance;
  }
  
  /* equation 14, Wigmosta et al [1994] */

  if (Vpd >= VpdThres) {
    Resistance = DHSVM_HUGE;
    return Resistance;
  }
  else
    VpdFactor = 1.0/(1 - Vpd/VpdThres);
  
  /* equation 15, Wigmosta et al [1994 */
  /* SRW 2-1-01 adds condition for when Rp is close to 0.  
  if (Rp < 0.001) {
    Resistance = DHSVM_HUGE;
    return Resistance;
  }
  else */
    RpFactor = 1.0/((RsMin/RsMax + Rp/Rpc)/(1 + Rp/Rpc));

  /* equation 16, Wigmosta et al [1994] */

  if (SoilMoisture <= WP) {
    Resistance = DHSVM_HUGE;
    return Resistance;
  }
  else if (SoilMoisture < MoistThres)
    MoistFactor = (MoistThres - WP)/(SoilMoisture - WP);
  else
    MoistFactor = 1.0;

  Resistance = TFactor * VpdFactor * RpFactor * MoistFactor * RsMin/LAI; 
/*
printf("\n 4 %f %f %f %f %f %f \n", TFactor, VpdFactor, RpFactor, MoistFactor, RsMin, LAI);
*/
  return Resistance;
  
}

