

      subroutine npzd_src (bioin, ntsb, tsb, gl, bct, impo, dzt
     &,                    dayfrac, wwd, rkw, nud, bioout, expoout
     &,                    grazout, morpout, morzout, graz_Det_out
     &,                    graz_Zout

     &,                    nppout, morptout, remiout, excrout




     &,                    npp_Dout, graz_Dout, morp_Dout, nfixout





     &,                    avej_out, avej_D_out, gmax_out, no3P_out
     &,                    po4P_out, po4_D_out
     &,                    jmax_out, jmax_D_out
     &,                    gd_out, gd_D_out
     &,                    bctz
     &,                    expofeout, impofe, remifeout
     &,                    feorgadsout, thetamaxout, deffeout
     &,                    feprimeout, fecolout
     &,                    o2
     &,                    chlout, thetachlout
     &,                    chl_D_out
     &                     ) 

!=======================================================================
!     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:4) = N,P,Z,D [mmol m-3]
!       bioin(5)   = nitrate [mmol m-3]
!       bioin(6)   = diazotrophs [mmol m-3]
!       bioin(7)   = dissolved iron [mmol m-3]
!       bioin(8)   = particulate iron [mmol m-3]
!{ wyao 20170710
! wyao 20170710}
!       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]
!{ wyao 20170710
! wyao 20170710}

!     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
!{ wyao 20170710
! new variables for including DON and DOP deducted by Ben
! wyao 20170710}
!{ wyao 20170712
! wyao 20170712}
!=======================================================================

      implicit none

      integer n, ntsb

      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 avej_out, avej_D_out, gmax_out, no3P_out, po4P_out, po4_D_out
      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, graz_Det_out, graz_Zout, thetaZ, ing_P, ing_D
      real ing_Det, ing_Z, g_Z, g_Det, graz_Z, graz_Det, gmax
      real no3P, po4P, po4_D, bctz
      real bct_out, jmax_out, jmax_D_out, gd_out, gd_D_out
      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, p1,p2,kfevar
      real expofeout, impofeout, feorgadsout, fecolout
      real remifeout, thetamaxout, deffeout, thetachlout
      real chlout, feprimeout, o2, chl_D_out
!{ wyao 20170710
! wyao 20170710}
!{ wyao 20170711
! wyao 20170711}

!======================= include file "size.h" =========================

!-----------------------------------------------------------------------
!     USER INPUT:
!-----------------------------------------------------------------------

!     imt    = number of grid points in the longitudinal direction
!              (calculated points are from 2 through imt-1. End points
!               are boundaries)
!     jmt    = number of grid points (latitude rows) in the latitudinal
!              direction (calculated points are from 2 through jmt-1.
!              End points are boundaries)
!     km     = number of grid points in the vertical direction
!              (calculated points are from 1 through km)
!     nt     = number of tracers (temperature, salinity, ...)
!     nsrc   = number of tracer with sources
!     kpzd   = depth for limited npzd model
!     jmz    = size for "unrotated" zonal averages
!     jmzm1  = jmz minus one
!     mnisle = maximum number of islands (unconnected land masses)
!     maxipp = maximum number of all island perimeter points
!-----------------------------------------------------------------------

      integer imt, jmt, km, nt, nsrc, kpzd, nat, jmz, jmzm1, mnisle
      integer maxipp, jmw, jsmw, jemw

      parameter (imt= 1, jmt= 1, km= 19)
      parameter (nt=2
     $             +1
     $             +1
     $             +1
     $             +4
     $             +2
     $             +2
     $               )
      parameter (nsrc=0
     $               +1
     $               +1
     $               +1
     $             +4
     $             +2
     $               +2 
     $                 )
      parameter (kpzd=km)

      parameter (nat=2
     $, jmz=jmt, jmzm1=jmz-1)
      parameter (mnisle=50, maxipp=5000)

      parameter (jmw=jmt)

!-----------------------------------------------------------------------
!     set first and last calculated row within the MW. other rows
!     are used as buffers
!-----------------------------------------------------------------------

