      subroutine npzd_src (bioin, ntsb, tsb, gl, bct, impo, dzt
     &,                    dayfrac, wwd, rkw, nud, bioout, expoout
     &,                    grazout, morpout, morzout, graz_Det_out
     &,                    graz_Zout, nudop, nudon
# if defined O_save_npzd
     &,                    nppout, npp_dopout, morptout, remiout
     &,                    excrout
#  if defined O_npzd_nitrogen
     &,                    npp_Dout, npp_D_dopout, graz_Dout, morp_Dout
     &,                    morpt_Dout, nfixout
#  endif
#  if defined O_npzd_extra_diagnostics
     &,                    avej_out, avej_D_out, gmax_out, no3P_out
     &,                    po4P_out, po4_D_out
#  end if
# endif
#  if defined O_npzd_fe_limitation
     &,                    felimit, felimit_D
#  endif
# if defined O_zoop_graz_upper_temp_limit
     &,                    bctz
# endif
# if defined O_npzd_nitrogen_15
     &,                    rn15impo, rn15expoout
# endif
# if defined O_carbon_13
     &,                    rc13impo, rc13expoout, ac13b
# endif
     &                     ) 

#if defined O_mom && defined O_npzd 
!=======================================================================
!     computes source terms of the NPZD model
!     initial version of code adapted from Xavier Giraud:
!     Giraud et al. 2000, J Mar Res, 58, 609-630
!     original model reference:
!     Oeschlies and Garcon 1999, Global Biogeochem. Cycles 13, 135-160
!     Schmittner et al. 2005,  Global Biogeochem. Cycles 19, GB3004,
!     doi:10.1029/2004GB002283.
!     Schmittner et al. 2008, Global Biogeochem. Cycles 22, GB1013
!
!     This version was modified by David Keller and corrects the zooplankton
!     grazing formulation.  Note that zooplankton are now allowed to graze
!     on themselves and detritus, in addition to phyt. and diazotrophs.
!     The calculation of light has also been corrected.
!
!     Note that nutrient now represents phosphate

!     input variables:

!       bioin(1:5) = N,DON,P,Z,D [mmol m-3]
!       bioin(6)   = nitrate [mmol m-3]
!       bioin(7)   = DON [mmol m-3]
!       bioin(8)   = diazotrophs [mmol m-3]

!       gl         = 2.*light at top of grid box
!       ntsb       = number of time steps
!       tsb        = time step [s]
!       bct        = bbio**(cbio*temperature)
!       impo       = import of detritus from above [mmol m-3]
!       dzt        = depth of grid box [cm]
!       dayfrac    = day length (fraction: 0 < dayfrac < 1)
!       wwd        = sinking speed of detritus/dzt
!       rkw        = reciprical of kw*dzt(k)
!       nud        = remineralisation rate of detritus [s-1]
!       nudop      = remineralisation rate of dop [s-1]
!       nudon      = remineralisation rate of don [s-1]
!     output variables:

!       bioout     = change from bioin [mmol m-3]
!       nppout     = net primary production [mmol m-3]
!       npp_dopout = net primary production from DOP uptake [mmol m-3]
!       grazout    = grazing [mmol m-3]
!       morpout    = quadratic mortality of phytoplankton [mmol m-3]
!       morptout   = specific mortality of phytoplankton [mmol m-3]
!       morzout    = mortality of zooplankton [mmol m-3]
!       remiout    = remineralisation [mmol m-3]
!       excrout    = excretion [mmol m-3]
!       expoout    = detrital export [mmol m-3]
!       npp_Dout   = NPP of diazotrophs
!       npp_D_dopout= NPP from DOP uptake of diazotrophs [mmol m-3]
!       graz_Dout  = grazing of diazotrophs
!       morp_Dout  = mortality of diazotrophs
!       nfixout    = rate of N2 fixation
!       graz_Det_out = grazing of detritus
!       graz_Zout   = grazing on othe zooplankton
!       avej_out    = light-depend phyt. growth rate
!       avej_D_out  = light-depend Diaz growth rate
!       gmax_out    = temp-depend. zoo growth rate
!       no3P_out    = no3 depend. phyt growth rate
!       po4P_out    = po4 depend. phyt growth rate
!       po4_D_out   = po4 depend. Diaz growth rate

!      New grazing formulation variables and parameters

