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