!     FILE hydrotop.f
!
!     SUBROUTINES IN THIS FILE CALLED FROM
!     subroutine hydrotop(j, jea, k, n) subbasin
!     subroutine init_mgt(j, jea) hydrotop
module hydrotope

  use utilities, only : dp, path_max_length, log_info, log_debug

  implicit none

  ! surface flow in HYDROTOPE, mm
  real(dp), save :: qi
  ! soil water index for hydrotope: swind = swe() / sumfc()
  real(dp), save :: swind
  ! preinf(j, je) - qd, mm
  real(dp), save :: rain
  !
  real(dp), save, dimension(30) :: smm
  !
  real(dp), save, dimension(30) :: smy
  !
  real(dp), save, dimension(30) :: sm
  !
  real(dp), save :: vl
  !
  real(dp), save :: vb
  !
  real(dp), save :: v1
  !
  real(dp), save :: v3
  !
  real(dp), save :: v5
  !
  real(dp), save :: v7
  ! SUM(syq) total SUM (whole period) surface runoff
  real(dp), save, dimension(:), allocatable :: sq
  ! SUM(sysq) total SUM (whole period) sub - surface runoff
  real(dp), save, dimension(:), allocatable :: ssq
  ! SUM(syy) total SUM (whole period) sed. yield
  real(dp), save, dimension(:), allocatable :: sy
  ! monthly SUM of sediment yield for subbasin
  real(dp), save, dimension(:), allocatable :: sym
  ! SUM(smq) annual SUM surface runoff for subbasin
  real(dp), save, dimension(:), allocatable :: syq
  ! SUM(smsq) annual SUM sub - surface runoff
  real(dp), save, dimension(:), allocatable :: sysq
  ! SUM(sym) annual SUM sed. yield for subbasin
  real(dp), save, dimension(:), allocatable :: syy
  !
  real(dp), save, dimension(:), allocatable :: tmpNsur
  !
  real(dp), save, dimension(:), allocatable :: tmpNsub
  !
  real(dp), save, dimension(:), allocatable :: tmpNgrw
  !
  real(dp), save, dimension(:), allocatable :: tmpPsur
  ! soil water content, mm
  real(dp), save, dimension(:, :), allocatable :: swe
  ! lag factor (residue and snow effect on temp)
  real(dp), save, dimension(:, :), allocatable :: bcv
  ! vegetation number (database)
  integer, save, dimension(:, :), allocatable :: nveg
  ! precipitation adjusted for canopy storage, mm
  real(dp), save, dimension(:, :), allocatable :: preinf
  ! sum of upper limit water content in soil, calc in subbasin, mm
  real(dp), save, dimension(:, :), allocatable :: hsumul
  ! sum of field capacity in soil, calc in subbasin, mm
  real(dp), save, dimension(:, :), allocatable :: hsumfc
  ! annual sums, analogue to dflow
  real(dp), save, dimension(:, :, :), allocatable :: dfloy
  ! total period average sums, analogue to dflow
  real(dp), save, dimension(:, :, :), allocatable :: dflav
  integer, save, dimension(4) :: k2 = (/ 135, 28, 203, 85/)
  integer, save, dimension(4) :: k3 = (/ 43, 54, 619, 33/)
  integer, save, dimension(4) :: k4 = (/ 645, 9, 948, 65/)
  integer, save, dimension(4) :: k5 = (/ 885, 41, 696, 62/)
  ! HRU structure vector to define k, n
  integer, save, dimension(:, :, :), allocatable :: mstruc
  ! max number of HRUs in subbasin (counted)
  integer, save :: meap
  ! fractional area of HYDROTOPE in SUBBASIN
  real(dp), save, dimension(:, :), allocatable :: frar


  ! Hydrotope data read from hydrotope.csv
  integer, dimension(:), allocatable :: hydrotope_ids
  integer, dimension(:), allocatable :: hydrotope_subbasin_id
  integer, dimension(:), allocatable :: landuse_ids
  integer, dimension(:), allocatable :: soil_ids
  integer, dimension(:), allocatable :: elevations
  integer, dimension(:), allocatable :: glaciers
  real(dp), dimension(:), allocatable :: hydrotope_area
  integer, dimension(:), allocatable :: irrigations
  integer, dimension(:), allocatable :: wetland
  integer, dimension(:), allocatable :: crop_management_id

  integer :: hydrotope_input_file_id

  ! Number of rows in `hydrotope.csv`; set by read_allocate_params
  integer, save :: hydrotope_csv_size

  ! Output variable ids
  integer :: &
    precipitation_output_id = 0, &
    tmin_output_id = 0, &
    tmax_output_id = 0, &
    water_yield_output_id = 0, &
    tmean_output_id = 0, &
    soil_water_content_output_id = 0

  character(len=path_max_length) :: hydrotope_input_file = "hydrotope.csv"

