! source file: /raid23/jmuglia/ironruns3/lig15/lgm-picsed/updates/npzd_src.F
      subroutine mobii

!     initialize MOBI parameters
      implicit none
      include "size.h"
      include "npzd.h"
      include "calendar.h"
      include "coord.h"
      include "stdunits.h"
      include "scalar.h"
      integer k, ioun
      real alpha, par
      include "param.h"
      include "pconst.h"
      integer ib(10), ic(10), iou, j, i
      character (120) :: fname, new_file_name
      logical exists
      real c1e3
      real tmpijkm(1:100,1:100,1:19)
      namelist /npzd/   alpha, kw, kc, abio, bbio, cbio, k1n, nup
     &,                 gamma1, gbio, epsbio, nuz, nud0, LFe
     &,                 wd0, par, dtnpzd, redctn, ki, redptn, capr
     &,                 dcaco3, nupt0, redotn, jdiar, kzoo, geZ
     &,                 zprefP, zprefZ, zprefDet, zprefD, kfe, kfe_D
     &,                 nupt0_D, sgbdfac, diazntp, nup_D, dfr, dfrt
     &,                 nudon0, nudop0, eps_assim, eps_excr, eps_nfix
     &,                 eps_wcdeni, eps_bdeni0, eps_recy, hdop
     &,                 mw, mwz