!     jsmw   = 1st calculated row within the MW
!     jemw   = last calculated row within the MW
      parameter (jsmw=1, jemw=1)
! Moses-Triffid land model

! POINTS = Maximum number of points in grid.
! STEPSM = Maximum number of timesteps in a day.
! klmax = maximum ocean depth levels over which the land model can exist

      integer POINTS, STEPSM, klmax
      parameter (POINTS=14300, STEPSM=24, klmax=0)

! NNVG  = Number of non-vegetation surface types.
! NPFT  = Number of plant functional types.
! NTYPE = Number of surface types.
! SOIL  = Index of the surface type 'Soil'
! Land surface types :
!     1 - Broadleaf Tree
!     2 - Needleleaf Tree
!     3 - C3 Grass
!     4 - C4 Grass
!     5 - Shrub
!     6 - Soil

      integer NNVG, NPFT, NTYPE, SOIL
      parameter (NNVG=4, NPFT=5, NTYPE=6, SOIL=6)
!======================= include file "param.h" ========================

!     main parameter file which sets ocean characteristics:

!     nvar   = number of prognostic variables
!     lseg   = maximum number of longitudinal stream function segments
!     nlatpr = maximum number of latitudes for matrix printouts
!              on diagnostic time steps
!     nhreg  = number of regions in the horizontal used for averaging
!              tracers.
!     nvreg  = number of regions in the vertical used for term balance
!              calculations. note "nvreg" is not used for tracer
!              averages
!     numreg = total number of regions ( = product of nhreg & nvreg)
!              used for term balance calculations

!     nvarbh = number of prognostic variables using biharmonic mixing

!     ncrows = number of calculated rows within the MW.
!              (the remaining rows are buffer rows).

!     rstd   = standard c14/c12 ratio

      integer lseg, nlatpr, nhreg, nvreg, numreg, nvar, nvarbh
      integer imtm1, kmm1, imtp1, imtm2, jmtp1, jmtm1, jmtm2, jscan
      integer kmp1, kmp2, imtkm, nwds, nkflds, nslab, ntmin2, ncrows

      real rstd

      parameter (lseg=5, nlatpr=10)
      parameter (nhreg=3, nvreg=1, numreg=nhreg*nvreg)
      parameter (nvar=nt+2)


      parameter (imtm1=imt-1, kmm1=km-1)
      parameter (imtp1=imt+1, imtm2=imt-2
     &,          jmtp1=jmt+1, jmtm1=jmt-1, jmtm2=jmt-2
     &,          jscan=jmtm2
     &,          kmp1=km+1, kmp2=km+2
     &,          imtkm=imt*km, nwds=imt*jmt, nkflds=2
     &,          nslab=imt*nvar*km, ntmin2=nt+1/nt)

      parameter (ncrows = jmw - 3 + jmw/jmt)
!====================== include file "pconst.h" ========================

!     rules for parameter constants

!     use prefix of "c" for whole real numbers (eg: c57 for 57.0)
!     use "m" after prefix to designate negative values (minus sign)
!       (eg: cm7 for -7.0)
!     use prefix of "p" for non repeating fractions (eg: p5 for 0.5)
!     use prefix of "r" for reciprocals (eg: r3 for 1/3.0)
!     combine use of prefix above and "e" for scientific notation, with
!       (eg: c5e4 for 5.0e4, c1em10 for 1.0e-10)

      real c0, c1, c2, c3, c4, c5, c7, c8, c14, c16, c360, p125, p25
      real p5, p75, epsln, c24, c60, c1440, r24, r60, r1440, secday

      parameter (c0=0.0, c1=1.0, c2=2.0, c3=3.0, c4=4.0, c5=5.0, c7=7.0)
      parameter (c8=8.0)
      parameter (c14=14.0, c16=16.0, c360=360.0)
      parameter (p125=0.125, p25=0.25, p5=0.5, p75=0.75)
      parameter (epsln=1.0e-20)
      parameter (c24=24.0, c60=60.0, c1440=1440.0)
      parameter (r24=c1/c24, r60=c1/c60, r1440=c1/c1440)
      parameter (secday=c1/(c60*c1440))
