subroutine nuse (bn1, bn2, bn3, bp1, bp2, bp4)

c     + + + PURPOSE + + +
c     This subroutine calculates the daily soil supply of N & P from each
c     soil layer in which there are roots.

c     + + + ARGUMENT DECLARATIONS + + +
      real bn1, bn2, bn3, bp1, bp2, bp4

c     + + + COMMON BLOCKS + + +
*$noereference
      include 'crop/cgrow.inc'
      include 'crop/csoil.inc'
      include 'crop/cfert.inc'
      include 'crop/cparm.inc'
      include 'crop/cenvr.inc'
      include 'crop/chumus.inc'
      include 'crop/p1crop.inc'
*$reference

c     + + + LOCAL VARIABLES + + +
      real vt, xx, clp, flu, xy
      integer j, k
c
c     + + + LOCAL VARIABLE DEFINITIONS + + +
c     Note:variable names in brackets are the names used in the EPIC manual
c     cnt(cnb) - optimal plant N concentration(kg/t) on day i
c     bn1,bn2,... - crop parameters for plant N concentration equation
c     hui         - heat unit index(0-1) on day i
c     un2(undef.) - optimal crop N concentration (kg/ha) on day i
c     un1(undef.) - sum of soil supplied plus fixed N (kg/ha) on day i ??
c     uno3(und)   - N demand rate for crop (kg/ha/d)
c     dm(b)       - accumulated plant biomass(tops+roots) (t/ha)
c     ddm(undef.) - daily plant biomass accumulation (t/ha/d)
c     cpt(cpb)    - optimal plant P concentration(kg/t) on day i
c     bp          - crop parameters for plant P concentration equation
c     hui         - heat unit index(0-1) on day i
c     up2(undef.) - optimal crop P concentration (kg/ha) on day i
c     up1(undef.) - actual crop P concentration  (kg/ha) on day i
c     upp(UPD)    - P demand rate for crop (kg/ha/d)
c     dm(b)       - accumulated plant biomass(tops+roots) (t/ha)
c     upp(upd)  = P demand rate (kg/ha/d)
c     rw(rwt)   = total root weight upto day i (t/ha)
c     ir        = deepest layer number to which roots have extended
c     un        = rate of N supplied by the soil from layer J (kg/ha/d)
c     wno3      = amount of N in a layer (kg/ha)
c     u         = water uptake from a layer by evpt (mm)
c     st(sw)    = soil water content of a layer (mm)
c     sunn(uns)  = sum of N supplied from all rooted layers (kg/ha)
c     clp(clp)  = P concentration in a layer(g/ton)
c     flu(lfu)  = labile P factor for crop uptake(0-1)
c     up        = soil supply of P from a layer (kg/ha)
c     sup(ups)  = soil supply of P from all rooted layers (kg/ha)
c     ap        = amount of labile P in a layer (kg/ha)
c     xy        =

c     + + + OUTPUT FORMATS + + +
c2000 format(1x,i3,1x,9(f7.3,1x))

c     + + + END OF SPECIFICATIONS + + +
c     This section is the EPIC subroutine NUP
c      dimension fr(10)
      un1=un2
c      sunn=0.
      sup=0.
c     calculate optimal N concentration for a crop using a modified version of
c     eq. 2.215 in the next 2 lines.
      cnt=bn1+bn2*exp(-bn3*hui)
      un2=cnt*dm*1000.
c      if (un2.lt.un1) un2=un1
c     allow positive N demand late in the season ?
c      uno3=amin1(4000.*.0023*ddm,un2-un1)
      uno3=un2-un1
      if (uno3.le.0.) uno3=0.
      vt=uno3
c     This section is the EPIC subroutine NPUP
c     Calculate P concentration for a crop using a modified form of eq.2.229.
      cpt=bp2+bp1*exp(-bp4*hui)
      up2=cpt*dm*1000.
      if (up2.lt.up1) up2=up1
      upp=up2-up1
c
c     This section is the EPIC subroutine NUSE
c     calculate parts of eq. 2.231 --- P demand rate
      if (rw.eq.0.) goto 45
      xx=1.5*upp/rw
c     loop for computing soil supply of N and P from each layer
      do 4 j=1,ir
c        compute soil supply of N in the next 2 lines
c        next line commented out unitl water use data is available
c        un(j)=wno3(j)*u(j)/(st(j)+.001)
c         un(j)=wno3(j)*0.05
c         sunn=sunn+un(j)
c        above 3 lines replaced by the following 5 lines
c         un(j)=uno3*rwt(j)/rw
          if (vt.le.0.) goto 75
          if (wno3(j).eq.0.) goto 75
          xy=wno3(j)-vt
          if (xy.ge.0.) then
             wno3(j)=wno3(j)-vt
             sunn=sunn+vt
             vt=0.
           endif
           if (xy.lt.0.) then
             sunn=sunn+wno3(j)
             vt=vt-wno3(j)
             wno3(j)=0.
           endif
c         xy=wno3(j)-un(j)
c         if (xy.le.0.) un(j)=wno3(j)
c         sunn=sunn+un(j)
c         wno3(j)=wno3(j)-un(j)
c        compute soil supply of P in the next 8 lines
c        F in the next line replaced by CLP
  75  continue
         clp=1000.*ap(j)/wt(j)
c        F in the next and subsequent lines replaced by FLU --- eq 2.232
c        *********** re-arranged to avoid the goto statement **************
         flu=clp/(clp+exp(a_s11-b_s11*clp))
         if (clp.gt.30.) flu=1.
c        use equation 2.231
c         rwt(j) = ???  if this is root mass by layer - change to bcmbgr(j)
c         bcmbgr should also be passed as an argument
         up(j)=xx*flu*rwt(j)
         if (up(j).ge.ap(j)) up(j)=.9*ap(j)
         sup=sup+up(j)
    4 continue
c     the following algorithms may be temporary depending whether subroutine
c     najn & najp are eliminated.
c     adjust soil supply and plant demand for N
c      sum=0.
c      rt=uno3/(sunn+1.e-20)
c      if (rt.lt.1.) then
c         do 2 k=1,ir
c            un(k)=un(k)*rt
c            sum=sum+un(k)
c   2     continue
c        sunn=sum
c      endif
cc     adjust soil supply and plant demand for P
c      sum=0.
c      rt=upp/(sup+1.e-20)
c      if (rt.lt.1.) then
c         do 3 k=1,ir
c            up(k)=up(k)*rt
c            sum=sum+up(k)
c   3     continue
c        sup=sum
c      endif
c     acummulate plant uptake of N and P
c      un1=un1+sunn
      up1=up1+sup
c     a new variable un3 added for debugging purposes
      suno3=suno3+uno3
c     update remaining no3 and labile p in the soil and ouput for debugging
      do 8 k=1,ir
c         wno3(k)=wno3(k)-un(k)
         ap(k)=ap(k)-un(k)
c        write(310,2133) jd,k,ap(k),rwt(k),wno3(k),un(k),up(k)
c2133 format (i3,1x,i3,1x,5(f8.4,1x))
   8  continue
      tno3=tno3-sunn+rmnr
      tap=tap-sup+wmp
      j=j-1
c     write(39,2000)jd,cnt,un1,un2,uno3,suno3,sunn,vt,dm,ddm
  45  continue
      return
      end