#include "definetab.h"
      subroutine ptend(k,km,solFe,remin,rkn,rks,rkf,bn,bf,cnratio,
     *  cchlratio,wsh,wshc,drate,rlamz,greff,dratez1,dratez2,
     *  regen,
#if PCOTAB
     *  pco2tab,
#endif
     *  P,H,atmfe,pnoice,
     *  tfac,fescav,rikd,rmuplsr,gcmax,
     *  tzoo,tirrq,atmco2,wspd,wspdvar,slp,T,S,
     *  P_tend,gronfix,ws,Hpst,pco2,cflx,pp2)
 
!  Computes tendencies of biological particles.  Outputs as 
!  units/h.
!  P(1)  = nitrate (uM)
!  P(2)  = ammonium (uM)
!  P(3)  = silica (uM)
!  P(4)  = iron (nM)
!  P(5)  = diatoms (mg chl m-3)
!  P(6)  = chlorophytes (mg chl m-3)
!  P(7)  = cyanobacteria (mg chl m-3)
!  P(8)  = coccolithophores (mg chl m-3)
!  P(9)  = herbivores (mg chl m-3)
!  P(10) = nitrate/carbon detritus (ugC l-1)
!  P(11) = silica detritus (uM)
!  P(12) = iron detritus (nM)
!  P(13) = DOC (uM)
!  P(14) = DIC (uM)
 
#include "definebio.h"
      real remin(ndet)      !detrital remineralization rate /hr
      real rkn(nchl),rks(nchl) !half-sat constants for N and S (uM)
      real rkf(nchl)        !half-sat constant for Fe (nM)
      real dphy(ntyp)       !phytoplankton death rate * P
      real zoo(ntyp)        !ingestion fraction by phyto group
      real rmu4(nchl),rmu3(nchl) !uptake of nutrients
      real wsh(ntyp)        !particle sinking rate m/h
      real P(ntyp)
      real H,tfac,fescav,wshc
      real atmfe,pnoice,tzoo
      real rikd(nchl)
      real rmuplsr(nchl)
      real gcmax,tirrq
      real T,S
      real P_tend(ntyp)
      real ws(ntyp)
      real Hpst
      real pp2(nchl)
      real gro(nchl)
      real gronfix
      real viscfac
      real, parameter :: Rm=1.20/24.0   !max zoopl. growth rate/hr
!               increase to account for excretion and respiration
      real, parameter :: excp=0.05      !excretion of DOC by phyto %growth
      real, parameter :: resp=0.05      !respiration of DIC by phyto %growth
      real, parameter :: phygross=1.0-(excp+resp)  !factor to derive net PP
      real, parameter :: uMtomgm3=12.0  !conversion uM to mg/m3
#if PCOTAB
      integer, parameter :: it0inc=1,nt0=80/it0inc
      integer, parameter :: isalinc=1,nsal=20
      integer, parameter :: idicinc=2,ndic=(650+idicinc)/idicinc
      integer, parameter :: itainc=2,nta=(500+itainc)/itainc
      real*4 pco2tab(nt0,nsal,ndic,nta)
#endif
      real mgchltouMC
 
!  Initialize
!      P_tend = 0.0
      gro = 0.0
      gronfix = 0.0
      pp2 = 0.0
      bs = 2.0*bn
      mgchltouMC = cchlratio/uMtomgm3
 
!  Start Model Space Loop
!   Iron + atm iron: disperse in layer and convert to nM
      if (k .eq. 1)then
       P_tend(4) = (atmFe*solFe)/H*0.001
      endif
!   Grazing/regeneration of ammonium
      ptot = 0.0
      do n = nnut+1,npe
       ptot = ptot + P(n)
      enddo
      ptot = max(ptot,1.0E-36)
      Pzoo = P(npe+1)
      gzoo = tzoo*Rm*(1.0-exp(-rlamz*ptot))*Pzoo
      dzoo1 = dratez1*Pzoo
      dzoo2 = dratez2*Pzoo*Pzoo
      P_tend(npe+1) = ((1.0-greff)*gzoo-dzoo1-dzoo2)
     *                      *pnoice
      do n = nnut+1,npe
!      fraction of grazing for this group
       zoo(n) = gzoo*P(n)/ptot
       dphy(n) = drate*P(n)
       P_tend(n) = (-zoo(n)-dphy(n))*pnoice
      enddo
      exc = greff*gzoo
      P_tend(1) = P_tend(1) + tfac*remin(1)*P(nds)/cnratio
     *                  *pnoice
      P_tend(2) = bn*exc*pnoice
     *                  + bn*regen*dzoo2*pnoice
      P_tend(3) = tfac*remin(2)*P(nds+1)*pnoice
      P_tend(4) = P_tend(4)
     *                  + bf*exc*pnoice
     *                  + bf*regen*dzoo2*pnoice
     *                  + tfac*remin(3)*P(nde)*pnoice
     *                  - fescav
      P_tend(ncs) = exc*mgchltouMC*pnoice
     *                  + regen*dzoo2*mgchltouMC*pnoice
      dphyt = 0.0
      do n = nnut+1,npe
       dphyt = dphyt + dphy(n)
      enddo