!====================== include file "stdunits.h" ======================

!     stdin  = unit number for standard input.
!     stdout = unit number for standard output.
!     stderr = unit number for standard error.

      integer stdin, stdout, stderr
      parameter (stdin = 5, stdout = 6, stderr = 6)
!====================== include file "calendar.h" ======================

!                       calendar specification arrays

!-----------------------------------------------------------------------
!     eqyear  = true to select a calendar in which each year
!               has the same number of days (i.e., no leap years)
!               false selects a Julian calendar

!     eqmon   = true to force all months to have the same number of days
!               false => the usual 31, 28, 31, 30, ..., days per month.
!               only used when eqyear = true

!     dayname = character names of days

!     monname = character names of months

!     monlen  = the length of each month (in days) when eqmon is true

!     yrlen   = the length of a typical (non-leap) year in days

!     daypm   = array of month lengths in days   (non-leap)

!     msum    = array of cumulative days preceding each month
!               (again, non-leap)

!     daylen  = the length of a day in seconds
!-----------------------------------------------------------------------

      character(10) :: dayname
      character(12) :: monname
      common /calen_c/ dayname(7), monname(12)

      logical eqyear, eqmon
      common /calen_l/ eqyear, eqmon

      integer daypm, msum, yrlen, monlen
      common /calen_i/ daypm(12), msum(12), yrlen, monlen

      real daylen
      common /calen_r/ daylen
!====================== include file "npzd.h" =========================

!   variables for npzd model