!       The following terms determine ingestion according to a
!       a Holling II curve (i.e. Michaelis Menten):
!
!       Ingestion = max_graz_rate * (Ft/(Ft + kzoo))
!
!       where Ft is the weighted measure of the total food available
!       and equals the sum of the different prey types times the
!       preference of Z for that type of prey 
!
!       zprefP   = Z preference for P
!       zprefD   = Z preference for Diaz
!       zprefDet = Z preference for detritus
!       zprefZ   = Z preference for other Z
!       kzoo = half saturation coefficienct for Z ingestion mmol N m-3
!       ing_P    = zooplankton ingestion of phytoplankon
!       ing_D    = zooplankton ingestion of diazotrophs
!       ing_Det  = zooplankton ingestion of detritus
!       ing_Z    = zooplankton ingestion of other zooplankton
!       thetaZ   = Michaelis-Menten denominator
!
!       felimit = Fe limitation parameter
!       felmit_D = Fe limitation parameter for diazotrophs
!
!=======================================================================

      implicit none

      integer n, ntsb

      real gl, f1, biopo4, biophyt, biozoop, biodetr, jmax, u_P, g_P
      real morp, morpt, morz, remi, excr, expo, impo, nppout, grazout
      real morpout, morptout, morzout, remiout, excrout, expoout, tsb
      real avej_out, avej_D_out, gmax_out, no3P_out, po4P_out, po4_D_out
      real dzt, po4flag, phytflag, zoopflag, detrflag, wwd, rkw, gd
      real nupt, nud, biono3, u_D,npp_D, npp_Dout, no3flag, biodiaz, bct
      real diazflag, g_D,graz_D, morpt_D, jmax_D, gd_D, avej_D, no3upt_D
      real morpt_Dout, graz_Dout, nfixout, biop2, u1, u2, phi1, phi2
      real avej, graz_Det_out, graz_Zout, thetaZ, ing_P, ing_D, dayfrac
      real ing_Det, ing_Z, g_Z, g_Det, graz_Z, graz_Det, gmax
      real no3P, po4P, po4_D, felimit, bctz, felimit_D, nudop, nudon
      real excrdiaz, biodop, biodon, dopflag, donflag, recy
      real recy_p, npp, graz, biodin15, biodon15, biophytn15, biozoopn15
      real biodetrn15, biodiazn15, bassim, fcassim, bexcr, fcexcr
      real bnfix, fcnfix, rn15impo, rn15expoout, dig, dig_P, dig_Z
      real dig_Det, dig_D, excr_P, excr_Z, excr_Det, excr_D, sf, sf_P
      real sf_Z, sf_Det, sf_D, nr_excr_D, uno3, rno3, rzoop
      real rtdin15, rtphytn15, rtzoopn15, rtdetrn15, rtdiazn15, rtdon15
      real din15flag, dopupt_flag, dopupt_D_flag, npp_dopout
      real don15flag, phytn15flag, zoopn15flag, detrn15flag, diazn15flag
      real limP, limP_D, uptdop_D, dopupt_D, morp_D, nupt_D, morp_Dout
      real udon, rdon, brecy, fcrecy, biodic, uptdop, dopupt
      real biodic13, biophytc13, biozoopc13, biodetrc13, biodiazc13
      real dic13flag, phytc13flag, zoopc13flag, detrc13flag, diazc13flag
      real biodoc13, doc13flag, biodoc, npp_D_dopout
      real rdic13, rtphytc13, rtzoopc13, rtdetrc13, rtdiazc13, ac13b
      real rc13impo, rc13expoout, dicflag, bc13npp, rtdic13, fcnpp
      real rtdoc13

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "calendar.h"
      include "npzd.h"

      real bioin(ntnpzd), bioout(ntnpzd)

      bioout(:) = 0.0
      biopo4 = bioin(1)
      biodop = bioin(2)
      biophyt = bioin(3)
      biozoop = bioin(4)
      biodetr = bioin(5)
# if defined O_npzd_nitrogen
      biodic = bioin(6)
      biono3 = bioin(7)
      biodon = bioin(8)
      biodiaz = bioin(9)
#  if defined O_npzd_nitrogen_15
      biodin15 = bioin(10)
      biodon15 = bioin(11)
      biophytn15 = bioin(12)
      biozoopn15 = bioin(13)
      biodetrn15 = bioin(14)
      biodiazn15 = bioin(15)
#  endif
# endif
# if defined O_carbon_13
      biodic13 = bioin(16)
      biodoc13 = bioin(17)
      biophytc13 = bioin(18)
      biozoopc13 = bioin(19)
      biodetrc13 = bioin(20)
#  if defined O_npzd_nitrogen
      biodiazc13 = bioin(21)
#  endif
# endif
!       flags prevent negative values by setting outgoing fluxes to
!       zero if tracers are lower than trcmin
      po4flag = 0.5 + sign(0.5,biopo4 - trcmin)
      dopflag = 0.5 + sign(0.5,biodop - trcmin)
      phytflag = 0.5 + sign(0.5,biophyt - trcmin)
      zoopflag = 0.5 + sign(0.5,biozoop - trcmin)
      detrflag = 0.5 + sign(0.5,biodetr - trcmin)
