! source file: /home/aschmitt/UVic/c13_lightbug_fixed/epsbio5/kvSO_low_rfac.5/zq101.9/felim/popp/updates/npzd_src.F
      subroutine npzd_src (bioin, ntsb, tsb, gl, temp, impo, dzt
     &,                    dayfrac, wwd, rkw, nud, bioout, expoout
     &,                    grazout, morpout, morzout
     &,                    i, j, mi
     &,                    nppout, morptout, remiout, excrout
     &,                    npp_Dout, graz_Dout, morp_Dout, nfixout
     &,                    rc13impo, rc13expoout, co2aq
     &                     )

!=======================================================================
!     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.
!     This is the optimized model version of Schartau & Oschlies, 2003,
!     J. Mar. Res. 62, 765-793 and Oschlies & Schartau, in prep.
!     as described in Schmittner and Oschlies, in prep.
!     note that nutrient now represents phosphate
!     Andreas Schmittner Sept 2005 (aschmitt@coas.oregonstate.edu)

!     input variables:

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

!       gl         = 2.*light at top of grid box
!       ntsb       = number of time steps
!       tsb        = time step [s]
!       temp       = temperature [deg C]
!       impo       = import of detritus from above [mmol m-3]
!       rc13impo   = import of 13C/12C of detritus from above
!       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]
!       co2aq      = aquatic CO2 concentration (mol/m^3)

!     output variables:

!       bioout     = change from bioin [mmol m-3]
!       nppout     = net primary production [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]
!       rc13expoout= 13C/12C ratio of detrital export
!       npp_Dout   = NPP of diazotrophs
!       graz_Dout  = grazing of diazotrophs
!       morp_Dout  = mortality of diazotrophs
!       nfixout    = rate of N2 fixation
!=======================================================================
      implicit none

      include "param.h"
      include "calendar.h"
      include "npzd.h"

      integer n, ntsb, i, j, mi

      real bioin(ntnpzd), bioout(ntnpzd+1)
      real gl, f1, bion, biop, bioz, biod, jmax, u_P, g_P, npp, graz
      real morp, morpt, morz, remi, excr, expo, impo, nppout, grazout
      real morpout, morptout, morzout, remiout, excrout, expoout, tsb
      real dzt, nflag, pflag, zflag, dflag, wwd, rkw, gd, dayfrac, bct
      real nupt, nud, biono3, u_D,npp_D, npp_Dout, no3flag, biodiaz
      real diazflag, g_D,graz_D, morp_D, jmax_D, gd_D, avej_D, no3upt_D
      real morp_Dout, graz_Dout, nfixout, biop2, u1, u2, phi1, phi2
      real avej
      real biodic13, biophytc13, biozoopc13, biodetrc13, biodiazc13
      real dic13flag, phytc13flag, zoopc13flag, detrc13flag, diazc13flag
      real rdic13, rc13phyt, rc13zoop, rc13detr, rc13diaz
      real rc13impo, rc13expoout, biodic, dicflag, ac13bp, ac13bd
      real ac13_aq_POC_P, ac13_aq_POC_D, ac13_DIC_aq, temp, co2aq
      real ktk, crau, kkrau, epsf, epsrau, mui

!     photosynthesis after Evans & Parslow (1985)
!     notation as in JGOFS report No. 23 p. 6
      f1 = exp((-kw - kc*max(bioin(2),trcmin))*dzt)
      bct = bbio**(cbio*temp)
      jmax = abio*bct*(0.3+0.7*felim(i,j,mi))
      gd = jmax*dayfrac
      u1 = max(gl/gd, 1.e-6)
      u2 = u1*f1
c      phi1 = u1*(0.555588 + 0.004926*u1)/(1. + 0.188721*u1)
c      phi2 = u2*(0.555588 + 0.004926*u2)/(1. + 0.188721*u2)
      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
      avej = gd*rkw*(phi1 - phi2)