!   ntnpzd   = number of npzd tracers
!   nbio     = number of npzd timesteps per ocean timestep
!   trcmin   = minimum tracer for flux calculations
!   alpha    = initial slope P-I curve [(W/m^2)^(-1)/day]
!   kw       = light attenuation due to water [1/m]
!   kc       = light attenuation by phytoplankton [1/(m*mmol m-3)]
!   ki       = light attenuation through sea ice & snow
!   abio     = maximum growth rate parameter [1/day]
!   bbio     = b
!   cbio     = [1/deg_C]
!   k1n      = half saturation constant for N uptake [mmol m-3]
!   gamma1   = assimilation efficiency (zpk)
!   gbio     = maximum grazing rate at 0 deg C [day-1]
!   nuz      = quadratic mortality (zpk)
!   nud0     = remineralization rate [day-1]
!{ wyao 20170707
! wyao 20170707}
!   LFe      = Iron limitation
!   wd       = sinking speed of detritus [m day-1]
!   ztt      = depth to top of grid cell [cm]
!   rkwz     = reciprical of light attenuation times grid depth
!   par      = fraction of photosythetically active radiation
!   dtnpzd   = time step of biology
!   capr     = carbonate to carbon production ratio
!   dcaco3   = remineralisation depth of calcite [cm]
!   rcak     = array used in calculating calcite remineralization
!   rcab     = array used in calculating bottom calcite remineralization
!{ wyao 20170707
!   nup    = quick recycling rate(Phytoplankton)[1/day]
!   nupt0    = quick recycling rate(Temprtr depnd.Phytoplankton)[1/day]
! wyao 20170707}
!   wd0      = sinking speed of detritus at surface [m/day]
!   k1p      = half saturation constant for P uptake
!   jdiar    = factor reducing the growth rate of diazotrophs
!   redctn   = C/N Redfield ratio (includes mol to mmol conversion)
!   redctp   = C/P Redfield ratio (includes mol to mmol conversion)
!   redptn   = P/N Redfield ratio
!   redntp   = N/P Redfield ratio
!   redotn   = O/N Redfield ratio (includes mol to mmol conversion)
!   redotp   = O/P Redfield ratio (includes mol to mmol conversion)
!   rnbio    = reciprical of nbio
!   rdtts    = reciprical of dtts [s-1]
!   dtbio    = npzd time step [s]
!   rnpp     = rate of net primary production [nmol cm-3 s-1]
!   rgraz    = rate of grazing [nmol cm-3 s-1]
!   rmorp    = rate of mortality of phytoplankton [nmol cm-3 s-1]
!{ wyao 20170707
!   rmorpt    = rate of quick recy of phytoplankton [nmol cm-3 s-1]
! wyao 20170707}
!   rmorz    = rate of mortality of zooplankton [nmol cm-3 s-1]
!   rremi    = rate of remineralization [nmol cm-3 s-1]
!   rexcr    = rate of excretion [nmol cm-3 s-1]
!   rexpo    = rate of export through the bottom [nmol cm-3 s-1]
!   rnpp_D   = npp for diazotraphs [nmol cm-3 s-1]
!   rgraz_D  = rgraz for diazotraphs [nmol cm-3 s-1]
!{ wyao 20170707
!   rmorp_D  = rate of mortality for diazotraphs [nmol cm-3 s-1]
! wyao 20170707}
!   rnfix    = rate of nitrogen fixation [nmol cm-3 s-1]
!   rdeni    = rate of denitrification [nmol cm-3 s-1]
!   kzoo     = half saturation constant for Z grazing
!   zprefP   = Z preference for grazing on P
!   zprefD   = Z preference for grazing on Diaz
!   zprefZ   = Z preference for grazing on other Z
!   zprefDet = Z preference for grazing on Detritus
!   rgraz_Det = rate of grazing on Detritus [nmol cm-3 s-1]
!   rgraz_Z   = rate of grazing on other Zooplankton [nmol cm-3 s-1]
!   geZ      = growth efficiency of zooplankton
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! New diagnostic output
!   ravej      = light-dependant growth rate of P
!   ravej_D    = light-dependant growth rate of Diaz
!   rgmax      = temp-dependant growth rate of zoo
!   rno3P      = nitrate-dependant growth rate of P
!   rpo4P       = phosphate-dependant growth rate of P
!   rpo4_D     = phosphate-dependant growth rate of D
!
!   fe_dissolved = dissolved iron concentration
!   kfe = Fe limitation half saturation parameter
!   kfe_D = Fe limitation half sat. param. for diaz.
!   kmfe = number of depth levels for the iron mask
!	
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Dynamic iron cycle
!   kfeleq = Fe-ligand stability constant [m^3/(mmol ligand)]
!   lig = Ligand concentration  [mmol/m^3]
!   thetamaxhi = Maximum Chl:C ratio, abundant iron [gChl/(gC)]
!   thetamaxlo = Maximum Chl:C ratio, extreme iron limitation [gChl/(gC)]
!   alphamax = Maximum initial slope in PI-curve [day^-1 (W/m^2)^-1 (gChl/(gC))^-1]
!   alphamin = Minimum intital slope in PI-curve [day^-1 (W/m^2)^-1 (gChl/(gC))^-1]
!   mc = Molar mass of carbon [g/mol]
!   fetopsed = Fe:P ratio for sedimentary iron source [molFe/molP]
!   o2min = Minimum O2 concentration for aerobic respiration [mmolO_2/m^3]
!   kfeorg = Organic-matter dependent scavenging rate [(m^3/(gC s))^0.58]
!   kfeorg_ps = Organic-matter dependent scavenging power lay scale parameter default 0.58
!   kfecol = Colloidal Fe production rate (inorganic scavenging) [s^-1]
!   fehydro = hydrothermal flux of iron [mol/m^3*s^-1]
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! DOP and DON guessing by Ben
!	option O_npzd_dop & O_npzd_don
!	dfr = refractory/semi-labile DOM fraction from mortality.
!	dfrt = refractory/semi-labile DOM fraction from quick recy.
!	redotc = o/c redfield ratio
!	redntc = n/c redfield ratio
!	redotc = b/c redfield ratio
!	redptc = p/c redfield ratio
!	hdop = DOP growth rate handicap
!	rnpp_dop = rate of net primary production from DOP [nmol cm-3 s-1]
!	rnpp_D_dop = rate of net primary production for diaz from DOP [nmol cm-3 s-1]
!	diazntp = n/p diazo ratio
!	diazptn = p/N diazo ratio
!	
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Sinking speed parameters
!	option O_npzd_caped_sinkspeed
!	wdi		the dependency on depth
!	wdit	the handicaped depth, fixed speed after certain depth. 
!
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! non-redfield iron (diaz and phyt have different uptake fe to N ratio)
!	option O_npzd_iron_diazratio
!
!	rfeton_D Uptake ratio of iron to nitrogen for Diaz


      integer ntnpzd, nbio, kmfe	  
      parameter (ntnpzd = 4
     $                  + 2
     $                  + 2
     $                    )

      common /npzd_i/ nbio(km)

      real trcmin
      parameter (trcmin=5e-12)

      real alpha, kw, kc, ki, abio, bbio, cbio, k1n, nup, gamma1, gbio
      real epsbio, nuz, gamma2, nud0, LFe, wd, ztt, rkwz, par, dtnpzd
      real capr, dcaco3, rcak, rcab, nupt0, wd0, k1p, jdiar, redctn
      real redctp, redptn, redntp, redotn, redotp, rnbio, rdtts, dtbio
      real rnpp, rgraz, rmorp, rmorpt, rmorz, rremi, rexcr, rexpo