# if defined O_npzd_nitrogen
      no3flag = 0.5 + sign(0.5,biono3 - trcmin)
      donflag = 0.5 + sign(0.5,biodon - trcmin)
      diazflag = 0.5 + sign(0.5,biodiaz - trcmin)
#  if defined O_npzd_nitrogen_15
      din15flag = 0.5 + sign(0.5, biodin15 - trcmin)
      don15flag = 0.5 + sign(0.5, biodon15 - trcmin)
      phytn15flag = 0.5 + sign(0.5, biophytn15 - trcmin)
      zoopn15flag = 0.5 + sign(0.5, biozoopn15 - trcmin)
      detrn15flag = 0.5 + sign(0.5, biodetrn15 - trcmin)
      diazn15flag = 0.5 + sign(0.5, biodiazn15 - trcmin)
#  else
      din15flag = 1.0
      don15flag = 1.0
      phytn15flag = 1.0
      zoopn15flag = 1.0
      detrn15flag = 1.0
      diazn15flag = 1.0
#  endif
# endif
# if defined O_carbon_13
      dic13flag = 0.5 + sign(0.5,biodic13 - trcmin)
      doc13flag = 0.5 + sign(0.5,biodoc13 - trcmin)
      phytc13flag = 0.5 + sign(0.5,biophytc13 - trcmin)
      zoopc13flag = 0.5 + sign(0.5,biozoopc13 - trcmin)
      detrc13flag = 0.5 + sign(0.5,biodetrc13 - trcmin)
      diazc13flag = 0.5 + sign(0.5,biodiazc13 - trcmin)
# endif

      bioin(:) = max(bioin(:), trcmin)
      biopo4 = max(biopo4, trcmin)
      biodop = max(biodop, trcmin)
      biophyt = max(biophyt, trcmin)
      biozoop = max(biozoop, trcmin)
      biodetr = max(biodetr, trcmin)
# if defined O_npzd_nitrogen
      biodic = max(biodic, trcmin)
      biono3 = max(biono3, trcmin)
      biodon = max(biodon, trcmin)
      biodiaz = max(biodiaz, trcmin)
#  if defined O_npzd_nitrogen_15
      biodin15 = max(biodin15, trcmin)
      biodon15 = max(biodon15, trcmin)
      biophytn15 = max(biophytn15, trcmin)
      biozoopn15 = max(biozoopn15, trcmin)
      biodetrn15 = max(biodetrn15, trcmin)
      biodiazn15 = max(biodiazn15, trcmin)
#  endif
# endif
# if defined O_carbon_13
      biodic13 = max(biodic13, trcmin)
      biodoc13 = max(biodoc13, trcmin)
      biophytc13 = max(biophytc13, trcmin)
      biozoopc13 = max(biozoopc13, trcmin)
      biodetrc13 = max(biodetrc13, trcmin)
      biodiazc13 = max(biodiazc13, trcmin)
# endif
!     photosynthesis after Evans & Parslow (1985)
!     notation as in JGOFS report No. 23 p. 6
# if defined O_npzd_cdom_attenuation
      f1 = exp(1.2*(-kw - kc*(biophyt+biodiaz))*dzt)
# else
      f1 = exp((-kw - kc*(biophyt+biodiaz))*dzt)
# endif
# if defined O_npzd_fe_limitation
! In the following "felimit" is determined by an iron mask and
! is used to limit phytoplankton growth in HNLC regions
      jmax = abio*bct*felimit
# else
      jmax = abio*bct
# endif
      gd = jmax*dayfrac
      u1 = max(gl/gd,1.e-6)
      u2 = u1*f1
!     for the following approximation ensure that u1 < 20
      phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1
      phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2
# if defined O_npzd_cdom_attenuation
      avej = gd*(phi1 - phi2)/(1.2*(kw+kc*(biophyt+biodiaz))*dzt)
# else
      avej = gd*(phi1 - phi2)/((kw+kc*(biophyt+biodiaz))*dzt)
# endif
! Make the max grazing rate a function of temperature
# if defined O_zoop_graz_upper_temp_limit
      gmax = gbio*bctz
! Note bctz, which sets an upper limit on the effects of temp on the
! grazing rate, is set in tracers.F
# else
      gmax = gbio*bct
# endif
# if defined O_npzd_nitrogen
# if defined O_npzd_fe_limitation
      jmax_D = max(0.,abio*bct*felimit_D)*jdiar
# else
      jmax_D = max(0.,abio*bct)*jdiar
# endif
      gd_D = max(1.e-14,jmax_D*dayfrac)
      u1 = max(gl/gd_D,1.e-6)
      u2 = u1*f1
!     for the following approximation ensure that u1 < 20
      phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1
      phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2
# if defined O_npzd_cdom_attenuation
      avej_D = gd_D*(phi1 - phi2)/(1.2*(kw+kc*(biophyt+biodiaz))*dzt)
