subroutine npzd_src (bioin, ntsb, tsb, gl, temp, impo, dzt &, dayfrac, wwd, rkw, nud, bioout, expoout &, grazout, morpout, morzout #if defined osu_felim &, i, j, mi #endif # if defined uvic_npzd_out &, nppout, morptout, remiout, excrout # if defined uvic_nitrogen &, npp_Dout, graz_Dout, morp_Dout, nfixout # endif # endif # if defined osu_c13 &, rc13impo, rc13expoout, co2aq # endif & ) !======================================================================= ! computes source terms of the NPZD model ! initial version of code adapted from Xavier Giraud: ! Giraud et al. 2000, J Mar Res, 58, 609-630 ! original model reference: ! Oeschlies and Garcon 1999, Global Biogeochem. Cycles 13, 135-160 ! Schmittner et al. 2005, Global Biogeochem. Cycles 19, GB3004, ! doi:10.1029/2004GB002283. ! This is the optimized model version of Schartau & Oschlies, 2003, ! J. Mar. Res. 62, 765-793 and Oschlies & Schartau, in prep. ! as described in Schmittner and Oschlies, in prep. ! note that nutrient now represents phosphate ! Andreas Schmittner Sept 2005 (aschmitt@coas.oregonstate.edu) ! input variables: ! bioin(1:4) = N,P,Z,D [mmol m-3] ! bioin(5) = nitrate [mmol m-3] ! bioin(6) = diazotrophs [mmol m-3] ! gl = 2.*light at top of grid box ! ntsb = number of time steps ! tsb = time step [s] ! temp = temperature [deg C] ! impo = import of detritus from above [mmol m-3] ! rc13impo = import of 13C/12C of detritus from above ! dzt = depth of grid box [cm] ! dayfrac = day length (fraction: 0 < dayfrac < 1) ! wwd = sinking speed of detritus/dzt ! rkw = reciprical of kw*dzt(k) ! nud = remineralisation rate of detritus [s-1] ! co2aq = aquatic CO2 concentration (mol/m^3) ! output variables: ! bioout = change from bioin [mmol m-3] ! nppout = net primary production [mmol m-3] ! grazout = grazing [mmol m-3] ! morpout = quadratic mortality of phytoplankton [mmol m-3] ! morptout = specific mortality of phytoplankton [mmol m-3] ! morzout = mortality of zooplankton [mmol m-3] ! remiout = remineralisation [mmol m-3] ! excrout = excretion [mmol m-3] ! expoout = detrital export [mmol m-3] ! rc13expoout= 13C/12C ratio of detrital export ! npp_Dout = NPP of diazotrophs ! graz_Dout = grazing of diazotrophs ! morp_Dout = mortality of diazotrophs ! nfixout = rate of N2 fixation !======================================================================= #if defined uvic_npzd implicit none # include "param.h" # include "calendar.h" # include "npzd.h" integer n, ntsb, i, j, mi real bioin(ntnpzd+1), bioout(ntnpzd) real gl, f1, bion, biop, bioz, biod, jmax, u_P, g_P, npp, graz real morp, morpt, morz, remi, excr, expo, impo, nppout, grazout real morpout, morptout, morzout, remiout, excrout, expoout, tsb real dzt, nflag, pflag, zflag, dflag, wwd, rkw, gd, dayfrac, bct real nupt, nud, biono3, u_D,npp_D, npp_Dout, no3flag, biodiaz real diazflag, g_D,graz_D, morp_D, jmax_D, gd_D, avej_D, no3upt_D real morp_Dout, graz_Dout, nfixout, biop2, u1, u2, phi1, phi2 real avej real biodic13, biophytc13, biozoopc13, biodetrc13, biodiazc13 real dic13flag, phytc13flag, zoopc13flag, detrc13flag, diazc13flag real rdic13, rc13phyt, rc13zoop, rc13detr, rc13diaz real rc13impo, rc13expoout, biodic, dicflag, ac13bp, ac13bd real ac13_aq_POC_P, ac13_aq_POC_D, ac13_DIC_aq, temp, co2aq real ktk, crau, kkrau, epsf, epsrau, mui ! photosynthesis after Evans & Parslow (1985) ! notation as in JGOFS report No. 23 p. 6 f1 = exp((-kw - kc*max(bioin(2),trcmin))*dzt) bct = bbio**(cbio*temp) jmax = abio*bct*(0.3+0.7*felim(i,j,mi)) gd = jmax*dayfrac u1 = max(gl/gd, 1.e-6) u2 = u1*f1 c phi1 = u1*(0.555588 + 0.004926*u1)/(1. + 0.188721*u1) c phi2 = u2*(0.555588 + 0.004926*u2)/(1. + 0.188721*u2) phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1 phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2 avej = gd*rkw*(phi1 - phi2) # if defined uvic_nitrogen c jmax_D = max(0.,abio_D*(bct-2.6))*(0.1+0.9*felim(i,j,mi)) jmax_D = max(0.,abio_D*bct)*(0.1+0.9*felim(i,j,mi)) gd_D = max(1.e-14,jmax_D*dayfrac) u1 = max(gl/gd_D, 1.e-6) u2 = u1*f1 c phi1 = u1*(0.555588 + 0.004926*u1)/(1. + 0.188721*u1) c phi2 = u2*(0.555588 + 0.004926*u2)/(1. + 0.188721*u2) phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1 phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2 avej_D = gd_D*rkw*(phi1 - phi2) # endif nupt = nupt0*bct # if defined osu_c13 kkrau = 0.3583 ktk = 5.019e-6*exp(-19510./(8.3143*(temp+273.15))) crau = 1.e4 + 1.82e-5/ktk epsf = 23. # endif bioout(:) = 0.0 bion = bioin(1) biop = bioin(2) bioz = bioin(3) biod = bioin(4) # if defined uvic_nitrogen biono3 = bioin(5) biodiaz = bioin(6) # endif # if defined osu_c13 biodic13 = bioin(7) biophytc13 = bioin(8) biozoopc13 = bioin(9) biodetrc13 = bioin(10) biodiazc13 = bioin(11) biodic = bioin(12) # endif 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) # if defined uvic_nitrogen no3flag = 0.5 + sign(0.5,biono3 - trcmin) diazflag = 0.5 + sign(0.5,biodiaz - trcmin) # endif # if defined osu_c13 dic13flag = 0.5 + sign(0.5,biodic13 - trcmin) phytc13flag = 0.5 + sign(0.5,biophytc13 - trcmin) zoopc13flag = 0.5 + sign(0.5,biozoopc13 - trcmin) detrc13flag = 0.5 + sign(0.5,biodetrc13 - trcmin) diazc13flag = 0.5 + sign(0.5,biodiazc13 - trcmin) dicflag = 0.5 + sign(0.5,biodic - trcmin) # endif bion = max(bion, trcmin) biop = max(biop, trcmin) bioz = max(bioz, trcmin) biod = max(biod, trcmin) # if defined uvic_nitrogen biono3 = max(biono3, trcmin) biodiaz = max(biodiaz, trcmin) # endif # if defined osu_c13 biodic13 = max(biodic13, trcmin*rc13std) biophytc13 = max(biophytc13, trcmin*rc13std) biozoopc13 = max(biozoopc13, trcmin*rc13std) biodetrc13 = max(biodetrc13, trcmin*rc13std) biodiazc13 = max(biodiazc13, trcmin*rc13std) biodic = max(biodic, trcmin) rc13expoout = 0.0 ac13_DIC_aq = -1.0512994e-4*temp+1.0117605 c ac13_DIC_aq = 1. # endif expoout = 0.0 grazout = 0.0 morpout = 0.0 morzout = 0.0 # if defined uvic_npzd_out nppout = 0.0 morptout = 0.0 remiout = 0.0 excrout = 0.0 # if defined uvic_nitrogen npp_Dout = 0.0 graz_Dout = 0.0 morp_Dout = 0.0 nfixout = 0.0 # endif # endif do n=1,ntsb ! growth rate of phytoplankton u_P = min(avej, jmax*bion/(k1p + bion)) # if defined uvic_nitrogen ! nitrate limitation u_P = min(u_P, jmax*biono3/(k1 + biono3)) ! growth rate of diazotrophs smaller than other phytoplankton and ! not nitrate limited u_D = min(avej_D, jmax_D*bion/(k1p + bion)) # endif npp = u_P*biop # if defined uvic_nitrogen npp_D = max(0.,u_D*biodiaz) ! grazing function for diazotrophs g_D = gbio*epsbio*biodiaz*biodiaz/(gbio+epsbio*biodiaz*biodiaz) graz_D = g_D*bioz*bct morp_D = nupt*biodiaz ! linear mortality no3upt_D = biono3/(k1 + biono3)*npp_D ! nitrate uptake # endif ! grazing function biop2 = biop*biop g_P = gbio*epsbio*biop2/(gbio+epsbio*biop2) graz = g_P*bioz*bct morp = nup*biop2 morpt = nupt*biop morz = nuz*bioz*bioz remi = nud*bct*biod excr = gamma2*bct*bioz expo = wwd*biod ! flags prevent negative values by setting outgoing fluxes to ! zero if tracers are lower than trcmin c nflag = 0.5 + sign(0.5,bion - trcmin) c pflag = 0.5 + sign(0.5,biop - trcmin) c zflag = 0.5 + sign(0.5,bioz - trcmin) c dflag = 0.5 + sign(0.5,biod - trcmin) c# if defined uvic_nitrogen c no3flag = 0.5 + sign(0.5,biono3 - trcmin) c diazflag = 0.5 + sign(0.5,biodiaz - trcmin) c# endif graz = graz*pflag*phytc13flag morp = morp*pflag*phytc13flag morpt = morpt*pflag*phytc13flag morz = morz*zflag*zoopc13flag remi = remi*dflag*detrc13flag excr = excr*zflag*zoopc13flag expo = expo*dflag*detrc13flag # if defined uvic_nitrogen npp = npp*nflag*no3flag*dic13flag npp_D = npp_D*nflag graz_D = graz_D*diazflag*diazc13flag morp_D = morp_D*diazflag*diazc13flag no3upt_D = no3upt_D*no3flag*dic13flag # else npp = npp*nflag # endif # if defined osu_c13 rdic13 = biodic13/biodic rdic13 = min(rdic13, 2.) rdic13 = max(rdic13, 0.5) c rdic13 = rc13std rc13phyt = biophytc13/(biop*redctn) rc13phyt = min(rc13phyt, 2.) rc13phyt = max(rc13phyt, .5) c rc13phyt = rc13std rc13zoop = biozoopc13/(bioz*redctn) rc13zoop = min(rc13zoop, 2.) rc13zoop = max(rc13zoop, 0.5) c rc13zoop = rc13std rc13detr = biodetrc13/(biod*redctn) rc13detr = min(rc13detr, 2.) rc13detr = max(rc13detr, 0.5) c rc13detr = rc13std rc13diaz = biodiazc13/(biodiaz*redctn) rc13diaz = min(rc13diaz, 2.) rc13diaz = max(rc13diaz, 0.5) c rc13diaz = rc13std # endif # if defined uvic_nitrogen ! !!!!!!!!!!!!!!!!! nutrients equation bion = bion + tsb*redptn*(remi + excr - (npp + npp_D) + morpt) ! !!!!!!!!!!!!!!!!! phytoplankton equation biop = biop + tsb*(npp - morp - graz - morpt) ! !!!!!!!!!!!!!!!!! zooplankton equation bioz = bioz + tsb*(gamma1*(graz + graz_D) - excr - morz) ! !!!!!!!!!!!!!!!!! detritus equation biod = biod + tsb*((1.-gamma1)*(graz + graz_D) + morp + morp_D & + morz - remi - expo + impo) ! !!!!!!!!!!!!!!!!! nitrate (NO3) equation biono3 = biono3 + tsb*(remi + excr - npp + morpt - no3upt_D) ! !!!!!!!!!!!!!!!!! diazotroph equation biodiaz = biodiaz + tsb*(npp_D - morp_D - graz_D) # else ! !!!!!!!!!!!!!!!!! nutrients equation bion = bion + tsb*redptn*(remi + excr - npp + morpt) ! !!!!!!!!!!!!!!!!! phytoplankton equation biop = biop + tsb*(npp - morp - graz - morpt) ! !!!!!!!!!!!!!!!!! zooplankton equation bioz = bioz + tsb*(gamma1*graz - excr - morz) ! !!!!!!!!!!!!!!!!! detritus equation biod = biod + tsb*((1.-gamma1)*graz + morp & + morz - remi - expo + impo) # endif # if defined osu_c13 biodic = biodic + tsb*redctn*(remi + excr - (npp+npp_D) + morpt) ! !!!!!!!!!!!!!!!!!! isotope equations c ac13b = 1. ac13_aq_POC_P = -0.017*log10(co2aq*1000.)+1.0034 c u_P is the daily averaged (24 h) growth rate c u_P*2.3 is the growth rate during the photoperiod (day time) c Laws et al. 1995 Geochim. et Cosm. Acta 59 No. 6, pp 1131 Fig. 2 c ac13_aq_POC_P = 1.-(24.7-66.7*u_P*86400*2.35/(co2aq*1000.))/1000. c Rau et al. 1996 with parameters from Hofmann et al. 2000 c ! mui = u_P*86400./dayfrac ! epsrau = kkrau*crau*mui/co2aq - epsf ! epsrau = max(min(epsrau,-15.),-40.) ! ac13_aq_POC_P = 1.+epsrau/1000. ac13bp=ac13_aq_POC_P/ac13_DIC_aq c ac13bp=0.975 c for diazotroph's use eps=15. (Carpenter et al. 1997 DSR I) ! ac13_aq_POC_D = 1.-15./1000. ! ac13bd=ac13_aq_POC_D/ac13_DIC_aq ac13bd=ac13_aq_POC_P/ac13_DIC_aq c ac13bd=0.975 biodic13 = biodic13 + tsb*redctn*(rc13detr*remi & + rc13zoop*excr - rdic13*(ac13bp*npp + ac13bd*npp_D) & + rc13phyt*morpt) biophytc13 = biophytc13 + tsb*redctn*(ac13bp*rdic13*npp & - rc13phyt*(morp + graz + morpt) ) biozoopc13 = biozoopc13 + tsb*redctn* & (gamma1*(rc13phyt*graz + rc13diaz*graz_D) & - rc13zoop*(excr + morz) ) biodetrc13 = biodetrc13 + tsb*redctn*((1.-gamma1)*(rc13phyt*graz & + rc13diaz*graz_D) + rc13phyt*morp + rc13diaz*morp_D & + rc13zoop*morz - rc13detr*(remi + expo) & + rc13impo*impo) biodiazc13 = biodiazc13 + tsb*redctn*(ac13bd*rdic13*npp_D & - rc13diaz*(morp_D + graz_D) ) rc13expoout = rc13expoout + rc13detr # endif expoout = expoout + expo grazout = grazout + graz morpout = morpout + morp morzout = morzout + morz # if defined uvic_npzd_out nppout = nppout + npp morptout = morptout + morpt remiout = remiout + remi excrout = excrout + excr # if defined uvic_nitrogen npp_Dout = npp_Dout + npp_D graz_Dout = graz_Dout + graz_D morp_Dout = morp_Dout + morp_D nfixout = nfixout + npp_D - no3upt_D # endif # endif if (nflag == 1) nflag = 0.5 + sign(0.5,bion - trcmin) if (pflag == 1) pflag = 0.5 + sign(0.5,biop - trcmin) if (zflag == 1) zflag = 0.5 + sign(0.5,bioz - trcmin) if (dflag == 1) dflag = 0.5 + sign(0.5,biod - trcmin) # if defined uvic_nitrogen if (no3flag == 1) no3flag = 0.5 + sign(0.5,biono3 - trcmin) if (diazflag == 1) diazflag = 0.5 + sign(0.5,biodiaz - trcmin) # endif # if defined osu_c13 dic13flag = 0.5 + sign(0.5,biodic13 - trcmin) phytc13flag = 0.5 + sign(0.5,biophytc13 - trcmin) zoopc13flag = 0.5 + sign(0.5,biozoopc13 - trcmin) detrc13flag = 0.5 + sign(0.5,biodetrc13 - trcmin) diazc13flag = 0.5 + sign(0.5,biodiazc13 - trcmin) dicflag = 0.5 + sign(0.5,biodic - trcmin) # endif enddo bioout(1) = bion - bioin(1) bioout(2) = biop - bioin(2) bioout(3) = bioz - bioin(3) bioout(4) = biod - bioin(4) # if defined uvic_nitrogen bioout(5) = biono3 - bioin(5) bioout(6) = biodiaz - bioin(6) # endif # if osu_c13 bioout(7) = biodic13 - bioin(7) bioout(8) = biophytc13 - bioin(8) bioout(9) = biozoopc13 - bioin(9) bioout(10) = biodetrc13 - bioin(10) bioout(11) = biodiazc13 - bioin(11) rc13expoout = rc13expoout/ntsb # endif #endif return end