!     FILE ncycle.f
!
!     SUBROUTINES IN THIS FILE CALLED FROM
!     subroutine ncycle(j, je, k, n)   hydrotop
!     subroutine nlch(j, je, k, n) h ydrotop
!     subroutine nuptake(j, je, k, n) crpmd(), veget()
!     subroutine npstress(u1, u2, uu) nuptake, puptake
!     subroutine fert(j, je, n, k, ii) hydrotop
module nutrient

  use utilities, only : dp

  implicit none

  ! Switch to enable/disable nutrient calculation
  logical, save :: nutrient_on = .false.
  ! subsurface flow in HYDROTOPE, mm
  real(dp), save :: ssf
  ! N - NO3 loss with surface flow, kg / ha
  real(dp), save :: yno3
  ! N uptake by the crop for a given day, kg / ha (SUPPLY)
  real(dp), save :: uno3
  ! optimal conc N in biomass, kg / kg
  real(dp), save :: cnb
  ! N stress for plant growth
  real(dp), save :: strsn
  ! P uptake in hydrotope, kg / ha
  real(dp), save :: uap
  ! P stress for plant growth
  real(dp), save :: strsp
  ! N loss with subsurface flow, kg / ha
  real(dp), save :: ssfn
  ! N - NO3 leaching to g - w, kg / ha
  real(dp), save :: percn
  ! soluble P leaching, kg / ha
  real(dp), save :: ysp
  ! N mineralization from humus mineralization, kg / ha
  real(dp), save :: humn
  ! combined water / temperature factor
  real(dp), save :: csf
  ! SUM(hump): humus P - mineralization for basin, kg / h
  real(dp), save :: shump
  ! SUM(fomp): fresh organic P mineralisation for basin, kg / ha
  real(dp), save :: sfomp
  ! SUM(aspflow): flow between active and stable P pool for basin, kg / ha
  real(dp), save :: saspf
  ! SUM(alpflow): flow between active and labile P pool for basin, kg / ha
  real(dp), save :: salpf
  !
  real(dp), save, dimension(20) :: xnflow
  real(dp), save :: shumn = 0. ! SUM(humn): humus N minerlisation for basin, kg / ha
  real(dp), save :: sasnf = 0. ! SUM(asnflow): flow between active and stable org. N for basin, kg / ha
  real(dp), save :: sfomn = 0. ! SUM(fomn): mineralisation from fresh org. N for basin, kg / ha
  real(dp), save :: sdnit = 0. ! SUM(denit): daily N - NO3 loss by denitrification for basin, kg / ha
  real(dp), save :: sbnup = 0. ! sum of N uptake for basin
  real(dp), save :: sbpup = 0. ! weighted average P uptake in the basin, kg / ha
  ! switch code to print from ncycle(), left to keep aggregation, should be purged
  integer, save :: inutr = 1, inusb = 1, inuhd = 1
  ! nitrate (NO3 - N) content in a layer
  real(dp), save, dimension(:, :, :), allocatable :: ano3
  ! active org. N content in a layer, kg / ha
  real(dp), save, dimension(:, :, :), allocatable :: anora
  ! stable organic N in sub j, HRU je, layer l
  real(dp), save, dimension(:, :, :), allocatable :: anors
  ! labile P content in a layer, kg / ha
  real(dp), save, dimension(:, :, :), allocatable :: plab
  ! organic P in a layer, kg / ha
  real(dp), save, dimension(:, :, :), allocatable :: porg
  ! active mineral P in a layer, kg / ha
  real(dp), save, dimension(:, :, :), allocatable :: pma
  ! stable mineral P in a layer, kg / ha
  real(dp), save, dimension(:, :, :), allocatable :: pms
  ! monthly flows for water and N (see writhru.f)
  real(dp), save, dimension(:, :, :), allocatable :: dflow
  ! normal fraction of N in plant biomass at emergence
  real(dp), save, dimension(:), allocatable :: bn1
  ! nitrogen uptake parameter #2: normal fraction of N
  real(dp), save, dimension(:), allocatable :: bn2
  ! nitrogen uptake parameter #3: normal fraction of N
  real(dp), save, dimension(:), allocatable :: bn3
  ! normal fraction of P in plant biomass at emergence
  real(dp), save, dimension(:), allocatable :: bp1
  ! phosphorus uptake parameter #2: normal fraction of P
  real(dp), save, dimension(:), allocatable :: bp2
  ! phosphorus uptake parameter #3: normal fraction of P
  real(dp), save, dimension(:), allocatable :: bp3
  ! coef used to calculate sp1 - S - curve parameter
  real(dp), save, dimension(:), allocatable :: bnu1
  ! coef used to calculate sp2 - S - curve parameter
  real(dp), save, dimension(:), allocatable :: bnu2
  ! shape parameter to calc optimal P fraction in crop biomass
  real(dp), save, dimension(:), allocatable :: bpu1
  ! shape parameter to calc optimal P fraction in crop biomass
  real(dp), save, dimension(:), allocatable :: bpu2
  ! org N loss with erosion, kg / ha
  real(dp), save :: yon
  ! P org. loss with erosion, kg / ha
  real(dp), save :: yph
  ! N retention in surface flow, days
  real(dp), save :: retNsur = 5.0
  ! N retention in subsurface flow, days
  real(dp), save :: retNsub = 365.0
  ! N retention in groundwater, days
  real(dp), save :: retNgrw = 15000.0
  ! P retention in surface flow, days
  real(dp), save :: retPsur = 20.0
  ! N degradation in surface flow, 1 / day
  real(dp), save :: degNsur = 0.02
  ! N degradation in subsurface flow, 1 / day
  real(dp), save :: degNsub = 0.3
  ! N degradation in groundwater, 1 / day
  real(dp), save :: degNgrw = 0.3
  ! P degradation in surface flow, 1 / day
  real(dp), save :: degPsur = 0.02

  real(dp), save :: cdn = - 1.4
  real(dp), save :: cmn = .0003
  real(dp), save :: rcn = 1.
  real(dp), save :: bk = .0006
  real(dp), save :: dkd = 175.

  namelist / NUTRIENT_PARAMETERS / &
    nutrient_on, &
    degNgrw, &
    degNsub, &
    degNsur, &
    degPsur, &
    retNgrw, &
    retNsub, &
    retNsur, &
    retPsur