c      jmax_D = max(0.,abio_D*(bct-2.6))*(0.1+0.9*felim(i,j,mi))
      jmax_D = max(0.,abio_D*bct)*(0.1+0.9*felim(i,j,mi))
      gd_D = max(1.e-14,jmax_D*dayfrac)
      u1 = max(gl/gd_D, 1.e-6)
      u2 = u1*f1
c      phi1 = u1*(0.555588 + 0.004926*u1)/(1. + 0.188721*u1)
c      phi2 = u2*(0.555588 + 0.004926*u2)/(1. + 0.188721*u2)
      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
      avej_D = gd_D*rkw*(phi1 - phi2)
      nupt = nupt0*bct
      kkrau = 0.3583
      ktk = 5.019e-6*exp(-19510./(8.3143*(temp+273.15)))
      crau = 1.e4 + 1.82e-5/ktk
      epsf = 23.

      bioout(:) = 0.0
      bion = bioin(1)
      biop = bioin(2)
      bioz = bioin(3)
      biod = bioin(4)
      biono3 = bioin(5)
      biodiaz = bioin(6)
      biodic13 = bioin(7)
      biophytc13 = bioin(8)
      biozoopc13 = bioin(9)
      biodetrc13 = bioin(10)
      biodiazc13 = bioin(11)
      biodic = bioin(12)
       nflag = 0.5 + sign(0.5,bion - trcmin)
       pflag = 0.5 + sign(0.5,biop - trcmin)
       zflag = 0.5 + sign(0.5,bioz - trcmin)
       dflag = 0.5 + sign(0.5,biod - trcmin)
       no3flag = 0.5 + sign(0.5,biono3 - trcmin)
       diazflag = 0.5 + sign(0.5,biodiaz - trcmin)
       dic13flag = 0.5 + sign(0.5,biodic13 - 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)
       dicflag = 0.5 + sign(0.5,biodic - trcmin)
       bion = max(bion, trcmin)
       biop = max(biop, trcmin)
       bioz = max(bioz, trcmin)
       biod = max(biod, trcmin)
       biono3 = max(biono3, trcmin)
       biodiaz = max(biodiaz, trcmin)
       biodic13 = max(biodic13, trcmin*rc13std)
       biophytc13 = max(biophytc13, trcmin*rc13std)
       biozoopc13 = max(biozoopc13, trcmin*rc13std)
       biodetrc13 = max(biodetrc13, trcmin*rc13std)
       biodiazc13 = max(biodiazc13, trcmin*rc13std)
       biodic = max(biodic, trcmin)
       rc13expoout = 0.0
       ac13_DIC_aq = -1.0512994e-4*temp+1.0117605
c       ac13_DIC_aq = 1.
      expoout = 0.0
      grazout = 0.0
      morpout = 0.0
      morzout = 0.0
      nppout = 0.0
      morptout = 0.0
      remiout = 0.0
      excrout = 0.0
      npp_Dout = 0.0
      graz_Dout = 0.0
      morp_Dout = 0.0
      nfixout = 0.0

      do n=1,ntsb

!       growth rate of phytoplankton
        u_P = min(avej, jmax*bion/(k1p + bion))
!       nitrate limitation
        u_P = min(u_P, jmax*biono3/(k1 + biono3))
!       growth rate of diazotrophs smaller than other phytoplankton and
!       not nitrate limited
        u_D = min(avej_D, jmax_D*bion/(k1p + bion))
        npp = u_P*biop
        npp_D = max(0.,u_D*biodiaz)
!       grazing function for diazotrophs
        g_D = gbio*epsbio*biodiaz*biodiaz/(gbio+epsbio*biodiaz*biodiaz)
        graz_D = g_D*bioz*bct
        morp_D = nupt*biodiaz ! linear mortality
        no3upt_D = biono3/(k1 + biono3)*npp_D ! nitrate uptake
