! source file: /Users/csomes/Research/Models/UVic_ESCM/2.9/sim/nonredfield_dom/clean_up_code/seaice_wind_bug/afs_high/updates/npzd_src.F subroutine npzd_src (bioin, ntsb, tsb, gl, bct, impo, dzt &, dayfrac, wwd, rkw, nud, bioout, expoout &, grazout, morpout, morzout, graz_Det_out &, graz_Zout, nudop, nudon &, nppout, npp_dopout, morptout, remiout &, excrout &, npp_Dout, npp_D_dopout, graz_Dout, morp_Dout &, morpt_Dout, nfixout &, felimit, felimit_D &, 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:5) = N,DON,P,Z,D [mmol m-3] ! bioin(6) = nitrate [mmol m-3] ! bioin(7) = DON [mmol m-3] ! bioin(8) = diazotrophs [mmol m-3] ! gl = 2.*light at top of grid box ! ntsb = number of time steps ! tsb = time step [s] ! bct = bbio**(cbio*temperature) ! impo = import of detritus from above [mmol m-3] ! dzt = depth of grid box [cm] ! dayfrac = day length (fraction: 0 < dayfrac < 1) ! wwd = sinking speed of detritus/dzt ! rkw = reciprical of kw*dzt(k) ! nud = remineralisation rate of detritus [s-1] ! nudop = remineralisation rate of dop [s-1] ! nudon = remineralisation rate of don [s-1] ! output variables: ! bioout = change from bioin [mmol m-3] ! nppout = net primary production [mmol m-3] ! npp_dopout = net primary production from DOP uptake [mmol m-3] ! grazout = grazing [mmol m-3] ! morpout = quadratic mortality of phytoplankton [mmol m-3] ! morptout = specific mortality of phytoplankton [mmol m-3] ! morzout = mortality of zooplankton [mmol m-3] ! remiout = remineralisation [mmol m-3] ! excrout = excretion [mmol m-3] ! expoout = detrital export [mmol m-3] ! npp_Dout = NPP of diazotrophs ! npp_D_dopout= NPP from DOP uptake of diazotrophs [mmol m-3] ! graz_Dout = grazing of diazotrophs ! morp_Dout = mortality of diazotrophs ! nfixout = rate of N2 fixation ! graz_Det_out = grazing of detritus ! graz_Zout = grazing on othe zooplankton ! avej_out = light-depend phyt. growth rate ! avej_D_out = light-depend Diaz growth rate ! gmax_out = temp-depend. zoo growth rate ! no3P_out = no3 depend. phyt growth rate ! po4P_out = po4 depend. phyt growth rate ! po4_D_out = po4 depend. Diaz growth rate ! New grazing formulation variables and parameters ! The following terms determine ingestion according to a ! a Holling II curve (i.e. Michaelis Menten): ! ! Ingestion = max_graz_rate * (Ft/(Ft + kzoo)) ! ! where Ft is the weighted measure of the total food available ! and equals the sum of the different prey types times the ! preference of Z for that type of prey ! ! zprefP = Z preference for P ! zprefD = Z preference for Diaz ! zprefDet = Z preference for detritus ! zprefZ = Z preference for other Z ! kzoo = half saturation coefficienct for Z ingestion mmol N m-3 ! ing_P = zooplankton ingestion of phytoplankon ! ing_D = zooplankton ingestion of diazotrophs ! ing_Det = zooplankton ingestion of detritus ! ing_Z = zooplankton ingestion of other zooplankton ! thetaZ = Michaelis-Menten denominator ! ! felimit = Fe limitation parameter ! felmit_D = Fe limitation parameter for diazotrophs ! !======================================================================= implicit none integer n, ntsb real gl, f1, biopo4, biophyt, biozoop, biodetr, jmax, u_P, g_P real morp, morpt, morz, remi, excr, expo, impo, nppout, grazout real morpout, morptout, morzout, remiout, excrout, expoout, tsb real avej_out, avej_D_out, gmax_out, no3P_out, po4P_out, po4_D_out real dzt, po4flag, phytflag, zoopflag, detrflag, wwd, rkw, gd real nupt, nud, biono3, u_D,npp_D, npp_Dout, no3flag, biodiaz, bct real diazflag, g_D,graz_D, morpt_D, jmax_D, gd_D, avej_D, no3upt_D real morpt_Dout, graz_Dout, nfixout, biop2, u1, u2, phi1, phi2 real avej, graz_Det_out, graz_Zout, thetaZ, ing_P, ing_D, dayfrac real ing_Det, ing_Z, g_Z, g_Det, graz_Z, graz_Det, gmax real no3P, po4P, po4_D, felimit, bctz, felimit_D, nudop, nudon real excrdiaz, biodop, biodon, dopflag, donflag, recy real recy_p, npp, graz, biodin15, biodon15, biophytn15, biozoopn15 real biodetrn15, biodiazn15, bassim, fcassim, bexcr, fcexcr real bnfix, fcnfix, rn15impo, rn15expoout, dig, dig_P, dig_Z real dig_Det, dig_D, excr_P, excr_Z, excr_Det, excr_D, sf, sf_P real sf_Z, sf_Det, sf_D, nr_excr_D, uno3, rno3, rzoop real rtdin15, rtphytn15, rtzoopn15, rtdetrn15, rtdiazn15, rtdon15 real din15flag, dopupt_flag, dopupt_D_flag, npp_dopout real don15flag, phytn15flag, zoopn15flag, detrn15flag, diazn15flag real limP, limP_D, uptdop_D, dopupt_D, morp_D, nupt_D, morp_Dout real udon, rdon, brecy, fcrecy, biodic, uptdop, dopupt real biodic13, biophytc13, biozoopc13, biodetrc13, biodiazc13 real dic13flag, phytc13flag, zoopc13flag, detrc13flag, diazc13flag real biodoc13, doc13flag, biodoc, npp_D_dopout real rdic13, rtphytc13, rtzoopc13, rtdetrc13, rtdiazc13, ac13b real rc13impo, rc13expoout, dicflag, bc13npp, rtdic13, fcnpp real rtdoc13 include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "calendar.h" include "npzd.h" real bioin(ntnpzd), bioout(ntnpzd) bioout(:) = 0.0 biopo4 = bioin(1) biodop = bioin(2) biophyt = bioin(3) biozoop = bioin(4) biodetr = bioin(5) biodic = bioin(6) biono3 = bioin(7) biodon = bioin(8) biodiaz = bioin(9) ! 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 = 1.0 don15flag = 1.0 phytn15flag = 1.0 zoopn15flag = 1.0 detrn15flag = 1.0 diazn15flag = 1.0 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) ! photosynthesis after Evans & Parslow (1985) ! notation as in JGOFS report No. 23 p. 6 f1 = exp((-kw - kc*(biophyt+biodiaz))*dzt) ! In the following "felimit" is determined by an iron mask and ! is used to limit phytoplankton growth in HNLC regions jmax = abio*bct*felimit 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*(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*felimit_D)*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*(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 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 do n=1,ntsb ! growth rate of phytoplankton if (hdop*biodop/(k1p + biodop) .gt. biopo4/(k1p + biopo4)) then dopupt_flag = 1. limP = biodop/(k1p + biodop) u_P = min(avej, jmax*limP*hdop) else dopupt_flag = 0. limP = biopo4/(k1p + biopo4) u_P = min(avej, jmax*limP) endif po4P = jmax*biopo4/(k1p + biopo4) ! 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 no3upt_D = biono3/(k1n + biono3)*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 = nudon*bct*biodon recy_p = nudop*bct*biodop morz = nuz*biozoop*biozoop remi = nud*bct*biodetr expo = wwd*biodetr graz = graz*phytflag*zoopflag*phytn15flag*zoopn15flag graz_Z = graz_Z*zoopflag*zoopn15flag graz_Det = graz_Det*detrflag*zoopflag*detrn15flag*zoopn15flag morp = morp*phytflag*phytn15flag morpt = morpt*phytflag*phytn15flag morz = morz*zoopflag*zoopn15flag remi = remi*detrflag*detrn15flag expo = expo*detrflag*detrn15flag recy_p = recy_p*dopflag if (dopupt_flag .eq. 1) then npp = npp*dopflag*no3flag*din15flag else npp = npp*po4flag*no3flag*din15flag endif if (dopupt_D_flag .eq. 1) then npp_D = npp_D*dopflag else npp_D = npp_D*po4flag endif graz_D = graz_D*diazflag*zoopflag*diazn15flag*zoopn15flag morpt_D = morpt_D*diazflag*diazn15flag morp_D = morp_D*diazflag*diazn15flag no3upt_D = no3upt_D*no3flag*din15flag recy = recy*donflag*don15flag ! Excretion is the difference between ! the zooplankton assimilation and growth efficiencies c excr_P = (gamma1-geZ)*graz c excr_Z = (gamma1-geZ)*graz_Z c excr_Det = (gamma1-geZ)*graz_Det c excr_D = (gamma1-geZ)*graz_D dig_P = gamma1*graz dig_Z = gamma1*graz_Z dig_Det = gamma1*graz_Det dig_D = gamma1*graz_D*(redntp/diazntp) dig = dig_P + dig_Z + dig_Det + dig_D excr_P = gamma1*(1-geZ)*graz excr_Z = gamma1*(1-geZ)*graz_Z excr_Det = gamma1*(1-geZ)*graz_Det excr_D = gamma1*(1-geZ)*graz_D*(redntp/diazntp) excr = excr_P + excr_Z + excr_Det + excr_D nr_excr_D = gamma1*graz_D*(1-(redntp/diazntp)) & + (1-gamma1)*graz_D*(1-(redntp/diazntp)) sf_P = (1-gamma1)*graz sf_Z = (1-gamma1)*graz_Z sf_Det = (1-gamma1)*graz_Det sf_D = (1-gamma1)*graz_D*(redntp/diazntp) sf = sf_P + sf_Z + sf_Det + sf_D ! nutrients equation biopo4 = biopo4 + tsb*(redptn*((1.-dfrt)*morpt + excr + remi & - (npp-dopupt)) + diazptn*(morpt_D - (npp_D-dopupt_D)) & + recy_p) ! DOP equation biodop = biodop + tsb*(redptn*dfr*morp + redptn*dfrt*morpt & - redptn*dopupt - diazptn*dopupt_D - recy_p) ! phytoplankton equation biophyt = biophyt + tsb*(npp - morp - graz - morpt) ! zooplankton equation biozoop = biozoop + tsb*(dig - morz - graz_Z - excr) ! detritus equation biodetr = biodetr + tsb*((1.-dfr)*morp + sf + morz - remi & - graz_Det - expo + impo + morp_D*(redntp/diazntp)) ! DIC equation biodic = biodic + tsb*redctn*((1.-dfrt)*morpt + excr + remi & - npp + morpt_D - npp_D + recy + nr_excr_D & + morp_D*(1-(redntp/diazntp))) ! nitrate (NO3) equation biono3 = biono3 + tsb*((1.-dfrt)*morpt + excr + morpt_D & + morp_D*(1-(redntp/diazntp)) + nr_excr_D + remi + recy & - npp - no3upt_D) ! DON equation biodon = biodon + tsb*(dfr*morp + dfrt*morpt - recy) ! diazotroph equation biodiaz = biodiaz + tsb*(npp_D - morp_D - morpt_D - graz_D) 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 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 ! 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) 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) return end