contains

  subroutine nutrient_initialise(mb, mcrdb, meap, ml)
    use input, only: get_config_fid
    integer, intent(in) :: mb, mcrdb, meap, ml

    read(get_config_fid(), NUTRIENT_PARAMETERS)
    call nutrient_allocate(mb, mcrdb, meap, ml)

  end subroutine nutrient_initialise

  subroutine nutrient_allocate(mb, mcrdb, meap, ml)
    integer, intent(in) :: mb, mcrdb, meap, ml
    allocate(ano3(mb, meap, ml))
    allocate(anora(mb, meap, ml))
    allocate(anors(mb, meap, ml))
    allocate(bn1(mcrdb))
    allocate(bn3(mcrdb))
    allocate(bnu1(mcrdb))
    allocate(bnu2(mcrdb))
    allocate(bp1(mcrdb))
    allocate(bp2(mcrdb))
    allocate(bn2(mcrdb))
    allocate(bp3(mcrdb))
    allocate(bpu1(mcrdb))
    allocate(bpu2(mcrdb))
    allocate(dflow(mb, meap, 20))
    allocate(plab(mb, meap, ml))
    allocate(pma(mb, meap, ml))
    allocate(pms(mb, meap, ml))
    allocate(porg(mb, meap, ml))
    ano3 = 0.
    anora = 0.
    anors = 0.
    bn1 = 0.
    bn3 = 0.
    bnu1 = 0.
    bnu2 = 0.
    bp1 = 0.
    bp3 = 0.
    bpu1 = 0.
    bpu2 = 0.
    bn2 = 0.
    bp2 = 0.
    dflow = 0.
    plab = 0.
    pma = 0.
    pms = 0.
    porg = 0.
  end subroutine nutrient_allocate

  subroutine dealloc_nutrient
   deallocate(ano3)
   deallocate(anora)
   deallocate(anors)
   deallocate(bn1)
   deallocate(bn3)
   deallocate(bnu1)
   deallocate(bnu2)
   deallocate(bp1)
   deallocate(bp3)
   deallocate(bpu1)
   deallocate(bpu2)
   deallocate(bn2)
   deallocate(bp2)
   deallocate(dflow)
   deallocate(plab)
   deallocate(pma)
   deallocate(pms)
   deallocate(porg)
  end subroutine dealloc_nutrient

  subroutine nutrient_nitrate_cycle(j, je, k, cbn, flu, fon, fop, frar, nn, rsd, rtn, ste, te, ul, yone)
    !**** PURPOSE: THIS SUBROUTINE ESTIMATES DAILY N AND P MINERALIZATION &
    !     IMMOBILIZATION CONSIDERING FRESH ORGANIC MATERIAL (CROP RESIDUE)
    !     AND ACTIVE AND STABLE HUMUS MATERIAL.
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      ano3(j, je, l) = nitrate (NO3-N) content in a layer, kg/ha
    !      anora(j, je, l) = active org. N content in a layer, kg/ha
    !      anors(j, je, l) = stable org. N content in a layer, kg/ha
    !      cbn(l, k) = organic carbon content in a layer, %
    !      csf = combined water/temperature factor
    !      dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
    !      flu(j) = fraction of subbasin area in the basin
    !      fon(j, je, l) = fresh organic N from residue in a layer, kg/ha
    !      fop(j, je, l) = fresh organic P from residue in a layer, kg/ha
    !      frar(j, je) = fractional areas of hydrotope in subbasin
    !      humn = N mineralization from humus mineralization, kg/ha
    !      ida = current day
    !      inuhd = number of hydrotope to print from ncycle(), if inutr=1
    !      inusb = number of subbasin to print from ncycle(), if inutr=1
    !      inutr = switch code to print from ncycle()
    !      nn = number of soil layers, from subbasin
    !      plab(j, je, l) = labile P content in a layer, kg/ha
    !      qd = daily surface runoff, mm
    !      rsd(j, je, 2) = crop residue, kg/ha
    !      rtn = active N pool fraction, = 0.15
    !      sasnf = SUM(asnflow): flow between active and stable org. N
    !                       for basin, kg/ha
    !      sdnit = SUM(denit): daily N-NO3 loss by denitrification
    !                       for basin, kg/ha
    !      sfomn = SUM(fomn): mineralisation from fresh org. N for basin,
    !                       kg/ha
    !      shumn = SUM(humn): humus N minerlisation for basin, kg/ha
    !      ste(j, je, l) = water storage in a layer, mm, recalc here
    !      te(j, je, l) = daily ave temp at the bottom of each layer, degree C
    !      ul(l, k) = upper limit water content in layer, mm
    !      xnflow(1:17) = N flows for a choosen hydrotope to write in nutr.prn
    !      xnflow(1) = N loss with surface flow calc nlch, kg/ha
    !      xnflow(2) = N loss with subsurface flow 	calc nlch, kg/ha
    !      xnflow(3) = N loss with percolation 	calc nlch, kg/ha
    !      xnflow(4) = N concentration in layer 2 calc nlch, kg/ha
    !      xnflow(5) = N input with precip calc nlch, kg/ha
    !      xnflow(6) = N loss from layer 2 calc nlch, kg/ha
    !      xnflow(7) = N fertilization         calc fert, kg/ha
    !      xnflow(8) = N uptake by plants calc nuptake, kg/ha
    !      xnflow(9) = N denit calc ncycle, kg/ha
    !      xnflow(10) = N miner from fresh org N calc ncycle, kg/ha
    !      xnflow(11) = N miner from humus calc ncycle, kg/ha
    !      xnflow(12) = xhumcdg/xhumn calc ncycle
    !      xnflow(13) = xhumsut/xhumn calc ncycle
    !      xnflow(14) = xhumcsf/xhumn calc ncycle
    !      xnflow(15) = xfomcdg 	 calc ncycle
    !      xnflow(16) = xfomsut 	 calc ncycle
    !      xnflow(17) = xfomcsf 	 calc ncycle
    !      yone(j) = org. N lost with erosion (calc in orgnsed), kg/ha
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      asnflow = flow between active and stable org. N Pools, kg/ha
    !      ca = minimum of CN and CP-ratio factor
    !      cdg = temperature factor for humus mineralization
    !      cdn = shape coefficient for combined temp-carbon factor
    !      cmn = humus rate constant for N (normally set to 0.0003)
    !      cnr = CN-ratio
    !      cnrf = CN-ratio facto
    !      cpr = CP-ratio
    !      cprf = CP-ratio factor
    !      decr = decompostion rate for residue
    !      denit = daily N-NO3 loss by denitrification
    !      deth = threshold of soil water content for denitrification
    !      fomn = mineralisation from fresh org. N, kg/ha
    !      ik = local par
    !      l = local par
    !      ll = local par
    !      nraz = local par
    !      r4 = CN- or CP-ratio
    !      resdc = decompostion rate for fresh org. material
    !      sut = soil water factor for humus mineralization
    !      sut4 = local par
    !      xden = accumulated denit
    !      xfomcdg = cdg (temperature factor)
    !      xfomcsf = csf (combined water/temperature factor)
    !      xfomn = accumulated fomn
    !      xfomsut = sut (water factor)
    !      xhumcdg = cdg*humn (temperature factor*humn)
    !      xhumcsf = csf*humn (combined water/temperature factor*humn)
    !      xhumn = accumulated humn
    !      xhumsut = sut*humn (water factor*humn)
    !      xx = local par
    !      xx1 = local par
    !      xx2 = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !**** Include common parameters

    real(dp), dimension(:, :), intent(in) :: cbn
    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:, :, :), intent(inout) :: fon
    real(dp), dimension(:, :, :), intent(in) :: fop
    real(dp), dimension(:, :), intent(in) :: frar
    integer, intent(in) :: nn
    real(dp), dimension(:, :, :), intent(inout) :: rsd
    real(dp), intent(in) :: rtn
    real(dp), dimension(:, :, :), intent(in) :: ste
    real(dp), dimension(:, :, :), intent(in) :: te
    real(dp), dimension(:, :), intent(in) :: ul
    real(dp), dimension(:), intent(inout) :: yone
    integer j, je, k, l, ll, nraz
    real(dp) asnflow, ca, cdg, cnr, cnrf, cpr, cprf, decr, denit
    real(dp) deth, fomn, r4, resdc, sut, sut4, xden, xx, xx1, xx2, xfomcdg
    real(dp) xfomcsf, xfomn, xfomsut, xhumcdg, xhumcsf, xhumn, xhumsut

    !**** INITIALIZATION
    xx1 = 0.
    xx2 = 0.
    xden = 0.
    xfomn = 0.
    xhumn = 0.
    xhumcdg = 0.
    xhumsut = 0.
    xhumcsf = 0.
    xfomcdg = 0.
    xfomsut = 0.
    xfomcsf = 0.
    nraz = 0

    !**** SUBTRACT org N lost with erosion
    anors(j, je, 1) = anors(j, je, 1) - yone(j)
    if (anors(j, je, 1) .le. 0.) anors(j, je, 1) = 0.
    yone(j) = 0.

    !*********************************************************** START OF CYCLE 10
    do l = 1, nn

      !**** CALC soil water factor for humus mineralization: sut
      ll = l
      if (l .eq. 1) ll = 2
      xx = te(j, je, ll)
      sut = ste(j, je, ll) / (ul(ll, k) + 1.e-10)
      if (sut .gt. 1.) sut = 1.
      if (sut .lt. 0.) sut = 1.e-10
      sut4 = .06 * exp(3. * sut)

      !**** CALC temperature factor for humus mineralization: cdg
      if (xx .gt. 0.) then
        cdg = xx / (xx + exp(6.82-.232 * xx) + 1.e-6)
      else
        cdg = 0.
      endif

      !**** CALC combined water & temperature factor: csf
      csf = sqrt(cdg * sut)

      !**** CALC daily NO3-N loss by denitrification: denit, xden
      deth = 0.9
      if (sut .ge. deth) then
        denit = sut4 * ano3(j, je, l) * (1. - exp(cdn * cdg * cbn(l, k)))
      else
        denit = 0.
      end if
      ano3(j, je, l) = ano3(j, je, l) - denit
      xden = xden + denit

      !**** CALC asnflow - N flow between act. & stab. org N; RECALC pools
      asnflow = .1e-4 * (anora(j, je, l) * (1. / rtn-1.)-anors(j, je, l))
      anors(j, je, l) = anors(j, je, l) + asnflow
      anora(j, je, l) = anora(j, je, l) - asnflow

      !**** CALC humus mineralization: humn; RECALC pools
      humn = cmn * csf * anora(j, je, l)
      xx = anora(j, je, l) + anors(j, je, l)
      xhumcdg = xhumcdg + cdg * humn
      xhumsut = xhumsut + sut * humn
      xhumcsf = xhumcsf + csf * humn
      if (humn .gt. 0.) nraz = nraz + 1
      anora(j, je, l) = anora(j, je, l) - humn
      ano3(j, je, l) = ano3(j, je, l) + humn
      xhumn = xhumn + humn

      !**** CALC mineralization of fresh organic matter: fomn; RECALC pools
      !       CALC residue decomposition (only here, not in pcycle!)
      if (l .le. 2) then
        r4 = .58 * rsd(j, je, l)
        cnr = r4 / (fon(j, je, l) + ano3(j, je, l) + 1.e-6)
        cpr = r4 / (fop(j, je, l) + plab(j, je, l) + 1.e-6)
        if (cnr .gt. 25.) then
          cnrf = exp(- .693 * (cnr - 25.) / 25.)
        else
          cnrf = 1.0_dp
        end if
        if (cpr .gt. 200.) then
          cprf = exp(- .693 * (cpr - 200.) / 200.)
        else
          cprf = 1.0_dp
        end if

        ca = amin1(real(cnrf, 4), real(cprf, 4))
        decr = .05 * ca * csf

        fomn = decr * fon(j, je, l)
        resdc = decr * rsd(j, je, l)

        rsd(j, je, l) = rsd(j, je, l) - resdc
        fon(j, je, l) = fon(j, je, l) - fomn
        xfomcdg = cdg
        xfomsut = sut
        xfomcsf = csf
      else
        fomn = 0.
      end if

      anora(j, je, l) = anora(j, je, l) + .2 * fomn
      ano3(j, je, l) = ano3(j, je, l) + .8 * fomn
      xfomn = xfomn + .8 * fomn

      !**** CALC SUMS for basin
      shumn = shumn + humn * flu(j) * frar(j, je)
      sasnf = sasnf + asnflow * flu(j) * frar(j, je)
      sfomn = sfomn + fomn * flu(j) * frar(j, je)
      sdnit = sdnit + denit * flu(j) * frar(j, je)

    end do
    !*********************************************************** END OF CYCLE 10

    !**** CALC N flows for a choosen hydrotope (output in ncycle)
    if (inutr .eq. 1 .and. j .eq. inusb .and. je .eq. inuhd) then
      xnflow(9) = xden
      xnflow(10) = xfomn
      xnflow(11) = xhumn
    endif

    !**** CALC monthly flows for selected HRUs (output in writhru.f)
    dflow(j, je, 16) = dflow(j, je, 16) + xfomn
    dflow(j, je, 17) = dflow(j, je, 17) + xhumn
    dflow(j, je, 9) = dflow(j, je, 9) + xden

    return
  end subroutine nutrient_nitrate_cycle



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



  subroutine nutrient_nitrate_leaching(j, je, k, flate, nn, poe, preinf, qd, ul)
    !**** PURPOSE: THIS SUBROUTINE CALCULATES NITRATE LEACHING
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      ano3(j, je, l) = nitrate (NO3-N) content in a layer, kg/ha
    !      anora(j, je, l) = active org. N content in a layer, kg/ha
    !      anors(j, je, l) = stable org. N content in a layer, kg/ha
    !      dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
    !      flate(j, je, l) = subsurface flow, mm, from purk
    !      inuhd = number of hydrotope to print from ncycle(), if inutr=1
    !      inusb = number of subbasin to print from ncycle(), if inutr=1
    !      nn = number of soil layers, from subbasin
    !      percn = N leaching to g-w, kg/ha
    !      poe(j, je, l) = percolation, mm, from purk
    !      precip = precipitation, mm
    !      qd = daily surface runoff, mm
    !      ssfn = N loss with subsurface flow, kg/ha
    !      ul(l, k) = upper limit water content in layer, mm
    !      xnflow(1:17) = N flow for a choosen hydrotope to write in nutr.prn
    !                      (see ncycle)
    !      yno3 = N loss with surface flow, kg/ha
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      co = average daily concentration of N-NO3 in the layer, kg/ha
    !      l = local par
    !      qip = input with precip, then - input in a layer from the layer above
    !      rcn = local par
    !      sro = local variable: surface runoff
    !      vno3 = amount of N-NO3 lost from the layer
    !      vv = total amount of water lost from the soil layer
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:, :, :), intent(in) :: flate
    integer, intent(in) :: nn
    real(dp), dimension(:, :, :), intent(in) :: poe
    real(dp), dimension(:, :), intent(in) :: preinf
    real(dp), intent(in) :: qd
    real(dp), dimension(:, :), intent(in) :: ul
    integer j, je, k, l
    real(dp) co, qip, sro, vno3, vv

    !**** INITIALIZATION
    qip = .01 * rcn * preinf(j, je)
    if (j .eq. inusb .and. je .eq. inuhd) xnflow(5) = qip
    sro = qd
    ssfn = 0.
    xnflow(6) = 0.
    vno3 = 0.
    co = 0.

    !**** CALC nitrate loss in surface and subsurface runoff: yno3, ssfn
    !     RECALC ano3()
    do l = 1, nn
      vv = poe(j, je, l) + sro + flate(j, je, l) + 1.e-10
      ano3(j, je, l) = ano3(j, je, l) + qip
      vno3 = ano3(j, je, l) * (1. - exp(- vv / ul(l, k)))
      co = vno3 / vv
      if (l .eq. 1.) yno3 = qd * co
      ano3(j, je, l) = ano3(j, je, l) - vno3
      qip = co * poe(j, je, l)
      sro = 0.
      ssfn = ssfn + co * flate(j, je, l)
    end do

    !**** CALC nitrate leaching into ground water: percn
    percn = qip

    !**** CALC N flows for a choosen hydrotope (output in ncycle)
    if (j .eq. inusb .and. je .eq. inuhd) then
      xnflow(1) = yno3
      xnflow(2) = ssfn
      xnflow(3) = percn
    endif
    if (j .eq. inusb .and. je .eq. inuhd .and. l .eq. 2) &
      xnflow(4) = co
    if (j .eq. inusb .and. je .eq. inuhd .and. l .eq. 2) &
      xnflow(6) = xnflow(6) + vno3 - qip

    !**** CALC monthly flows for selected HRUs (output in writhru.f)
    dflow(j, je, 5) = dflow(j, je, 5) + yno3
    dflow(j, je, 6) = dflow(j, je, 6) + ssfn
    dflow(j, je, 7) = dflow(j, je, 7) + percn
    dflow(j, je, 10) = dflow(j, je, 10) + ano3(j, je, 1) + &
        ano3(j, je, 2) + ano3(j, je, 3) + ano3(j, je, 4) + ano3(j, je, 5)
    dflow(j, je, 11) = dflow(j, je, 11) + anora(j, je, 1) + anora(j, je, 2) &
      + anora(j, je, 3) + anora(j, je, 4) + anora(j, je, 5)
    dflow(j, je, 12) = dflow(j, je, 12) + anors(j, je, 1) + anors(j, je, 2) &
      + anors(j, je, 3) + anors(j, je, 4) + anors(j, je, 5)

    dflow(j, je, 13) = dflow(j, je, 13) + ano3(j, je, 3)
    dflow(j, je, 14) = dflow(j, je, 14) + ano3(j, je, 4)
    dflow(j, je, 15) = dflow(j, je, 15) + ano3(j, je, 5)

    return
  end subroutine nutrient_nitrate_leaching



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



  subroutine nutrient_nitrate_uptake(j, je, nv, dm, flu, frar, g, ida, nn, snup)
    !**** PURPOSE: CALCULATES N UPTAKE by PLANTS, calls npstress()
    !**** CALLED IN: GROWTH
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      ano3(j, je, l) = nitrate (NO3-N) content in a layer , kg/ha
    !      bn1(nv) = normal fraction of N in plant biomass at emergence
    !      bn3(nv) = normal fraction of N in plant biomass at maturity
    !      bnu1(nv) = coef used to calculate sp1 - S-curve parameter
    !      bnu2(nv) = coef used to calculate sp2 - S-curve parameter
    !      cnb = optimal conc N in biomass, kg/kg
    !      dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
    !      dm(j, je) = total biomass, kg/ha
    !      flu(j) = fraction of subbasin area in the basin
    !      frar(j, je) = fractional areas of hydrotope in subbasin
    !      g(j, je) = fraction of heat units to maturity accumulated
    !      ida = current day
    !      inuhd = number of hydrotope to print from ncycle(), if inutr=1
    !      inusb = number of subbasin to print from ncycle(), if inutr=1
    !      nn = number of soil layers, from subbasin
    !      sbnup = sum of N uptake for basin
    !      snup(j, je) = N uptake accumulated, kg/ha
    !      strsn = N stress factor
    !      uno3 = N uptake by the crop for a given day, kg/ha (SUPPLY)
    !      xnflow(1:17) = N flows for a choosen hydrotope to write in nutr.prn
    !                      (see ncycle)
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      l = local par
    !      sp1 = local par
    !      sp2 = local par
    !      uno3pot = optimal N uptake by the crop until given day. kg/ha (DEMAND)
    !      uu = nutrient stress
    !      xx = local par
    !      yy = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:, :), intent(in) :: dm
    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:, :), intent(in) :: frar
    real(dp), dimension(:, :), intent(in) :: g
    integer, intent(in) :: ida
    integer, intent(in) :: nn
    real(dp), dimension(:, :), intent(inout) :: snup
    integer j, je, nv, l
    real(dp) sp1, sp2, uno3pot, uu, xx, yy
    logical update_uno3

    update_uno3 = .true.

    !**** INITIALIZATION
    xx = g(j, je)
    sp1 = bnu1(nv)
    sp2 = bnu2(nv)

    !**** CALC N uptake by plants & RECALC pools
    !#### CALL nutrient stress
    cnb = bn3(nv) + &
      (bn1(nv) - bn3(nv)) * (1 - xx / (xx + exp(sp1 - sp2 * xx)))
    uno3 = cnb * dm(j, je) - snup(j, je)
    if (uno3 .ge. 0) then
      if (ida .eq. 1) uno3 = 0.
      uno3pot = uno3
      xx = uno3
      do l = 1, nn
        yy = ano3(j, je, l) - xx
        if (yy .gt. 0.) then
          ano3(j, je, l) = yy
          update_uno3 = .false.
          exit
        else
          xx = xx - ano3(j, je, l)
          ano3(j, je, l) = 0.
        end if
      end do
      if (update_uno3) uno3 = uno3 - xx
      if (uno3 .lt. 0.) uno3 = 0.
      snup(j, je) = snup(j, je) + uno3
      sbnup = sbnup + uno3 * flu(j) * frar(j, je)
      call vegetation_nutrient_stress(uno3, uno3pot, uu)
      strsn = uu
    else
      uno3 = 0.
      strsn = 1.
    end if

    !**** CALC N flows for a choosen hydrotope (output in ncycle)
    if (j .eq. inusb .and. je .eq. inuhd) xnflow(8) = uno3

    !**** CALC monthly flows for selected HRUs (output in writhru.f)
    dflow(j, je, 8) = dflow(j, je, 8) + uno3

    return
  end subroutine nutrient_nitrate_uptake



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



  subroutine nutrient_fertilisation(j, je, ii, fen, feno, fep)
    !**** PURPOSE: this subroutine applies N and P fertilizers as
    !              specified by date and amount in initcrop()
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      ano3(j, je, l) = nitrate (NO3-N) content in a layer
    !      anora(j, je, l) = active org. N content in a layer
    !      dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
    !      fen(ii) = amount of min N fertilizers applied, kg N/ha
    !      feno(ii) = amount of org N fertilizers applied, kg N/ha
    !      fep(ii) = amount of P fertilizers applied, kg P/ha
    !      ida = current day
    !      inuhd = number of hydrotope to print from ncycle(), if inutr=1
    !      inusb = number of subbasin to print from ncycle(), if inutr=1
    !      plab(j, je, l) = labile P content in a layer, kg/ha
    !      xnflow(7) = NO3-N fertilization in a choosen hydrotope, kg N/ha
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:), intent(in) :: fen
    real(dp), dimension(:), intent(in) :: feno
    real(dp), dimension(:), intent(in) :: fep
    integer j, je, ii

    !**** ADD FERTILIZERS
    ano3(j, je, 2) = ano3(j, je, 2) + fen(ii)
    anora(j, je, 2) = anora(j, je, 2) + feno(ii)
    plab(j, je, 1) = plab(j, je, 1) + fep(ii)

    !**** CALC N flows for a choosen hydrotope (output in ncycle)
    if (j .eq. inusb .and. je .eq. inuhd) xnflow(7) = fen(ii)

    !**** CALC monthly flows for selected HRUs (output in writhru.f)
    dflow(j, je, 20) = dflow(j, je, 20) + fen(ii) + feno(ii)

    return
  end subroutine nutrient_fertilisation


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


  subroutine nutrient_phosphorus_cycle(j, je, flu, fon, fop, frar, nn, psp, rsd, yphe)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES P FLUX BETWEEN
    !              THE LABILE, ACTIVE MINERAL AND STABLE MINERAL P POOLS
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      ano3(j, je, l) = nitrate (NO3-N) content in a layer, kg/ha
    !      anora(j, je, l) = active org. N content in a layer, kg/ha
    !      anors(j, je, l) = stable org. N content in a layer, kg/ha
    !      csf = combined water/temperature factor
    !      flu(j) = fraction of subbasin area in the basin
    !      fon(j, je, l) = fresh organic N from residue in a layer, kg/ha
    !      fop(j, je, l) = fresh organic P from residue in a layer, kg/ha
    !      frar(j, je) = fractional areas of hydrotope in subbasin
    !      humn = N mineralization from humus, kg/ha
    !      nn = number of soil layers, from subbasin
    !      plab(j, je, l) = labile P content in a layer, kg/ha
    !      pma(j, je, l) = active mineral P content, kg/ha
    !      pms(j, je, l) = stable mineral P content, kg/ha
    !      porg(j, je, l) = organic P content, kg/ha
    !      psp = phosphorus availability index
    !      rsd(j, je, 2) = crop residue, kg/ha
    !      salpf = SUM(alpflow): flow between active and labile P pool
    !                       for basin, kg/ha
    !      saspf = SUM(aspflow): flow between active and stable P pool
    !                       for basin, kg/ha
    !      sfomp = SUM(fomp): fresh organic P mineralisation for basin,
    !                       kg/ha
    !      shump = SUM(hump): humus P-mineralization for basin, kg/ha
    !      yphe(j)       = org. P loss with erosion, kg/ha
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      alpflow = flow between active and labile P pool
    !      aspflow = flow between active and stable P pool
    !      bk = rate constant set to 0.0006
    !      ca = minimum of CN and CP-ratio factor
    !      cnr = CN-ratio
    !      cnrf = CN-ratio factor
    !      cpr = CP-ratio
    !      cprf = CP-ratio factor
    !      decr = decompostion rate for residue
    !      fomp = mineralization from fresh org. P pool
    !      hump = humus P-mineralization
    !      l = layer
    !      r4 = CN- or CP-ratio
    !      resdc = decompostion rate for fresh org. material
    !      rto = equilibrium constant (normally set to 1)
    !      xx = locale variable
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:, :, :), intent(in) :: fon
    real(dp), dimension(:, :, :), intent(inout) :: fop
    real(dp), dimension(:, :), intent(in) :: frar
    integer, intent(in) :: nn
    real(dp), intent(in) :: psp
    real(dp), dimension(:, :, :), intent(in) :: rsd
    real(dp), dimension(:), intent(inout) :: yphe
    integer j, je, l
    real(dp) alpflow, aspflow, ca, cnr, cnrf, cpr, cprf, decr
    real(dp) fomp, hump, r4, resdc, rto, xx

    fomp = 0.
    hump = 0.

    !**** SUBTRACT org P lost with erosion
    porg(j, je, 1) = porg(j, je, 1) - yphe(j)
    if (porg(j, je, 1) .le. 0.) porg(j, je, 1) = 0.
    yphe(j) = 0.

    !*********************************************************** START OF CYCLE 10
    do l = 1, nn

      !**** CALC humus mineralization for P: hump
      xx = anora(j, je, l) + anors(j, je, l)
      if (xx .gt. 0.) then
        hump = 1.4 * humn * porg(j, je, l) / xx
      else
        hump = 0.
      endif

      !**** RECALC pools
      porg(j, je, l) = porg(j, je, l) - hump
      plab(j, je, l) = plab(j, je, l) + hump

      !**** CALC fresh organic matter mineralization: fomp & RECALC pools
      if (l .le. 2) then
        r4 = .58 * rsd(j, je, l)
        cnr = r4 / (fon(j, je, l) + ano3(j, je, l) + 1.e-6)
        cpr = r4 / (fop(j, je, l) + plab(j, je, l) + 1.e-6)
        if (cnr .gt. 25.) then
          cnrf = exp(- .693 * (cnr - 25.) / 25.)
        else
          cnrf = 1.
        end if
        if (cpr .gt. 200.) then
          cprf = exp(- .693 * (cpr - 200.) / 200.)
        else
          cprf = 1.
        end if
        ca = amin1(real(cnrf, 4), real(cprf, 4))
        decr = .05 * ca * csf
        fomp = decr * fop(j, je, l)
        resdc = decr * rsd(j, je, l)
        fop(j, je, l) = fop(j, je, l) - fomp
      end if

      if (l .ge. 2) then
        fomp = 0.
      end if

      porg(j, je, l) = porg(j, je, l) - hump + .2 * fomp
      plab(j, je, l) = plab(j, je, l) + hump + .8 * fomp

    end do
    !*********************************************************** END OF CYCLE 10

    rto = psp / (1. - psp)
    do l = 1, nn
      alpflow = (plab(j, je, l) - pma(j, je, l) * rto)
      if (alpflow .lt. 0.) alpflow = alpflow * .1
      aspflow = bk * (4. * pma(j, je, l) - pms(j, je, l))
      if (aspflow .lt. 0.) aspflow = aspflow * .1
      pms(j, je, l) = pms(j, je, l) + aspflow
      pma(j, je, l) = pma(j, je, l) - aspflow + alpflow
      plab(j, je, l) = plab(j, je, l) - alpflow

      shump = shump + hump * flu(j) * frar(j, je)
      sfomp = sfomp + fomp * flu(j) * frar(j, je)
      saspf = saspf + aspflow * flu(j) * frar(j, je)
      salpf = salpf + alpflow * flu(j) * frar(j, je)
    end do
    return
  end subroutine nutrient_phosphorus_cycle



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




  subroutine nutrient_phosphorus_loss(j, je, k, bd, qd)
    !**** PURPOSE: COMPUTES soluble P loss with surface runoff
    !**** CALLED IN: HYDROTOP
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      bd(k) = bulk density of the upper soil layer, g/cm3
    !      plab(j, je, 1) = P labile in the I soil layer, kg/ha
    !                     plab(j, je, 1) * (10. / bd(k)), g/t
    !      qd = surface runoff, mm
    !      ysp = soluble P leaching, kg/ha
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      dkd = P conc in sediment divided by that of water in m3/t
    !      xx = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:), intent(in) :: bd
    real(dp), intent(in) :: qd
    integer j, je, k
    real(dp) xx

    ysp = .01 * plab(j, je, 1) * (10. / bd(k)) * qd / dkd
    xx = plab(j, je, 1) - ysp
    if (xx .le. 0.) then
      ysp = 0.
    endif
    plab(j, je, 1) = plab(j, je, 1) - ysp
    return
  end subroutine nutrient_phosphorus_loss



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




  subroutine nutrient_phosphorus_uptake(j, je, nv, dm, flu, frar, g, ida, nn, spup)
    !**** PURPOSE: CALCULATES P UPTAKE by PLANTS, calls npstress
    !**** CALLED IN: GROWTH
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      bp1(nv) = normal fraction of P in plant biomass at emergence
    !      bp3(nv) = normal fraction of P in plant biomass at maturity
    !      bpu1(nv) = used to calculate sp1 - S-curve parameter
    !      bpu2(nv) = used to calculate sp2 - S-curve parameter
    !      dm(j, je) = total biomass, kg/ha
    !      flu(j) = fraction of subbasin area in the basin
    !      frar(j, je) = fractional areas of hydrotope in subbasin
    !      g(j, je) = fraction of heat units to maturity accumulated
    !      ida = current day
    !      nn = number of soil layers, from subbasin
    !      plab(j, je, l) = P content in soil layer, recalc here, kg/ha
    !      sbpup = weighted average P uptake in the basin, kg/ha
    !      spup(j, je) = P uptake in hydrotop, kg/ha
    !      strsp = P stress factor for plants
    !      uap = P uptake by plants for a given day (SUPPLY), kg/ha
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      cpt = optimal conc N in biomass
    !      l = local par
    !      sp1 = local par
    !      sp2 = local par
    !      uapot = optimal N uptake by the crop until given day (DEMAND)
    !      uu = P stress
    !      xx = local par
    !      yy = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), dimension(:, :), intent(in) :: dm
    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:, :), intent(in) :: frar
    real(dp), dimension(:, :), intent(in) :: g
    integer, intent(in) :: ida
    integer, intent(in) :: nn
    real(dp), dimension(:, :), intent(inout) :: spup
    integer j, je, l, nv
    real(dp) cpt, sp1, sp2, uapot, uu, xx, yy
    logical update_uap

    !**** INITIALIZATION
    xx = g(j, je)
    sp1 = bpu1(nv)
    sp2 = bpu2(nv)
    update_uap = .true.

    !**** CALC P uptake by plants & RECALC pools
    !#### CALL nutrient stress
    cpt = bp3(nv) + &
      (bp1(nv) - bp3(nv)) * (1 - xx / (xx + exp(sp1 - sp2 * xx)))
    uap = cpt * dm(j, je) - spup(j, je)
    if (uap .ge. 0) then
      if (ida .eq. 1) uap = 0.
      uapot = uap
      xx = uap
      do l = 1, nn
        yy = plab(j, je, l) - xx
        if (yy .gt. 0.) then
          plab(j, je, l) = yy
          update_uap = .false.
          exit
        else
          xx = xx - plab(j, je, l)
          plab(j, je, l) = 0.
        end if
      end do
      if (update_uap) uap = uap - xx
      if (uap .lt. 0.) uap = 0.
      call vegetation_nutrient_stress(uap, uapot, uu)
      strsp = uu
      spup(j, je) = spup(j, je) + uap
      sbpup = sbpup + uap * flu(j) * frar(j, je)
    else
      uap = 0.
      strsp = 1.
    end if

    return
  end subroutine nutrient_phosphorus_uptake



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


  subroutine vegetation_nutrient_stress(u1, u2, uu)
    !**** PURPOSE: THIS FUNCTION CALCULATES THE NUTRIENT STRESS FACTOR uu
    !              CAUSED BY LIMITED SUPPLY of N or P
    !**** CALLED IN: NUPTAKE, PUPTAKE
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      u1 = N or P uptake in current day (SUPPLY), kg/ha
    !      u2 = optimal uptake of N or P until current day (DEMAND), kg/ha
    !      uu = nutrient stress factor
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


    real(dp) u1, u2, uu

    uu = 200. * (u1 / (u2 + .0001) - .5)
    if (uu .le. 0.) then
      uu = 0.
    else
      uu = uu / (uu + exp(3.535 - .02597 * uu))
    end if

    return
  end subroutine vegetation_nutrient_stress


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

  subroutine nutrient_leaching(j, jea, k, bd, cbn, flate, flu, fon, fop, frar, nn, poe, preinf, psp, qd, rsd, rtn, ste, te, ul, yone, yphe)

    real(dp), dimension(:), intent(in) :: bd
    real(dp), dimension(:, :), intent(in) :: cbn
    real(dp), dimension(:, :, :), intent(in) :: flate
    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:, :, :), intent(inout) :: fon
    real(dp), dimension(:, :, :), intent(inout) :: fop
    real(dp), dimension(:, :), intent(in) :: frar
    integer, intent(in) :: nn
    real(dp), dimension(:, :, :), intent(in) :: poe
    real(dp), dimension(:, :), intent(in) :: preinf
    real(dp), intent(in) :: psp
    real(dp), intent(in) :: qd
    real(dp), dimension(:, :, :), intent(inout) :: rsd
    real(dp), intent(in) :: rtn
    real(dp), dimension(:, :, :), intent(in) :: ste
    real(dp), dimension(:, :, :), intent(in) :: te
    real(dp), dimension(:, :), intent(in) :: ul
    real(dp), dimension(:), intent(inout) :: yone
    real(dp), dimension(:), intent(inout) :: yphe
    integer, intent(in) :: j, jea, k
    !#### CALC NITRATE LEACHING, SOLUBLE P LEACING, N & P CYCLES
    call nutrient_nitrate_leaching(j, jea, k, flate, nn, poe, preinf, qd, ul)
    call nutrient_phosphorus_loss(j, jea, k, bd, qd)
    call nutrient_nitrate_cycle(j, jea, k, cbn, flu, fon, fop, frar, nn, rsd, rtn, ste, te, ul, yone)
    call nutrient_phosphorus_cycle(j, jea, flu, fon, fop, frar, nn, psp, rsd, yphe)
  end subroutine nutrient_leaching

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


end module nutrient
