!     FILE vegfun.f95
!
!     SUBROUTINES IN THIS FILE CALLED FROM
!     subroutine wstress(j, je, k, n) crpmd, vegmd
!     subroutine tstress(tgx, j, je, n) crpmd, vegmd
!     subroutine scurve(x1, x2, x3, x4) crpmd, adjustbe
!     subroutine ascrv(x1, x2, x3, x4, x5, x6) readcrp
!     subroutine adjustbe(j, je, xtem, beadj) crpmd
module vegetation

  use utilities, only : dp

  implicit none

  ! initial Leaf Area Index for land use type
  real(dp), save, dimension(:), allocatable :: alai0
  integer, parameter :: nih = 7
  ! water use rate - depth parameter, used in wstress
  real(dp) :: ub = 3.065
  ! alpha factor for vegetation
  real(dp), save :: alfa
  ! plant transpiration, mm
  real(dp), save :: ep
  ! temperature stress factor
  real(dp), save :: ts
  ! plant evaporation, potential, mm
  real(dp), save :: potentl
  real(dp), save :: sdt = 0. ! sum temp stress days
  real(dp), save :: sdw = 0. ! sum water stress days
  ! |hour |time threshold used to define dormant , dimension(subtot)
  real(dp), save, dimension(:), allocatable :: dormhr
  ! dormancy threshold for day length
  logical, save :: bDormancy = .false.
  real(dp), save :: epco = 1. ! limit plant water demand compensation from lower layers.
  ! pap calc in readbas from ub, ub=3.065 blockdata
  real(dp), save :: uob
  ! switch parameter: to calc CO2 effect on net photosynthesis?
  integer, save :: ialpha
  ! switch parameter: to calc CO2 effect on transpiration?
  integer, save :: ibeta
  ! switch parameter: 3 / 4 - C3 or C4 crop?
  real(dp), save :: ic3c4
  ! atm CO2 in the reference period, ppm
  real(dp), save :: co2ref
  ! atm CO2 in the scenario period, ppm
  real(dp), save :: co2sce
  real(dp), save, dimension(:), allocatable :: gwqLastday
  ! maximum additional additional uptake rate in wetlands, mm
  real(dp), save :: rzmaxup = 0.
  ! min day length, h, calc. in readwet
  real(dp), save, dimension(:), allocatable :: daylmn
  ! day length in subbasin, h, calc in readcli
  real(dp), save, dimension(:), allocatable :: daylen
  ! index for dormant period
  integer, save, dimension(:, :), allocatable :: idorm
  ! crop residue in two upper soil layers, kg / ha
  real(dp), save, dimension(:, :, :), allocatable :: rsd
  ! total biomass, kg / ha
  real(dp), save, dimension(:, :), allocatable :: dm
  ! fraction heat units accumulated, 1 / 1
  real(dp), save, dimension(:, :), allocatable :: g
  ! alai(j, je) - leaf area index
  real(dp), save, dimension(:, :), allocatable :: olai
  ! leaf area index
  real(dp), save, dimension(:, :), allocatable :: alai
  ! root depth, mm
  real(dp), save, dimension(:, :), allocatable :: rd
  ! water stress factor
  real(dp), save, dimension(:, :), allocatable :: ws
  ! water stress, sum for the growth period
  real(dp), save, dimension(:, :), allocatable :: wsav
  ! temp stress, sum for the growth period
  real(dp), save, dimension(:, :), allocatable :: tsav
  ! base temp. for plant growth, degree C
  real(dp), save, dimension(:), allocatable :: tb
  ! max LAI for crop
  real(dp), save, dimension(:), allocatable :: blai
  ! fraction of growing season, when LAI declines
  real(dp), save, dimension(:), allocatable :: dlai
  ! minimum Leaf Area Index (for forest and natural vegetation)
  real(dp), save, dimension(:), allocatable :: almn
  ! specific leaf area, m2 / kg, LAI / SLA in kg / m2
  real(dp), save, dimension(:), allocatable :: sla
  ! land cover category for vegetation iv
  integer, save, dimension(:), allocatable :: ilcc
  ! shape parameter for the LAI developement equation
  real(dp), save, dimension(:), allocatable :: ald1
  ! shape parameter for the LAI developement equation
  real(dp), save, dimension(:), allocatable :: ald2
  ! optimal temperature for plant growth, degrees C
  real(dp), save, dimension(:), allocatable :: to

  ! output ids
  integer :: &
    transpiration_output_id = 0, &
    vegetation_water_stress_oid = 0, &
    vegetation_temp_stress_oid = 0, &
    heat_unit_fract_output_id = 0, &
    biomass_total_output_id = 0, &
    leaf_area_index_output_id = 0, &
    root_depth_output_id = 0

  namelist / VEGETATION_PARAMETERS / &
    co2ref, &
    co2sce, &
    epco, &
    ialpha, &
    ibeta, &
    ic3c4, &
    rzmaxup, &
    ub, &
    bDormancy