!       grazing function
        biop2 = biop*biop
        g_P = gbio*epsbio*biop2/(gbio+epsbio*biop2)
        graz = g_P*bioz*bct
        morp = nup*biop2
        morpt = nupt*biop
        morz = nuz*bioz*bioz
        remi = nud*bct*biod
        excr = gamma2*bct*bioz
        expo = wwd*biod
!       flags prevent negative values by setting outgoing fluxes to
!       zero if tracers are lower than trcmin
c        nflag = 0.5 + sign(0.5,bion - trcmin)
c        pflag = 0.5 + sign(0.5,biop - trcmin)
c        zflag = 0.5 + sign(0.5,bioz - trcmin)
c        dflag = 0.5 + sign(0.5,biod - trcmin)
c# if defined 1
c        no3flag = 0.5 + sign(0.5,biono3 - trcmin)
c        diazflag = 0.5 + sign(0.5,biodiaz - trcmin)
c# endif
        graz = graz*pflag*phytc13flag
        morp = morp*pflag*phytc13flag
        morpt = morpt*pflag*phytc13flag
        morz = morz*zflag*zoopc13flag
        remi = remi*dflag*detrc13flag
        excr = excr*zflag*zoopc13flag
        expo = expo*dflag*detrc13flag
        npp = npp*nflag*no3flag*dic13flag
        npp_D = npp_D*nflag
        graz_D = graz_D*diazflag*diazc13flag
        morp_D = morp_D*diazflag*diazc13flag
        no3upt_D = no3upt_D*no3flag*dic13flag
        rdic13 = biodic13/biodic
        rdic13 = min(rdic13, 2.)
        rdic13 = max(rdic13, 0.5)
c        rdic13 = rc13std

        rc13phyt = biophytc13/(biop*redctn)
        rc13phyt = min(rc13phyt, 2.)
        rc13phyt = max(rc13phyt, .5)
c        rc13phyt = rc13std

        rc13zoop = biozoopc13/(bioz*redctn)
        rc13zoop = min(rc13zoop, 2.)
        rc13zoop = max(rc13zoop, 0.5)
c        rc13zoop = rc13std

        rc13detr = biodetrc13/(biod*redctn)
        rc13detr = min(rc13detr, 2.)
        rc13detr = max(rc13detr, 0.5)
c        rc13detr = rc13std

        rc13diaz = biodiazc13/(biodiaz*redctn)
        rc13diaz = min(rc13diaz, 2.)
        rc13diaz = max(rc13diaz, 0.5)
c        rc13diaz = rc13std
!  !!!!!!!!!!!!!!!!!  nutrients equation
        bion = bion + tsb*redptn*(remi + excr - (npp + npp_D) + morpt)
!  !!!!!!!!!!!!!!!!!  phytoplankton equation
        biop = biop + tsb*(npp - morp - graz - morpt)
!  !!!!!!!!!!!!!!!!!  zooplankton equation
        bioz = bioz + tsb*(gamma1*(graz + graz_D) - excr - morz)
!  !!!!!!!!!!!!!!!!!  detritus equation
        biod = biod + tsb*((1.-gamma1)*(graz + graz_D) + morp + morp_D
     &       + morz - remi - expo + impo)
!  !!!!!!!!!!!!!!!!!  nitrate (NO3) equation
        biono3 = biono3 + tsb*(remi + excr - npp + morpt - no3upt_D)
!  !!!!!!!!!!!!!!!!!  diazotroph equation
        biodiaz = biodiaz + tsb*(npp_D - morp_D - graz_D)
        biodic = biodic + tsb*redctn*(remi + excr - (npp+npp_D) + morpt)
! !!!!!!!!!!!!!!!!!! isotope equations
c        ac13b = 1.
        ac13_aq_POC_P = -0.017*log10(co2aq*1000.)+1.0034