contains

  subroutine hydrotope_initialise(mb, neap, sbar)
    use input, only : get_config_fid
    use output, only: output_register_hydrotope_var, output_register_subbasin_var
    integer, intent(in) :: mb
    integer, intent(inout) :: neap(:)
    real(dp), intent(inout) :: sbar(:)

    precipitation_output_id = output_register_hydrotope_var("precipitation", .false.)
    tmin_output_id = output_register_subbasin_var("tmin")
    tmax_output_id = output_register_subbasin_var("tmax")
    water_yield_output_id = output_register_subbasin_var("water_yield", .false.)
    tmean_output_id = output_register_subbasin_var("tmean")
    soil_water_content_output_id = output_register_hydrotope_var("soil_water_content")

    call hydrotope_read_input
    call hydrotope_subbasin_division(mb, neap, sbar)
    call hydrotope_allocate(mb)

  end subroutine hydrotope_initialise

  subroutine hydrotope_subbasin_division(mb, neap, sbar)
    use utilities, only : log_error
    integer, intent(in) :: mb
    integer, intent(inout) :: neap(:)
    real(dp), intent(inout) :: sbar(:)
    integer i, j, jea

    ! Count max number of hydrotopes per subbasin
    meap = 0
    do i = 1, hydrotope_csv_size
      j = hydrotope_subbasin_id(i)
      if (j > 0) then
        neap(j) = neap(j) + 1 ! number of HRUs per subbasin
        sbar(j) = sbar(j) + hydrotope_area(i) ! total subbasin area in m^2
        if (j.gt.mb) call log_error("hydrotope_subbasin_division", &
          "subbasin_id in hydrotope input out of range.")
        if (neap(j) .gt. meap) meap = neap(j)
      else
        exit
      end if ! (j > 0)
    end do ! i = 1, mbeap

    ! Fill frar and mstruct sparse arrays
    allocate(frar(mb, meap))
    frar = 0.
    allocate(mstruc(mb, meap, 7))
    mstruc = 0

    ! fraction of hydrotop area of subbasin area
    i = 1
    do j = 1, mb
      frar(j, :neap(j)) = hydrotope_area(i:i+neap(j)-1) ! area in m^2 !!! recalculated below as fraction of subbasin
      mstruc(j, :neap(j), 1) = landuse_ids(i:i+neap(j)-1) ! land use ID
      mstruc(j, :neap(j), 2) = soil_ids(i:i+neap(j)-1) ! soil type ID
      mstruc(j, :neap(j), 3) = wetland(i:i+neap(j)-1) ! wetland
      mstruc(j, :neap(j), 4) = crop_management_id(i:i+neap(j)-1) ! LU management (this parameter has disappeared from hydrotope.csv)
      mstruc(j, :neap(j), 5) = elevations(i:i+neap(j)-1) ! HRU mean elevation
      mstruc(j, :neap(j), 6) = glaciers(i:i+neap(j)-1) ! initial glacier depth
      mstruc(j, :neap(j), 7) = irrigations(i:i+neap(j)-1) ! initial glacier depth
      i = i + neap(j)
      do jea = 1, neap(j)
        frar(j, jea) = frar(j, jea) / sbar(j)
      end do
    end do

    call log_info("hydrotope_subbasin_division", &
      "Max. number HRUs in subbasin:", int=meap)
  end subroutine hydrotope_subbasin_division

  subroutine hydrotope_read_input
    use input, only : read_integer_column, read_real_column, input_open_file, input_count_rows

    hydrotope_input_file_id = input_open_file(hydrotope_input_file)
    hydrotope_csv_size = input_count_rows(hydrotope_input_file_id, .true.)
    call log_info('hydrotope_read_input', 'Number of hydrotopes:', &
      int=hydrotope_csv_size)

    allocate(hydrotope_ids(hydrotope_csv_size))
    allocate(hydrotope_subbasin_id(hydrotope_csv_size))
    allocate(landuse_ids(hydrotope_csv_size))
    allocate(soil_ids(hydrotope_csv_size))
    allocate(elevations(hydrotope_csv_size))
    allocate(glaciers(hydrotope_csv_size))
    allocate(hydrotope_area(hydrotope_csv_size))
    allocate(irrigations(hydrotope_csv_size))
    allocate(wetland(hydrotope_csv_size))
    allocate(crop_management_id(hydrotope_csv_size))

    call read_integer_column(hydrotope_input_file_id, "subbasin_id", hydrotope_subbasin_id)
    call read_integer_column(hydrotope_input_file_id, "hydrotope_id", hydrotope_ids)
    call read_integer_column(hydrotope_input_file_id, "landuse_id", landuse_ids)
    call read_integer_column(hydrotope_input_file_id, "soil_id", soil_ids)
    call read_integer_column(hydrotope_input_file_id, "elevation", elevations)
    call read_integer_column(hydrotope_input_file_id, "glacier_thickness", glaciers, 0)
    call read_real_column(hydrotope_input_file_id, "area", hydrotope_area)
    call read_integer_column(hydrotope_input_file_id, "irrigation_user_id", irrigations, 0)
    call read_integer_column(hydrotope_input_file_id, "wetland", wetland, 0)
    call read_integer_column(hydrotope_input_file_id, "crop_management_id", crop_management_id, 1)

  end subroutine hydrotope_read_input

  subroutine hydrotope_allocate(mb)
    use utilities, only : random_n
    integer, intent(in) :: mb

    ! Allocate arrays for hydrotope.csv
    allocate(bcv(mb, meap))
    allocate(dflav(mb, meap, 20))
    allocate(dfloy(mb, meap, 20))
    allocate(hsumfc(mb, meap))
    allocate(hsumul(mb, meap))
    allocate(nveg(mb, meap))
    allocate(preinf(mb, meap))
    allocate(sq(mb))
    allocate(ssq(mb))
    allocate(swe(mb, meap))
    allocate(sy(mb))
    allocate(sym(mb))
    allocate(syq(mb))
    allocate(sysq(mb))
    allocate(syy(mb))
    allocate(tmpNgrw(mb))
    allocate(tmpNsub(mb))
    allocate(tmpNsur(mb))
    allocate(tmpPsur(mb))
    bcv = 0.
    dflav = 0.
    dfloy = 0.
    hsumfc = 0.
    hsumul = 0.
    nveg = 0
    preinf = 0.
    sq = 0.
    ssq = 0.
    ssq = 0.
    swe = 0.
    sy = 0.
    sym = 0.
    syq = 0.
    sysq = 0.
    sysq = 0.
    syy = 0.
    tmpNgrw = 0.005
    tmpNsub = 0.005
    tmpNsur = 0.005
    tmpPsur = 0.005
    sym = 0.
    syq = 0.
    sysq = 0.
    syy = 0.
    sq = 0.
    ssq = 0.
    sy = 0.
    smm = 0.
    smy = 0.
    sm = 0.
    preinf = 0.
    hsumul = 0.
    hsumfc = 0.
    dfloy = 0.
    dflav = 0.
    vl = 100.
    vb = 0.
    v1 = random_n(k2)
    v3 = random_n(k3)
    v5 = random_n(k4)
    v7 = random_n(k5)
  end subroutine hydrotope_allocate

  subroutine dealloc_hydrotope
    deallocate(bcv)
    deallocate(dflav)
    deallocate(dfloy)
    deallocate(hsumfc)
    deallocate(hsumul)
    deallocate(mstruc)
    deallocate(nveg)
    deallocate(preinf)
    deallocate(sq)
    deallocate(ssq)
    deallocate(swe)
    deallocate(sy)
    deallocate(sym)
    deallocate(syq)
    deallocate(sysq)
    deallocate(syy)
    deallocate(tmpNgrw)
    deallocate(tmpNsub)
    deallocate(tmpNsur)
    deallocate(tmpPsur)
    deallocate(frar)
  end subroutine dealloc_hydrotope

  subroutine hydrotope_process(j, jea, k, n, wet, bSubcatch, dart, daycounter, flu, frar, ida, iy, iyr, mo, nbyr, precip, sbar, tx)
    !**** PURPOSE: THIS SUBROUTINE CALCULATES ALL PROCESSES in HYDROTOPES
    !**** CALLED IN: SUBBASIN
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !   >>>>> COMMON PARAMETERS & VARIABLES
    !   bcv(j, je) = lag factor (residue and snow effect on temp)
    !   canstor(j, je) = canopy water storage, mm
    !   cn = current CN
    !   cn2(k, n, j) = Curve Numbers for soil k, land use n, and subbasin
    !   cva(j, je) = vegetation cover, kg/ha
    !   dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
    !   eo = potential evapotranspiration, mm
    !   ep = plant transpiration, mm
    !   es = soil evaporation, mm
    !   icc = index for cover crop corr. number in crop database
    !   ida = current day
    !   idfe(n, if) = day of fertilization, if - number of fertilisation
    !   FORD = index for deciduous forest corr. number in crop database
    !   FORE = index for evergreen forest corr. number in crop database
    !   FORM = index for mixed forest corr. number in crop database
    !   RNGB = index for heather corr. number in crop database
    !   imea = index for meadows corr. number in crop database
    !   ipas = index for pasture corr. number in crop database
    !   iwet = index for wetland corr. number in crop database
    !   iwetf = index for forested wetland corr. number in crop database
    !   nn = number of soil layers
    !   nveg(j, je) = number of vegetation from crop database
    !   percn = N-NO3 leaching to g-w, kg/ha
    !   precip = precipitation, mm
    !   preinf(j, je) = precipitation adjusted for canopy storage, mm
    !   prk = lateral subsurface flow, mm, calc in perc for 4mm layers
    !   qd = surface flow in HYDROTOPE, mm
    !   qi = surface flow in HYDROTOPE, mm
    !   rain = preinf(j, je) - qd, mm
    !   sep = percolation, mm
    !   sml = snow melt, calc in snom(), mm
    !   snoa(j, jea) = snow water content in HYDROTOPE, mm
    !   snoev = snow evap. in HYDROTOPE, mm
    !   ssf = subsurface flow in HYDROTOPE, mm
    !   ssfn = N-NO3 in subsurface flow, kg/ha
    !   ste(j, je, l) = water storage, recalc here, mm
    !   strsn = N stress for plant growth
    !   strsp = P stress for plant growth
    !   sumfc(k) = sum of field capacity in soil, mm
    !   swe(j, je) = soil water content, mm
    !   swind = soil water index for hydrotope: swind=swe()/sumfc()
    !   ts = temp. stress
    !   tx(j) = daily average temperature in the subbasin, degree C
    !   uno3 = N uptake by plants, kg/ha
    !   vb = CN max - for output in main
    !   vl = CN min - for output in main
    !   ws(j, je) = water stress factor
    !   yno3 = N-NO3 loss with surface flow, kg/ha
    !   ysp = soluble P leaching, kg/ha
    !   >>>>>

    !   >>>>> STATI! PARAMETERS
    !   ii = local par
    !   l = local par
    !   xx = local par
    !   zz = local par
    !   >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    use crop, only : &
      be, &
      crop_initialise_hydrotope, &
      crop_process, &
      cva, &
      fen, &
      feno, &
      fep, &
      fon, &
      fop, &
      huharv, &
      hun, &
      icc, &
      idfe, &
      igro, &
      mfe, &
      nucr, &
      rdmx, &
      rwt, &
      snup, &
      spup
    use erosion, only : erosion_cklsp_factor, yone, yphe
    use evapotranspiration, only : &
      canstor, &
      eo, &
      es, &
      evapotranspiration_process, &
      humi, &
      pit, &
      ra, &
      snoev, &
      tmn, &
      ylc, &
      yls
    use groundwater, only : additionalGwUptake
    use landuse, only : &
      LULC, &
      landuse_index, &
      landuse_is_cropland, &
      landuse_is_forest, &
      landuse_is_natural_vegetation
    use management, only : &
      TSubbasin, &
      TWaterUser, &
      bWAM_Module, &
      management_is_active_period, &
      management_is_transfer_subbasin, &
      management_subbasin_pointer, &
      management_user_pointer
    use nutrient, only : &
      dflow, &
      nutrient_fertilisation, &
      nutrient_leaching, &
      percn, &
      ssf, &
      ssfn, &
      strsn, &
      strsp, &
      uap, &
      uno3, &
      yno3, &
      ysp
    use snow, only : &
      bSnowModule, &
      precipe, &
      rnew, &
      rsn, &
      sml, &
      snoa, &
      tmax, &
      tmax_tmp, &
      tmin, &
      tmin_tmp, &
      tmit, &
      tmx, &
      tx_tmp
    use soil, only : &
      avt, &
      bd, &
      cbn, &
      cn, &
      cn2, &
      cncor, &
      ek, &
      fc, &
      flate, &
      nn, &
      poe, &
      prk, &
      psp, &
      qd, &
      rtn, &
      sep, &
      soil_curve_number, &
      soil_curve_number_runoff, &
      soil_process, &
      soil_temperature, &
      ste, &
      sumfc, &
      te, &
      ul, &
      z
    use vegetation, only : &
      alai, &
      blai, &
      ep, &
      olai, &
      rsd, &
      to, &
      ts, &
      vegetation_process, &
      vegetation_store_output, &
      ws

    logical, intent(in) :: bSubcatch
    real(dp), dimension(:), intent(in) :: dart
    integer, intent(in) :: daycounter
    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:, :), intent(in) :: frar
    integer, intent(in) :: ida
    integer, intent(in) :: iy
    integer, intent(in) :: iyr
    integer, intent(in) :: mo
    integer, intent(in) :: nbyr
    real(dp), intent(in) :: precip
    real(dp), dimension(:), intent(in) :: sbar
    real(dp), dimension(:), intent(in) :: tx

    ! subbasin number
    integer :: j
    ! number of hydrotope in subbasin
    integer :: jea
    ! land use number from LULC map and * .lut file
    integer :: n
    ! soil number
    integer :: k
    ! wetland code 0 / 1
    integer :: wet
    integer :: ii
    real(dp) :: xx, zz
    ! the curve number values passed to function curno
    real(dp) :: cn2curno
    !#### WATER MANAGEMENT MODULE ####
    real(dp) :: area ! hydrotop area, m2
    !#### WATER MANAGEMENT MODULE ####
    TYPE (TSubbasin), POINTER :: pS ! pointer on subbasin
    !#### WATER MANAGEMENT MODULE ####
    TYPE (TWaterUser), POINTER :: pWU ! pointer on actual TWU

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

    !**** INITIALIZATION
    qd = 0.
    qi = 0.
    ssf = 0.
    sep = 0.
    ep = 0.
    es = 0.
    snoev = 0.
    ysp = 0.
    yno3 = 0.
    uno3 = 0.
    ssfn = 0.
    percn = 0.
    ws(j, jea) = 1.
    ts = 1.
    strsn = 1.
    strsp = 1.

    !**** INITIZALIZE MANAGEMENT OPERATIONS
    if (landuse_is_cropland(n) ) then
      call crop_initialise_hydrotope(j, jea, iy, mstruc)
    end if

    !#### CALL CURNO: to set Curve Number parameters
    if (bSubcatch) then
      cn2curno = 0.
      cn2curno = cn2(j, jea) * cncor(j)
      if (cn2curno * cncor(j) > 100.) cn2curno = 100.
      if (cn2curno * cncor(j) < 1.) cn2curno = 1.
      call soil_curve_number(cn2curno, j, jea, hsumfc, hsumul)
    else
      call soil_curve_number(cn2(j, jea), j, k, hsumfc, hsumul)
    end if

    !**** CALC bcv() - lag factor for soil temperature
    bcv(j, jea) = cva(j, jea) / (cva(j, jea) + exp(7.563-1.297e-4 * cva(j, jea)))
    if (snoa(j, jea) .gt. 0.) then
      if (snoa(j, jea) .le. 120.) then
        xx = snoa(j, jea) / (snoa(j, jea) + exp(6.055 - .3002 * snoa(j, jea)))
      else
        xx = 1.
      end if
      bcv(j, jea) = amax1(real(xx, 4), real(bcv(j, jea), 4))
    end if

    if (bSnowModule) then
      tx_tmp = tmit
      tmax_tmp = tmax
      tmin_tmp = tmin
    else
      tx_tmp = tx(j)
      tmax_tmp = tmx(j)
      tmin_tmp = tmn(j)
    end if
    !#### CALL SOLT - COMPUTE TEMP. of SOIL LAYERS
    call soil_temperature(zz, j, jea, k, bcv, ida, mo, pit, preinf, swe, tmax_tmp, tmin_tmp, tx_tmp)


    !#### CALL VOLQ TO CALC RUNOFF VOLUME,
    !#### CALL ECKLSP TO CALC COMBINED CKLSP factor for hydrotope
      ! preinf = precipitation adjusted for canopy storage
      !        = precipitation + snow melt + canopy water storage

    !###########################
    !#### SNOW MODULE ####
    !###########################
    if (bSnowModule) then
      !preinf(j, jea) = precipe + canstor(j, jea)
      !preinf(j, jea) = precipe + canstor(j, jea) + vsn
      preinf(j, jea) = precipe + canstor(j, jea)
      !if ( vsn > 0. ) preinf(j, jea) = vsn + canstor(j, jea)
    else
      preinf(j, jea) = precip + sml + canstor(j, jea)
    end if
    !preinf(j, jea) = precip + sml + canstor(j, jea)
    !###########################

    !#################################
    !#### WATER MANAGEMENT MODULE ####
    !#################################
    if (bWAM_Module .AND. daycounter > 1) then
      if (management_is_transfer_subbasin(j) ) then
        if (landuse_is_cropland(n) .AND. mstruc(j, jea, 7) >= 1 ) then ! if hydrotope is irrigated cropland
          area = frar(j, jea) * sbar(j)
          pS => management_subbasin_pointer(j)
          !pWU => TWU(pS%pos_irr)
          pWU => management_user_pointer(pS % pos_irr) ! pointer on current water user
          ! if sprinkler irrigation, add available irrigation volume to precipitation
          ! supply in m3/s, precipitation in mm
          ! pWU%supplied is the total amount of water for all irrigation hydrotopes within this subbasin.
          ! By converting to mm, the area weighted volume is maintained.
          if (pWU % wu_opt == 4 .AND. pWU % irr_practice == 1 .AND. management_is_active_period(iyr, ida, pWU) ) &
            preinf(j, jea) = preinf(j, jea) + pWU % supplied(daycounter - 1) * 1000. * 86400. / area !pWU % area
        end if
      end if
    end if
    !#################################

    if (tx_tmp > 0.) then
      if (preinf(j, jea) > 0. ) then
        call soil_curve_number_runoff(j, jea, alai, blai, canstor, igro, nucr, preinf, LULC%canmx(landuse_index(n)))
        qi = qd
        call erosion_cklsp_factor(j, jea, k, cva, ek, igro, nucr, landuse_is_cropland(n), landuse_is_natural_vegetation(n), landuse_is_forest(n))
        if (cn .ge. vl) then
          if (cn .ge. vb) &
            vb = cn
        else
          vl = cn
        end if
      end if
    end if

    !#### CALL PURK TO PERFORM SOIL WATER ROUTING
    rain = preinf(j, jea) - qd

    !#################################
    !#### WATER MANAGEMENT MODULE ####
    !#################################
    if (bWAM_Module .AND. daycounter > 1) then
      if (management_is_transfer_subbasin(j) ) then
        if (landuse_is_cropland(n) .AND. mstruc(j, jea, 7) >= 1 ) then ! if hydrotope is irrigated cropland
          area = frar(j, jea) * sbar(j)
          pS => management_subbasin_pointer(j)
          pWU => management_user_pointer(pS % pos_irr) ! pointer on current water user
          ! if drip irrigation, add available irrigation volume to precipitation, mm
          ! supply in m3/s, precipitation in mm
          ! pWU%supplied is the total amount of water for all irrigation hydrotopes within this subbasin.
          ! By converting to mm, the area weighted volume is maintained.
          if (pWU % wu_opt == 4 .AND. pWU % irr_practice == 2 .AND. management_is_active_period(iyr, ida, pWU) ) &
            rain = rain + pWU % supplied(daycounter - 1) * 1000. * 86400. / area !pWU % area
        end if
      end if
    end if
    !#################################

    call soil_process(j, jea, k, rain, swe)
    ssf = prk

    !#### CALC swe = soil water content [mm]
    ! swe is re-calculated after vegetation growth below
    ! ste = water storage per layer [mm], calculated in purk
    swe(j, jea) = 0.
    swe(j, jea) = sum(ste(j, jea, :))

    call evapotranspiration_process(j, jea, k, alai, cva, ep, ida, mo, nn, preinf, qd, snoa, ste, tx, z, bSnowModule, rnew, tmit, tx_tmp, rsn, LULC%ETcor(landuse_index(n)))
    ! re-calc soil water content
    swe(j, jea) = 0.
    swe(j, jea) = sum(ste(j, jea, :))

    !#### CALL FERT, CRPMD & VEGMD to calculate vegetation growth & fertilisation
    select case(LULC % lutype(landuse_index(n)))
      case (0) !#### ! NO VEGETATION - - > plant transpiration = 0.
        ep = 0.

      case (1) !#### CROPLAND, managed
        ! Fertilizer application
        do ii = 1, mfe
          if (ida == idfe(ii) ) call nutrient_fertilisation(j, jea, ii, fen, feno, fep)
        end do
        call crop_process(j, jea, k, n, wet, additionalGwUptake, avt, bWAM_Module, dart, daycounter, es, fc, flu, frar, humi, ida, iy, iyr, mstruc, nbyr, nn, nveg, pit, ra, sbar, sep, ste, tmn, tx, uap, ylc, yls, z, bSnowModule, tmit)
        call nutrient_leaching(j, jea, k, bd, cbn, flate, flu, fon, fop, frar, nn, poe, preinf, psp, qd, rsd, rtn, ste, te, ul, yone, yphe)

      case(2) !#### NATURAL VEGETATION
        nveg(j, jea) = LULC % veg_code(landuse_index(n))
        call vegetation_process(j, jea, 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)
        call nutrient_leaching(j, jea, k, bd, cbn, flate, flu, fon, fop, frar, nn, poe, preinf, psp, qd, rsd, rtn, ste, te, ul, yone, yphe)

      case(3) !#### WATER
        es = eo

      case(4) !#### FORESTS
        nveg(j, jea) = LULC % veg_code(landuse_index(n))
        call vegetation_process(j, jea, 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)
        call nutrient_leaching(j, jea, k, bd, cbn, flate, flu, fon, fop, frar, nn, poe, preinf, psp, qd, rsd, rtn, ste, te, ul, yone, yphe)

    end select

    call vegetation_store_output(j, jea)

    !**** RE-CALC SOIL WATER CONTENT swe() & SOIL WATER INDEX swind
    swe(j, jea) = 0.
    swe(j, jea) = sum(ste(j, jea, :))

    swind = swe(j, jea) / sumfc(k)

    ! dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
    dflow(j, jea, 1) = dflow(j, jea, 1) + qd
    dflow(j, jea, 2) = dflow(j, jea, 2) + ssf
    dflow(j, jea, 3) = dflow(j, jea, 3) + sep
    dflow(j, jea, 4) = dflow(j, jea, 4) + es + ep

  end subroutine hydrotope_process

end module hydrotope