!     set defaults for namelist npzd

      alpha = 0.16  ! Initial slope P-I curve [(W/m^2)^(-1)/day]
      kw = 0.04  ! Light attenuation due to water [1/m]
      kc = 0.047  ! Light atten. by phytoplankton [1/(m*mmol/m^3)]
      ki = 5.0  ! Light attenuation through sea ice & snow
      abio = 0.6  ! a; Maximum growth rate parameter [1/day]
      bbio = 1.066  ! b
      cbio = 1.0  ! c [1/deg_C]
      k1n = 0.7  ! Half saturation constant for N uptake [mmol/m^3]
      nup = 0.03  ! Specific mortality rate (Phytoplankton) [1/day]
      nup_D = 0.0001 !  Specific mortaility rate (Diazotrophs [1/day]
      nupt0 = 0.015  ! Fast-recycling  mortality rate (Phytoplankton) [1/day]
      nupt0_D = 0.001 ! Fast-recyling mortality rate (Diazotrophs) [1/day]
      gamma1 = 0.70  ! gama1; Assimilation efficiency (zpk)
      gbio = 0.38  ! Maximum grazing rate [1/day]
      epsbio = 1.6  ! Prey capture rate [(mmol/m^3)^(-2)day^(-1)]
      nuz = 0.06  ! Quadratic mortality (zpk) [(mmol/m^3)^(-2)day^(-1)]
      nud0 = 0.07  ! detritus remineralization rate [1/day]
      nudon0 = 2.33e-5 ! DON remineralization rate [1/day]
      nudop0 = 7.e-5  ! DOP remineralization rate [1/day]
      LFe = 1.0  ! Iron limitation
      wd0 = 16.0  ! Sinking speed of detritus at surface [m/day]
      mwz = 100000. !Depth where sinking below remains constant (cm)
      mw = 0.045 !Sinking speed increase with depth (d-1)
      par = 0.43  ! fraction of photosythetically active radiation
      dtnpzd = dtts/4.  ! time step of biology [s]
      redctn = 7.1  ! C/N Redfield ratio (Schneider et al., 2003, GBC)
      redptn = 1./16.  ! P/N Redfield ratio
      capr = 0.022  ! carbonate to carbon production ratio
      dcaco3 = 650000.0  ! remineralisation depth of calcite [cm]
      redotn = 10.6  ! O2/N Redfield ratio (Anderson & Sarmiento, 1994, GBC)
      jdiar = 0.08  ! factor reducing the growth rate of diazotrophs
      kzoo = 0.15  ! half sat. constant for Z grazing [mmol N m^3]
      geZ = 0.6   ! Zooplankton growth efficiency
      sgbdfac = 1.0 ! sub-grid benthic denitrification rate factor
      diazntp = 28. !diazotroph N:P ratio
      dfr = 0.08   ! phyt mortality refractory/semi-labile DOM fraction
      dfrt = 0.01  ! phyt fast-recy refracotyr/semi-labile DOM fraction
      hdop = 0.4 ! DOP growth rate handicap
      zprefP = 0.3       ! Zooplankton preference for P
      zprefZ = 0.3       ! Zooplankton preference for other Z
      zprefDet = 0.3     ! Zooplankton preference for Detritus
      zprefD = 0.1       ! Zooplankton preference for diazotrophs
      kfemin = 0.04e-3    ! Minimum half saturation constant for Fe limitation [mmol Fe / m-3]
      kfemax = 0.2e-3    ! Maximum half saturation constant for Fe limitation [mmol Fe / m-3]
      pmax = 0.15    ! Phytoplankton biomass above which kfe increases [mmol N / m-3]
      kfe_D = 0.1e-3 ! Half saturation constant for Diaz Fe limitation [mmol Fe / m-3]
      kfeleq = 10.**5.5   !  Fe-ligand stability constant [m^3/(mmol ligand)]
      lig = 1.5e-3      !  Ligand concentration  [mmol/m^3]
      thetamaxhi = 0.04    !  Maximum Chl:C ratio, abundant iron [gChl/(gC)]
      thetamaxlo = 0.01    !  Maximum Chl:C ratio, extreme iron limitation [gChl/(gC)]
      alphamax = 73.6e-6*86400  !   Maximum initial slope in PI-curve [day^-1 (W/m^2)^-1 (gChl/(gC))^-1]
      alphamin = 18.4e-6*86400   !   Minimum initial slope in PI-curve [day^-1 (W/m^2)^-1 (gChl/(gC))^-1]
      mc = 12.011           ! Molar mass of carbon [g/mol]
      fetopsed = 0.005  !  Fe:P for sedimentary iron source [molFe/molP]
      o2min = 5.      !  Minimum O2 concentration for aerobic respiration [mmolO_2/m^3]
      kfeorg = 1.0*(1./86400.)     !  Organic-matter dependent scavenging rate [(m^3/(gC s))^0.58]
      rfeton = 10.e-6*6.625        ! Uptake ratio of iron to nitrogen [mol Fe/mol N] = 10 micromol Fe /mol C
      kfecol = 0.016/86400.     ! Colloidal production and precipitation rate [s^-1]
!     read namelist
      call getunit (ioun, 'control.in', 'f s r')
      read  (ioun, npzd, end=108)
108   continue
      write (stdout, npzd)
      call relunit (ioun)

!     convert units of NPZD parameters to MOM units
      redctn = redctn*1.e-3
      redotn = redotn*1.e-3
      redotp = redotn/redptn
      redctp = redctn/redptn
      redotc = redotn/redctn
      redntp = 1./redptn
      redntc = 1./redctn
      diazptn = 1./diazntp
      k1p   = k1n*redptn
      kw = kw*1.e-2
      kc = kc*1.e-2
      ki = ki*1.e-2
      wd0 = wd0*1.e2
      alpha = alpha/daylen
      abio = abio/daylen
      nup = nup/daylen
      nup_D = nup_D/daylen
      nupt0 = nupt0/daylen
      nupt0_D = nupt0_D/daylen
      gbio = gbio/daylen
      epsbio = epsbio/daylen
      nuz = nuz/daylen
      nud0 = nud0/daylen
      nudop0 = nudop0/daylen
      nudon0 = nudon0/daylen
      tap = 2.*alpha*par
!     calculate sinking speed of detritus divided by grid width
      do k=1,km
!     linear increase wd0-200m with depth
         if (zt(k) .lt. mwz) then
            wd(k) = (wd0+mw*zt(k))/daylen/dzt(k) ! [s-1]
         else
            wd(k) = (wd0+mw*mwz)/daylen/dzt(k) ! [s-1]
         endif
         rkwz(k) = 1./(kw*dzt(k))
      enddo
      ztt(1)=0.0
      do k=1,km
         ztt(k+1)=(-1)*zw(k)
      enddo
! read the hydrothermal Fe input from O_fe_hydr.nc
      fe_hydr(:,:,:) = 0.
      ib(:) = 1
      ic(:) = imtm2
      ic(2) = jmtm2
      ic(3) = km

      fname = new_file_name ("O_fe_hydr.nc")
      inquire (file=trim(fname), exist=exists)
      if (exists) then
         c1e3 = 1000
         call openfile (trim(fname), iou)
         call getvara ('O_fe_hydr', iou, ic(1)*ic(2)*ic(3)
     &,                 ib, ic, tmpijkm, c1e3, c0)
         fe_hydr(2:imtm1,2:jmtm1,:) = tmpijkm(1:imtm2
     &,                                            1:jmtm2,:)

            do k=1,km
               do j=1,jmt
                  fe_hydr(1,j,k) = fe_hydr(imtm1,j,k)
                  fe_hydr(imt,j,k) = fe_hydr(2,j,k)
               enddo
               do i=1,imt
                  fe_hydr(i,1,k) = fe_hydr(i,2,k)
                  fe_hydr(i,jmt,k) = fe_hydr(2,j,k)
               enddo
            enddo
      else
         print*,"Warning => Cannot find", trim(fname)
      endif

!---------------------------------------------------------------------
!     calculate variables used in calcite remineralization
!---------------------------------------------------------------------

      rcak(1) = -(exp(-zw(1)/dcaco3)-1.0)/dzt(1)
      rcab(1) = -1./dzt(1)
      do k=2,km
        rcak(k) = -(exp(-zw(k)/dcaco3))/dzt(k)
     &          + (exp(-zw(k-1)/dcaco3))/dzt(k)
        rcab(k) = (exp(-zw(k-1)/dcaco3))/dzt(k)
      enddo
      return
      end
!     END mobii

      subroutine mobi (kmx, twodt, rctheta, dayfrac, swr, tnpzd
     &,                t_in, po4_in
     &,                o2_in
     &,                s_in, dic_in, alk_in, co2_in, dic13_in
     &,                no3_in, sgb_in
     &,                din15_in
     &,                src)

!     main driver for Model of Ocean Biogeochemistry and Isotopes (MOBI)
!     input: kmx, twodt, rctheta, dayfrac, swr, tnpzd, t_in,
!            po4_in, (o2_in, s_in, dic_in, alk_in, co2_in, dic13_in, no3_in,
!            sgb_in, din15_in, felimit_in, felimit_D_in)
!     output: src, (sedcorgflx, sedcalflx)
!             if 1 is switch on additional output (rnpp etc.)
!                is in npzd.h
!     copied from tracer.F (Andreas Schmittner, Oct 2, 2013)
!     modified to a vertical column model
      implicit none
      integer k,kmx
      include "size.h"   ! km maximum # of levels
      include "coord.h"  ! zt, dzt
      include "grdvar.h" ! dztr = 1/dzt
      include "mw.h"     ! ispo4 etc.
      include "npzd.h"
      real snpzd(ntnpzd), tnpzd(km,ntnpzd), src(km,nsrc)
      real dic_npzd_sms(km), twodt, rctheta, gl, impo, expo
      real npp, remi, excr, graz, morp, morpt, morz, temp, swr, dayfrac
      real graz_Det, graz_Z, phin, dz, prca, dprca, nud, bct, bctz
      real t_in(km), po4_in(km), sedrr, nudop, nudon, npp_dop
      real impofe, expofe, remife
      real fesed, bfe, sgb
      real o2_in(km), fo2, so2
      real s_in(km), dic_in(km), alk_in(km), co2_in, dic13_in(km)
      real rc13impo, rc13expo, ac13_DIC_aq
      real ac13_aq_POC, ac13b, prca13, rdic13, rtdic13
      real pH, co2star, dco2star, pCO2, dpco2, CO3
      real Omega_c, Omega_a, pCO2_old, pCO2_new
      real atmpres, depth
      real no3_in(km), sgb_in(km), lno3, sg_bdeni
      real morp_D, lntp, npp_D_dop
      real npp_D, graz_D, morpt_D, no3flag, wcdeni, bdeni(km), nfix(km)
      real din15_in(km), din15flag, eps_bdeni
      real rn15impo, rn15expo, uno3, rno3, bwcdeni, bbdeni

      expo = 0.0
      impo = 0.0
      phin = 0.0                ! integrated phytoplankton
      prca = 0.0                ! integrated production of calcite

      sedrr = 0.0
      rsedrr = 0.0
      rprca = 0.0
      rn15impo = 0.0
      rn15expo = 0.0
      rc13impo = 0.0
      rc13expo = 0.0
      prca13 = 0.0

            expofe = 0.0
            impofe = 0.0

!1111 start of main k-loop
      do k=1,kmx
         rn15impo = rn15expo
!        c13 biological fractionation
         atmpres = 1.0          !atm
         depth = zt(k)/100.     !m
         call co2calc_SWS (t_in(k), s_in(k), dic_in(k), alk_in(k)
     &,                    co2_in, atmpres, depth, pH, co2star, dco2star
     &,                    pCO2, dpco2
     &,                    CO3, Omega_c, Omega_a)

         ac13_DIC_aq = -1.0512994e-4*t_in(k)+1.011765
!        Popp et al. (1989) Am. J. Sci.
         ac13_aq_POC = -0.017*log10(min(max(co2star*1000.,2.0),74.0))+1.0034
         ac13b = ac13_aq_POC/ac13_DIC_aq

         rc13impo = rc13expo
         swr = swr*exp(-kc*phin)
!     fix phin bug
!     previous phyt layers already saved in swr term
c         phin = phin + max(tnpzd(k,3), trcmin)*dzt(k)
         phin = max(tnpzd(k,3), trcmin)*dzt(k)
     &                + max(tnpzd(k,9), trcmin)*dzt(k)
         gl = swr*exp(ztt(k)*rctheta)
         impo = expo*dztr(k)
              impofe = expofe*dztr(k)

         bct = bbio**(cbio*t_in(k))
         bctz = (0.5*(tanh(o2_in(k) - 8.)+1))
     &           *bbio**(cbio*min(t_in(k),20.))
!        decrease remineralisation rate in oxygen minimum zone
         nud = nud0*(0.65+0.35*tanh(o2_in(k) - 3.0))
         nudon = nudon0
         nudop = nudop0

!-----------------------------------------------------------------------
!        call the npzd ecosystem dynamics model
!-----------------------------------------------------------------------

         call npzd_src (tnpzd(k,:), gl, bct, impo, dzt(k)
     &,                 dayfrac, wd(k), rkwz(k), nud, nudop, nudon
     &,                 snpzd, expo, graz, morp, morz, graz_Det, graz_Z
     &,                 npp, npp_dop, morpt, remi, excr
     &,                 npp_D, npp_D_dop, graz_D, morp_D, morpt_D
     &,                 nfix(k)
     &,                 bctz
     &,                 rn15impo, rn15expo
     &,                 rc13impo, rc13expo, ac13b
     &,                    expofe, impofe, remife
cjuan: use o2_in(k) instead of t(i,k,j,io2,taum1)
     &,                    o2_in(k)
!
     &                      )
! These are source/sink terms
         snpzd = snpzd*rdtts
              expofe=expofe*rnbio
         expo = expo*rnbio
         rn15expo = rn15expo*rnbio
         rc13expo = rc13expo*rnbio
         rexpo(k) = expo
         rgraz(k) = graz*rnbio
         rgraz_Det(k) = graz_Det*rnbio
         rgraz_Z(k) = graz_Z*rnbio
         rmorp(k) = morp*rnbio
         rmorz(k) = morz*rnbio
         rnpp(k) = npp*rnbio
         rnpp_dop(k) = npp_dop*rnbio
         rmorpt(k) = morpt*rnbio
         rremi(k) = remi*rnbio
         rexcr(k) = excr*rnbio
         rnpp_D(k) = npp_D*rnbio
         rnpp_D_dop(k) = npp_D_dop*rnbio
         rgraz_D(k) = graz_D*rnbio
         rmorp_D(k) = morp_D*rnbio
         rmorpt_D(k) = morpt_D*rnbio
         rnfix(k) = nfix(k)*rnbio
              rexpofe(k) = expofe
              rremife(k) = remife*rnbio
!-----------------------------------------------------------------------
!             calculate detritus at the bottom and remineralize
!-----------------------------------------------------------------------

         sedrr = sgb_in(k)*expo*dzt(k)
         snpzd(1) = snpzd(1) + sgb_in(k)*expo*redptn
         snpzd(6) = snpzd(6) + sgb_in(k)*expo*redctn
c juan: CHANGED sgb FOR sgb_in(k)
                fesed=fetopsed*bct*redptn*expo*sgb_in(k)
                bfe=fesed
                snpzd(22) = snpzd(22) + fesed
                if (o2_in(k) .lt. o2min) then
                  snpzd(22) = snpzd(22) + expofe*sgb_in(k)
                  bfe=bfe+expofe*sgb_in(k)
                  expofe = expofe - sgb_in(k)*expofe
                endif
                rremife(k) = rremife(k) + bfe
         snpzd(16) = snpzd(16) + rc13expo*sgb_in(k)*expo*redctn
!-----------------------------------------------------------------------
!        benthic denitrification model of Bohlen et al., 2012, GBC (eq. 10)
!        NO3 is removed out of bottom water nitrate.
!        See Somes et al., 2012, BGS for additional details/results
!-----------------------------------------------------------------------
!        limit denitrification as nitrate approaches 0 uM
         no3flag = 0.5+sign(0.5,no3_in(k)-trcmin)
         din15flag = 0.5+sign(0.5,din15_in(k)-trcmin)
         lno3 = 0.5*tanh(no3_in(k)*10 - 5.0)

         sg_bdeni = (0.06 + 0.19*0.99**(max(o2_in(k),trcmin)
     &                 - max(no3_in(k),trcmin)))
     &                   *max(expo*sgb_in(k),trcmin)*redctn*1.e3

         sg_bdeni = min(sg_bdeni, sgb_in(k)*expo)
         sg_bdeni = max(sg_bdeni, 0.)
c remove ?
c         if (zt(k) .le. 200000.) sg_bdeni = sgbdfac*sg_bdeni

         sg_bdeni = sg_bdeni*(0.5 + lno3)*no3flag
     &                   *din15flag
         snpzd(7) = snpzd(7) + sgb_in(k)*expo - sg_bdeni
         rno3 = max(din15_in(k), trcmin*rn15std/(1+rn15std))
     &            /max(no3_in(k)-din15_in(k),trcmin*rn15std/(1+rn15std))
         rno3 = min(rno3, 2.*rn15std)
         rno3 = max(rno3, rn15std/2.)
         eps_bdeni = 6.0*exp(-2.5e-6*(zt(k)))
         bbdeni = rno3 - eps_bdeni*rno3/1000.

         snpzd(10) = snpzd(10) + rn15expo*sgb_in(k)*expo
     &                     - bbdeni/(1+bbdeni)*sg_bdeni

         expo = expo - sgb_in(k)*expo
         rremi(k) = rremi(k) + sgb_in(k)*expo
         rsedrr = rsedrr + sedrr
         bdeni(k) = sg_bdeni
         rbdeni(k) = bdeni(k)

!-----------------------------------------------------------------------
!             set source/sink terms
!-----------------------------------------------------------------------

         src(k,ispo4) = snpzd(1)
         src(k,isdop) = snpzd(2)
         src(k,isphyt) = snpzd(3)
         src(k,iszoop) = snpzd(4)
         src(k,isdetr) = snpzd(5)
         src(k,isdic) = snpzd(6)
         src(k,isno3) = snpzd(7)
         src(k,isdon) = snpzd(8)
         src(k,isdiaz) = snpzd(9)
         src(k,isdin15) = snpzd(10)
         src(k,isdon15) = snpzd(11)
         src(k,isphytn15) = snpzd(12)
         src(k,iszoopn15) = snpzd(13)
         src(k,isdetrn15) = snpzd(14)
         src(k,isdiazn15) = snpzd(15)
         src(k,isdic13) = snpzd(16)
         src(k,isdoc13) = snpzd(17)
         src(k,isphytc13) = snpzd(18)
         src(k,iszoopc13) = snpzd(19)
         src(k,isdetrc13) = snpzd(20)
         src(k,isdiazc13) = snpzd(21)
cjuan: MOVED THE COMMENTED PART TO TRACER.F
c              if ( k .eq. 1 ) then
c                snpzd(22) = snpzd(22) +    ! comes in mol/m2/s gets converted to mmol/m3/s
c     &                   sbc(i,j,idfeadep)*1000./(dzt(1)/100.)  !dzt is in cm
c              endif
         src(k,isdfe) = snpzd(22)
         src(k,isdetrfe) = snpzd(23)
!
!        production of calcite
         dprca = (morp+morz+(graz+graz_Z)*(1.-gamma1))
     &              *capr*redctn*rnbio
         prca = prca + dprca*dzt(k)

!        These are sources and sinks of DIC (i.e. remin - pp)
!        all are based on po4 uptake and remineralization
!        dprca is a correction term
         dic_npzd_sms(k) = snpzd(6)
         src(k,isdic) = (snpzd(6) - dprca)
         rtdic13 = max(dic13_in(k),trcmin*rc13std/(1+rc13std))
     &             / max(dic_in(k),trcmin)
         rtdic13 = min(rtdic13, 2.*rc13std/(1+rc13std))
         rtdic13 = max(rtdic13, 0.5*rc13std/(1+rc13std))

         prca13 = prca13 + dprca*dzt(k)*rtdic13

         src(k,isdic13) = src(k,isdic13) - rtdic13*dprca
         src(k,isalk) = (-snpzd(6)*redntc*1.e-3 - 2.*dprca)
!        calculate total export to get total import for next layer
         expo = expo*dzt(k)
              expofe = expofe*dzt(k)

      enddo
!1111 end of main k-loop

      rprca = prca

!2222 start of second k-loop
      do k=1,kmx
!        limit oxygen consumption below concentrations of
!        5umol/kg as recommended in OCMIP
         fo2 = 0.5*tanh(o2_in(k) - 2.5)
!        sink of oxygen
! O2 is needed to generate the equivalent of NO3 from N2 during N2 fixation
! 0.5 H2O + 0.5 N2+1.25O2 -> HNO3
! note that so2 is -dO2/dt
         so2 = dic_npzd_sms(k)*redotc + nfix(k)*rnbio*1.25e-3
c         so2 = dic_npzd_sms(k)*redotc

         src(k,iso2) = -so2*(0.5 + fo2)
!        add denitrification as source term for NO3
         no3flag = 0.5+sign(0.5,no3_in(k)-trcmin)
         din15flag = 0.5+sign(0.5,din15_in(k)-trcmin)
         lno3 = 0.5*tanh(no3_in(k) - 2.5)
         lntp = 0.5*tanh(no3_in(k)/(redntp*po4_in(k))*100 - 60.)
!        800 = 0.8*1000 = (elec/mol O2)/(elec/mol NO3)*(mmol/mol)
         wcdeni = 800.*no3flag*so2*(0.5 - fo2)*(0.5 + lno3)*din15flag
         wcdeni = max(wcdeni, 0.)
         src(k,isno3) = src(k,isno3) - wcdeni
!        calculate isotope effect of water column denitrification
         uno3 = wcdeni*twodt/no3_in(k)
         uno3 = min(uno3, 0.999)
         uno3 = max(uno3, trcmin)
         rno3 = max(din15_in(k),trcmin*rn15std/(1+rn15std))
     &             / max(no3_in(k)-din15_in(k),
     &                   trcmin*rn15std/(1+rn15std))
         rno3 = min(rno3, 2.*rn15std)
         rno3 = max(rno3, rn15std/2.)
         bwcdeni = rno3 + eps_wcdeni*(1-uno3)/uno3*log(1-uno3)*rno3/1000.

         src(k,isdin15) = src(k,isdin15) - (bwcdeni/(1+bwcdeni))*wcdeni
! Correct the ALK stoichiometry to account for N2 fixation
         src(k,isalk) = src(k,isalk) + wcdeni*1.e-3
!bdeni
         src(k,isalk) = src(k,isalk) + bdeni(k)*1.e-3
! Now add account for N2 fixation (ALK production is tied to PO4 change and
! thus in the case of N2 fixation is not correct).
         src(k,isalk) = src(k,isalk) - nfix(k)*rnbio*1.e-3
         rwcdeni(k) = wcdeni
      enddo
!2222 end of second k-loop

!-----------------------------------------------------------------------
!           remineralize calcite
!-----------------------------------------------------------------------
!3333 start of third k-loop
      do k=1,kmx-1
         src(k,isdic) = src(k,isdic) + prca*rcak(k)
         src(k,isdic13) = src(k,isdic13) + prca13*rcak(k)
         src(k,isalk) = src(k,isalk) + 2.*prca*rcak(k)
      enddo
!3333 end of third k-loop
      src(kmx,isdic) = src(kmx,isdic) + prca*rcab(kmx)
      src(kmx,isdic13) = src(kmx,isdic13) + prca13*rcab(kmx)
      src(kmx,isalk) = src(kmx,isalk) + 2.*prca*rcab(kmx)
      return
      end
!     END mobi

      subroutine npzd_src (bioin, gl, bct, impo, dzt
     &,                    dayfrac, wwd, rkw, nud, nudop, nudon
     &,                    bioout, expoout, grazout, morpout
     &,                    morzout, graz_Det_out, graz_Zout
     &,                    nppout, npp_dopout, morptout, remiout
     &,                    excrout
     &,                    npp_Dout, npp_D_dopout, graz_Dout, morp_Dout
     &,                    morpt_Dout, nfixout
     &,                    bctz
     &,                    rn15impo, rn15expoout
     &,                    rc13impo, rc13expoout, ac13b
     &,                    expofeout, impofe, remifeout
     &,                    o2
     &                     )

!=======================================================================
!     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.
!     Keller et al. 2012, Geosci. Model Dev., 5, 1195-1220
!
!     The nitrogen isotope model (nitrogen_15) has been developed by
!     Chris Somes and is described in
!     Somes et al. 2010, Global Biogeochem. Cycles 24, GB4019
!     Somes et al. 2010, Geophys. Res. Lett. 37, L23605 and
!     Somes et al. 2013, Biogeosc. 10, 5889-5910
!
!     The carbon isotope model (carbon_13) has been written by Andreas
!     Schmittner and is documented in
!     Schmittner et al. 2013 Biogeosc. 10, 5793-5816
!     Chris Somes modified the code converting from alpha to beta formulation
!
!     Note that nutrient (N) 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
!       nbio       = number of time steps
!       dtbio        = 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]
!
!     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]
!       npp_Dout   = NPP of diazotrophs
!       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

      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
      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
      real excrdiaz, biodop, biodon, dopflag, donflag, recy_don
      real recy_dop, npp, graz, biodin15, biodon15, biophytn15
      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, biozoopn15, morp_Dout
      real don15flag, phytn15flag, zoopn15flag, detrn15flag, diazn15flag
      real limP, limP_D, dopupt_D_flag, dopupt_D, morp_D, nupt_D
      real udon, rdon, brecy, fcrecy, biodic, dopupt_flag, dopupt
      real biodic13, biophytc13, biozoopc13, biodetrc13, biodiazc13
      real dic13flag, phytc13flag, zoopc13flag, detrc13flag, diazc13flag
      real biodoc13, doc13flag, biodoc
      real rdic13, rtphytc13, rtzoopc13, rtdetrc13, rtdiazc13, ac13b
      real rc13impo, rc13expoout, dicflag, bc13npp, rtdic13, fcnpp
      real rtdoc13, nudop, nudon, npp_dopout, npp_D_dopout

      real biodfe, biodetrfe, dfeflag, detrfeflag
      real expofe, impofe, feorgads, remife, thetamax, deffe, fepa
      real thetamax_D, deffe_D, thetachl, thetachl_D, chl, feprime
      real fecol, irrtop, kirr, aveirr, fesed, gl_O, gl_D, alpha_O
      real alpha_D, par

      real expofeout, impofeout, feorgadsout, fecolout
      real remifeout, thetamaxout, deffeout, thetachlout
      real chlout, feprimeout, o2, chl_D_out

      real p1,p2,kfevar

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

      real bioin(ntnpzd), bioout(ntnpzd)

      p1=min(bioin(3),pmax)
      p2=max(0.0,bioin(3)-pmax)
      kfevar=(kfemin*p1+kfemax*p2)/(p1+p2)
      deffe=bioin(22)/(kfevar+bioin(22))
      thetamax=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe
      alpha_O=alphamin+(alphamax-alphamin)*deffe
      gl_O=gl*thetamax*alpha_O
      deffe_D=bioin(22)/(kfe_D+bioin(22))
      thetamax_D=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe_D
      alpha_D=alphamin+(alphamax-alphamin)*deffe_D
      gl_D=gl*thetamax_D*alpha_D

      bioout(:) = 0.0
      biopo4 = bioin(1)
      biodop = bioin(2)
      biophyt = bioin(3)
      biozoop = bioin(4)
      biodetr = bioin(5)
      biodic = bioin(6)
      biono3 = bioin(7)
      biodon = bioin(8)
      biodiaz = bioin(9)
      biodin15 = bioin(10)
      biodon15 = bioin(11)
      biophytn15 = bioin(12)
      biozoopn15 = bioin(13)
      biodetrn15 = bioin(14)
      biodiazn15 = bioin(15)
      biodic13 = bioin(16)
      biodoc13 = bioin(17)
      biophytc13 = bioin(18)
      biozoopc13 = bioin(19)
      biodetrc13 = bioin(20)
      biodiazc13 = bioin(21)
      biodfe = bioin(22)
      biodetrfe = bioin(23)

