evapotranspiration_process Subroutine

  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