vegetation_water_stress Subroutine

  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