!       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)
      no3flag = 0.5 + sign(0.5,biono3 - trcmin)
      donflag = 0.5 + sign(0.5,biodon - trcmin)
      diazflag = 0.5 + sign(0.5,biodiaz - trcmin)
      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)
      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)
        dfeflag = 0.5 + sign(0.5,biodfe - trcmin)
        detrfeflag = 0.5 + sign(0.5,biodetrfe - trcmin)

!     limit tracers to positive values
      bioin(:) = max(bioin(:), trcmin)
      biopo4 = max(biopo4, trcmin)
      biodop = max(biodop, trcmin)
      biophyt = max(biophyt, trcmin)
      biozoop = max(biozoop, trcmin)
      biodetr = max(biodetr, trcmin)
      biodic = max(biodic, trcmin)
      biono3 = max(biono3, trcmin)
      biodon = max(biodon, trcmin)
      biodiaz = max(biodiaz, trcmin)
      biodin15 = max(biodin15, trcmin)
      biodon15 = max(biodon15, trcmin)
      biophytn15 = max(biophytn15, trcmin)
      biozoopn15 = max(biozoopn15, trcmin)
      biodetrn15 = max(biodetrn15, trcmin)
      biodiazn15 = max(biodiazn15, trcmin)
      biodic13 = max(biodic13, trcmin)
      biodoc13 = max(biodoc13, trcmin)
      biophytc13 = max(biophytc13, trcmin)
      biozoopc13 = max(biozoopc13, trcmin)
      biodetrc13 = max(biodetrc13, trcmin)
      biodiazc13 = max(biodiazc13, trcmin)
      biodfe = max(biodfe, trcmin)
      biodetrfe = max(biodetrfe, trcmin)