#else
      avej_D = gd_D*(phi1 - phi2)/((kw+kc*(biophyt+biodiaz))*dzt)
# endif
# endif
# if defined O_npzd_nitrogen
!     check grazing preferences = 1 for N case
        IF ((zprefP + zprefDet + zprefZ + zprefD).ne.1) THEN
           zprefP = 0.3
           zprefZ = 0.3
           zprefDet = 0.3
           zprefD = 0.1
        END IF
# else
!      check that grazing preferences = 1 for no N case
        IF ((zprefP + zprefDet + zprefZ).ne.1) THEN
           zprefP = 0.35
           zprefZ = 0.35
           zprefDet = 0.30
        END IF
# endif
      nupt = nupt0*bct
      nupt_D = nupt0_D*bct

      expoout = 0.0
      grazout = 0.0
      morpout = 0.0
      morzout = 0.0
      graz_Det_out = 0.0
      graz_Zout = 0.0
# if defined O_npzd_nitrogen_15
      rn15expoout = 0.0
# endif
# if defined O_carbon_13
      rc13expoout = 0.0
# endif
# if defined O_save_npzd
      nppout = 0.0
      npp_dopout = 0.0
      morptout = 0.0
      remiout = 0.0
      excrout = 0.0
#  if defined O_npzd_nitrogen
      npp_Dout = 0.0
      npp_D_dopout = 0.0
      graz_Dout = 0.0
      morpt_Dout = 0.0
      morp_Dout = 0.0
      nfixout = 0.0
#  endif
# if defined O_npzd_extra_diagnostics
      avej_out = 0.0
      avej_D_out = 0.0
      gmax_out = 0.0
      no3P_out = 0.0
      po4P_out = 0.0
      po4_D_out = 0.0
# endif
# endif

      do n=1,ntsb

!       growth rate of phytoplankton
        if (hdop*biodop/(k1p + biodop) .gt. biopo4/(k1p + biopo4)) then
           dopupt_flag = 1.
           limP = biodop/(k1p + biodop)
           u_P = min(avej, jmax*limP*hdop)
        else
           dopupt_flag = 0.
           limP = biopo4/(k1p + biopo4)
           u_P = min(avej, jmax*limP)
        endif
        po4P = jmax*biopo4/(k1p + biopo4) 

# if defined O_npzd_nitrogen
!       nitrate limitation
        u_P = min(u_P, jmax*biono3/(k1n + biono3))
        no3P = jmax*biono3/(k1n + biono3)
!       growth rate of diazotrophs smaller than other phytoplankton and
!       not nitrate limited
        if (hdop*biodop/(k1p + biodop) .gt. biopo4/(k1p + biopo4)) then
           dopupt_D_flag = 1.
           limP_D = biodop/(k1p + biodop)
           u_D = min(avej_D, jmax_D*limp_D*hdop)
        else
           dopupt_D_flag = 0.
           limP_D = biopo4/(k1p + biopo4)
           u_D = min(avej_D, jmax_D*limp_D)
        endif
        po4_D = jmax_D*biopo4/(k1p + biopo4)
!       Set the grazing coefficients for the N case
        thetaZ = zprefP*biophyt + zprefDet*biodetr + zprefZ*biozoop
     &       + zprefD*biodiaz + kzoo
        ing_P = zprefP/thetaZ
        ing_Det = zprefDet/thetaZ
        ing_Z = zprefZ/thetaZ
        ing_D = zprefD/thetaZ
# else
!       If "else" then set the grazing coefficients for the no N case
!       note kzoo is in terms of N so convert to P
        thetaZ = zprefP*biophyt + zprefDet*biodetr + zprefZ*biozoop 
     &       +kzoo*redptn
        ing_P = zprefP/thetaZ
        ing_Det = zprefDet/thetaZ
        ing_Z = zprefZ/thetaZ
# endif
        npp = u_P*biophyt
        dopupt = npp*dopupt_flag
# if defined O_npzd_nitrogen
        npp_D = max(0.,u_D*biodiaz)
!       grazing on diazotrophs
        g_D = gmax*ing_D*biodiaz
        graz_D = g_D*biozoop
        morpt_D = nupt_D*biodiaz ! linear mortality
        morp_D = nup_D*biodiaz
        no3upt_D = biono3/(k1n + biono3)*npp_D ! nitrate uptake
        dopupt_D = npp_D*dopupt_D_flag
# endif
!       grazing on P
        g_P = gmax*ing_P*biophyt
        graz = g_P*biozoop
!       grazing on Z
        g_Z = gmax*ing_Z*biozoop
        graz_Z = g_Z*biozoop
!       grazing on Detritus
        g_Det = gmax*ing_Det*biodetr
        graz_Det = g_Det*biozoop