!  1st detrital fraction is carbon
      P_tend(nds) = (dphyt+dzoo1)*cchlratio*pnoice
     *                  + (1.0-regen)*dzoo2*cchlratio*pnoice
     *                  - tfac*remin(1)*P(nds)*pnoice
!  2nd detrital fraction is silica
      P_tend(nds+1) = bs*dphy(nnut+1)*pnoice
     *                  + bs*zoo(nnut+1)*pnoice
     *                  - tfac*remin(2)*P(nds+1)*pnoice
!  3rd detrital fraction is iron
      P_tend(nde) = bf*(dphyt+dzoo1)*pnoice
     *                  + bf*(1.0-regen)*dzoo2*pnoice
     *                  - tfac*remin(3)*P(nde)*pnoice
     *                  + 0.1*fescav
      
!   Day: Grow
      if (tirrq .gt. 0.0)then
       tirrqice = tirrq*0.01  !reduce light in ice by half
!   Light-regulated growth
#if NCHL_DEFINED > 0
!   Diatoms
       n = 1
!    Nutrient-regulated growth; Michaelis-Menton uptake kinetics
       rnut2 = P(2)/(rkn(n)+P(2))     !ammonium
       tnit = P(1)/(rkn(n)+P(1))      !nitrate
       tmp = 1.0 - rnut2
!    Enforce preferential utilization of ammonium
       rnut1 = min(tmp,tnit)
       rmmn = rnut1 + rnut2
       framm = rnut2/rmmn
       rmml = tirrq/(tirrq+0.5*rikd(n))
       rmmlice = tirrqice/(tirrqice+0.5*rikd(n))
!    Silica
       rmms = P(3)/(rks(n)+P(3)) 
       rmmf = P(4)/(rkf(n)+P(4))      !iron
       rlim = min(rmml,rmmn,rmms,rmmf)
       rlimice = min(rmmlice,rmmn,rmms,rmmf)
       grate = rmuplsr(n)*rlim*pnoice
     *          + rmuplsr(n)*rlimice*(1.0-pnoice)
       rmu4(n) = grate*framm
       rmu3(n) = grate*(1.0-framm)
       gro(n) = grate*P(n+nnut)
       P_tend(n+nnut) = P_tend(n+nnut) + gro(n)
!    Net primary production
       pp2(n) = gro(n)*phygross*H*cchlratio
#endif
#if NCHL_DEFINED > 1
!   Chlorophytes
       n = 2
!    Nutrient-regulated growth; Michaelis-Menton uptake kinetics
       rnut2 = P(2)/(rkn(n)+P(2))     !ammonium
       tnit = P(1)/(rkn(n)+P(1))      !nitrate
       tmp = 1.0 - rnut2
!    Enforce preferential utilization of ammonium
       rnut1 = min(tmp,tnit)
       rmmn = rnut1 + rnut2
       framm = rnut2/rmmn
       rmml = tirrq/(tirrq+0.5*rikd(n))
       rmmlice = tirrqice/(tirrqice+0.5*rikd(n))
       rmmf = P(4)/(rkf(n)+P(4))      !iron
       rlim = min(rmml,rmmn,rmmf)
       rlimice = min(rmmlice,rmmn,rmmf)
       grate = rmuplsr(n)*rlim*pnoice
     *          + rmuplsr(n)*rlimice*(1.0-pnoice)
       rmu4(n) = grate*framm
       rmu3(n) = grate*(1.0-framm)
       gro(n) = grate*P(n+nnut)
       P_tend(n+nnut) = P_tend(n+nnut) + gro(n)
!    Net primary production
       pp2(n) = gro(n)*phygross*H*cchlratio
#endif
#if NCHL_DEFINED > 2
!   Cyanobacteria
       n = 3
!    Nutrient-regulated growth; Michaelis-Menton uptake kinetics
       rnut2 = P(2)/(rkn(n)+P(2))     !ammonium
       tnit = P(1)/(rkn(n)+P(1))      !nitrate
       tmp = 1.0 - rnut2
!    Enforce preferential utilization of ammonium
       rnut1 = min(tmp,tnit)
       rmmn = rnut1 + rnut2
       framm = rnut2/rmmn
       rmml = tirrq/(tirrq+0.5*rikd(n))
       rmmf = P(4)/(rkf(n)+P(4))      !iron
       rlim = min(rmml,rmmn,rmmf)
       rlimnfix = min(rmml,rmmf)         !limitation for N2 fixation
       rlimrkn = min(rmml,rkn(n),rmmf)   !limitation at kn
       grate = rmuplsr(n)*rlim*pnoice
       rmu4(n) = grate*framm
       rmu3(n) = grate*(1.0-framm)
       rfix = 0.25*exp(-(75.0*P(n+nnut)))
       rfix = max(rfix,0.0)