!{ wyao 20170707
	  real wdi, wdit
! wyao 20170707}
      real rnpp_D, rgraz_D, rmorp_D, rnfix, rdeni, kzoo, zprefP
      real zprefD, zprefZ, zprefDet, rgraz_Det, rgraz_Z, geZ, kfe
      real ravej, ravej_D, rgmax, rno3P, rpo4P, rpo4_D, kfe_D
      real kfemax, kfemin, pmax


      real kfeleq, lig, thetamaxhi, alphamax, alphamin
      real thetamaxlo, mc, fetopsed, o2min, kfeorg, rfeton
      real kfecol
      real kfeorg_ps
      real rfeorgads, rdeffe, rremife, rexpofe, rfeprime
      real rfesed, rbfe, rfecol
      real fe_hydr
      real rchl, rchl_D 
	  real fehydro, rfehydro

      common /npzd_r/ alpha, kw, kc, ki, abio, bbio, cbio, k1n, nup
      common /npzd_r/ gamma1, gbio, epsbio, nuz, nud0, LFe
      common /npzd_r/ wd(km), ztt(km), rkwz(km), par, dtnpzd, capr
      common /npzd_r/ dcaco3, rcak(km), rcab(km), nupt0, wd0, k1p
      common /npzd_r/ jdiar, redctn, redctp, redptn, redntp, redotn
      common /npzd_r/ redotp, rnbio(km), rdtts(km), dtbio(km), geZ
      common /npzd_r/ kzoo, zprefP, zprefD, zprefZ, zprefDet, kfe, kfe_D
      common /npzd_r/ kfeleq, alphamax, alphamin
      common /npzd_r/ thetamaxhi, thetamaxlo, lig, fetopsed, o2min
      common /npzd_r/ mc, kfeorg, rfeton, kfecol, kfeorg_ps
      common /npzd_r/ kfemax, kfemin, pmax, fe_hydr(imt,jmt,km)