!
        morp = nup*biophyt
        morpt = nupt*biophyt
        recy = nudon*bct*biodon
        recy_p = nudop*bct*biodop
        morz = nuz*biozoop*biozoop
        remi = nud*bct*biodetr
        expo = wwd*biodetr

        graz = graz*phytflag*zoopflag*phytn15flag*zoopn15flag
        graz_Z = graz_Z*zoopflag*zoopn15flag
        graz_Det = graz_Det*detrflag*zoopflag*detrn15flag*zoopn15flag
        morp = morp*phytflag*phytn15flag
        morpt = morpt*phytflag*phytn15flag
        morz = morz*zoopflag*zoopn15flag
        remi = remi*detrflag*detrn15flag
        expo = expo*detrflag*detrn15flag
        recy_p = recy_p*dopflag
# if defined O_npzd_nitrogen
        if (dopupt_flag .eq. 1) then
           npp = npp*dopflag*no3flag*din15flag
        else
           npp = npp*po4flag*no3flag*din15flag
        endif
        if (dopupt_D_flag .eq. 1) then
           npp_D = npp_D*dopflag
        else
           npp_D = npp_D*po4flag
        endif
        graz_D = graz_D*diazflag*zoopflag*diazn15flag*zoopn15flag
        morpt_D = morpt_D*diazflag*diazn15flag
        morp_D = morp_D*diazflag*diazn15flag
        no3upt_D = no3upt_D*no3flag*din15flag
        recy = recy*donflag*don15flag
# else
        npp = npp*po4flag
# endif

!   Excretion is the difference between
!   the zooplankton assimilation and growth efficiencies
# if defined O_npzd_nitrogen
c        excr_P = (gamma1-geZ)*graz
c        excr_Z = (gamma1-geZ)*graz_Z
c        excr_Det = (gamma1-geZ)*graz_Det
c        excr_D = (gamma1-geZ)*graz_D

        dig_P = gamma1*graz
        dig_Z = gamma1*graz_Z
        dig_Det = gamma1*graz_Det
        dig_D = gamma1*graz_D*(redntp/diazntp)
        dig = dig_P + dig_Z + dig_Det + dig_D

        excr_P = gamma1*(1-geZ)*graz
        excr_Z = gamma1*(1-geZ)*graz_Z
        excr_Det = gamma1*(1-geZ)*graz_Det
        excr_D = gamma1*(1-geZ)*graz_D*(redntp/diazntp)
        excr = excr_P + excr_Z + excr_Det + excr_D
        nr_excr_D = gamma1*graz_D*(1-(redntp/diazntp)) 
     &       + (1-gamma1)*graz_D*(1-(redntp/diazntp))

        sf_P = (1-gamma1)*graz
        sf_Z = (1-gamma1)*graz_Z
        sf_Det = (1-gamma1)*graz_Det
        sf_D = (1-gamma1)*graz_D*(redntp/diazntp)
        sf = sf_P + sf_Z + sf_Det + sf_D
# else
c        excr = (gamma1-geZ)*(graz+graz_Z+graz_Det)
#endif
# if defined O_npzd_nitrogen_15
!       calculate isotope parameters
!       See Somes et al., 2010, GBC for details/results
        uno3 = npp*tsb/biono3
        uno3 = min(uno3, 0.999)
        uno3 = max(uno3, trcmin)
        rno3 = biodin15/(biono3-biodin15)
        rno3 = min(rno3, 2*rn15std)
        rno3 = max(rno3, rn15std/2.)
        bassim = rno3 + 6.*(1-uno3)/uno3*log(1-uno3)*rno3/1000.
        fcassim = bassim/(1+bassim)

        udon = recy*tsb/biodon
        udon = min(udon, 0.999)
        udon = max(udon, trcmin)
        rdon = biodon15/(biodon-biodon15)
        rdon = min(rdon, 2*rn15std)
        rdon = max(rdon, rn15std/2.)
        brecy = rdon + 1.3*(1-udon)/udon*log(1-udon)*rdon/1000.
        fcrecy = brecy/(1+brecy)

        rzoop = biozoopn15/(biozoop-biozoopn15)
        rzoop = min(rzoop, 2.*rn15std)
        rzoop = max(rzoop, rn15std/2.)
        bexcr = rzoop - 4.*rzoop/1000.
        fcexcr = bexcr/(1+bexcr)

        bnfix = rn15std - 1.*rn15std/1000.
        fcnfix = bnfix/(1+bnfix)

        rtdin15 = biodin15/biono3
        rtdin15 = min(rtdin15, 2*rn15std/(1+rn15std))
        rtdin15 = max(rtdin15, rn15std/(1+rn15std)/2.)

        rtdon15 = biodon15/biodon
        rtdon15 = min(rtdon15, 2*rn15std/(1+rn15std))
        rtdon15 = max(rtdon15, rn15std/(1+rn15std)/2.)

        rtphytn15 = biophytn15/biophyt
        rtphytn15 = min(rtphytn15, 2.*rn15std/(1+rn15std))
        rtphytn15 = max(rtphytn15, rn15std/(1+rn15std)/2.)

        rtzoopn15 = biozoopn15/biozoop
        rtzoopn15 = min(rtzoopn15, 2.*rn15std/(1+rn15std))
        rtzoopn15 = max(rtzoopn15, rn15std/(1+rn15std)/2.)

        rtdetrn15 = biodetrn15/biodetr
        rtdetrn15 = min(rtdetrn15, 2.*rn15std/(1+rn15std))
        rtdetrn15 = max(rtdetrn15, rn15std/(1+rn15std)/2.)

        rtdiazn15 = biodiazn15/biodiaz
        rtdiazn15 = min(rtdiazn15, 2.*rn15std/(1+rn15std))
        rtdiazn15 = max(rtdiazn15, rn15std/(1+rn15std)/2.)