!     photosynthesis after Evans & Parslow (1985)
!     notation as in JGOFS report No. 23 p. 6
      f1 = exp((-kw - kc*(biophyt+biodiaz))*dzt)
      jmax = abio*bct*deffe
      gd = jmax*dayfrac
      u1 = max(gl_O/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
      avej = gd*(phi1 - phi2)/((kw+kc*(biophyt+biodiaz))*dzt)
! Make the max grazing rate a function of temperature
      gmax = gbio*bctz
! Note bctz, which sets an upper limit on the effects of temp on the
! grazing rate, is set in tracers.F
      jmax_D = max(0.,abio*(bct - 2.6)*deffe_D)*jdiar
!
      gd_D = max(1.e-14,jmax_D*dayfrac)
      u1 = max(gl_D/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
      avej_D = gd_D*(phi1 - phi2)/((kw+kc*(biophyt+biodiaz))*dzt)
!     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
      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
      rn15expoout = 0.0
      rc13expoout = 0.0
      nppout = 0.0
      npp_dopout = 0.0
      morptout = 0.0
      remiout = 0.0
      excrout = 0.0
      npp_Dout = 0.0
      npp_D_dopout = 0.0
      graz_Dout = 0.0
      morpt_Dout = 0.0
      morp_Dout = 0.0
      nfixout = 0.0
      expofeout = 0.0
      remifeout = 0.0

      do n=1,nbio

        p1 = min(biophyt,pmax)
        p2 = max(0.0,biophyt - pmax)
        kfevar = (kfemin*p1+ kfemax*p2)/(p1 + p2)
        deffe = biodfe/(kfevar + biodfe)
        jmax = abio*bct*deffe
        deffe_D = biodfe/(kfe_D + biodfe)
        jmax_D = max(0.,abio*(bct - 2.6)*deffe_D)*jdiar

!       growth rate of phytoplankton
!       consume DOP when it is more efficient
        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)

!       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
        npp = u_P*biophyt
        dopupt = npp*dopupt_flag
        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*biodiaz
c        no3upt_D = biono3/(k1n + biono3)*npp_D ! nitrate uptake
        no3upt_D = (0.5+0.5*tanh(biono3-5.))*npp_D ! nitrate uptake
        dopupt_D = npp_D*dopupt_D_flag
!       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_don = nudon*bct*biodon
        recy_dop = nudop*bct*biodop
        morz = nuz*biozoop*biozoop
        remi = nud*bct*biodetr
        expo = wwd*biodetr

!   Calculation of average light in mixed layer for calculation of Chl diagnostic
        par = 0.43  ! fraction of photosythetically active radiation
        irrtop=gl/2/par
        kirr=-kw - kc*(bioin(2)+bioin(6))
        aveirr=-1/dzt/kirr*(irrtop-irrtop*exp(kirr*dzt))
!  remineralization of iron from organic matter
        remife=nud*bct*biodetrfe

!   Scavenging of dissolved iron is based on Honeymoon et al. (1988)
!   and Parekh et al. (2004).

        if (o2*1e3 .gt. o2min) then

! KFeL proportional to biomass

          fepa = 1.0 + kfeleq * (lig - biodfe)

          feprime = (-fepa +(fepa
     &         * fepa + 4.0 * kfeleq
     &         * biodfe)**(0.5)) /(2.0 * kfeleq)

          feorgads = kfeorg*((biodetr*mc*redctn)**0.58)*feprime
          fecol = kfecol*feprime
       else
!   These are default values in case of O2 < O2min, no scavenging is taking place
          feorgads = 0.0
          fecol = 0.0
          feprime = 0.0
          fepa = 0.0
       endif

        expofe = wwd * biodetrfe

!

        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_dop = recy_dop*dopflag
        if (dopupt_flag .eq. 1) then
           npp = npp*dopflag*no3flag*din15flag
        else
           npp = npp*po4flag*no3flag*din15flag
        endif
        if (dopupt_D .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_don = recy_don*donflag*don15flag

        remife=remife*detrfeflag
        feorgads=feorgads*dfeflag
        expofe=expofe*detrfeflag
        fecol=fecol*dfeflag

!     digestion of grazed material
        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

!     excretion based on growth efficiency
        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
!       excrete "extra" non-Redfield diazotroph N
        nr_excr_D = gamma1*graz_D*(1-(redntp/diazntp))
     &       + (1-gamma1)*graz_D*(1-(redntp/diazntp))

!     sloppy feeding
        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
!       calculate isotope parameters
!       See Somes et al., 2010, GBC for details/results
        uno3 = npp*dtbio/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 + eps_assim*(1-uno3)/uno3*log(1-uno3)*rno3/1000.
        fcassim = bassim/(1+bassim)

        udon = recy_don*dtbio/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 + eps_recy*(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 - eps_excr*rzoop/1000.
        fcexcr = bexcr/(1+bexcr)

        bnfix = rn15std - eps_nfix*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.)
        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))