!       rfix = min(rfix,0.2)
       gratenfix1 = rmuplsr(n)*rlimnfix*rfix  !N fix
       graterkn = rmuplsr(n)*rlimrkn
       gratenlim = graterkn - grate
       gratenfix = min(gratenlim,gratenfix1)  !N fix cannot exceed kn
       gratenfix = max(gratenfix,0.0)*pnoice
       gron = grate*P(n+nnut)
       gronfix = gratenfix*P(n+nnut)
       gro(n) = gron + gronfix
       P_tend(n+nnut) = P_tend(n+nnut) + gro(n)
!    Net primary production
       pp2(n) = gro(n)*phygross*H*cchlratio
#endif
#if NCHL_DEFINED > 3
!   Coccolithophores
       n = 4
!    Nutrient-regulated growth; Michaelis-Menton uptake kinetics
       rnut2 = P(2)/(rkn(n)+P(2))     !ammonium
       tnit = P(1)/(rkn(n)+P(1))      !nitrate
       tmp = 1.0 - rnut2
!    Enforce preferential utilization of ammonium
       rnut1 = min(tmp,tnit)
       rmmn = rnut1 + rnut2
       framm = rnut2/rmmn
       rmml = tirrq/(tirrq+0.5*rikd(n))
       rmmf = P(4)/(rkf(n)+P(4))      !iron
       rlim = min(rmml,rmmn,rmmf)
       grate = rmuplsr(n)*rlim*pnoice
       rmu4(n) = grate*framm
       rmu3(n) = grate*(1.0-framm)
       gro(n) = grate*P(n+nnut)
       P_tend(n+nnut) = P_tend(n+nnut) + gro(n)
       gcmax = max(gcmax,grate)
!    Net primary production
       pp2(n) = gro(n)*phygross*H*cchlratio
#endif
!#if NCHL_DEFINED > 4
!!   Dinoflagellates
!       n = 5
!!    Nutrient-regulated growth; Michaelis-Menton uptake kinetics
!       rnut2 = P(k,2)/(rkn(n)+P(k,2))     !ammonium
!       tnit = P(k,1)/(rkn(n)+P(k,1))      !nitrate
!       tmp = 1.0 - rnut2
!!    Enforce preferential utilization of ammonium
!       rnut1 = min(tmp,tnit)
!       rmmn = rnut1 + rnut2
!       framm = rnut2/rmmn
!       rmml = tirrq(k)/(tirrq(k)+0.5*rikd(k,n))
!       rmmf = P(k,4)/(rkf(n)+P(k,4))      !iron
!       rlim = min(rmml,rmmn,rmmf)
!       grate = rmuplsr(k,n)*rlim*pnoice
!       rmu4(n) = grate*framm
!       rmu3(n) = grate*(1.0-framm)
!       gro(k,n) = grate*P(k,n+nnut)
!       P_tend(k,n+nnut) = P_tend(k,n+nnut) + gro(k,n)
!    Net primary production
!        pp2(n) = gro(n)*phygross*H*cchlratio
!#endif
!   Nutrient uptake
       upn = 0.0
       upa = 0.0
       upf = 0.0
       do n = 1,nchl
        upn = upn + bn*(rmu3(n)*P(nnut+n))
        upa = upa + bn*(rmu4(n)*P(nnut+n))
        upf = upf + bf*gro(n)
       enddo
       P_tend(1) = P_tend(1) - upn
       P_tend(2) = P_tend(2) - upa
       ups = bs*gro(1)
       P_tend(3) = P_tend(3) - ups
       P_tend(4) = P_tend(4) - upf
      endif
 
      call carbon(k,cnratio,cchlratio,bn,bf,remin,P(:),tfac,
     * pnoice,tzoo,
#if PCOTAB
     * pco2tab,
#endif
     * tirrq,gro(:),atmco2,wspd,wspdvar,slp,T,S,H,
     * P_tend(:),pco2,cflx)
 
!  Sinking rate temperature (viscosity) dependence (also convert
!   to /hr)
      viscfac = 0.451 + 0.0178*T
      do n = nnut+1,npe
       ws(n) = wsh(n)*viscfac*pnoice
      enddo
      n = nnut+4    !cocco's
      ws(n) = wshc*viscfac*pnoice
      do n = nds,nde
       ws(n) = wsh(n)*viscfac*pnoice
      enddo
 
!  Save method for hard boundary condition (no flux)
!      srate = 0.0 - wsh(n)*P(i,k-1,m,n)
 
!  Save array for "past" layer thickness
      Hpst = H
 
      return
      end