# endif
# if defined O_carbon_13
        rdic13 = biodic13/(biodic-biodic13)
        rdic13 = min(rdic13, 2.*rc13std)
        rdic13 = max(rdic13, 0.5*rc13std)
        bc13npp = ac13b*rdic13
        fcnpp = bc13npp/(1+bc13npp)
        
        rtdic13 = biodic13/biodic
        rtdic13 = min(rtdic13, 2*rc13std/(1+rc13std))
        rtdic13 = max(rtdic13, 0.5*rc13std/(1+rc13std))

        rtdoc13 = biodoc13/(biodon*redctn)
        rtdoc13 = min(rtdoc13, 2*rc13std/(1+rc13std))
        rtdoc13 = max(rtdoc13, 0.5*rc13std/(1+rc13std))

        rtphytc13 = biophytc13/(biophyt*redctn)
        rtphytc13 = min(rtphytc13, 2.*rc13std/(1+rc13std))
        rtphytc13 = max(rtphytc13, 0.5*rc13std/(1+rc13std))

        rtzoopc13 = biozoopc13/(biozoop*redctn)
        rtzoopc13 = min(rtzoopc13, 2.*rc13std/(1+rc13std))
        rtzoopc13 = max(rtzoopc13, 0.5*rc13std/(1+rc13std))

        rtdetrc13 = biodetrc13/(biodetr*redctn)
        rtdetrc13 = min(rtdetrc13, 2.*rc13std/(1+rc13std))
        rtdetrc13 = max(rtdetrc13, 0.5*rc13std/(1+rc13std))

        rtdiazc13 = biodiazc13/(biodiaz*redctn)
        rtdiazc13 = min(rtdiazc13, 2.*rc13std/(1+rc13std))
        rtdiazc13 = max(rtdiazc13, 0.5*rc13std/(1+rc13std))
# endif
# if defined O_npzd_nitrogen
!       nutrients equation
        biopo4 = biopo4 + tsb*(redptn*((1.-dfrt)*morpt + excr + remi 
     &       - (npp-dopupt)) + diazptn*(morpt_D - (npp_D-dopupt_D)) 
     &       + recy_p)
!       DOP equation
        biodop = biodop + tsb*(redptn*dfr*morp + redptn*dfrt*morpt 
     &       - redptn*dopupt - diazptn*dopupt_D - recy_p)
!       phytoplankton equation
        biophyt = biophyt + tsb*(npp - morp - graz - morpt)
!       zooplankton equation
        biozoop = biozoop + tsb*(dig - morz - graz_Z - excr)
!       detritus equation
        biodetr = biodetr + tsb*((1.-dfr)*morp + sf + morz - remi 
     &       - graz_Det - expo + impo + morp_D*(redntp/diazntp))
!       DIC equation
        biodic = biodic + tsb*redctn*((1.-dfrt)*morpt + excr + remi 
     &       - npp + morpt_D - npp_D + recy + nr_excr_D
     &       + morp_D*(1-(redntp/diazntp)))
!       nitrate (NO3) equation
        biono3 = biono3 + tsb*((1.-dfrt)*morpt + excr + morpt_D  
     &       + morp_D*(1-(redntp/diazntp)) + nr_excr_D + remi + recy 
     &       - npp - no3upt_D)
!       DON equation
        biodon = biodon + tsb*(dfr*morp + dfrt*morpt - recy)
!       diazotroph equation
        biodiaz = biodiaz + tsb*(npp_D - morp_D - morpt_D - graz_D)
# else
!       nutrients equation
        biopo4 = biopo4 + tsb*redptn*(remi + excr - npp + morpt)
!       phytoplankton equation
        biophyt = biophyt + tsb*(npp - morp - graz - morpt)