!       nutrients equation
        biopo4 = biopo4 + dtbio*(redptn*((1.-dfrt)*morpt + excr + remi
     &       - (npp-dopupt)) + diazptn*(morpt_D - (npp_D-dopupt_D))
     &       + recy_dop)
!       DOP equation
        biodop = biodop + dtbio*(redptn*dfr*morp + redptn*dfrt*morpt
     &       - redptn*dopupt - diazptn*dopupt_D - recy_dop)
!       phytoplankton equation
        biophyt = biophyt + dtbio*(npp - morp - graz - morpt)
!       zooplankton equation
        biozoop = biozoop + dtbio*(dig - morz - graz_Z - excr)
!       detritus equation
        biodetr = biodetr + dtbio*((1.-dfr)*morp + sf + morz - remi
     &       - graz_Det - expo + impo + morp_D*(redntp/diazntp))
!       DIC equation
        biodic = biodic + dtbio*redctn*((1.-dfrt)*morpt + excr + remi
     &       - npp + morpt_D - npp_D + recy_don + nr_excr_D
     &       + morp_D*(1-(redntp/diazntp)))
!       nitrate (NO3) equation
        biono3 = biono3 + dtbio*((1.-dfrt)*morpt + excr + morpt_D
     &       + morp_D*(1-(redntp/diazntp)) + nr_excr_D + remi + recy_don
     &       - npp - no3upt_D)
