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 & ) !======================================================================= ! 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] !{ 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 !{ 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 $ ) parameter (nsrc=0 $ +1 $ +1 $ +1 $ +4 $ +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 $ ) 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 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 !{ 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) 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) !{ 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 ! 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 gd = jmax*dayfrac u1 = max(gl/gd,1.e-6) u2 = u1*f1 ! for the following approximation ensure that u1 < 20 phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1 phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2 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))*jdiar gd_D = max(1.e-14,jmax_D*dayfrac) u1 = max(gl/gd_D,1.e-6) u2 = u1*f1 ! for the following approximation ensure that u1 < 20 phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1 phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2 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} 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 ! 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 ! 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} 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 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} 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 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) !{ wyao 20170710 ! wyao 20170710} return end