!       zooplankton equation
        biozoop = biozoop + tsb*(geZ*(graz + graz_Det + graz_Z) - morz
     &       - graz_Z)
!       detritus equation
        biodetr = biodetr + tsb*((1.-gamma1)*(graz + graz_Det + graz_Z) 
     &       + morp + morz - remi - graz_Det - expo + impo)
# endif
# if defined O_npzd_nitrogen_15
!       isotope equations
        biodin15 = biodin15 + tsb*(rtphytn15*(1.-dfrt)*morpt 
     &       + fcexcr*excr + rtdiazn15*morpt_D + rtdiazn15*nr_excr_D 
     &       + rtdiazn15*morp_D*(1-(redntp/diazntp)) + rtdetrn15*remi 
     &       + fcrecy*recy - fcassim*npp - fcassim*no3upt_D)

        biodon15 = biodon15 + tsb*(dfr*rtphytn15*morp 
     &       + dfrt*rtphytn15*morpt - fcrecy*recy)

        biophytn15 = biophytn15 + tsb*(fcassim*npp - rtphytn15*morp 
     &       - rtphytn15*graz - rtphytn15*morpt)

        biozoopn15 = biozoopn15 + tsb*(rtphytn15*dig_P + rtzoopn15*dig_Z 
     &       + rtdetrn15*dig_Det + rtdiazn15*dig_D - rtzoopn15*morz 
     &       - rtzoopn15*graz_Z - fcexcr*excr)

        biodetrn15 = biodetrn15 + tsb*(rtphytn15*(1.-dfr)*morp 
     &       + rtphytn15*sf_P + rtzoopn15*sf_Z + rtdetrn15*sf_Det 
     &       + rtdiazn15*sf_D + rtzoopn15*morz - rtdetrn15*remi 
     &       - rtdetrn15*graz_Det - rtdetrn15*expo + rn15impo*impo 
     &       + rtdiazn15*morp_D*(redntp/diazntp))

        biodiazn15 = biodiazn15 + tsb*(fcnfix*(npp_D-no3upt_D)
     &       + fcassim*no3upt_D - rtdiazn15*morp_D - rtdiazn15*graz_D 
     &       - rtdiazn15*morpt_D)
# endif
# if defined O_carbon_13
! !!!!!!!!!!!!!!!!!! isotope equations
        biodic13 = biodic13 + tsb*redctn*(rtphytc13*(1.-dfrt)*morpt 
     &       + rtzoopc13*excr + rtdiazc13*morpt_D + rtdiazc13*nr_excr_D 
     &       + rtdiazc13*morp_D*(1-(redntp/diazntp)) + rtdetrc13*remi 
     &       + rtdoc13*recy - fcnpp*npp - fcnpp*npp_D)

        biodoc13 = biodoc13 + tsb*redctn*(dfr*rtphytc13*morp 
     &       + rtphytc13*dfrt*morpt - rtdoc13*recy)

        biophytc13 = biophytc13 + tsb*redctn*(fcnpp*npp - rtphytc13*morp 
     &       - rtphytc13*graz - rtphytc13*morpt)

        biozoopc13 = biozoopc13 + tsb*redctn*(rtphytc13*dig_P 
     &       + rtzoopc13*dig_Z + rtdetrc13*dig_Det + rtdiazc13*dig_D 
     &       - rtzoopc13*morz - rtzoopc13*graz_Z - rtzoopc13*excr)

        biodetrc13 = biodetrc13 + tsb*redctn*(rtphytc13*(1.-dfr)*morp 
     &       + rtphytc13*sf_P + rtzoopc13*sf_Z + rtdetrc13*sf_Det 
     &       + rtdiazc13*sf_D + rtzoopc13*morz - rtdetrc13*remi 
     &       - rtdetrc13*graz_Det - rtdetrc13*expo + rc13impo*impo 
     &       + rtdiazc13*morp_D*(redntp/diazntp))

        biodiazc13 = biodiazc13 + tsb*redctn*(fcnpp*npp_D 
     &       - rtdiazc13*morp_D - rtdiazc13*graz_D - rtdiazc13*morpt_D)
# endif
        expoout = expoout + expo
# if defined O_npzd_nitrogen_15
        rn15expoout = rn15expoout + rtdetrn15
# endif
# if defined O_carbon_13
        rc13expoout = rc13expoout + rtdetrc13
#  endif
        grazout = grazout + graz
        morpout = morpout + morp
        morzout = morzout + morz
        graz_Det_out = graz_Det_out + graz_Det
        graz_Zout = graz_Zout + graz_Z
# if defined O_save_npzd
        nppout = nppout + npp
        npp_dopout = npp_dopout + npp*dopupt_flag
        morptout = morptout + morpt
        remiout = remiout + remi
        excrout = excrout + excr
