!************************************************************************* ! ! Using Penman-Monteith Method to Calculate the evapotranspirationation ! !************************************************************************* subroutine evappm(nowcrp,itype,nsl,j2,j1, & & elevm,radly,tave,tdpt,tmax,tmin,radpot,vwind,eo, & & lai, rtd, kcb, canhgt, kcbcon, solthk, thetfc, thetdr, & & st, fin, rawp, resint, cv, es, ep, et) implicit none include 'constants.inc' !************************************************************************* ! ! nowcrp - ? ! itype - plant type, plant growth scenario index(?) ! nsl - number of soil layers ! j2 - soil layers(?) ! j1 - soil layers subjected to soil evaporation ! elevm - elevation of the climate station in meters ! radly - daily solar radiation (Langleys/day) ! tave - average daily temperature (degrees C) ! tdpt - mean daily dew point temperature (degrees C) ! tmax - maximum daily temperature (degrees C) ! tmin - minimum daily temperature (degrees C) ! radpot - potential radiation at slope (cal/cm^2/day) ! vwind - mean daily wind speed (m/s) ! eo - potential daily evapotranspiration (m/day) ! lai - leaf area index ! rtd - root depth (m) ! kcb - mid-season crop coefficient for Penman-Monteith ! canhgt - canopy height (m) ! kcbcon - basal crop coefficient ! solthk - cumulative thickness of soil layer (m) ! thetfc - 1/3-bar soil water content (field capacity) ! thetdr - 15-bar soil water content (wilting point) ! st - current available water content per soil layer (m) ! fin - infiltrated water amount ! rawp - coefficient p for readily available root zone ! resint - interception of rain water by flat residue(m) ! cv - residue amount ! es - soil evaporation (m/day) ! ep - plant transpiration (m/day) ! et - ? ! !************************************************************************** integer, intent(in) :: nowcrp,itype(mxcrop),nsl real, intent(in) :: elevm, radly, tave,tdpt,tmax,tmin,radpot,vwind real, intent(in) :: lai, rtd, kcb, canhgt real, intent(in) :: solthk(mxnsl), thetfc(mxnsl), thetdr(mxnsl) real, intent(in) :: fin, rawp(mxcrop), cv integer, intent(inout) :: j1,j2 real, intent(inout) :: eo, kcbcon, es, ep, et, resint, st(mxnsl) ! ! + + + ARGUMENT DEFINITIONS + + + ! elevm - elevation of the climate station in meters ! nowcrp - current crop ! ! ! + + + LOCAL VARIABLES + + + real gma,alb,xx,yy,rto,dlt,ed,ee,fwv,pb, & & ra, ralb1, rbo, rhd, rn, rso, xl, & & emaxt,emint,kcbadj,TEW,REW,wfevp, & & etke,etkr,etks,TAW,RAW,wftrp,etorc,etcsc,rawpaj, & & epdp,tpdp,etcadj,kcmax,eaj,kecon,potes,bpotes integer xitflg,crpindx,i ! ! + + + LOCAL DEFINITIONS + + + ! eaj - soil cover index (actually, the fraction UN-covered.) ! tk - daily average air temperature (degree Kelvin) ! dlt - slope of the saturated vapor pressure ! gma - the second part of Priestly Taylor equation ! alb - Albedo (fraction) ! xx - soil evaporation (potential of stage 1 and 2) ! yy - depth of soil layer which provide water for soil ! evaporation 7.3.1 ! rto - available water for soil evaporation/m ! ho - net radiation, used in the Priestly Taylor equation (unused) ! aph - is 1.28, conversion factor ETp/ETeq (not used) ! xitflg - flag. 1=exit do-140 loop. ! fphu - ratio of total growing degree days received so far, to ! total growing degree days expected at senescence (0-1). ! rhd - Relative humidity. ! ed - saturated water vapor pressure at dew point ! ee - daily average water vapor pressure. ! emaxt - saturated water vapor pressure at maximum temperature. ! emint - saturated water vapor pressure at minimum temperature. ! kcbadj - adjust basal crop coefficient. ! strstg - current crop development stage start date. ! stgend - current crop development stage end date. ! TEW - maximum water that can be evaporated from soil (mm) ! REW - cumulative evaporation at the end of stage 1 (mm) ! wfevp - water available at the begining of the day for evaporation(mm) ! wetting events are assumed occour before evaporation. ! etke - soil evaporation coefficient ! etkr - evaporation reduction coefficient ! etfew - exposed and wetted soil fraction ! etfw - soil fraction wetted by wetting events. ! TAW - total root zone soil water that can be extracted by crop (mm). ! RAW - raidly soil water in the root zone for crop to extract (mm). ! wftrp - water available for traspiration at the begining of the day ! soil water recharge is countered next day. ! etorc - crop ET under refference condition (mm). ! etcsc - crop ET under standard condition (mm). ! rawpaj - adjusted water stress coefficient. ! etcadj - adjusted crop ET (mm) ! kcbcon - basal crop coefficient. ! kcmax - maximum value of Kc following rain or irrigation ! ! + + + DATA INITIALIZATIONS + + + ! ! + + + END SPECIFICATIONS + + + ! ! FAO PENMAN-MONTIETH reference ET0 ! kcmax = 0.0 eaj = 0.0 kecon = 0.0 potes = 0.0 bpotes = 0.0 ! alb = 0.23 ! albedo of hypothetical grass surface ra = radly / 23.9 xl = 2.501 - 0.002361 * tave ralb1 = ra * (1-alb) ! ! FAO 56 page 40. ed = 0.6108 * exp(17.27*tdpt/(tdpt+237.3)) emaxt = 0.6108 * exp(17.27*tmax/(tmax+237.3)) emint = 0.6108 * exp(17.27*tmin/(tmin+237.3)) ee = 0.5*(emaxt+emint) ! ! radpot is calculated in winter routines only in the winter ! it needs to be calculated every day therefore the calculation ! in sunmap must be called from contin not in winter ! as in this version JEF 3/20/91 rso = radpot / 23.9 ! ! -------- estimate the net outgoing longwave radiation (RBO) rbo = (0.34-0.14*sqrt(ed)) * 4.9e-9 * ((tmax+273.2) ** 4 + & & (tmin+273.2)**4)/2*(1.35*ra/rso-0.35) ! ! -------- calculate net radiation (RN) rn = ralb1 - rbo ! ! NOTE - The wind generated ! by CLIGEN is for a 10 meter height ! fwv = vwind * 4.87/(alog(67.8*10.0-5.42)) dlt = 4098.0/(tave+237.3)**2 * 0.6108 * exp(17.27*tave/(tave+237.3)) pb = 101.3*(1.0- 0.0065*elevm/293.0)**5.26 ! ------ compute psychrometric constant kpa/c (GMA) ! gma = 0.000665 * pb ! ! the unit of eo is (m/d) etorc = (0.408 * dlt * rn + & & gma*(900.0/(tave+273))*(ee-ed)*fwv) & & /(dlt+ gma *(1.+ 0.34* fwv)) ! eo = etorc*0.001 ! Dual crop coefficient approach to crop evapotranspiration ! under standard condition ! crpindx = itype(nowcrp) rhd = ed/emaxt*100 ! if((lai.gt.0.0).and.(rtd.gt.0.0)) then ! Adjustment of middle season crop coefficient ! for climate condition kcbadj = kcb(crpindx) +(0.04*(fwv-2.)-0.004*(rhd-45))*(canhgt/3.)**0.3 else kcbadj = 0 endif !d According to Stöckle et al ! basal crop coefficient kcbcon = kcbadj*(1-exp(-0.45*lai)) ! soil evaporation coefficient if(kcbadj.gt.0.0) then etke = kcbadj*exp(-0.45*lai) else etke = 1.2 endif ! ! calculation of soil evaporation reduction coefficient ! In the TEW formula 0.1m is the depth of the surface soil layer ! that is subject to drying by way of evaporation 0.10~0.15m. ! REW is a formula from linear regression of Table 19 (FAO56 p144) ! Notes for previous day water depletion calculation: ! st(mxnsl,mxplan) is available water content per soil layer(m) ! at the end of previous day. This value is from the difference ! between soil water content and wilting point soil water ! content. Because the available water for evaporation is from ! the difference between water content and half of wilting ! point water content. So we add half of wilting point water ! content to calculate available water for using st value. ! TEW = 0.0 REW = 0.0 wfevp = 0.0 ! epdp = 0.1 ! do 10 i = 1, nsl if (i.eq.1) then yy = 0.0 else yy = solthk(i-1) endif if (solthk(i).lt.epdp) then TEW = TEW+1000*(thetfc(i)- 0.3*thetdr(i))*(solthk(i)-yy) REW = REW+(57.856*(thetfc(i)-thetdr(i))+0.280)*(solthk(i)-yy)/epdp wfevp = wfevp+1000*(st(i) + 0.7*thetdr(i)*(solthk(i)-yy)) else TEW =TEW + 1000*(thetfc(i)- 0.3*thetdr(i))*(epdp-yy) REW = REW+(57.856*(thetfc(i)-thetdr(i))+0.280)*(epdp-yy)/epdp wfevp = wfevp+1000*(st(i)* & & (epdp-yy)/(solthk(i)-yy)+ & & 0.7*thetdr(i)*(epdp-yy)) goto 20 endif 10 continue 20 continue ! ! adding in wetting event influence before evaporation wfevp = wfevp + fin*1000 ! if ((TEW-wfevp).le.REW) then etkr = 1.0 else etkr = (wfevp/(TEW - REW))**2 endif ! ! Soil water stress on crop ! Double count the avaiable water for evaporation and transpiration ! for the top 0.1~0.15m. But I think it get correction when water ! is actually substract from its availvabel water. TAW = 0.0 RAW = 0.0 wftrp = 0.0 ! tpdp = rtd if (solthk(nsl).lt.tpdp) then tpdp = solthk(nsl) endif ! do 30 i = 1, nsl if (i.eq.1) then yy = 0.0 else yy = solthk(i-1) endif if (solthk(i).lt.tpdp) then TAW = TAW+1000*(thetfc(i)-thetdr(i))*(solthk(i)-yy) wftrp = wftrp+st(i)*1000 else TAW =TAW + 1000*(thetfc(i)-thetdr(i))*(tpdp-yy) wftrp = wfevp+st(i)*1000*(tpdp-yy)/(solthk(i)-yy) goto 40 endif 30 continue 40 continue ! etcsc = kcbadj*etorc rawpaj = rawp(crpindx) +0.04*(5.-etcsc) RAW = rawpaj*TAW ! if ((TAW-wftrp).le.RAW) then etks = 1.0 else etks = wftrp/(TAW - RAW) endif ! ! ************************************************************************ ! ** This section of the code computes soil evaporation with residue ** ! ************************************************************************ ! potes = etorc*etke*0.001 if (potes.gt.resint) then bpotes = potes - resint ! ---- compute the fraction of soil that is uncovered (EAJ) ! (WEPP Equation 7.2.3) cv : residue amount eaj = exp(-0.5*(cv+.1)) kcmax = 1.2 +(0.04*(fwv-2.)-0.004*(rhd-45))*(canhgt/3.)**0.3 ! kecon = min(etke*etkr, eaj*kcmax) es = kecon*bpotes/etke + resint else es = potes endif ! ep = etorc*etks*kcbcon*0.001 et = es + ep ! ! water redistribution in the soil xx = es if (resint.ge.0.0) then xx = es - resint resint = 0.0 end if if (xx.lt.0.0) then st(1) = st(1)-xx xx = 0.0 goto 60 endif ! Evaporate maximum potential water from each soil layer. xitflg = 0 j2 = 0 50 continue j2 = j2 + 1 ! ------ if cum. soil depth exceeds "soil evaporated depth"... if (solthk(j2).gt.epdp) then j1 = j2 - 1 yy = 0. if (j1.gt.0) yy = solthk(j1) rto = st(j2) * (epdp-yy) / (solthk(j2)-yy) ! if (rto.gt.xx) then st(j2) = st(j2) - xx if (st(j2).lt.1e-10) st(j2) = 0.00 else xx = xx - rto st(j2) = st(j2) - rto if (st(j2).lt.1e-10) st(j2) = 0.00 es = es - xx et = et - xx end if xitflg = 1 ! ! ------ if water available in soil layer exceeds potential soil evap... else if (st(j2).gt.xx) then st(j2) = st(j2) - xx if (st(j2).lt.1e-10) st(j2) = 0.00 xitflg = 1 ! else xx = xx - st(j2) st(j2) = 0. end if ! ------ (force immediate exit from loop.) if (xitflg.ne.0) j2 = nsl if (j2.lt.nsl.and.xitflg.eq.0) go to 50 ! if (xitflg.eq.0) then es = es - xx et = et - xx end if 60 continue return end subroutine ! end module