!       DON equation
        biodon = biodon + dtbio*(dfr*morp + dfrt*morpt - recy_don)
!       diazotroph equation
        biodiaz = biodiaz + dtbio*(npp_D - morp_D - morpt_D - graz_D)
!       dissolved iron equation
        biodfe = biodfe + dtbio*(rfeton*(excr + (1.-dfrt)*morpt
     &            - npp + morpt_D - npp_D + recy_don + nr_excr_D
     &            + morp_D*(1-(redntp/diazntp))) - feorgads
     &            + remife - fecol)
!       particulate iron equation
        biodetrfe = biodetrfe + dtbio*(rfeton*(sf + (1.-dfr)*morp
     &            + morp_D*(redntp/diazntp) + morz - graz_Det)
     &            + feorgads + fecol - remife - expofe + impofe)
!       isotope equations
        biodin15 = biodin15 + dtbio*(rtphytn15*(1.-dfrt)*morpt
     &       + fcexcr*excr + rtdiazn15*morpt_D + rtdiazn15*nr_excr_D
     &       + rtdiazn15*morp_D*(1-(redntp/diazntp)) + rtdetrn15*remi
     &       + fcrecy*recy_don - fcassim*npp - fcassim*no3upt_D)

        biodon15 = biodon15 + dtbio*(dfr*rtphytn15*morp
     &       + dfrt*rtphytn15*morpt - fcrecy*recy_don)

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

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

        biodetrn15 = biodetrn15 + dtbio*(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 + dtbio*(fcnfix*(npp_D-no3upt_D)
     &       + fcassim*no3upt_D - rtdiazn15*morp_D - rtdiazn15*graz_D
     &       - rtdiazn15*morpt_D)
! !!!!!!!!!!!!!!!!!! isotope equations
        biodic13 = biodic13 + dtbio*redctn*(rtphytc13*(1.-dfrt)*morpt
     &       + rtzoopc13*excr + rtdiazc13*morpt_D + rtdiazc13*nr_excr_D
     &       + rtdiazc13*morp_D*(1-(redntp/diazntp)) + rtdetrc13*remi
     &       + rtdoc13*recy_don - fcnpp*npp - fcnpp*npp_D)

        biodoc13 = biodoc13 + dtbio*redctn*(dfr*rtphytc13*morp
     &       + rtphytc13*dfrt*morpt - rtdoc13*recy_don)

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

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

        biodetrc13 = biodetrc13 + dtbio*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 + dtbio*redctn*(fcnpp*npp_D
     &       - rtdiazc13*morp_D - rtdiazc13*graz_D - rtdiazc13*morpt_D)
        expoout = expoout + expo
        rn15expoout = rn15expoout + rtdetrn15
        rc13expoout = rc13expoout + rtdetrc13
        grazout = grazout + graz
        morpout = morpout + morp
        morzout = morzout + morz
        graz_Det_out = graz_Det_out + graz_Det
        graz_Zout = graz_Zout + graz_Z
        nppout = nppout + npp
        npp_dopout = npp_dopout + npp*dopupt_flag
        morptout = morptout + morpt
        remiout = remiout + remi
        excrout = excrout + excr
        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
        expofeout = expofeout + expofe
        remifeout = remifeout + remife
!     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 (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 (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)
        if (dfeflag .eq. 1) dfeflag = 0.5 + sign(0.5,biodfe - trcmin)
        if (detrfeflag .eq. 1) detrfeflag = 0.5
     &                                 + sign(0.5,biodetrfe - trcmin)
        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 (diazc13flag .eq. 1)
     &       diazc13flag = 0.5 + sign(0.5,biodiazc13 - trcmin)
      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)
      bioout(6) = biodic - bioin(6)
      bioout(7) = biono3 - bioin(7)
      bioout(8) = biodon - bioin(8)
      bioout(9) = biodiaz - bioin(9)
      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)
      bioout(22) = biodfe - bioin(22)
      bioout(23) = biodetrfe - bioin(23)
      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)
      return
      end

