! 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 ! 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 / & 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