contains

  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

  subroutine vegetation_initialise(mb, nlut, mcrdb, meap, lat)
    use input, only : get_config_fid
    use output, only : output_register_hydrotope_var
    integer, intent(in) :: mb, nlut, mcrdb, meap
    real(dp), intent(in) :: lat(:)

    read(get_config_fid(), VEGETATION_PARAMETERS)

    call vegetation_allocate(mb, nlut, mcrdb, meap)
    call vegetation_prepare_parameters(lat)

    ! output variables
    transpiration_output_id = output_register_hydrotope_var("transpiration", .false.)
    vegetation_water_stress_oid = &
      output_register_hydrotope_var("vegetation_water_stress")
    vegetation_temp_stress_oid = &
      output_register_hydrotope_var("vegetation_temperature_stress")
    heat_unit_fract_output_id = output_register_hydrotope_var("heat_unit_fraction")
    biomass_total_output_id = output_register_hydrotope_var("biomass_total")
    leaf_area_index_output_id = output_register_hydrotope_var("leaf_area_index")
    root_depth_output_id = output_register_hydrotope_var("root_depth")

  end subroutine vegetation_initialise

  subroutine vegetation_prepare_parameters(lat)
    real(dp), intent(in) :: lat(:)
    real(dp) :: ch(size(lat)), hrs(size(lat))

    ch = .439 * abs(tan(lat / 57.296))
    hrs = merge(0., acos(ch), ch >= 1.)
    daylmn = 7.72 * hrs
    if (bDormancy) then
      dormhr = merge(1., (abs(lat) - 20.) / 20., abs(lat) > 40.)
      dormhr = merge(-0.1, dormhr, abs(lat) < 20.)
    else
      dormhr = - 1.
    end if
  end subroutine vegetation_prepare_parameters

  subroutine vegetation_allocate(mb, nlut, mcrdb, meap)
    integer, intent(in) :: mb, nlut, mcrdb, meap

    allocate(alai(mb, meap))
    allocate(alai0(nlut))
    allocate(ald1(mcrdb))
    allocate(ald2(mcrdb))
    allocate(almn(mcrdb))
    allocate(blai(mcrdb))
    allocate(daylen(mb))
    allocate(daylmn(mb))
    allocate(dlai(mcrdb))
    allocate(dm(mb, meap))
    allocate(dormhr(mb))
    allocate(g(mb, meap))
    allocate(gwqLastday(mb))
    allocate(idorm(mb, meap))
    allocate(ilcc(mcrdb))
    allocate(olai(mb, meap))
    allocate(rd(mb, meap))
    allocate(rsd(mb, meap, 2))
    allocate(sla(mcrdb))
    allocate(tb(mcrdb))
    allocate(tsav(mb, meap))
    allocate(ws(mb, meap))
    allocate(wsav(mb, meap))
    allocate(to(mcrdb))

    alai(:, :) = 0.
    alai0 = 0.
    ald1 = 0.
    ald2 = 0.
    almn = 0.
    blai = 0.
    daylen = 0.
    daylmn = 0.
    dlai = 0.
    dm = 0.
    dormhr = 0.
    g = 0.
    gwqLastday = 0.
    idorm = 0
    ilcc = 0
    olai(: , :) = 0.
    rd = 0.
    rsd = 0.
    sla = 0.
    tb = 0.
    tsav = 0.
    uob = 1. - exp(- ub)
    ws = 0.
    wsav = 0.
    to = 0.

  end subroutine vegetation_allocate


  subroutine dealloc_vegetation
    deallocate(alai)
    deallocate(alai0)
    deallocate(ald1)
    deallocate(ald2)
    deallocate(almn)
    deallocate(blai)
    deallocate(daylen)
    deallocate(daylmn)
    deallocate(dlai)
    deallocate(dm)
    deallocate(g)
    deallocate(gwqLastday)
    deallocate(idorm)
    deallocate(ilcc)
    deallocate(olai)
    deallocate(rd)
    deallocate(rsd)
    deallocate(sla)
    deallocate(tb)
    deallocate(tsav)
    deallocate(ws)
    deallocate(wsav)
    deallocate(to)

  end subroutine dealloc_vegetation

  subroutine vegetation_water_stress(j, je, k, n, wet, additionalGwUptake, bWAM_Module, dart, daycounter, fc, frar, humi, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, rdmx, sbar, sep, ste, z)
    !**** PURPOSE: THIS SUBROUTINE ESTIMATES WATER STRESS FACTOR.
    !              IT DISTRIBUTES POTENTIAL PLANT EVAPORATION THROUGH
    !              THE ROOT ZONE AND CALCULATES ACTUAL PLANT WATER USE
    !              BASED ON SOIL WATER AVAILABILITY
    !**** CALLED IN: CRPMD, VEGMD
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      alfa = calculated in adjustbe
    !      co2ref = CO2 in reference period, ppm
    !      co2sce = CO2 in simulated period, ppm
    !      ep = plant evaporation, mm
    !      g(j, je) = fraction heat units accumulated, 1/1
    !      humi(j) = air humidity, %
    !      ialpha = switch parameter: to calc CO2 effect on net photosynthesis?
    !      ibeta = switch parameter: to calc CO2 effect on transpiration?
    !      ida = current day
    !      ilcc(iv) = land cover category for vegetation iv
    !      iwshd = number of hydrotope to print from wstress(), if iwstr = 1
    !      iwssb = number of subbasin to print from wstress(), if iwstr = 1
    !      iwstr = switch code to print from wstress()
    !      nn = number of soil layers
    !      nucr(j, je) = crop number (database)
    !      nveg(j, je) = vegetation number (database)
    !      potentl = plant evaporation, potential, mm
    !      precip = precipitation, mm
    !      rd(j, je) = root depth, mm
    !      rdmx(iv) = max root depth, mm
    !      sep = water percolation
    !      ste(j, je, l) = water storage, recalc here, mm
    !block ub = 3.065
    !      ul(l, k) = upper limit water content in layer, mm
    !      uob = pap calc in readbas from ub, ub=3.065 blockdata
    !      ws(j, je) = water stress factor
    !      z(l, k) = depth of layer, mm
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      beta = beta factor for plant transpiration under higher CO2
    !      ccfi = local par
    !      epin = local par
    !      gx = local par
    !      ir = local par
    !      l = local par
    !      num = local par
    !      sum = local par
    !      u = local par |mm H2O |water uptake by plants in each soil layer
    !      ul4 = local par
    !      xx = local par |mm H2O |water uptake by plants from all layers
    !      xx1 = local par
    !      xx2 = local par
    !      xx3 = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters
    use management, only : &
        TSubbasin, TWaterUser, management_is_transfer_subbasin, management_subbasin_pointer, &
        management_user_pointer, management_is_active_period, wam_correct_irrigationdemand, &
        wam_plantdemand_summary, wam_plantdemand_summary, wam_irrigationdemand_summary
    use landuse, only : landuse_is_cropland

    real(dp), dimension(:), intent(inout) :: additionalGwUptake
    logical, intent(in) :: bWAM_Module
    real(dp), dimension(:), intent(in) :: dart
    integer, intent(in) :: daycounter
    real(dp), dimension(:, :), intent(in) :: fc
    real(dp), dimension(:, :), intent(in) :: frar
    real(dp), dimension(:), intent(in) :: humi
    integer, intent(in) :: icc
    integer, intent(in) :: ida
    integer, intent(in) :: iy
    integer, intent(in) :: iyr
    integer, dimension(:, :, :), intent(in) :: mstruc
    integer, intent(in) :: nbyr
    integer, intent(in) :: nn
    integer, dimension(:, :), intent(in) :: nucr
    integer, dimension(:, :), intent(in) :: nveg
    real(dp), dimension(:), intent(in) :: rdmx
    real(dp), dimension(:), intent(in) :: sbar
    real(dp), intent(inout) :: sep
    real(dp), dimension(:, :, :), intent(inout) :: ste
    real(dp), dimension(:, :), intent(in) :: z
    integer j, je, n, k, ir, l, num, wet
    real(dp) beta, ccfi, epin, gx, sum, ul4, xx, xx1, xx2, xx3
    real(dp) areahyd, additionalUptake, co2sce_interpolated
    real(dp), dimension(10) :: u
    ! epco as a global calibration parameter in * .bsn
    real(dp) zu, sump
    real(dp) rdmx_temp
    !#### WATER MANAGEMENT MODULE ####
    real(dp) :: area ! hydrotop area, m2
    !#### WATER MANAGEMENT MODULE ####
    TYPE (TSubbasin), POINTER :: pS ! pointer on subbasins
    !#### WATER MANAGEMENT MODULE ####
    TYPE (TWaterUser), POINTER :: pWU ! pointer on actual TWU

    u = 0.
    potentl = ep

    ! The factor "epco" is introduced to limit plant water demand
    ! compensation from lower layers.
    ! Without epco (or epco=1) water demand which could not be met
    ! by an individual layer was moved to the next lower layer.
    ! epco=0 means no compensation. epco=0.5 means 50% compensation
    ! from next lower layer.
    ! It is basically a calibration parameter and reflects
    ! the plants flexibility to respond to water limitations
    ! with increased root growth. -- Pia Gottschalk
    !**** SL
          !epco=0.25
          ! epco is now a calibration paramter in *.bsn
    !**** SL

    !**** Including factor beta for plant transpiration under higher CO2
    !               as a constant beta coupled to alfa
    if (ialpha .eq. 1 .and. ibeta .eq. 1) then
      ! Calculation of transient CO2 increase = simple interpolation
      ! from start year CO2 concentration (co2ref)
      ! to end year CO2 concentration (co2sce):
      co2sce_interpolated = ((co2sce - co2ref) / (nbyr - 1)) * (iy - 1) + co2ref
      !        beta = alfa * co2sce / co2ref
      ! Correction of calculation of beta, Pia Gottschalk, 23/08/2011
      beta = alfa * co2ref / co2sce_interpolated
      ep = ep * beta
    end if

    epin = ep
    xx = 0.

    if (landuse_is_cropland(n) ) then
      num = nucr(j, je)
    else
      num = nveg(j, je)
    end if
    !**** CALC plant transpiration ep and water stress ws()
    !     Correction from F. Badeck: ul4 = ul()*0.6, 60% instead of 25% before
    !     Correction from P. Gottschalk, 19/07/11:
    !     ul4 = ul()*0.5, 50% instead of 25% (SWIM original) before
    !     Correction from P. Gottschalk, 19/07/11, is based on
    !     "Leitfaden zur Beregnung landwirtschaftlicher Kulturen",
    !     LVLF, September 2005, Tabelle 3
    if (ep .le. 0.) then
      ws(j, je) = 1.
    else
      l = 1
      ir = 0
      sum = 0.
      sump = 0.
      !       select case (ilcc(num))
      !          case (1, 2, 4, 5)
      !             rd(j, je) = 2.5 * g(j, je) * rdmx(num)
      !             if ( rd(j, je) > rdmx(num) ) rd(j, je)= rdmx(num)
      !             !if (sol_rd < 10.) sol_rd = 10.
      !          case default
      !             rd(j, je) = rdmx(num)
      !       end select
      ! Pia: rdmx (max. rooting depth) is adjusted to the depth
      ! of the arable layer or max rooting depth (from crop.dat)
      ! if the depth of the arable layer is greater than the
      ! crop specific max root depth.
      ! To adjust hydrotop specific water uptake the variable
      ! "rdmx_temp" is introduced here locally.
      ! This prevents that the loop over the arable layer to calculate
      ! water uptake generates a water stress artefact.
      ! Before, water demand was calculated along the whole root profile
      ! but the loop only ran over the arable layer.
      ! Note, this effectively reduces plant water stress
      ! because max. root depth for winter wheat is set to 2 m
      ! whereas the arable layer is usually set to a shallower layer.
      rdmx_temp = min(rdmx(num), z(nn, k))
      if (ilcc(num) .le. 2 .or. ilcc(num) .eq. 4 .or. ilcc(num) .eq. 5 &
          .OR. nucr(j, je) == icc ) then ! why not include cover crop?
        rd(j, je) = 2.5 * g(j, je) * rdmx_temp
        if (rd(j, je) .gt. rdmx_temp) then
          rd(j, je) = rdmx_temp
        end if
      else
        rd(j, je) = rdmx_temp
      end if

      xx = 0.
      zu = 0.
      do l = 1, nn
        if (ir .gt. 0) exit
        if (rd(j, je) .le. z(l, k)) then
          gx = rd(j, je)
          ir = l
        else
          gx = z(l, k)
        end if
        if (rd(j, je) .le. 0.) then
          sum = ep / uob
        else
          if (num .eq. 69) then
            sum = ep * (1. - exp(- ub)) / uob
          else
            sum = ep * (1. - exp(- ub * gx / rd(j, je))) / uob
          end if
        end if

        sump = sum
        !          u(l) = sum - xx
          u(l) = sum - sump + (sump - xx) * epco
        ! Pia: Commented out as the threshold from which plant water
        ! uptake is reduced due to increasing matrix potential (here 0.5)
        ! refers to the plant available water (=fc) and not to saturation (=ul)
        ! Note, this reduces plant water stress.
        !         ul4 = ul(l, k) * 0.5
        ul4 = fc(l, k) * 0.75

        ! Pia: Commented out: linear decline of plant water uptake below fc*0.5
        !            u(l) = u(l) * ste(j, je, l) / ul4
        ! Pia: Now using SWAT approach which features a exponential decline
        ! of plant water uptake below fc*0.5
        if (ste(j, je, l) .lt. ul4 ) u(l) = u(l) * exp(5 * (ste(j, je, l) / ul4 - 1))

        if (ste(j, je, l) .lt. u(l) ) u(l) = ste(j, je, l)
        ste(j, je, l) = ste(j, je, l) - u(l)
        xx = xx + u(l)
        zu = z(l, k)
      end do ! l = 1, nn

      ws(j, je) = xx / ep

      !#################################
      !#### WATER MANAGEMENT MODULE ####
      !#################################
      ! calculate plant water demand for irrigated crops
      if (bWAM_Module) then
        if (management_is_transfer_subbasin(j) ) then
          if (landuse_is_cropland(n) .AND. mstruc(j, je, 7) > 0 & ! if HRU is cropland .AND. if value in column IRR in * .str file is > 0
                              .AND. nucr(j, je) /= icc ) then ! if hydrotope is irrigated cropland, no cover crop (icc)
            pS => management_subbasin_pointer(j) ! pointer on current subbasin (could be source or destination)
            pWU => management_user_pointer(pS % pos_irr) ! pointer on current water user
            area = frar(j, je) * sbar(j) ! HRU area in m^2
            !-----------------------------------------------------------------
            ! Plant and total irrigation requirements are calculated here.
            ! If irr_opt == 1, these requirements are written to the output files, but
            ! the amount of water used for irrigation will be overwritten by input time series in: wam_withdraw_Transfer_Out() !!!
            !-----------------------------------------------------------------
            if (pWU % wu_opt == 4 .AND. management_is_active_period(iyr, ida, pWU) ) then ! if source is not external but another subbasin

              ! calculate plant demand but only for output purposes because under this option,
              ! water is provided (input time series) independently from requirements.
              ! xx = water uptake by plants from all layers [mm]
              ! ep = plant transpiration [mm]
              pWU % plantDemand(daycounter) = real(pWU % plantDemand(daycounter) + (xx - ep) / 1000. * area / 86400. * (- 1.)) ! mm - > m3 / s

              ! calculate actual irrigation demand including losses depending on irrigation practices
              pWU % irrigationDemand(daycounter) = wam_correct_irrigationDemand(pWU, pWU % plantDemand(daycounter))

              ! summarise total plant demand and irrigation demand for summary output file
              wam_plantDemand_summary(daycounter) = &
                  wam_plantDemand_summary(daycounter) + pWU % plantDemand(daycounter)
              wam_irrigationDemand_summary(daycounter) = &
                  wam_irrigationDemand_summary(daycounter) + pWU % irrigationDemand(daycounter)
              end if ! ( pWU % wu_opt == 4 )
          end if ! ( landuse_is_cropland(n) .AND. mstrucx(j, je, 1) >= 1 )
        end if ! ( wam_is_transfer_subbasin(j) )
      end if ! ( wam_bTransfer )
        !#################################

      ep = xx
    end if ! if (ep .le. 0.)

    !**** Correction on air humidity
    if (ida .gt. 100 .and. humi(j) .gt. 75.) then
        xx2 = potentl - ep
      ! ccfi = (humi(j) - 75.)/10.
      ! ccfi = 2.*(humi(j) - 75.)/10.
      ccfi = 1.
      if (ccfi .gt. 1.) ccfi = 1.
      if (xx2 .gt. 0.) then
        ! xx3=min(sep, ccfi*xx2)
        xx3 = ccfi * xx2
        if (xx3 .gt. sep) xx3 = sep
        ep = ep + xx3
        sep = sep - xx3
        if (sep .lt. 0.) sep = 0.
        ws(j, je) = ep / potentl
      end if
    end if

    !       new Riperian zone
    if (wet == 1 .and. gwqLastday(j) > 0) then

      ! calculate transpiration deficite xx1
      xx1 = potentl - ep
      if (xx1 .gt. 0.) then
        ! uptake is limitited to rzmaxup
        if (xx1 > rzmaxup) xx1 = rzmaxup
        ep = ep + xx1
        ! calc hydrotop area
        areahyd = frar(j, je) * dart(j)
        ! calculate additional uptake in m3
        additionalUptake = xx1 * areahyd * 1000

        ! to avoid negative flow
        if ((gwqLastday(j) - additionalUptake) < 0) then
          ep = ep - xx1
          additionalUptake = gwqLastday(j)
          ! increase Et by additional uptake
          ep = ep + additionalUptake / (areahyd * 1000)
        end if

        gwqLastday(j) = gwqLastday(j) - additionalUptake
        additionalGwUptake(j) = additionalGwUptake(j) + additionalUptake
        ws(j, je) = ep / potentl
      end if
    end if

    return
  end subroutine vegetation_water_stress


  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


  subroutine vegetation_temperature_stress(tgx, j, je, n, avt, nucr, nveg, tmn, to, tx)
    !**** PURPOSE: COMPUTES TEMPERATURE STRESS FOR CROP/VEGETATION GROWTH
    !**** CALLED IN: GROWTH, VEGMD
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      avt(j) = average annual daily temerature., degree C
    !      nucr(j, je) = crop number
    !      nveg(j, je) = vegetation number
    !      tb(iv) = base temp for plant growth, degree C
    !      tmn(j) = daily min temp, from readcli, degree C
    !      to(iv) = opt temp for growth, degree C
    !      ts = temp stress factor
    !      tx(j) = daily aver temp, degree C
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      num = local par
    !      rto = local par
    !      rto1 = local par
    !      rto2 = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters
    use landuse, only : landuse_is_cropland
    use snow, only : bSnowModule, tmin, tx_tmp, tmit, tmin_tmp

    real(dp), dimension(:), intent(in) :: avt
    integer, dimension(:, :), intent(in) :: nucr
    integer, dimension(:, :), intent(in) :: nveg
    real(dp), dimension(:), intent(in) :: tmn
    real(dp), dimension(:), intent(in) :: to
    real(dp), dimension(:), intent(in) :: tx
    integer j, je, n, num
    real(dp) tgx, rto, rto1, rto2
    !###########################
    !#### SNOW MODULE ####
    !###########################
    if (bSnowModule) then
      tx_tmp = tmit
      tmin_tmp = tmin
    else
      tx_tmp = tx(j)
      tmin_tmp = tmn(j)
    end if
    !###########################

    if (landuse_is_cropland(n) ) then
        num = nucr(j, je)
    else
        num = nveg(j, je)
    endif

    rto = 0.
    rto1 = 0.
    rto2 = 0.

    !**** CALC temp stress factor ts
    if (tx_tmp .le. to(num)) then
      tgx = (to(num) + tb(num)) / (to(num) - tb(num))
      rto1 = tgx * (to(num) - tx_tmp) ** 2
      rto2 = rto1 / (tx_tmp + 1.e-6)
      if (rto2 .le. 200.) then
        ts = exp(- 0.1054 * rto2)
      else
        ts = 0.
      end if
    else
      tgx = 2. * to(num) - tb(num) - tx_tmp
      rto = ((to(num)-tx_tmp) / (tgx + 1.e-6)) ** 2
      if (rto .le. 200.) then
        ts = exp(- 0.1054 * rto)
      else
        ts = 0.
      end if
    endif

    if (tmin_tmp .le. (avt(j) - 15.)) ts = 0.

    return
  end subroutine vegetation_temperature_stress


  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



  subroutine vegetation_s_curve(x1, x2, x3, x4)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES S CURVE PARAMETERS x1 & x2
    !              GIVEN 2 (X, Y) POINTS
    !**** CALLED IN: CRPMD
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    real(dp) x1, x2, x3, x4, xx

    xx =log(x3 / x1 - x3)
    x2 = (xx -log(x4 / x2 - x4)) / (x4 - x3)
    x1 = xx + x3 * x2
    return
  end subroutine vegetation_s_curve



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


  subroutine vegetation_s_curve_parameters(x1, x2, x3, x4, x5, x6)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES S CURVE PARAMETERS x5 & x6
    !              GIVEN 2 (X, Y) POINTS
    !**** CALLED IN: READCRP
    !~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


    real(dp) xx, x1, x2, x3, x4, x5, x6

    xx = 0.0
    x5 = 0.0
    x6 = 0.0
    xx = Log(x3 / x1 - x3)
    x6 = (xx - Log(x4 / x2 - x4)) / (x4 - x3)
    x5 = xx + (x3 * x6)

    return
  end subroutine vegetation_s_curve_parameters


  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


  subroutine vegetation_adjust_energy_ratio(j, je, beadj, be, nucr, tx)
    !**** PURPOSE: this subroutine adjusts biomass-energy ratio for crop
    !              for higher atmosphere CO2 concentrations
    !              F. Wechsung method (temperature-dependent)
    !**** CALLED IN: CRPMD
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      alfa = alpha factor for vegetation
    !      be(icr) = biomass-energy ratio for crop
    !      co2ref = atm CO2 in the reference period, ppm
    !      co2sce = atm CO2 in the scenario period, ppm
    !      ic3c4 = switch parameter: 3/4 - C3 or C4 crop?
    !      nucr(j, je) = crop number (database)
    !      tx(j) = average daily temp., degree C
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      aa = coef
    !      bb = coef
    !      c1 = co2ref
    !      c2 = co2sce
    !      cc = coef
    !      ci1 = CO2 internal corr co2ref
    !      ci2 = CO2 internal corr co2sce
    !      deg = local par
    !      del = local par
    !      delq = local par
    !      xlnalfa = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters
    use snow, only : bSnowModule, tx_tmp, tmit

    real(dp), dimension(:), intent(in) :: be
    integer, dimension(:, :), intent(in) :: nucr
    real(dp), dimension(:), intent(in) :: tx
    integer j, je
    real(dp) beadj, aa, bb, c1, c2, cc, ci1, ci2, deg, del, delq, xlnalfa
    !###########################
    !#### SNOW MODULE ####
    !###########################
    if (bSnowModule) then
      tx_tmp = tmit
    else
      tx_tmp = tx(j)
    end if
    !###########################
    !**** CALC xlnalfa = f(tx)
    aa = 0.3898E-2
    bb = 0.3769E-5
    cc = 0.3697E-4
    c1 = co2ref
    c2 = co2sce
    ci1 = c1 * 0.7
    ci2 = c2 * 0.7
    del = ci2 - ci1
    delq = ci2 * ci2 - ci1 * ci1
    xlnalfa = aa * del - bb * delq + cc * del * tx_tmp

    !**** CALC degree - different for C3 and C4 crops
    if (ic3c4 .eq. 4) then
      deg = 0.36
    else
      deg = 0.6
    endif

    !**** CALC alfa and beadj
    alfa = exp(xlnalfa) ** deg
    if (alfa .lt. 1) alfa = 1.
    beadj = be(nucr(j, je)) * alfa

    return
  end subroutine vegetation_adjust_energy_ratio



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



  subroutine vegetation_process(j, je, k, n, wet, additionalGwUptake, avt, bWAM_Module, be, cva, dart, daycounter, fc, flu, frar, huharv, humi, hun, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, ra, rdmx, rwt, sbar, sep, snup, spup, ste, sumfc, swe, tmn, to, tx, uap, z, alai, olai)
    !**** PURPOSE: CALC daily growth of plant biomass for natural vegetation
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      alai(j, je) = leaf area index
    !      ald1(iv) = shape parameter for the LAI developement equation
    !                       for veget iv
    !      ald2(iv) = shape parameter for the LAI developement equation
    !                       for veget iv
    !      almn(iv) = LAI minimum (for forest and natural vegetation)
    !      be(iv) = biomass-energy ratio for crop, kg m2 MJ-1 ha-1 d-1
    !      blai(iv) = max LAI for crop
    !      cva(j, je) = vegetation cover, kg/ha
    !      daylen(j) = day length in subbasin, h, calc in readcli
    !      daylmn(j) = min day length, h, calc. in readwet
    !      dlai(iv) = fraction of season, when LAI declines
    !      dm(j, je) = total biomass, kg/ha
    !      flu(j) = fraction of subbasin area in the basin
    !      g(j, je) = fraction of heat units to maturity accumulated
    !      huharv(j, je) = harvest index heat unit
    !      hun(iv) = potential heat units required for maturity of crop
    !      ida = current day
    !      idorm(j, je) = index for dormant period
    !      ih1, ih2, ih3 = hydrotopes for HYDROTOPE PRINTOUT
    !      isb1, isb2, isb3 = subbasins for HYDROTOPE PRINTOUT
    !      iy = current year as counter (1, ..., nbyr)
    !      iyr = current year
    !      nveg(j, je) = number of vegetation from crop database
    !      olai(j, je) = alai(j, je) - leaf area index
    !      ra(j) = solar radiation in subbasin j, J/cm^2
    !      rd(j, je) = root depth, mm
    !      rdmx(iv) = maximum root depth, mm
    !      rsd(j, je, 2) = residue, kg/ha
    !      rwt(j, je) = fraction of root weight
    !      sdt = sum temp stress days
    !      sdw = sum water stress days
    !      sla(iv) = specific leaf area, m2/kg, LAI/SLA in kg/m2
    !      tb(iv) = base temperature for plant growth, degrees C
    !      ts = temp. stress
    !      tsav(j, je) = temp. stress, accumulated
    !      tx(j) = daily aver temp in the subbasin, degrees C
    !      uap = P uptake, kg/ha
    !      ws(j, je) = water stress
    !      wsav(j, je) = water stress, accumulated
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      ddm = delta dm
    !      delg = delta g
    !      deltalai = delta LAI
    !      f = fraction of plant's maximum leaf area index corresponding
    !                   to a given fraction of potential heat units for plant
    !      ff = delta f for the day
    !      par = local par
    !      reg = local par
    !      resnew = new residue
    !      tgx = local par
    !      xx = local par
    !      yy = local par
    !      zz = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters
    use snow, only: bSnowModule, tmit, tx_tmp
    use nutrient, only : nutrient_phosphorus_uptake, nutrient_nitrate_uptake, nutrient_on

    real(dp), dimension(:, :), intent(inout) :: alai
    real(dp), dimension(:, :), intent(inout) :: olai
    real(dp), dimension(:), intent(inout) :: additionalGwUptake
    real(dp), dimension(:), intent(in) :: avt
    logical, intent(in) :: bWAM_Module
    real(dp), dimension(:), intent(in) :: be
    real(dp), dimension(:, :), intent(inout) :: cva
    real(dp), dimension(:), intent(in) :: dart
    integer, intent(in) :: daycounter
    real(dp), dimension(:, :), intent(in) :: fc
    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:, :), intent(in) :: frar
    real(dp), dimension(:, :), intent(inout) :: huharv
    real(dp), dimension(:), intent(in) :: humi
    integer, dimension(:), intent(in) :: hun
    integer, intent(in) :: icc
    integer, intent(in) :: ida
    integer, intent(in) :: iy
    integer, intent(in) :: iyr
    integer, dimension(:, :, :), intent(in) :: mstruc
    integer, intent(in) :: nbyr
    integer, intent(in) :: nn
    integer, dimension(:, :), intent(in) :: nucr
    integer, dimension(:, :), intent(in) :: nveg
    real(dp), dimension(:), intent(in) :: ra
    real(dp), dimension(:), intent(in) :: rdmx
    real(dp), dimension(:, :), intent(inout) :: rwt
    real(dp), dimension(:), intent(in) :: sbar
    real(dp), intent(inout) :: sep
    real(dp), dimension(:, :), intent(inout) :: snup
    real(dp), dimension(:, :), intent(inout) :: spup
    real(dp), dimension(:, :, :), intent(inout) :: ste
    real(dp), dimension(:), intent(in) :: sumfc
    real(dp), dimension(:, :), intent(in) :: swe
    real(dp), dimension(:), intent(in) :: tmn
    real(dp), dimension(:), intent(in) :: to
    real(dp), dimension(:), intent(in) :: tx
    real(dp), intent(out) :: uap
    real(dp), dimension(:, :), intent(in) :: z
    integer j, je, k, n, wet
    real(dp) ddm, delg, deltalai, f, ff, par, reg, resnew, tgx, xx, yy, zz
    integer :: veg_end
    real(dp) :: veg_wstress, veg_dieoff, dieoff

    ! ToDo: veg_end (end of vegetation period to restart growth)
    ! this parameter is to be defined somewhere as input parameter, not fixed like here!
    veg_end = 365
    veg_wstress = .1
    veg_dieoff = .99
    dieoff = 0.

    !###########################
    !#### SNOW MODULE ####
    !###########################
    if (bSnowModule) then
      tx_tmp = tmit
    else
      tx_tmp = tx(j)
    end if
    !###########################

    uap = 0.
    ts = 0.

    !**** CALC land cover for natural vegetation:
    cva(j, je) = 1000. * alai(j, je) + rsd(j, je, 1)

    !**** CHECK if start of dormant period
    !     CALC rsd(), dm(), alai(), g()
    !     Residue allocation formula is changed: to check!
    !      if (idorm(j, je).eq.0) then
    !      if (daylen(j)-1..lt.daylmn(j)) then
    !      if (iy.ne.1.or.ida.ge.180) then
    if (bDormancy) then
      if (idorm(j, je) == 0 ) then
        if (daylen(j) + dormhr(j) < daylmn(j) ) then
          if (iy > 1.OR.ida >= veg_end) then
            idorm(j, je) = 1
            resnew = (olai(j, je) - alai(j, je)) * 10000. / sla(nveg(j, je)) ! sla = specific leaf area, m2 / kg, LAI / SLA in kg / m2
            rsd(j, je, 1) = rsd(j, je, 1) + resnew * 0.5
            rsd(j, je, 2) = rsd(j, je, 2) + resnew * 0.5
            dm(j, je) = dm(j, je) - resnew
            alai(j, je) = almn(nveg(j, je))
            g(j, je) = 0.
          else
            idorm(j, je) = 1
            alai(j, je) = almn(nveg(j, je))
            g(j, je) = 0.
          end if
        end if
      end if
    else ! ( bDormancy )
      if (idorm(j, je) .eq. 0) then
        if (daylen(j) - 1. .lt. daylmn(j)) then
          if (iy .ne. 1 .or. ida .ge. 180) then
            idorm(j, je) = 1
            resnew = (olai(j, je) - alai(j, je)) * 10000. / sla(nveg(j, je)) ! sla = specific leaf area, m2 / kg, LAI / SLA in kg / m2
            rsd(j, je, 1) = rsd(j, je, 1) + resnew * 0.5
            rsd(j, je, 2) = rsd(j, je, 2) + resnew * 0.5
            dm(j, je) = dm(j, je) - resnew
            alai(j, je) = almn(nveg(j, je))
            g(j, je) = 0.
          else
            idorm(j, je) = 1
            alai(j, je) = almn(nveg(j, je))
            g(j, je) = 0.
          end if
        end if
      end if
    end if ! ( bDormancy )


    !**** check if end of dormant period
    if (bDormancy) then
      if ((idorm(j, je) >= 1).AND.(daylen(j) + dormhr(j) >= daylmn(j)).AND.ida < veg_end ) then
        idorm(j, je) = 0
      end if
    else
      if (idorm(j, je) .ge. 1) then
        if (daylen(j) - 1. .ge. daylmn(j)) then
        idorm(j, je) = 0
        end if
      end if
    end if ! ( bDormancy )

    !**** Assuming that the root depth of natural vegetation is not allowed to decrease
    rd(j, je) = rdmx(nveg(j, je))

    !#### CALL WSTRESS
    call vegetation_water_stress(j, je, k, n, wet, additionalGwUptake, bWAM_Module, dart, daycounter, fc, frar, humi, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, rdmx, sbar, sep, ste, z)


    !sl begin
    !###########################################################
    !**** This code snippet has been included to better account for vegetation growth
    !     controlled by water stress rather than temperature or daylength.
    !     Vegetation growth is limited to water availability using the soil water index.
    !     This part is only active if parameter nat_veg is 1 (*.bsn)
    if (bDormancy) then
      ! swe(j, jea)/sumfc(k) = soil water index (0-1)
      ! swe = soil water content [mm]
      ! sumfc = sum of field capacity in soil [mm]
      if (swe(j, je) / sumfc(k) <= veg_wstress ) then
        idorm(j, je) = 1 ! vegetation stops growing (is dormant due to water stress)
        olai(j, je) = alai(j, je) ! olai = alai of previous day
        alai(j, je) = alai(j, je) * veg_dieoff ! assuming dying vegetation up to a certain amount per time step
        alai(j, je) = max(alai(j, je), almn(nveg(j, je)))

        ! calcuate residue from difference between acutal lai and lai of previous day
        resnew = (olai(j, je) - alai(j, je)) * 10000. / sla(nveg(j, je)) ! sla = specific leaf area, m2 / kg, LAI / SLA in kg / m2
        resnew = max(0.0_dp, resnew)

        rsd(j, je, 1) = rsd(j, je, 1) + resnew * .5 ! half of this residue is contributing to residue of upper storage
        rsd(j, je, 2) = rsd(j, je, 2) + resnew * .5 ! the other half contributes to second storage

        dm(j, je) = dm(j, je) - resnew ! reduce biomass by value of residue
        dm = max(0.0_dp, dm)
      else
        idorm(j, je) = 0
      end if

      ! fraction of heat units and alai need to be initialized after growing season
      ! ToDo: the definition of the growing season needs re-thinking!
      if (idorm(j, je) == 0 ) then
        if (g(j, je) >= 1. ) then
          g(j, je) = 0.
          !alai(j, je) = almn(nveg(j, je))
        end if
      end if

      if (ida == veg_end) then
        g(j, je) = 0.
      end if
    end if ! if (bDormancy)


    !###########################################################
    !sl end

    !       if (idorm(j, je).ge.1) goto 10
    if (idorm(j, je) == 0) then
      !**** COMPUTE DAILY INCREASE IN HEAT UNITS delg
      delg = (tx_tmp - tb(nveg(j, je))) / hun(nveg(j, je))
      if (delg .lt. 0.) delg = 0.
      g(j, je) = g(j, je) + delg
      if (g(j, je) .gt. 1.) g(j, je) = 1.

      !*********************************************************** START IF (G<=1)
      !**** GROWTH SEASON
      if (g(j, je) .le. 1) then
        !#### CALL TSTRESS: COMPUTE TEMPERATURE STRESS
        tgx = tx_tmp - tb(nveg(j, je))
        if (tgx .le. 0.) then
          ts = 0.
        else
          call vegetation_temperature_stress(tgx, j, je, n, avt, nucr, nveg, tmn, to, tx)
        end if

        !**** CALC daily biomass increase: ddm
        par = .005 * ra(j) * (1. - exp(- .65 * (alai(j, je) + .05)))
        ddm = be(nveg(j, je)) * par
        if (ddm .lt. 0.) ddm = 0.

        !#### CALL NUPTAKE, PUPTAKE: CALCULATE N AND P UPTAKE
        if (nutrient_on) then
            call nutrient_nitrate_uptake(j, je, nveg(j, je), dm, flu, frar, g, ida, nn, snup)
            call nutrient_phosphorus_uptake(j, je, nveg(j, je), dm, flu, frar, g, ida, nn, spup)
        end if

        !**** CALC BIOMASS dm(), ROOT WEIGHT rwt(9
        xx = dm(j, je) + ddm
        reg = amin1(real(ws(j, je), 4), real(ts, 4))
        if (reg .lt. 0.) reg = 0.
        if (reg .gt. 1.) reg = 1.
        dm(j, je) = dm(j, je) + ddm * reg
        rwt(j, je) = (.4 - .2 * g(j, je))

        tsav(j, je) = tsav(j, je) + ts
        wsav(j, je) = wsav(j, je) + ws(j, je)

        !**** CALC f, ff, huharv()
        f = g(j, je) / (g(j, je) + exp(ald1(nveg(j, je)) - ald2(nveg(j, je)) * g(j, je)))
        ff = f - huharv(j, je)
        huharv(j, je) = f

        !**** CALC SUM STRESS DAYS
        sdw = sdw + (1. - ws(j, je)) * flu(j)
        sdt = sdt + (1. - ts) * flu(j)

        !**** CALC LAI and adjust for lower limit of LAI for forest alnm()
        if (g(j, je) .le. dlai(nveg(j, je)) ) then
          if (alai(j, je) .gt. blai(nveg(j, je)) ) alai(j, je) = blai(nveg(j, je))
          deltalai = ff * blai(nveg(j, je)) * (1. - exp(5. * (alai(j, je) - blai(nveg(j, je))))) * sqrt(reg)
          alai(j, je) = alai(j, je) + deltalai
          if (alai(j, je) .gt. blai(nveg(j, je)) ) alai(j, je) = blai(nveg(j, je))
          olai(j, je) = alai(j, je)
          if (alai(j, je) .lt. almn(nveg(j, je)) ) alai(j, je) = almn(nveg(j, je))
        else
          yy = sqrt(1. - g(j, je))
          zz = 1. / sqrt(1. - dlai(nveg(j, je)))
          alai(j, je) = zz * olai(j, je) * yy
          if (alai(j, je) .lt. almn(nveg(j, je)) ) alai(j, je) = almn(nveg(j, je))
        end if ! ( g(j, je) .le. dlai(nveg(j, je)) )
      end if ! (g(j, je) .le. 1)
      !*********************************************************** END IF (G<=1)
    end if ! ( idorm(j, je) == 0)

  end subroutine vegetation_process

  subroutine vegetation_store_output(j, je)
    use output, only : output_store_hydrotope_value
    integer, intent(in) :: j, je

    call output_store_hydrotope_value(vegetation_water_stress_oid, j, je, ws(j, je))
    call output_store_hydrotope_value(vegetation_temp_stress_oid, j, je, ts)
    call output_store_hydrotope_value(heat_unit_fract_output_id, j, je, g(j, je))
    call output_store_hydrotope_value(biomass_total_output_id, j, je, dm(j, je))
    call output_store_hydrotope_value(leaf_area_index_output_id, j, je, alai(j, je))
    call output_store_hydrotope_value(root_depth_output_id, j, je, rd(j, je))
  end subroutine vegetation_store_output

end module vegetation