!{ wyao 20170711
! wyao 20170711}
!{ wyao 20170707
	  common /npzd_r/ wdi
! wyao 20170707}
      common /npzd_r/ rnpp(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rgraz(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rmorp(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rmorpt(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rmorz(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rexcr(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rremi(imt,km,jsmw:jemw)
      common /npzd_r/ rexpo(imt,km,jsmw:jemw)
      common /npzd_r/ rgraz_Det(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rgraz_Z(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rnpp_D(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rgraz_D(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rmorp_D(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rnfix(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rdeni(imt,km,jsmw:jemw)
      common /npzd_r/ rremife(imt,km,jsmw:jemw)
      common /npzd_r/ rexpofe(imt,km,jsmw:jemw)
! wyao 20170427{
	  common /npzd_r/ fehydro(imt,jmt,km)
	  common /npzd_r/ rfehydro(imt,km,jsmw:jemw)	  
! wyao 20170427}
      common /npzd_r/ rfeorgads(imt,km,jsmw:jemw)
      common /npzd_r/ rdeffe(imt,km,jsmw:jemw)
      common /npzd_r/ rfeprime(imt,km,jsmw:jemw)
      common /npzd_r/ rfesed(imt,km,jsmw:jemw)
      common /npzd_r/ rbfe(imt,km,jsmw:jemw)
      common /npzd_r/ rfecol(imt,km,jsmw:jemw)
      common /npzd_r/ rchl(imt,km,jsmw:jemw)
      common /npzd_r/ rchl_D(imt,km,jsmw:jemw)
	  real rbct, rjmax, rjmax_D, rgd, rgd_D
	  common /npzd_r/ rbct(imt,kpzd,jsmw:jemw)
	  common /npzd_r/ rjmax(imt,kpzd,jsmw:jemw)
	  common /npzd_r/ rjmax_D(imt,kpzd,jsmw:jemw)
	  common /npzd_r/ rgd(imt,kpzd,jsmw:jemw)
	  common /npzd_r/ rgd_D(imt,kpzd,jsmw:jemw)
      common /npzd_r/ ravej(imt,kpzd,jsmw:jemw)
      common /npzd_r/ ravej_D(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rgmax(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rno3P(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rpo4P(imt,kpzd,jsmw:jemw)
      common /npzd_r/ rpo4_D(imt,kpzd,jsmw:jemw)



      real bioin(ntnpzd), bioout(ntnpzd)


!{ wyao 20170710
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!here the tracers' start values are assigned to the private 
!variables of npzd_src(). 
!
!
!
! wyao 20170710}
   
      bioout(:) = 0.0
      bion = bioin(1)
      biop = bioin(2)
      bioz = bioin(3)
      biod = bioin(4)
      biono3 = bioin(5)
      biodiaz = bioin(6)
      biodfe = bioin(7)
      biodetrfe = bioin(8)
!{ wyao 20170710
! wyao 20170710}



      
!{ wyao 20170710
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!here starts the code block for paramter checking and calculation of 
!non-npzd_time-dependent parameters
!
!
! wyao 20170710}

!{ wyao 20170710
      p1=min(max(trcmin,biop),pmax)
      p2=max(0.0,max(trcmin,biop)-pmax)
      kfevar=(kfemin*p1+kfemax*p2)/(p1+p2) 
      deffe=max(trcmin,biodfe)/(kfevar+max(trcmin,biodfe))
      thetamax=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe
      alpha_O=alphamin+(alphamax-alphamin)*deffe
      gl_O=gl*thetamax*alpha_O
      deffe_D=max(trcmin,biodfe)/(kfe_D+max(trcmin,biodfe))
      thetamax_D=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe_D
      alpha_D=alphamin+(alphamax-alphamin)*deffe_D
      gl_D=gl*thetamax_D*alpha_D
! wyao 20170710} this part is for the calculation of avej and avej_D, 
! which are ocean timestep dependent, not npzd timestep dependent.


!     photosynthesis after Evans & Parslow (1985)
!     notation as in JGOFS report No. 23 p. 6
      f1 = exp((-kw - kc*(max(trcmin,biop)
     &     +max(trcmin,biodiaz)))*dzt)
! wyao 20170710} not sure if the attenuation formular is right, since 
! it differs between levin's and chris'. f1 is not used anyway.


      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*(max(trcmin,biop)
     &       + max(trcmin,biodiaz)))*dzt)
! Make the max grazing rate a function of temperature
      gmax = gbio*bctz
! Note that bctz 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*(max(trcmin,biop)
     &         +max(trcmin,biodiaz)))*dzt)
!     check grazing preferences = 1 for N case
        IF ((zprefP + zprefDet + zprefZ + zprefD).ne.1) THEN
           zprefP = 0.30
           zprefZ = 0.30
           zprefDet = 0.30
           zprefD = 0.10
        END IF
      nupt = nupt0*bct
!{ wyao 20170710
! wyao 20170710}


!{ wyao 20170710
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!here starts the code block for zeroing out all the output variables
!
!
! wyao 20170710}
      expoout = 0.0
      grazout = 0.0
      morpout = 0.0
      morzout = 0.0
      graz_Det_out = 0.0
      graz_Zout = 0.0
      nppout = 0.0
      morptout = 0.0
      remiout = 0.0
      excrout = 0.0
!{ wyao 20170710
! wyao 20170710}
      npp_Dout = 0.0
      graz_Dout = 0.0
      morp_Dout = 0.0
      nfixout = 0.0
!{ wyao 20170710
! wyao 20170710}
      expofeout = 0.0
      remifeout = 0.0
      feorgadsout = 0.0
      deffeout = 0.0
      feprimeout = 0.0
      fecolout = 0.0
      chlout = 0.0
      chl_D_out = 0.0
      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
      bct_out = 0.0
      jmax_out = 0.0
      jmax_D_out = 0.0
      gd_out = 0.0
      gd_D_out = 0.0


!{ wyao 20170710
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!here starts the code block for calculation of source terms in each 
!npzd timestep.
!
!
! wyao 20170710}
      do n=1,ntsb

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

!       growth rate of phytoplankton
!{ wyao 20170710
        u_P = min(avej, jmax*max(trcmin,bion)/(k1p + max(trcmin,bion)))
! wyao 20170710}
        po4P = max(trcmin,bion)/(k1p + max(trcmin,bion)) 
        

!       nitrate limitation
        u_P = min(u_P, jmax*max(trcmin,biono3)/(k1n
     &        + max(trcmin,biono3)))
        no3P = max(trcmin,biono3)/(k1n + max(trcmin,biono3))


!       growth rate of diazotrophs smaller than other phytoplankton and
!       not nitrate limited
!{ wyao 20170710
        u_D = min(avej_D, jmax_D*max(trcmin,bion)/(k1p
     &        + max(trcmin,bion)))
! wyao 20170710}
        po4_D = jmax_D*max(trcmin,bion)/(k1p + max(trcmin,bion))
        
        
!       Set the grazing coefficients for the N case
        thetaZ = zprefP*max(trcmin,biop)+zprefDet*max(trcmin,biod)
     &           +zprefZ*max(trcmin,bioz)+zprefD*max(trcmin,biodiaz)
     &         + kzoo
        ing_P = zprefP/thetaZ
        ing_Det = zprefDet/thetaZ
        ing_Z = zprefZ/thetaZ
        ing_D = zprefD/thetaZ
        npp = u_P*max(trcmin,biop) 
!{ wyao 20170710
! wyao 20170710}        
        npp_D = max(0.,u_D*biodiaz) 
!       grazing on diazotrophs
        g_D = gmax*ing_D*max(trcmin,biodiaz)
        graz_D = g_D*max(trcmin,bioz)

!{ wyao 20170710
        morp_D = nupt*biodiaz ! linear mortality. in levin's code, he use quick recycling rate for the diaz mortality, not sure why.
! wyao 20170710}
        no3upt_D = max(trcmin,biono3)/(k1n + max(trcmin,biono3))*npp_D ! nitrate uptake 
!       grazing on P
        g_P = gmax*ing_P*biop
        graz = g_P*bioz
!       grazing on Z
        g_Z = gmax*ing_Z*bioz
        graz_Z = g_Z*bioz
!       grazing on Detritus
        g_Det = gmax*ing_Det*biod
        graz_Det = g_Det*bioz
!
        morp = nup*biop
        morpt = nupt*biop
!{ wyao 20170710
! wyao 20170710}
        morz = nuz*bioz*bioz
        remi = nud*bct*biod
        expo = wwd*biod
!  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

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

          feorgads = kfeorg*((max(trcmin,biod)*mc*redctn)**kfeorg_ps)
     &         *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

!   Calculation of average light in mixed layer for calculation of Chl diagnostic
!        irrtop=gl/2/par
!        kirr=-kw - kc*(bioin(2)+bioin(6))
!        aveirr=-1/dzt/kirr*(irrtop-irrtop*exp(kirr*dzt))
!   Diagnostic calcalation of Chl:C ratio to calculate Chl concentration
!        thetachl=thetamax/(1+thetamax*aveirr/2/u_P)      !Chl:C after Geider et al. (1997), note that gl includes alpha and PAR already
!{ wyao20180508 change o_chl to the nickelsen 2015
         thetachl=thetamax
!        thetachl_D=thetamax_D/(1+thetamax_D*aveirr/2/u_D)
         thetachl_D=thetamax_D
!wyao20180508 }
!       flags prevent negative values by setting outgoing fluxes to
!       zero if tracers are lower than trcmin
        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)
!{ wyao 20170712        
! wyao 20170712}
        dfeflag = 0.5 + sign(0.5,biodfe - trcmin)
        detrfeflag = 0.5 + sign(0.5,biodetrfe - trcmin)
        no3flag = 0.5 + sign(0.5,biono3 - trcmin)
        diazflag = 0.5 + sign(0.5,biodiaz - trcmin)
        graz = graz*pflag*zflag
        graz_Z = graz_Z*zflag
        graz_Det = graz_Det*dflag*zflag
        morp = morp*pflag
        morpt = morpt*pflag
        morz = morz*zflag
        remi = remi*dflag
        expo = expo*dflag
!{ wyao 20170710        
! wyao 20170710}
!{ wyao 20170710        
! wyao 20170710}
        graz_D = graz_D*diazflag*zflag
        morp_D = morp_D*diazflag
        no3upt_D = no3upt_D*no3flag*diazflag
        remife=remife*detrfeflag
        feorgads=feorgads*dfeflag
        expofe=expofe*detrfeflag
        fecol=fecol*dfeflag
        avej = avej*pflag
        avej_D = avej_D*diazflag
        gmax = gmax*zflag
        no3P = no3P*no3flag
        po4P = po4P*nflag
        po4_D = po4_D*nflag
        bct = bct
        jmax=jmax*pflag
        jmax_D=jmax_D*diazflag
        gd = gd*pflag
        gd_D=gd_D*diazflag

!   Excretion is the difference between
!   the zooplankton assimilation and growth efficiencies
!{ wyao 20170710        
        excr = (gamma1-geZ)*(graz+graz_Z+graz_Det+graz_D)
! wyao 20170710}

!{ wyao 20170710        

!       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*(geZ*(graz + graz_D + graz_Det + graz_Z) 
     &       - morz - graz_Z)
!       detritus equation
        biod = biod + tsb*((1.-gamma1)*(graz + graz_D + graz_Det 
     &       + graz_Z) + morp + morp_D + morz - remi - graz_Det 
     &       - 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)
! wyao 20170710}







!{ wyao 20170710        
!       dissolved iron equation
        biodfe = biodfe + tsb*(rfeton*(excr - (npp + npp_D) + morpt)
     &         - feorgads + remife - fecol)
!       particulate iron equation
        biodetrfe = biodetrfe + tsb*(rfeton*((1.-gamma1)*(graz + graz_D
     &            + graz_Det
     &            + graz_Z) + morp + morp_D + morz - graz_Det)
     &            + feorgads + fecol - remife - expofe + impofe)

! wyao 20170710}

        expoout = expoout + expo
        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
!{ wyao 20170710        
! wyao 20170710}
        morptout = morptout + morpt
        remiout = remiout + remi
        excrout = excrout + excr
        npp_Dout = npp_Dout + npp_D
!{ wyao 20170710        
! wyao 20170710}
        graz_Dout = graz_Dout + graz_D
        morp_Dout = morp_Dout + morp_D
        nfixout = nfixout + npp_D - no3upt_D
        expofeout = expofeout + expofe
        remifeout = remifeout + remife
        feorgadsout = feorgadsout + feorgads
        deffeout = deffeout + deffe
        feprimeout = feprimeout + feprime
        fecolout = fecolout + fecol
        chlout = chlout + biop * redctn * mc * thetachl
        chl_D_out = chl_D_out + biodiaz * redctn * mc * thetachl_D
        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
        bct_out = bct_out + bct
        jmax_out = jmax_out + jmax
        jmax_D_out = jmax_D_out + jmax_D
        gd_out = gd_out + gd
        gd_D_out = gd_D_out + gd_D
      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) = biodfe - bioin(7)
      bioout(8) = biodetrfe - bioin(8)
!{ wyao 20170710        
! wyao 20170710}   
      return
      end