#  if defined O_npzd_nitrogen
        npp_Dout = npp_Dout + npp_D
        npp_D_dopout = npp_D_dopout + npp_D*dopupt_D_flag
        graz_Dout = graz_Dout + graz_D
        morpt_Dout = morpt_Dout + morpt_D
        morp_Dout = morp_Dout + morp_D
        nfixout = nfixout + npp_D - no3upt_D
#  endif
# if defined O_npzd_extra_diagnostics
        avej_out = avej_out + avej
        avej_D_out = avej_D_out + avej_D
        gmax_out = gmax_out + gmax
        no3P_out = no3P_out + no3P
        po4P_out = po4P_out + po4P
        po4_D_out = po4_D_out + po4_D
# endif
# endif
!       flags prevent negative values by setting outgoing fluxes to
!       zero if tracers are lower than trcmin
      if (po4flag .eq. 1) po4flag = 0.5 + sign(0.5,biopo4 - trcmin)
      if (dopflag .eq. 1) dopflag = 0.5 + sign(0.5,biodop - trcmin)
      if (phytflag .eq. 1) phytflag = 0.5 + sign(0.5,biophyt - trcmin)
      if (zoopflag .eq. 1) zoopflag = 0.5 + sign(0.5,biozoop - trcmin)
      if (detrflag .eq. 1) detrflag = 0.5 + sign(0.5,biodetr - trcmin)
# if defined O_npzd_nitrogen
      if (no3flag .eq. 1) no3flag = 0.5 + sign(0.5,biono3 - trcmin)
      if (donflag .eq. 1) donflag = 0.5 + sign(0.5,biodon - trcmin)
      if (diazflag .eq. 1) diazflag = 0.5 + sign(0.5,biodiaz - trcmin)
#  if defined O_npzd_nitrogen_15
      if (din15flag .eq. 1) 
     &     din15flag = 0.5 + sign(0.5, biodin15 - trcmin)
      if (don15flag .eq. 1)
     &     don15flag = 0.5 + sign(0.5, biodon15 - trcmin)
      if (phytn15flag .eq. 1) 
     &     phytn15flag = 0.5 + sign(0.5, biophytn15 - trcmin)
      if (zoopn15flag .eq. 1)
     &     zoopn15flag = 0.5 + sign(0.5, biozoopn15 - trcmin)
      if (detrn15flag .eq. 1)
     &     detrn15flag = 0.5 + sign(0.5, biodetrn15 - trcmin)
      if (diazn15flag .eq. 1)
     &     diazn15flag = 0.5 + sign(0.5, biodiazn15 - trcmin)
#  endif
# endif
# if defined O_carbon_13
        if (dic13flag .eq. 1)
     &       dic13flag = 0.5 + sign(0.5,biodic13 - trcmin)
        if (doc13flag .eq. 1)
     &       doc13flag = 0.5 + sign(0.5,biodoc13 - trcmin)
        if (phytc13flag .eq. 1)
     &       phytc13flag = 0.5 + sign(0.5,biophytc13 - trcmin)
        if (zoopc13flag .eq. 1)
     &       zoopc13flag = 0.5 + sign(0.5,biozoopc13 - trcmin)
        if (detrc13flag .eq. 1)
     &       detrc13flag = 0.5 + sign(0.5,biodetrc13 - trcmin)
#  if defined O_npzd_nitrogen
        if (diazc13flag .eq. 1)
     &       diazc13flag = 0.5 + sign(0.5,biodiazc13 - trcmin)
#  endif
# endif
      enddo

      bioout(1) = biopo4 - bioin(1)
      bioout(2) = biodop - bioin(2)
      bioout(3) = biophyt - bioin(3)
      bioout(4) = biozoop - bioin(4)
      bioout(5) = biodetr - bioin(5)
# if defined O_npzd_nitrogen
      bioout(6) = biodic - bioin(6)
      bioout(7) = biono3 - bioin(7)
      bioout(8) = biodon - bioin(8)
      bioout(9) = biodiaz - bioin(9)
#   if defined O_npzd_nitrogen_15
      bioout(10) = biodin15 - bioin(10)
      bioout(11) = biodon15 - bioin(11)
      bioout(12) = biophytn15 - bioin(12)
      bioout(13) = biozoopn15 - bioin(13)
      bioout(14) = biodetrn15 - bioin(14)
      bioout(15) = biodiazn15 - bioin(15)
#   endif
# endif
# if defined O_carbon_13
      bioout(16) = biodic13 - bioin(16)
      bioout(17) = biodoc13 - bioin(17)
      bioout(18) = biophytc13 - bioin(18)
      bioout(19) = biozoopc13 - bioin(19)
      bioout(20) = biodetrc13 - bioin(20)
      bioout(21) = biodiazc13 - bioin(21)
# endif
#endif
      return
      end
