!     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