nutrient_phosphorus_uptake Subroutine

  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