! FILE evap.f ! ! SUBROUTINES IN THIS FILE CALLED FROM ! subroutine evap(j, je, k, n) hydrotop module evapotranspiration use utilities, only : dp, log_warn, log_debug implicit none ! general potential evap calibration factor real(dp), save, dimension(:), allocatable :: bsn_ecal ! correction factor for potential evapotranspiration; range for thc: (0.8 - 1.0), value 1. - from R. Muttiah real(dp), save, dimension(:), allocatable :: bsn_thc ! parameter used to calc day length in crop real(dp) :: pit = 58.13 real(dp), save, dimension(12) :: turc_ivanov = & (/ 1.15, 1.15, 1.15, 1.15, 1.15, .85, .85, .85, .85, .85, 1.15, 1.15/) ! potential evapotranspiration, mm real(dp), save :: eo ! soil evaporation, mm real(dp), save :: es ! es + ep + canev real(dp), save :: et real(dp), save, dimension(12) :: omega = (/ .7, .85, .95, 1.05, 1.25, 1.15, 1.05, .95, .9, .8, .75, .7/) ! month factors for Turc(Ivanov) evap (Glugla 1989) ! snow evaporation, mm real(dp), save :: snoev ! canopy evaporation, mm / INTERCEPTION real(dp), save :: canev ! potential evapotranspiration, mm / INTERCEPTION real(dp), save :: eopot ! integer, save :: iemeth ! integer, save :: idvwk ! coefs for each subbasin and for entire basin real(dp), save, dimension(:), allocatable :: ecal real(dp) :: ecal0 = 1. ! coefs for each subbasin and for entire basin real(dp), save, dimension(:), allocatable :: thc real(dp) :: thc0 = 1. ! empirical constant in Hargreaves (1985) equation, readcli.f real(dp), save :: ec1 = 0.135 ! soil albedo real(dp), save, dimension(:), allocatable :: salb real(dp), save, dimension(:), allocatable :: yls real(dp), save, dimension(:), allocatable :: ylc ! daily min temp. for subbasin, readcli, degree C real(dp), save, dimension(:), allocatable :: tmn ! solar radiation in subbasin, J / cm^2 real(dp), save, dimension(:), allocatable :: ra ! air humidity in the subbasin, % real(dp), save, dimension(:), allocatable :: humi ! latitude of subbasin centroid real(dp), save, dimension(:), allocatable :: lat ! internal func. for Richie's method to estimate es real(dp), save, dimension(:, :), allocatable :: s1 ! internal func. for Richie's method to estimate es real(dp), save, dimension(:, :), allocatable :: s2 ! internal func. for Richie's method to estimate es real(dp), save, dimension(:, :), allocatable :: tv ! canopy water storage, mm real(dp), save, dimension(:, :), allocatable :: canstor ! 0 = read radiation data from clim1; 1 = calculate radiation after Hargreaves (latitude in degrees required) integer, save :: radiation_switch = 0 ! Output variable ids integer :: & eta_output_id = 0, & etp_output_id = 0, & soil_evaporation_output_id namelist / EVAPOTRANSPIRATION_PARAMETERS / & ec1, & ecal0, & thc0, & idvwk, & iemeth, & pit, & radiation_switch contains subroutine evapotranspiration_initialise(mb, meap, subbasin_input_file_id) use input, only: get_config_fid use output, only: output_register_hydrotope_var integer, intent(in) :: mb, meap, subbasin_input_file_id read(get_config_fid(), EVAPOTRANSPIRATION_PARAMETERS) if (iemeth.eq.0) call log_debug('evapotranspiration_initialise', & 'Evaporation - Priestley-Taylor, iemeth =', int=iemeth) if (iemeth.eq.1) call log_debug('evapotranspiration_initialise', & 'Evaporation - Turc-Ivanov, iemeth =', int=iemeth) if (idvwk.eq.1) call log_debug('evapotranspiration_initialise', & 'Evaporation variables after DVWK-M 238, idvwk =', int=idvwk) if (idvwk.eq.0) call log_debug('evapotranspiration_initialise', & 'Evaporation variables from legacy method, idvwk =', int=idvwk) eta_output_id = output_register_hydrotope_var("eta", .false.) etp_output_id = output_register_hydrotope_var("etp", .false.) soil_evaporation_output_id = & output_register_hydrotope_var("soil_evaporation", .false.) call evapotranspiration_allocate(mb, meap) call evapotranspiration_read_input(subbasin_input_file_id) yls = sin(lat / 57.296) ylc = cos(lat / 57.296) end subroutine evapotranspiration_initialise subroutine evapotranspiration_read_input(subbasin_input_file_id) use input, only : read_real_column integer, intent(in) :: subbasin_input_file_id call read_real_column(subbasin_input_file_id, "lat", lat) call read_real_column(subbasin_input_file_id, "salb", salb, 0.15_dp) end subroutine evapotranspiration_read_input subroutine evapotranspiration_allocate(mb, meap) integer, intent(in) :: mb, meap allocate(canstor(mb, meap)) allocate(ecal(mb)) allocate(humi(mb)) allocate(lat(mb)) ! subbasins latitude from stat - outdat.csv allocate(ra(mb)) allocate(s1(mb, meap)) allocate(s2(mb, meap)) allocate(salb(mb)) allocate(thc(mb)) allocate(tmn(mb)) allocate(tv(mb, meap)) allocate(ylc(mb)) allocate(yls(mb)) canstor = 0. ecal = ecal0 humi = 0. lat = 0. ra = 0. s1 = 0. s2 = 0. salb = 0.15 thc = thc0 tmn = 0. tv = 0. ylc = 0. yls = 0. end subroutine evapotranspiration_allocate subroutine evapotranspiration_allocate_sc(n_subcatch) integer, intent(in) :: n_subcatch allocate(bsn_ecal(n_subcatch)) bsn_ecal = 0. allocate(bsn_thc(n_subcatch)) bsn_thc = 0. end subroutine evapotranspiration_allocate_sc subroutine dealloc_evapotranspiration deallocate(canstor) deallocate(humi) deallocate(lat) deallocate(ra) deallocate(s1) deallocate(s2) deallocate(salb) deallocate(tmn) deallocate(tv) deallocate(ylc) deallocate(yls) end subroutine dealloc_evapotranspiration subroutine evapotranspiration_process(j, je, k, alai, cva, ep, ida, mo, nn, preinf, qd, snoa, ste, tx, z, bSnowModule, rnew, tmit, tx_tmp, rsn, ETcor) !**** PURPOSE: THIS SUBROUTINE COMPUTES THE AMOUNT OF SOIL EVAPORATION ! & THE POTENTIAL PLANT TRANSPIRATION USING RITCHIE'S METHOD !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! alai(j, je) = leaf area index ! canev = canopy evaporation, mm ! canstor(j, je) = canopy water storage, mm ! cva(j, je) = vegetation cover, kg/h ! ecal = general potential evap calibration factor ! eo = potential evapotranspiration, mm ! eopot = potential evapotranspiration, mm ! ep = plant transpiration, mm ! es = soil evaporation, mm ! et = es + ep, mm ! humi(j) = air humidity in the subbasin, % ! ida = current day ! ievap = switch code to print from evap() ! ievhd = number of hydrotope to print from evap(), if ievap = 1 ! ievsb = number of subbasin to print from evap(), if ievap = 1 ! nn = number of soil layers ! omega = month factor for Turc(Ivanov) evap (Glugla 1989) !block pit = parameter for estimation of the day length ! preinf(j, je) = precipitation adjusted for canopy storage, mm ! qd = daily surface runoff, mm ! ra(j) = solar radiation, J/cm^2 ! s1(j, je) = internal func. for Richie's method to estimate es ! s2(j, je) = internal func. for Richie's method to estimate es ! salb(j) = soil albedo ! snoa(j, je) = water content in snow, mm ! snoev = snow evaporation, mm ! ste(j, je, l) = water storage in a layer, , mm, calc in hydrotop & purk ! thc = correction factor for potential evapotranspiration ! range for thc: (0.8-1.0), value 1. - from R. Muttiah ! tv(j, je) = internal func. for Richie's method to estimate es ! tx(j) = average daily temperatue, degree C ! ylc(j) = cos(lat()/clt), lat() - lat, clt=57.296, to calc rmx ! yls(j) = sin(lat()/clt), lat() - lat, clt=57.296, to calc rmx ! (convert degrees to radians (2pi/360=1/57.296) ) ! z(l, k) = soil depth, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! alb = soil albedo ! aph = local par ! cej = local par ! ch = local par ! corn = correction on land use for Turc-Ivanov evaporation ! d = inclination of saturation vapor pressure curve ! dayl = day length ! dd = local par ! eaj = local par ! eos = local par ! esd = local par ! esx = local par ! ff = local par ! gma = local par ! h = local par ! hh = local par ! ho = local par ! jk1 = local par ! jk2 = local par ! p = local par ! rmx = max solar radiation ! rto = local par ! sb = local par ! sd = local par ! skyem = sky emissivity ! sp = local par ! suu = local par ! thrad = add. radiation due to sky emissivity ! tk = local par ! tkk = local par ! turc = potential evapotranspiration according to Turc ! u = local par ! vp = vapor pressure ! vps = saturated vapor pressure ! xi = current day ! xx = local par ! yc = local par ! ys = local par ! yy = local par ! zz = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ real(dp), dimension(:, :), intent(in) :: alai real(dp), dimension(:, :), intent(in) :: cva real(dp), intent(inout) :: ep integer, intent(in) :: ida integer, intent(in) :: mo integer, intent(in) :: nn real(dp), dimension(:, :), intent(in) :: preinf real(dp), intent(in) :: qd real(dp), dimension(:, :), intent(inout) :: snoa real(dp), dimension(:, :, :), intent(inout) :: ste real(dp), dimension(:), intent(in) :: tx real(dp), dimension(:, :), intent(in) :: z logical, intent(in) :: bSnowModule real(dp), intent(inout) :: rnew real(dp), intent(inout) :: tmit real(dp), intent(inout) :: tx_tmp real(dp), intent(inout) :: rsn(:, :) real(dp), intent(inout) :: ETcor integer j, je, k, jk1, jk2 logical set_suu, update_es_et real(dp) alb, aph, cej, ch, corn, d, dayl, dd, eaj, eos, esd, esx, ff, gma real(dp) h, hh, ho, p, rmx, rto, sb, sd, skyem, sp, suu, thrad, tk, tkk real(dp) u, vp, vps, xi, xx, yc, ys, yy, zz real(dp) turc esd = 300. cej = -5.e-5 u = 6. vps = 0.0 corn = 1.0 !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then tx_tmp = tmit else tx_tmp = tx(j) end if !########################### ! sl begin ! If radiation = 0.0 there will be caused an error while computing evap. ! This was my experience in a catchment in northern Sweden, hence I include here following correction: if (ra(j) <= 0. ) ra(j) = 1.e-6 ! sl end !**** CALC albedo p = preinf(j, je) - qd eaj = exp(cej * (cva(j, je) + .1)) tk = tx_tmp + 273. tkk = tk * tk if (idvwk .eq. 1) then if (tx_tmp .ge. 0.) then vps = 6.11 * exp((17.62 * tx_tmp) / (243.12 + tx_tmp)) d = vps * 4284. / (243.12 + tx_tmp) ** 2. else vps = 6.11 * exp((22.46 * tx_tmp) / (272.62 + tx_tmp)) d = vps * 6123. / (272.62 + tx_tmp) ** 2. end if gma = d / (d + .655) else d = exp(21.255 - 5304. / tk) * 5304. / tkk gma = d / (d + .68) end if if (snoa(j, je) .le. 5.) then alb = salb(j) if (alai(j, je) .gt. 0.) alb = .23 * (1. - eaj) + salb(j) * eaj else alb = .6 !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then if (rsn(j, je) .eq. rnew ) alb = .8 end if !########################### end if !**** CALC max solar radiation rmx, coef 711 changed to 916 ! 916 (ly) ==> 3847.2 (J/cm^2) xi = ida sd = .4102 * sin((xi - 80.25) / pit) ch = - yls(j) * tan(sd) / ylc(j) if (ch .lt. 1.) then if (ch .le. - 1.) then h = 3.1416 else h = acos(ch) end if else h = 0. end if dayl = 7.72 * h dd = 1. + .0335 * sin((xi + 88.2) / pit) ys = yls(j) * sin(sd) yc = ylc(j) * cos(sd) !!! rmx = 711. * (h*ys+yc*sin(h)) rmx = 3847.2 * (h * ys + yc * sin(h)) !**** CORRECTION of radiation on sky emissivity (Ranjan) ! ho1 - old, ho - from Ranjan, divisor 58.3 for cal ==> 244.86 for J ! thrad correction, thc - calib. coef.(0.8-1.) ! ho1 = ra(j) * (1.-alb) / 58.3 ! ho = (ra(j) * (1.-alb)+thrad) / 58.3 if (idvwk .ne. 1) then ! Idso & Jackson (1969) J. Geophys. Res. 74(23):5397--5403 skyem = 1. - 0.261 * exp(-7.77e-4 * (tk-273.15) ** 2) thrad = (skyem-0.96) * 4.914e-7 * tk ** 4 * (.2 + .8 * (ra(j) / rmx)) ho = (ra(j) * (1. - alb) + thc(j) * thrad) / 244.86 else ! DVWK (1995) Merkblatt Nr. 238 / Brunt (1932) if (humi(j) < 0. ) then vp = vps * .65 else vp = vps * (humi(j) / 100.) end if !vp = vps * (humi(j)/100.) skyem = 0.34 - 0.044 * sqrt(vp) thrad = -5.67 * 8.64e-8 * tk ** 4. * skyem * (.1 + .9 * (ra(j) / rmx)) ho = (ra(j) * (1. - alb) + thrad) / (249.8 - .242 * tx_tmp) end if if (ho .lt. 0.) ho = 0.001 !**** CORRECTION of aph=1.28 on humidity: ! aph = f(humidity): aph(60) = 1.74, aph(90)=1.28 if (humi(j) < 0. ) then hh = .65 else hh = humi(j) / 100. end if !hh = humi(j)/100. zz = 40. * hh - 29. if (zz .gt. 10.) zz = 10. if (zz .lt. - 10.) zz = - 10. ff = hh / (hh + exp(zz)) aph = 1.28 + 0.46 * ff !**** CALC POTENTIAL ET if (iemeth .eq. 0) then if (idvwk .eq. 1) then eo = ecal(j) * 1.26 * ho * gma else eo = ecal(j) * aph * ho * gma end if else ! TURC-IVANOV POTENTIAL EVAP if (tx_tmp .ge. 5.) then ! DEFINE MONTH FACTOR FOR TURC-IVANOV EVAP (assigned in common.f90) if (iemeth == 1) omega(mo) = turc_ivanov(mo) if (iemeth == 2) omega(mo) = 1. ! No monthly correction turc = .0031 * omega(mo) * (ra(j) + 209.4) * (tx_tmp / (tx_tmp + 15.)) else if (humi(j) < 0. ) then turc = 0.000036 * (25. + tx_tmp) ** 2. * 35. else turc = 0.000036 * (25. + tx_tmp) ** 2. * (100. - humi(j)) end if !turc = 0.000036 * (25.+tx(j))**2 * (100.-humi(j)) end if ! CORRECTION ON LAND USE eo = ecal(j) * ETcor * turc end if !**** CORRECTION of pot. evap. on interception ! eopot is used later as Pot Evap for write() ! eo is corrected for interception and used to calc es, ep eopot = eo canev = 0. if (canstor(j, je) .gt. 0.) then eo = eo - canstor(j, je) if (eo .lt. 0.) then canstor(j, je) = - eo canev = eopot eo = 0. else canev = canstor(j, je) canstor(j, je) = 0. endif endif !**** CALC soil evaporation es set_suu = .true. eos = eo * eaj if (s1(j, je) .ge. u) then sb = p - s2(j, je) if (sb .ge. 0.) then p = sb s1(j, je) = u - p tv(j, je) = 0. if (s1(j, je) .lt. 0.) s1(j, je) = 0. else tv(j, je) = tv(j, je) + 1. es = 3.5 * sqrt(tv(j, je)) - s2(j, je) if (p .le. 0.) then if (es .gt. eos) es = eos else esx = 0.8 * p if (esx .le. es) esx = es + p if (esx .gt. eos) esx = eos es = esx end if s2(j, je) = s2(j, je) + es - p tv(j, je) = (s2(j, je) / 3.5) ** 2. set_suu = .false. end if else sp = s1(j, je) - p if (sp .gt. 0.) then s1(j, je) = sp else s1(j, je) = 0. end if end if if (set_suu) then s1(j, je) = s1(j, je) + eos suu = s1(j, je) - u if (suu .le. 0.) then es = eos else es = eos - .4 * suu s2(j, je) = .6 * suu tv(j, je) = (s2(j, je) / 3.5) ** 2. end if end if if (es .le. 0.) es = 0. !**** CALC plant transpiration ep: preliminary estimation, ! later recalc in wstress (vegfun.f) if (alai(j, je) .le. 3.0) then ep = alai(j, je) * eo / 3. else ep = eo end if et = es + ep + canev if (eopot .lt. et) then et = eopot es = et - ep - canev end if xx = es !**** RECALC ste(), sno() if (snoa(j, je) .lt. es) then xx = xx - snoa(j, je) snoa(j, je) = 0. update_es_et = .true. do jk2 = 1, nn if (z(jk2, k) .gt. esd) then jk1 = jk2 - 1 yy = 0. if (jk1 .gt. 0) yy = z(jk1, k) rto = ste(j, je, jk2) * (esd - yy) / (z(jk2, k) - yy) if (rto .gt. xx) then ste(j, je, jk2) = ste(j, je, jk2) - xx else xx = xx - rto ste(j, je, jk2) = ste(j, je, jk2) - rto es = es - xx et = et - xx update_es_et = .false. exit end if else if (ste(j, je, jk2) .gt. xx) then ste(j, je, jk2) = - xx + ste(j, je, jk2) update_es_et = .false. exit else xx = xx - ste(j, je, jk2) ste(j, je, jk2) = 0. end if end if end do if (update_es_et) then es = es - xx et = et - xx end if else snoa(j, je) = snoa(j, je) - es snoev = es end if return end subroutine evapotranspiration_process !------------------------------------------------------------------------------- ! Author : aich@pik-potsdam.de ! Date : 2013-08-04 ! modified: 2013-08-04 ! ! PURPOSE : generating radiation with latitude, tmin and tmax ! ! CALLED : from program readcli !------------------------------------------------------------------------------- subroutine evapotranspiration_radiation(ida, mb, tmx) integer, intent(in) :: ida integer, intent(in) :: mb real(dp), dimension(:), intent(in) :: tmx integer :: k real(dp) :: gb, dr, ds, ws2, ram_1, radi, ek2 gb = 0. ! ek1 = empirical constant !ek1=0.16 ! used in West Africa (Niger Basin) !ek1 = 0.115 ! used in East Africa (Blue Nile Basin), old value !ec1 = 0.145 ! used in East Africa (Blue Nile Basin), for new SWIM trunk version ek2 = 0. !empirische Konstante 1 dr = 1 + (.033 * cos(((2. * 3.141592) / 365.) * ida)) !Inverse relative Distanz Erde - Sonne ds = .409 * sin(((2. * 3.141592) / 365.) * ida - 1.39) ! Deklination der Sonne do k = 1, mb ! Note, if latitude is larger 67 degrees North, this function does not work! if (lat(k) > 66.5 ) then lat(k) = 66.5 call log_warn("evapotranspiration_radiation", "Latitude is > 66.5, "// & "in subbasin (N) set to 66.5", i1=k) end if if (lat(k) < - 66.5 ) then lat(k) = - 66.5 call log_warn("evapotranspiration_radiation", "Latitude is < -66.5, "// & "in subbasin (N) set to -66.5", i1=k) end if gb = (3.1415 / 180.) * lat(k) !Breite in Bogenwinkel ws2 = acos((- tan(gb)) * tan(ds)) !Stundenwinkel beim Sonnenuntergang ram_1 = (ws2 * sin(gb) * sin(ds)) + (cos(gb) * cos(ds) * sin(ws2)) radi = ((24. * 60.) / 3.141592) * .0820 * dr * ram_1 radi = (radi * 100.) ! Formel in MJ / m2 / d ra(k) = ec1 * radi * ((tmx(k) - tmn(k)) ** .5) + ek2 ! Formel von Hargreaves(1985) end do end subroutine evapotranspiration_radiation end module evapotranspiration