!************************************************************************** ! ! Calculate the evaporation from bare soil. First, compute the ! amount of potential evaporation as if all the soil were bare ! (EO); second, prorate EO by the fraction of soil that actually ! is bare (EAJ) to estimate water demand (EOS); third, determine ! how much water is actually available to satisfy this demand. ! Calculations are performed using Ritchie's model. ! !************************************************************************** subroutine evap(iwind,mon,j2,j1,nsl,elevm,cv,tave,lai, & & salb,radly,obmaxt,obmint,eo,tdpt,radpot,vwind, & & canhgt,resint,s1,tu,s2,fin,su,es,tv,ep,et,st,solthk) !************************************************************************ ! iwind - wind data flag in climate input file ( 0 - wind data, 1 - no wind data ) ! mon - month of year ! j2 - soil layers(?) ! j1 - soil layers subjected to soil evaporation ! nsl - number of soil layers ! elevm - elevation of the climate station in meters ! cv - residue amount ! tave - average daily temperature (degrees C) ! lai - leaf area index ! salb - soil albedo (0-1) ! radly - daily solar radiation (Langleys/day) ! obmaxt - observed average monthly maximum temperature (degrees C) ! obmint - observed average monthly minimum temperature (degrees C) ! eo - potential daily evapotranspiration (m/day) ! tdpt - mean daily dew point temperature (degrees C) ! radpot - potential radiation at slope (cal/cm^2/day) ! vwind - mean daily wind speed (m/s) ! canhgt - canopy height (m) ! resint - interception of rain water by flat residue(m) ! s1 - stage 1, soil evaporation ! tu - upper limit soil evaporation, stage 1 ! s2 - stage 2, soil evaporation ! fin - infiltrated water amount ! su - soil water available for evaporation ! es - soil evaporation (m/day) ! tv - ? ! ep - plant transpiration (m/day) ! et - ? ! st - current available water content per soil layer (m) ! solthk - cumulative thickness of soil layer (m) !*********************************************************************** implicit none include 'constants.inc' integer, intent(in) :: iwind, mon, nsl real, intent(in) :: elevm,cv,tave, lai, salb, radly, obmaxt(12) real, intent(in) :: obmint(12), tdpt, radpot, canhgt real, intent(in) :: tu real, intent(in) :: solthk(mxnsl) integer, intent(inout) :: j2, j1 real, intent(inout) :: eo, vwind, s2, s1, su, es, tv real, intent(inout) :: fin, ep, et, resint, st(mxnsl) ! ! + + + LOCAL VARIABLES + + + real eaj, tk, delta, gma, alb, eos, sp, sb, esx, xx, yy, rto, & & d1, datd, dlt, ed, ee, eon, etoh, fwv, harrad, hc, pb, & & ra, ra1, ralb1, rbo, rc, rhd, rn, rso, ttdd, xl, xxll, zm, zom double precision tmpvr1 integer xitflg, ievap ! ! + + + LOCAL DEFINITIONS + + + ! eaj - soil cover index (actually, the fraction UN-covered.) ! tk - daily average air temperature (degree Kelvin) ! delta - slope of the saturated vapor pressure ! gma - the second part of Priestly Taylor equation ! alb - Albedo (fraction) ! eos - potential daily evapotranspiration adjusted for ! cover m/day ! sp - soil evaporation - infiltrated water (stage 1) ! sb - adjusted soil evaporation for water infiltration ! esx - evaporation from the infiltrated water ! 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. ! ! + + + DATA INITIALIZATIONS + + + ! ! + + + END SPECIFICATIONS + + + ! ! ************************************************************************ ! ** This section of the code computes evapotranspiration (EO). ** ! ************************************************************************ ! ! ---- compute the fraction of soil that is uncovered (EAJ) ! (WEPP Equation 7.2.3) eaj = exp(-0.5*(cv+.1)) ! ! ---- convert average daily temp to degrees Kelvin tk = tave + 273. ! ! ... Compute vapor pressure slope, SVP ! (WEPP Equation 7.2.4) ! (EPIC Equation 2.53) ! delta = exp(21.255-5304.0/tk) * 5304.0 / tk ** 2 gma = delta / (delta+0.68) ! ! ---- Compute corrected soil albedo (ALB), using albedo input by user ! (SALB) and soil cover index (EAJ). ! (WEPP Eq. 7.2.2) ! if (lai.gt.0.0) then alb = 0.23 * (1.0-eaj) + salb * eaj else alb = salb end if ! ! ... Compute potential evaporation from bare soil (EO) ! ! ievap - flag that determines which equation is used ! for PET calculations. ! ievap = 1 HARGRAVES ! ievap = 2 PRIESTLEY-TAYLOR ! ievap = 3 PENMAN ! ievap = 4 PENMAN-MONTIETH ! ! if wind data not available use HARGRAVES equation ! ---- (Select the HARGRAVES Equation) ! if(iwind .eq. 1)ievap = 1 ! ! ---- (Select the Priestley-Taylor Equation) if (iwind.eq.1) ievap = 2 ! ! ---- (Select the Penman-Montieth Equation) ! ievap = 4 ! ! ---- (Select the Penman Equation) ! if wind data available use PENMAN equation if (iwind.eq.0) ievap = 3 ! ! ! ! ! ** L0 IF ** ! PRIESTLEY-TAYLOR EQUATION ! ievap = 1 HARGRAVES ! ievap = 2 PRIESTLEY-TAYLOR ! ievap = 3 PENMAN ! ievap = 4 PENMAN-MONTIETH ! ! (WEPP Eq. 7.2.1) ! !-------------------------------------------------------------- ! HARGRAVES equation * ! iwind=1 * !-------------------------------------------------------------- if (ievap.eq.1) then ! print*,'Hargraves' xxll = 2.501 - 0.002361 * tave if (radly.le.0.0) then ! ! ! EOHA not used 12-16-93 02:04pm sjl ! eoha = 0.0 else harrad = radly / xxll ! ---- Note that harrad is extraterrestral radiation mm/day ! ttdd = sqrt(obmaxt(mon)-obmint(mon)) datd = tave + 17.8 if (datd.lt.0.0) datd = 0.0 etoh = 0.0023 * harrad * datd * ttdd ! ------ etoh is daily potential et in mm/day calculated by !------- Hargraves 1985 equation ! eo = etoh / 1000. ! eo = daily potential et in m/day end if ! ! ** L0 ELSE IF** else if (ievap.eq.2) then !------------------------------------------------------------- ! PRIESTLEY-TAYLOR !------------------------------------------------------------- ! print*,'Priestley-Taylor' eo = 0.00128 * (radly*(1.0-alb)/58.3) * gma ! else !------------------------------------------------------------- ! calculate PENMAN and PENMAN-MONTIETH common portions !------------------------------------------------------------- ra = radly / 23.9 ! ! ------ estimate latent heat of vaporization (XL) ! (EPIC Equation 2.38) xl = 2.501 - 0.0022 * tave ! ! Estimate Saturation vapor pressure (EE), using the average ! temperature in degrees Kelvin (TK). ! (EPIC Equation 2.39) ee = 0.1 * exp(54.879-5.029*alog(tk)-6790.5/tk) ! ! ------ calculate relative humidity (RHD) [fraction] from dewpoint ! temperature (TDPT) and average temperature (TAVE). ! (Rosenberg, Equations 5.17 & 5.18) rhd = exp((17.27*tdpt)/(237.3+tdpt)) / exp((17.27*tave)/(237.3+tave)) if (rhd.gt.1.0) rhd = 1.0 ! ! Calculate the vapor pressure (ED) from the saturation vapor ! pressure (EE) and the relative humidity (RHD). ! (EPIC Equation 2.40) ed = ee * rhd ! ralb1 = ra * (1-alb) ! ! Estimate the slope of the saturation vapor pressure curve (DLT). ! (EPIC Equation 2.41) dlt = ee * (6790.5/tk-5.029) / tk ! ! ------ compute barometric pressure in kpa (PB) ! (Approximation of EPIC Equation 2.43) pb = 101.3 - (0.01055*elevm) ! ! ------ compute psychrometric constant kpa/c (GMA) ! (EPIC Equation 2.42) gma = 0.00066 * pb ! xx = dlt + gma ! end if ! if (ievap.eq.3) then !---------------------------------------------------------------------- ! PENMAN calculation !---------------------------------------------------------------------- ! print*,'Penman' ! 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 ! ! -------- calculate mean daily wind speed (vwind) ! (Modified EPIC Equation 2.50) ! fwv = 2.7 + 0.537 * vwind ! ! NOTE - The value 1.63 is used below because the wind generated ! by CLIGEN is for a 10 meter height (not at 2 meters which would ! use the 0.537 value above). dcf 6/2/94 ! fwv = 2.7 + 1.63 * vwind ! ! -------- estimate the net outgoing longwave radiation (RBO) ! (EPIC Equation 2.46) rbo = (0.34-0.14*sqrt(ed)) * 4.9e-9 * tk ** 4 ! ! -------- calculate net radiation (RN) ! (EPIC Equation 2.45) rn = ralb1 - rbo * (0.9*ra/rso+0.1) if (rn.le.0.0) rn = 0.0 ! ! -------- Compute potential evaporation (EO) ! (Approximation of EPIC Equation 2.37, assuming the soil ! flux is close to zero, and thus negligible.) eo = ((dlt*rn)/xl+gma*fwv*(ee-ed)) * 0.001 / xx ! ! ** L0 ELSE-IF ** else if (ievap.eq.4) then !------------------------------------------------------------------ ! PENMAN-MONTEITH EQUATION ! (This is the only one of the four equations that takes plants ! (and therefore, evapotranspiration, into consideration.) !------------------------------------------------------------------ if (lai.lt.1.0) then rc = 1.0 else rc = 200.0 / lai end if ! if (vwind.le.0.0) vwind = 1.0e-10 hc = canhgt if (hc.lt.1) hc = 0.01 zm = 2.0 zom = 0.123 * hc d1 = 0.67 * hc tmpvr1 = alog((zm-d1)/zom) ra1 = tmpvr1 * (tmpvr1+1.0) / (0.1681*vwind) ! eon = (dlt*ralb1+(gma*(1710.-6.85*tave)*(ee-ed))/ra1) / (dlt+gma*(1.+rc/ra1)) eo = 0.00001712 * eon ! ! ! ** L0 ENDIF ** end if ! ! ! ... Compute potential evaporation from bare soil (EOS), bare soil ! water demand (EO), and fraction of soil bare (EAJ). eos = eo * eaj eos = eos - resint if (eos.le.0.0) eos = 0.0 ! ! ************************************************************************ ! ** The following section computes how much soil water is available to ** ! ** satisfy bare soil evaporation potential (EOS). The physical pro- ** ! ** cess is assumed to occur in two stages. The first (stage 1) is ** ! ** limited by the amount of energy stiking the soil surface. During ** ! ** stage 1 the amount of water evaporated is linearly related to the ** ! ** amount of energy striking the soil surface. During the second ** ! ** stage (stage 2), the relationship becomes non-linear. Increas- ** ! ** ingly more energy is required per unit of water evaporated; ie, ** ! ** evaporation is limited by the soil water content (the amount of ** ! ** water available to be evaporated). The variable TU represents the ** ! ** value on the current soil, at which this transition between stages ** ! ** occurs. Presently TU is set to 6 mm for all soil types. The var- ** ! ** iable S1 represents cumulative evaporation. It is reset whenever ** ! ** precipitation occurs. ** ! ************************************************************************ ! ! ! ** M0 IF ** ! If the cumulative water evaporated doesn't go beyond stage 1.... if (s1.lt.tu) then ! ------ compute adjusted stage 1 evaporation by subtracting infiltration. s2 = 0.0 sp = s1 - fin ! ! ------ if adjusted evaporation is greater than infiltration, set ! stage 1 evaporation (S1) equal to the difference. if (sp.gt.0.0) then s1 = sp else s1 = 0.0 end if ! ! ------ add evapotranspiration (EOS) to stage 1 evaporation (S1) s1 = s1 + eos ! ! ------ compute the water deficit (SU) by subtracting the upper limit ! for soil evaporation (TU) from the stage 1 evaporation (S1). su = s1 - tu ! ! If there is a water deficit, compute stage 2 evaporation. if (su.gt.0) then ! -------- set soil evaporation (ES) equal to evapotranspiration (EOS) ! minus 0.4 of the water deficit (SU). es = eos - .4 * su ! -------- set stage 2 evaporation (S2) equal to 0.6 of the water deficit. s2 = .6 * su ! -- XXX -- What does TV represent? -- CRM -- 9/2/92. tv = (s2/.0035) ** 2 ! ! If there is NOT a water deficit.... else ! ------- set soil evaporation (ES) equal to adjusted evapotranspiration es = eos end if ! ! ! ** M0 ELSE ** ! If the cumulative water evaporated DOES go beyond stage 1.... else ! sb = fin - s2 ! ! ** M1 IF ** ! Is stage 2 evap > 0? ! if (sb.lt.0) then tv = tv + 1 es = .0035 * sqrt(tv) - s2 ! ! If there is infiltration.... if (fin.gt.0) then esx = 0.8 * fin ! ---------- if soil evaporation (ES) exceeds 0.8 infiltration (FIN) if (es.gt.esx) esx = es + fin if (esx.gt.eos) esx = eos es = esx ! ! Otherwise, if soil evaporation (ES) is greater than the ! adjusted evapotranspiration (EOS).... else if (es.gt.eos) then ! ---------- set soil evap. equal to evapo-trans. es = eos end if ! s2 = s2 + es - fin tv = (s2/.0035) ** 2 ! ! ** M1 ELSE ** ! Stage 2 evap <= 0 ! else fin = sb s1 = tu - fin tv = 0.0 s2 = 0.0 if (s1.lt.0) s1 = 0.0 s1 = s1 + eos su = s1 - tu ! ! Soil water avail for evap > 0? ! if (su.gt.0) then es = eos - .4 * su s2 = .6 * su tv = (s2/.0035) ** 2 else es = eos end if ! ! ** M1 ENDIF ** end if ! ! ** M0 ENDIF ** end if ! ! Add water intercepted by residue back onto soil evaporation ! es = es + resint ! if (es.le.0) es = 0.0 ! ! Compute plant transpiration ! (WEPP Equation 7.2.10) ! (EPIC Equation 2.56 & 2.57) ! if (lai.gt.3.0) then ep = eo else ep = lai * eo / 3.0 end if ! et = es + ep if (eo.lt.et) then et = eo es = et - ep end if 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 10 continue j2 = j2 + 1 ! ------ if cum. soil depth exceeds "soil evaporated depth"... if (solthk(j2).gt.0.10) then j1 = j2 - 1 yy = 0. if (j1.gt.0) yy = solthk(j1) rto = st(j2) * (0.10-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 10 ! if (xitflg.eq.0) then es = es - xx et = et - xx end if ! 60 continue return end subroutine !end module