PIHMgis: pihm/et_is.cpp Source File

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