PIHM: et_is.c Source File

00001 /*******************************************************************************
00002  * File        : et_is.c                                                       *
00003  * Function    : for calculation of evapotranspiration and interception        *
00004  * Programmers : Yizhong Qu   @ Pennsylvania State Univeristy                  *
00005  *               Mukesh Kumar @ Pennsylvania State Univeristy                  *
00006  *               Gopal Bhatt  @ Pennsylvania State Univeristy                  *
00007  * Version     : 2.0 (July 10, 2007)                                           *
00008  *-----------------------------------------------------------------------------*
00009  *                                                                             *
00010  * Interception process and evapotranspiration are considered as weakly coupled*
00011  * processes compared with overland, channel routing and groundwater flow pro- *
00012  * cesses. Therefore, before each time step, interception is calculated and    *
00013  * deducted. After each time step, envapotranspiration is calculated and all   *
00014  * state variables are adjusted accordingly. In this file, ET is calculated    *
00015  * using Penman-Monteith equation. Most parameters are calculated based on     *
00016  * Crop Evapotranspiration (FAO No.56).                                        *
00017  *                                                                             *
00018  * This code is free for users with research purpose only, if appropriate      *
00019  * citation is refered. However, there is no warranty in any format for this   *
00020  * product.                                                                    *
00021  *                                                                             *
00022  * For questions or comments, please contact the authors of the reference.     *
00023  * One who want to use it for other consideration may also contact Dr.Duffy    *
00024  * at cxd11@psu.edu.                                                           *
00025  *******************************************************************************/
00026 
00028 
00029  /* TODO */
00030  //eliminate 'ret' variable by using EleTF
00031 
00032 /* C Header Files */
00033 #include <stdio.h>
00034 #include <stdlib.h>
00035 #include <math.h>
00036 #include <string.h>
00037 
00038 /* SUNDIALS Header Files */
00039 #include "nvector_serial.h"
00040 #include "sundials_types.h"
00041 
00042 /* PIHM Header Files */
00043 #include "pihm.h"
00044 #include "calib.h"
00045 //#include "et_is.h"
00046 
00047 #define EPSILON 0.05
00048 
00049 /* Calibration Parameters */
00050 realtype is_CALIB;                
00051 realtype et0_CALIB;               
00052 realtype mf_CALIB;                
00053 realtype tf_CALIB;                
00056 realtype Interpolation(TSD *Data, realtype t);        /* Data Value at time=t from a TimeSeries  */
00057 
00058 
00059 /********************************************************************
00060     Calculates Evaporation from Canopy & Interception Storage
00061 *********************************************************************/
00062 void calET_IS(realtype t, realtype stepsize, Model_Data MD, N_Vector VY)
00064 
00069 {
00070       int i;
00071       realtype totEvap;
00072       realtype Delta, Gamma;
00073       realtype Rn, G, T, Vel, RH, VP,P,LAI,zero_dh,cnpy_h,rl,r_a;
00074       realtype isval=0,etval=0;
00075       realtype fracSnow,snowRate,MeltRate,MF,Ts=-3.0,Tr=1.0,To=0.0,massMelt=0,ret;
00076 
00077       //Model_Data MD;
00078 
00079       is_CALIB = setis_CALIB();
00080       et0_CALIB = setet0_CALIB();
00081       mf_CALIB=setmf_CALIB();
00082       tf_CALIB = settf_CALIB();
00083 
00084       //MD = (Model_Data)DS;
00085 
00086       stepsize=stepsize/(24.0*60.0);
00087       for(i=0; i<MD->NumEle; i++)
00088       {
00089         MD->ElePrep[i] = Interpolation(&MD->TSD_Prep[MD->Ele[i].prep-1], t);
00090         Rn = Interpolation(&MD->TSD_Rn[MD->Ele[i].Rn-1], t);
00091         T = Interpolation(&MD->TSD_Temp[MD->Ele[i].temp-1], t);
00092         Vel = Interpolation(&MD->TSD_WindVel[MD->Ele[i].WindVel-1], t);
00093         RH = Interpolation(&MD->TSD_Humidity[MD->Ele[i].humidity-1], t);
00094         VP = Interpolation(&MD->TSD_Pressure[MD->Ele[i].pressure-1], t);
00095         P = 101.325*pow(10,3)*pow((293-0.0065*MD->Ele[i].zmax)/293,5.26);
00096         LAI = Interpolation(&MD->TSD_LAI[MD->Ele[i].LC-1], t);
00097         MF = Interpolation(&MD->TSD_MeltF[0], t);
00098         MF=mf_CALIB*MF;
00099 
00100         /*****************************************Snow Calculation ****************************************************/
00101         fracSnow = T<Ts?1.0:T>Tr?0:(Tr-T)/(Tr-Ts);
00102         snowRate = fracSnow*MD->ElePrep[i];
00103         MD->EleSnow[i]=MD->EleSnow[i]+snowRate*stepsize;
00104         MeltRate=(T>To?(T-To)*MF:0);
00105         //MeltRate=(T>To?Rn*(T-To)*MF:0);
00106         if(MD->EleSnow[i]>MeltRate*stepsize)
00107         {
00108             MD->EleSnow[i]=MD->EleSnow[i]-MeltRate*stepsize;
00109             //MeltRate=(T>To?Rn*(T-To)*MF/60.0:0);
00110         }
00111         else
00112         {
00113             MeltRate=MD->EleSnow[i]/stepsize;
00114             MD->EleSnow[i]=0;
00115         }
00116         /**************************************************************************************************************/
00117 
00118 
00119 
00120         /**************************************Evaporation from canopy*************************************************/
00121         MD->EleISmax[i] = MD->SIFactor[MD->Ele[i].LC-1]*Interpolation(&MD->TSD_LAI[MD->Ele[i].LC-1], t);
00122         //MD->EleIS[i] = is_CALIB*MD->EleISmax[i];
00124         MD->EleISmax[i] = is_CALIB*MD->EleISmax[i];
00125         if(LAI>0.0)
00126         {
00127             Delta = 2503*pow(10,3)*exp(17.27*T/(T+237.3))/(pow(237.3 + T, 2));
00128             Gamma = P*1.0035*0.92/(0.622*2441);
00129             zero_dh=Interpolation(&MD->TSD_DH[MD->Ele[i].LC-1], t);
00130             //zero_dh=0;
00131             cnpy_h = zero_dh/(1.1*(0.0000001+log(1+pow(0.007*LAI,0.25))));
00132             /*if(LAI<2.85)
00133             {
00134                 rl= 0.0002 + 0.3*cnpy_h*pow(0.07*LAI,0.5);
00135             }
00136             else
00137             {
00138                 rl= 0.3*cnpy_h*(1-(zero_dh/cnpy_h));
00139             }
00140             */
00141             rl=Interpolation(&MD->TSD_DH[MD->Ele[i].LC-1], t);
00142             r_a = log(MD->Ele[i].windH/rl)*log(10*MD->Ele[i].windH/rl)/(Vel*0.16);
00143 
00144             /* $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ DELETE 0.01 BELOW $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ */
00145 
00146             MD->EleET[i][0] = (LAI/MD->Ele[i].LAImax)*(pow((MD->EleIS[i]<0?0:MD->EleIS[i])/MD->EleISmax[i],2.0/3.0))*(Rn*(1-MD->Ele[i].Albedo)*Delta+(1.2*1003.5*((VP/RH)-VP)/r_a))/(1000*2441000.0*(Delta+Gamma));
00147             MD->EleTF[i]=tf_CALIB*5.65*pow(10,-2)*MD->EleISmax[i]*exp(3.89*(MD->EleIS[i]<0?0:MD->EleIS[i])/MD->EleISmax[i]);
00148             //printf("\n %f %f %f %f %f %f %f %f %f",MD->EleIS[i],LAI,MD->Ele[i].LAImax,r_a,rl,cnpy_h,Delta,Gamma,MD->EleET[i][0]);
00149 
00150         }
00151          else
00152          {
00153              MD->EleET[i][0]=0.0;
00154             MD->EleTF[i]=0.0;
00155          }
00156 
00157         /**********************************************************************************************************************/
00158 
00159         MD->EleET[i][0]=et0_CALIB*MD->EleET[i][0];
00160         //printf("\n%f %f",MD->EleIS[i],MD->EleISmax[i]); getchar();
00161         if(MD->EleIS[i] >= MD->EleISmax[i])
00162         {
00163               if((1-fracSnow)*MD->ElePrep[i]>=MD->EleET[i][0]+MD->EleTF[i])
00164               {
00165                   MD->Ele2IS[i] = MD->EleET[i][0];
00166                 ret = MD->EleTF[i];
00167               }
00168              else if(((1-fracSnow)*MD->ElePrep[i]<MD->EleET[i][0]+MD->EleTF[i])&&(MD->EleIS[i]+stepsize*((1-fracSnow)*MD->ElePrep[i]-MD->EleET[i][0]-MD->EleTF[i])<=0))
00169              {
00170                  MD->Ele2IS[i] =(1-fracSnow)*MD->ElePrep[i];
00171                  MD->EleET[i][0]=(MD->EleET[i][0]/(MD->EleET[i][0]+MD->EleTF[i]))*(MD->EleIS[i]/stepsize+(1-fracSnow)*MD->ElePrep[i]);
00172                 MD->EleIS[i] = 0;
00173                 ret =(MD->EleTF[i]/(MD->EleET[i][0]+MD->EleTF[i]))*(MD->EleIS[i]/stepsize+(1-fracSnow)*MD->ElePrep[i]);
00174              }
00175              else
00176              {
00177                  MD->Ele2IS[i] =(1-fracSnow)*MD->ElePrep[i];
00178                  MD->EleIS[i] = MD->EleIS[i]+stepsize*((1-fracSnow)*MD->ElePrep[i]-MD->EleET[i][0]-MD->EleTF[i]);
00179                 ret =  MD->EleTF[i];
00180               }
00181 
00182             /*MD->Ele2IS[i] = MD->ElePrep[i]>=MD->EleET[i][0]?MD->EleET[i][0]:MD->ElePrep[i];
00183               MD->EleIS[i] = MD->ElePrep[i]>=MD->EleET[i][0]?MD->EleIS[i]:(MD->EleIS[i]+MD->ElePrep[i]-MD->EleET[i][0]<0?0:MD->EleIS[i]+MD->ElePrep[i]-MD->EleET[i][0]);
00184               MD->EleET[i][0]=MD->EleIS[i]/stepsize+MD->ElePrep[i];
00185             */
00186         }
00187         else if(((MD->EleIS[i] + ((1-fracSnow)*MD->ElePrep[i]-MD->EleET[i][0]-MD->EleTF[i])*stepsize) >= MD->EleISmax[i]))
00188         {
00189               MD->Ele2IS[i] =  (MD->EleISmax[i]/stepsize - (MD->EleIS[i]/stepsize-MD->EleET[i][0]-MD->EleTF[i]));
00190               MD->EleIS[i] = MD->EleISmax[i];
00191               ret =  MD->EleTF[i];
00192         }
00193         else if(((MD->EleIS[i] + ((1-fracSnow)*MD->ElePrep[i]-MD->EleET[i][0]-MD->EleTF[i])*stepsize) <=0))
00194         {
00195               MD->Ele2IS[i] =  (1-fracSnow)*MD->ElePrep[i];
00196               if((MD->EleET[i][0]>0)||(MD->EleTF[i]>0))
00197             {
00198                 MD->EleET[i][0] = (MD->EleET[i][0]/(MD->EleET[i][0]+MD->EleTF[i]))*(MD->EleIS[i]/stepsize+(1-fracSnow)*MD->ElePrep[i]);
00199                 ret =(MD->EleTF[i]/(MD->EleET[i][0]+MD->EleTF[i]))*(MD->EleIS[i]/stepsize+(1-fracSnow)*MD->ElePrep[i]);
00200             }
00201               else
00202             {
00203                 MD->EleET[i][0] = 0;
00204                 ret = 0;
00205             }
00206               MD->EleIS[i] = 0;
00207         }
00208         else
00209         {
00210               MD->Ele2IS[i] = (1-fracSnow)*MD->ElePrep[i];
00211               MD->EleIS[i] = MD->EleIS[i] + ((1-fracSnow)*MD->ElePrep[i]-MD->EleET[i][0]-MD->EleTF[i])*stepsize;
00212                ret =  MD->EleTF[i];
00213         }
00214         massMelt=massMelt+MeltRate*stepsize;
00215         MD->EleNetPrep[i] = (1-MD->Ele[i].VegFrac)*(1-fracSnow)*MD->ElePrep[i] + ((1-fracSnow)*MD->ElePrep[i]+ret - MD->Ele2IS[i])*MD->Ele[i].VegFrac+MeltRate;
00216         MD->EleTF[i] = ret;
00217 
00218         //MD->EleNetPrep[i] =MD->ElePrep[i];
00219 
00220 
00221         //printf("\n%d %f %f %f %f %f %f %e %e %e",i,MD->EleIS[i],MD->ElePrep[i],MD->Ele2IS[i],et_Y[i + MD->NumEle],LAI,zero_dh,MD->EleET[i][0],MD->EleET[i][1],(Rn*(1-MD->Ele[i].Albedo)*Delta+(1.2*1003.5*((VP/RH)-VP)/r_a))/(1000*2441000.0*(Delta+Gamma)));
00222         //printf("\n%d %f %f %f %f %f",i,MD->ElePrep[i],ret,MD->EleTF[i],MD->Ele2IS[i],MD->EleNetPrep[i]);
00223         //gtchar();
00224       }
00225     //printf("\t%e",massMelt);
00226     //getchar();
00227 }