c       u_P is the daily averaged (24 h) growth rate
c       u_P*2.3 is the growth rate during the photoperiod (day time)
c       Laws et al. 1995 Geochim. et Cosm. Acta 59 No. 6, pp 1131 Fig. 2
c        ac13_aq_POC_P = 1.-(24.7-66.7*u_P*86400*2.35/(co2aq*1000.))/1000.
c       Rau et al. 1996 with parameters from Hofmann et al. 2000
c
!        mui = u_P*86400./dayfrac
!        epsrau = kkrau*crau*mui/co2aq - epsf
!        epsrau = max(min(epsrau,-15.),-40.)
!        ac13_aq_POC_P = 1.+epsrau/1000.
        ac13bp=ac13_aq_POC_P/ac13_DIC_aq
c        ac13bp=0.975
c       for diazotroph's use eps=15. (Carpenter et al. 1997 DSR I)
!        ac13_aq_POC_D = 1.-15./1000.
!        ac13bd=ac13_aq_POC_D/ac13_DIC_aq
        ac13bd=ac13_aq_POC_P/ac13_DIC_aq
c        ac13bd=0.975

        biodic13 = biodic13 + tsb*redctn*(rc13detr*remi
     &       + rc13zoop*excr - rdic13*(ac13bp*npp + ac13bd*npp_D)
     &       + rc13phyt*morpt)

        biophytc13 = biophytc13 + tsb*redctn*(ac13bp*rdic13*npp
     &       - rc13phyt*(morp + graz + morpt) )

        biozoopc13 = biozoopc13 + tsb*redctn*
     &         (gamma1*(rc13phyt*graz + rc13diaz*graz_D)
     &       - rc13zoop*(excr + morz) )

        biodetrc13 = biodetrc13 + tsb*redctn*((1.-gamma1)*(rc13phyt*graz
     &       + rc13diaz*graz_D) + rc13phyt*morp + rc13diaz*morp_D
     &       + rc13zoop*morz - rc13detr*(remi + expo)
     &       + rc13impo*impo)

        biodiazc13 = biodiazc13 + tsb*redctn*(ac13bd*rdic13*npp_D
     &       - rc13diaz*(morp_D + graz_D) )

        rc13expoout = rc13expoout + rc13detr
        expoout = expoout + expo
        grazout = grazout + graz
        morpout = morpout + morp
        morzout = morzout + morz
        nppout = nppout + npp
        morptout = morptout + morpt
        remiout = remiout + remi
        excrout = excrout + excr
        npp_Dout = npp_Dout + npp_D
        graz_Dout = graz_Dout + graz_D
        morp_Dout = morp_Dout + morp_D
        nfixout = nfixout + npp_D - no3upt_D
        if (nflag == 1) nflag = 0.5 + sign(0.5,bion - trcmin)
        if (pflag == 1) pflag = 0.5 + sign(0.5,biop - trcmin)
        if (zflag == 1) zflag = 0.5 + sign(0.5,bioz - trcmin)
        if (dflag == 1) dflag = 0.5 + sign(0.5,biod - trcmin)
        if (no3flag == 1) no3flag = 0.5 + sign(0.5,biono3 - trcmin)
        if (diazflag == 1) diazflag = 0.5 + sign(0.5,biodiaz - trcmin)
        dic13flag = 0.5 + sign(0.5,biodic13 - 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)
        dicflag = 0.5 + sign(0.5,biodic - trcmin)

      enddo

      bioout(1) = bion - bioin(1)
      bioout(2) = biop - bioin(2)
      bioout(3) = bioz - bioin(3)
      bioout(4) = biod - bioin(4)
      bioout(5) = biono3 - bioin(5)
      bioout(6) = biodiaz - bioin(6)
      bioout(7) = biodic13 - bioin(7)
      bioout(8) = biophytc13 - bioin(8)
      bioout(9) = biozoopc13 - bioin(9)
      bioout(10) = biodetrc13 - bioin(10)
      bioout(11) = biodiazc13 - bioin(11)
      rc13expoout = rc13expoout/ntsb
      return
      end