! source file: /gxfs_work1/fs1/work-geomar6/smomw258/Kvale_etal_2020_GMD/model_output_for_paper/updates/npzd_src.F subroutine npzd_ini include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "calendar.h" include "coord.h" include "grdvar.h" include "index.h" include "iounit.h" include "levind.h" include "mw.h" include "scalar.h" include "stab.h" include "state.h" include "switch.h" include "tmngr.h" include "vmixc.h" include "npzd.h" real tmpij(imtm2,jmtm2) real tmpijkm(1:100,1:100,1:3,1:12) real tmpijk1(1:100,1:100,1:19) character (120) :: fname, new_file_name, logfile integer ie, is, je, js, i, ioun, j, jrow, n, m, indp, ip, k, kr integer jq, mask, kz, nv, ke, ks, iou, isle integer ib(10), ic(10) logical error, vmixset, hmixset, exists, inqvardef real snapint, snapls, snaple, snapde, trajint, cksum, checksum real ocnp, boxat, boxau, dvolt, dvolu, sum, pctocn, diffa, amix real runstep, dtatms, ahbkg, spnep, senep, c1e3, c1e9 c1e3 = 1000. c1e9 = 1.e9 ! convert units of NPZD parameters to MOM units redctn = redctn*1.e-3 redotn = redotn*1.e-3 redotp = redotn/redptn redctp = redctn/redptn redntp = 1./redptn 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 alpha_D = alpha_D/daylen kfeorg = kfeorg/daylen kfeorg_ca = kfeorg_ca/daylen kfecol = kfecol/daylen alphamax = alphamax/daylen alphamin = alphamin/daylen k1p_C = k1n_C*redptn alpha_C = alpha_C/daylen abioc = abioc/daylen nuc = nuc/daylen nuct0 = nuct0/daylen wc0 = wc0*1.e2 kc_c = kc_c*1.e-2 dissk0 = dissk0/daylen k1p_Diat = k1n_Diat*redptn alpha_Diat = alpha_Diat/daylen abiodiat = abiodiat/daylen nu_diat = nu_diat/daylen nudt0 = nudt0/daylen ws0 = (ws0*1.e2)/daylen abio = abio/daylen nup = nup/daylen nupt0 = nupt0/daylen gbio0 = gbio0/daylen nuz = nuz/daylen nud0 = nud0/daylen ! calculate sinking speed of detritus divided by grid width do k=1,km ! calcite sinking speed wc(k) = (wc0+wdc*zt(k))/daylen/dzt(k) ! [s-1] ! linear increase wd0-200m with depth wd(k) = (wd0+wdd*zt(k))/daylen/dzt(k) ! [s-1] ! Note that the above equation has been modified from the original ! one which was wd(k) = (wd0+4.0e-2*zt(k))/daylen/dzt(k). rkwz(k) = 1./(kw*dzt(k)) enddo ! Note that the following is a bug correction that fixes the ! level of light that the upper ocean recieves. The original ! code was ztt(k) = -zt(k)+ dzt(k)/2. and resulted in too ! much light at the surface because ztt(1) did not equal 0. ztt(1)=0.0 do k=1,km-1 ztt(k+1)=(-1)*zw(k) enddo ! read in iron hydrothermal input file as annual average mask ! Note the conversion from mol/m^3/s to mmol/m^3/s 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 call openfile (trim(fname), iou) call getvara ('O_fe_hydr', iou, ic(1)*ic(2)*ic(3) &, ib, ic, tmpijk1, c1e3, c0) fe_hydr(2:imtm1,2:jmtm1,:) = tmpijk1(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(i,jmtm1,k) enddo enddo else print*,"Warning => Cannot find" , trim(fname) endif ! read in Silicate hydrothermal input file as annual average mask ! read in as mol Si per m^3 per sec si_hydr(:,:,:) = 0. ib(:) = 1 ic(:) = imtm2 ic(2) = jmtm2 ic(3) = km fname = new_file_name ("O_si_hydr.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) call getvara ('O_SIL_HYDR', iou, ic(1)*ic(2)*ic(3) &, ib, ic, tmpijk1, c1, c0) si_hydr(2:imtm1,2:jmtm1,:) = tmpijk1(1:imtm2,1:jmtm2,:) where (si_hydr(:,:,:) .lt. -1.e30) si_hydr(:,:,:) = 0.0 endwhere do k=1,km do j=1,jmt si_hydr(1,j,k) = si_hydr(imtm1,j,k) si_hydr(imt,j,k) = si_hydr(2,j,k) enddo do i=1,imt si_hydr(i,1,k) = si_hydr(i,2,k) si_hydr(i,jmt,k) = si_hydr(i,jmtm1,k) enddo enddo else print*,"Warning => Cannot find" , trim(fname) endif !--------------------------------------------------------------------- ! calculate variables used in Opal remineralization !--------------------------------------------------------------------- ! Test a few options. All ropk(k) unitless. From ! HAMOCC5 (Aumont et al. 2003 from Maier-Reimer et al. 1993) ropk(1) = -(exp(-zw(1)/dopal)-1.0) do k=2,km ropk(k) = -(exp(-zw(k)/dopal)) + (exp(-zw(k-1)/dopal)) enddo !--------------------------------------------------------------------- ! 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 end subroutine npzd_src (bioin, ntsb, tsb, gl, bct, impo, dzt &, dayfrac, wwd, gbio &, rkw, nud, bioout, expoout &, grazout, morpout, morzout, graz_Det_out &, graz_Zout &, nppout, morptout, remiout, excrout &, impo_B, expoout_B, remiout_B, graz_Det_out_B &, npp_Dout, graz_Dout, morp_Dout, nfixout &, npp_Cout, morpt_Cout,graz_Cout &, morp_Cout &, thetamaxout_C, deffeout_C &, npp_Diatout, morpt_Diatout,graz_Diatout &, morp_Diatout &, thetamaxout_Diat, deffeout_Diat &, biosi &, impocaco3, wwc, expocaco3out, disslout &, calproout &, calattout,dissk1 &, avej_out, avej_D_out, gmax_out, no3P_out &, po4P_out, po4_D_out &, bctz &, 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. ! ! 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] ! bioin(7) = coccolithophores [mmol m-3] ! bioin(8) = calcite [mmol C 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 ! wwdz = sinking speed of poo/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 ! zprefMP = Z preference for MP ! zprefDiaz = 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 !======================================================================= 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, 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, gbio 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 real p1,p2,kfevar real expofeout, impofeout, feorgadsout, fecolout real remifeout, thetamaxout, deffeout, thetachlout real chlout, feprimeout, o2, chl_D_out real feorgads_ca, feorgadsout_ca real expo_B, impo_B, remi_B, remiout_B, expoout_B, dflag_B real graz_Det_out_B, graz_Det_B, g_Det_B, biod_B real graz_Z_B, graz_C_B, morp_C_B, morz_B real biodiat,jmax_Diat,u_Diat,avej_Diat,gd_Diat,npp_Diat real npp_Diatout,morp_Diat,morpt_Diat, biosi, gl_Diat real nudt,g_Diat,diatflag,graz_Diat,graz_Diatout,morp_Diatout real morpt_Diatout, ing_Diat, felimit_Diat real p1_Diat, p2_Diat, thetamaxout_Diat real thetamax_Diat, deffe_Diat real deffeout_Diat, thetachl_Diat real bioc,jmax_C,u_C,avej_C,gd_C,npp_C,npp_Cout,morp_C,morpt_C real nuct,g_C,bioc2,cflag,graz_C,graz_Cout,morp_Cout real morpt_Cout, ing_C, felimit_C, gl_C real p1_C, p2_C, thetamaxout_C, thetamax_C, deffe_C real deffeout_C, thetachl_C real caco3flag, wwc,expocaco3 real expocaco3out,calpro,calproout real biocaco3, impocaco3, dissk1, dissl, disslout real calatt,calattout 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(max(trcmin,bioin(ibiop)),pmax) p2=max(0.0,max(trcmin,bioin(ibiop))-pmax) kfevar = kfe deffe=max(trcmin,bioin(ibiodfe)) &/(kfevar+max(trcmin,bioin(ibiodfe))) thetamax=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe alpha=alphamin+(alphamax-alphamin)*deffe gl_O=gl*thetamax*alpha deffe_D=max(trcmin,bioin(ibiodfe))/(kfe_D & +max(trcmin,bioin(ibiodfe))) thetamax_D=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe_D alpha_D=alphamin+(alphamax-alphamin)*deffe_D gl_D=gl*thetamax_D*alpha_D deffe_C=max(trcmin,bioin(ibiodfe)) &/(kfe_C+max(trcmin,bioin(ibiodfe))) thetamax_C=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe_C alpha_C=alphamin+(alphamax-alphamin)*deffe_C gl_C=gl*thetamax_C*alpha_C deffe_Diat=max(trcmin,bioin(ibiodfe)) &/(kfe_Diat+max(trcmin,bioin(ibiodfe))) thetamax_Diat=thetamaxlo+(thetamaxhi-thetamaxlo)*deffe_Diat alpha_Diat=alphamin+(alphamax-alphamin)*deffe_Diat gl_Diat=gl*thetamax_Diat*alpha_Diat ! photosynthesis after Evans & Parslow (1985) ! notation as in JGOFS report No. 23 p. 6 f1 = exp(((-kw & ) - kc*(max(trcmin,bioin(ibiop)) & +max(trcmin,bioin(ibiodiaz)) & +max(trcmin,bioin(ibioc)) & +max(trcmin,bioin(ibiodiat)) & ) & - kc_c*max(trcmin,bioin(ibiocaco3)) & )*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*(max(trcmin,bioin(ibiop)) & + max(trcmin,bioin(ibiodiaz)) & + max(trcmin,bioin(ibioc)) & + max(trcmin,bioin(ibiodiat)) & ) & + kc_c*max(trcmin,bioin(ibiocaco3)) & )*dzt) ! Make the max grazing rate a function of temperature gmax = gbio*bctz ! Note that bctz is set in tracers.F jmax_C = abioc*bct*deffe_C gd_C = jmax_C*dayfrac u1 = max(gl_C/gd_C,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_C = gd_C*(phi1 - phi2)/(( & kw + kc*(max(trcmin,bioin(ibiop)) & + max(trcmin,bioin(ibiodiaz)) & + max(trcmin,bioin(ibioc)) & + max(trcmin,bioin(ibiodiat)) & ) & + kc_c*max(trcmin,bioin(ibiocaco3)) & )*dzt) jmax_Diat = abiodiat*bct*deffe_Diat gd_Diat = jmax_Diat*dayfrac u1 = max(gl_Diat/gd_Diat,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_Diat = gd_Diat*(phi1 - phi2)/(( & kw + kc*(max(trcmin,bioin(ibiop)) & + max(trcmin,bioin(ibiodiaz)) & + max(trcmin,bioin(ibioc)) & + max(trcmin,bioin(ibiodiat)) & ) & + kc_c*max(trcmin,bioin(ibiocaco3)) & )*dzt) 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,bioin(ibiop)) & + max(trcmin,bioin(ibiodiaz)) & + max(trcmin,bioin(ibioc)) & + max(trcmin,bioin(ibiodiat)) & ) & + kc_c*max(trcmin,bioin(ibiocaco3)) & )*dzt) ! check grazing preferences = 1 for N case IF ((zprefP + zprefDet + zprefZ + zprefDiaz & + zprefDiat + zprefMP + zprefC).ne.1) THEN print*, ' ' print*, '==> Error: grazing preference values must be' print*, ' set in control.in' print*, 'zprefP=',zprefP print*, 'zprefZ=',zprefZ print*, 'zprefC=',zprefC print*, 'zprefDiat=',zprefDiat print*, 'zprefDiaz=',zprefDiaz print*, 'zprefDet=',zprefDet print*, 'zprefMP=',zprefMP stop '=>npzd_src' END IF nupt = nupt0*bct bioout(:) = 0.0 bion = bioin(ibion) biop = bioin(ibiop) bioz = bioin(ibioz) biod = bioin(ibiod) bioc = bioin(ibioc) biodiat = bioin(ibiodiat) biono3 = bioin(ibiono3) biodiaz = bioin(ibiodiaz) biodfe = bioin(ibiodfe) biodetrfe = bioin(ibiodetrfe) nuct = nuct0*bct npp_Cout = 0.0 graz_Cout = 0.0 morp_Cout = 0.0 morpt_Cout = 0.0 nudt = nudt0*bct npp_Diatout = 0.0 graz_Diatout = 0.0 morp_Diatout = 0.0 morpt_Diatout = 0.0 biocaco3 = bioin(ibiocaco3) calproout = 0.0 calattout = 0.0 disslout = 0.0 expocaco3out = 0.0 biod_B = bioin(ibiod_B) expoout_B = 0.0 graz_Det_out_B = 0.0 remiout_B = 0.0 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 npp_Dout = 0.0 graz_Dout = 0.0 morp_Dout = 0.0 nfixout = 0.0 expofeout = 0.0 remifeout = 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 do n=1,ntsb p1 = min(max(biop,trcmin),pmax) p2 = max(0.0,max(biop,trcmin) - pmax) kfevar = kfe deffe = max(trcmin,biodfe)/(kfevar + max(trcmin,biodfe)) jmax = abio*bct*deffe deffe_C = max(trcmin,biodfe)/(kfe_C + max(trcmin,biodfe)) jmax_C = abioc*bct*deffe_C deffe_Diat = max(trcmin,biodfe) & /(kfe_Diat + max(trcmin,biodfe)) jmax_Diat = abiodiat*bct*deffe_Diat 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 u_P = min(avej, jmax*max(trcmin,bion)/(k1p + max(trcmin,bion))) po4P = jmax*bion/(k1p + bion) u_C = min(avej_C, jmax_C*max(trcmin,bion)/(k1p_C & + max(trcmin,bion))) ! Si limitation here is a mask unless Si tracer also used u_Diat = min(avej_Diat, jmax_Diat*max(trcmin,biosi) & /(k1si + max(trcmin,biosi))) u_Diat = min(u_Diat, jmax_Diat*max(trcmin,bion) & /(k1p_Diat + max(trcmin,bion))) ! nitrate limitation u_P = min(u_P, jmax*max(trcmin,biono3)/(k1n & + max(trcmin,biono3))) u_C = min(u_C, jmax_C*max(trcmin,biono3)/(k1n_C & + max(trcmin,biono3))) u_Diat = min(u_Diat, jmax_Diat*max(trcmin,biono3) & /(k1n_Diat + max(trcmin,biono3))) no3P = jmax*max(trcmin,biono3)/(k1n + max(trcmin,biono3)) ! growth rate of diazotrophs smaller than other phytoplankton and ! not nitrate limited u_D = min(avej_D, jmax_D*max(trcmin,bion)/(k1p & + max(trcmin,bion))) po4_D = jmax_D*max(trcmin,bion)/(k1p + max(trcmin,bion)) ! Set the grazing coefficients thetaZ = zprefP*max(trcmin,biop)+zprefDet*max(trcmin,biod & ) & +zprefZ*max(trcmin,bioz) + kzoo & +zprefDiaz*max(trcmin,biodiaz) & + zprefC*max(trcmin,bioc) & + zprefDiat*max(trcmin,biodiat) & + zprefDet*max(trcmin,biod_B) ing_P = zprefP/thetaZ ing_Det = zprefDet/thetaZ ing_Z = zprefZ/thetaZ ing_D = zprefDiaz/thetaZ ing_C = zprefC/thetaZ ing_Diat = zprefDiat/thetaZ npp = u_P*biop npp_C = u_C*bioc npp_Diat = u_Diat*biodiat npp_D = u_D*biodiaz ! grazing on diazotrophs g_D = gmax*ing_D*biodiaz graz_D = g_D*bioz morp_D = nupt*biodiaz ! linear mortality 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 g_C = gmax*ing_C*bioc graz_C = g_C*bioz morp_C = nuc*bioc morpt_C = nuct*bioc g_Diat = gmax*ing_Diat*biodiat graz_Diat = g_Diat*bioz morp_Diat = nu_diat*biodiat morpt_Diat = nudt*biodiat expocaco3 = wwc*biocaco3 !mmol/m3/s dissl = biocaco3*dissk1 !mmol/m3 g_Det_B = gmax*ing_Det*biod_B graz_Det_B = g_Det_B*bioz remi_B = bapr*dissl/(capr*redctn*1.e3) !mmol N/m3 if(remi_B.lt.0) remi_B = 0. if(remi_B.ge.biod_B) remi_B = biod_B expo_B = wwc*biod_B !ballasted by CaCO3 morp = nup*biop morpt = nupt*biop 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)**0.58)*feprime feorgads_ca = kfeorg_ca*(((max(trcmin,biocaco3*1.e-3) & )*60.01)**0.58)*feprime !molar mass caco3 fecol = kfecol*feprime else ! These are default values. In case of O2 < O2min, no scavenging is taking place feorgads = 0.0 feorgads_ca = 0.0 fecol = 0.0 feprime = 0.0 fepa = 0.0 endif expofe = wwd * biodetrfe ! 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) 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) dflag_B = 0.5 + sign(0.5,biod_B - trcmin) cflag = 0.5 + sign(0.5,bioc - trcmin) diatflag = 0.5 + sign(0.5,biodiat - trcmin) caco3flag = 0.5 + sign(0.5,biocaco3 - 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 npp = npp*nflag*no3flag*pflag & *dfeflag npp_D = npp_D*nflag*diazflag graz_D = graz_D*diazflag*zflag morp_D = morp_D*diazflag no3upt_D = no3upt_D*no3flag*diazflag remife=remife*detrfeflag feorgads=feorgads*dfeflag feorgads_ca=feorgads_ca*dfeflag expofe=expofe*detrfeflag fecol=fecol*dfeflag graz_Det_B = graz_Det_B*dflag_B*zflag dissl = dissl*caco3flag expocaco3 = expocaco3*caco3flag calatt = calatt*caco3flag npp_C = npp_C*nflag*cflag & *no3flag & *dfeflag graz_C = graz_C*cflag*zflag morp_C = morp_C*cflag morpt_C = morpt_C*cflag npp_Diat = npp_Diat*nflag*diatflag & *no3flag & *dfeflag graz_Diat = graz_Diat*diatflag*zflag morp_Diat = morp_Diat*diatflag morpt_Diat = morpt_Diat*diatflag avej = avej*pflag avej_D = avej_D*diazflag gmax = gmax*zflag no3P = no3P*no3flag po4P = po4P*nflag po4_D = po4_D*nflag ! Excretion is the difference between ! the zooplankton assimilation and growth efficiencies excr = (gamma1-geZ)*(graz+graz_Z+graz_Det & + graz_D & + graz_C & + graz_Diat & + graz_Det_B & ) ! nutrients equation bion = bion + tsb*redptn*(remi + excr - npp + morpt & - npp_D & - npp_C + morpt_C & - npp_Diat + morpt_Diat & ) ! phytoplankton equation biop = biop + tsb*(npp - morp - graz - morpt) ! coccolithophores equation bioc = bioc + tsb*(npp_C - morp_C - graz_C & - morpt_C) ! diatoms equation biodiat = biodiat + tsb*(npp_Diat - morp_Diat - graz_Diat & - morpt_Diat) ! zooplankton equation bioz = bioz + tsb*(geZ*(graz + graz_Det + graz_Z & + graz_D & + graz_C & + graz_Diat & + graz_Det_B & ) & - morz - graz_Z) ! detritus equations biod = biod + tsb*((1.-gamma1)* & (graz + graz_Det & + (graz_Z & - bapr*graz_Z & ) & + graz_D & + graz_Diat & + (graz_C & - bapr*graz_C & ) & ) & + morp_D & + morp_Diat & + (morp_C & - bapr*morp_C & ) & + (morz & - bapr*morz & ) & + morp - remi - graz_Det - expo + impo) & + remi_B biod_B = biod_B + tsb*((1.-gamma1)*(graz_Det_B & + bapr*graz_Z & + bapr*graz_C & ) & + bapr*morp_C & + bapr*morz & - graz_Det_B - expo_B + impo_B) & - remi_B ! calcite equations ! net formation of attached (living) tests ! total primary production by coccs calatt = (capr*npp_C & - (morpt_C + morp_C + graz_C)*capr & + capr*geZ*( & graz_Det_B & + graz_Diat & + graz + graz_D + graz_Det & + graz_Z + graz_C) - (graz_Z + morz)*capr & )*redctn*1.e3 !convert to carbon units mmol C ! formation of detached (dead) tests, or PIC calpro = ((1.-gamma1)*(capr*graz_C + capr*graz_Z) & + capr*morp_C + capr*morz)*redctn*1.e3 !stay in mmol biocaco3 = biocaco3 -dissl + tsb*(calpro-expocaco3+impocaco3) ! nitrate (NO3) equation biono3 = biono3 + tsb*(remi + excr - npp + morpt - no3upt_D & - npp_C + morpt_C & - npp_Diat + morpt_Diat & ) ! diazotroph equation biodiaz = biodiaz + tsb*(npp_D - morp_D - graz_D) ! dissolved iron equation biodfe = biodfe + tsb*(rfeton*(excr - (npp & + npp_D & + npp_C & + npp_Diat & ) + morpt & + morpt_C & + morpt_Diat & )- feorgads & - feorgads_ca & + remife - fecol) ! particulate iron equation ! KK ignoring ballasting of iron for now biodetrfe = biodetrfe + tsb*(rfeton*((1.-gamma1)*(graz+graz_Det & + graz_D & + graz_C & + graz_Det_B & + graz_Diat & + graz_Z) + morp & + morp_D & + morp_C & + morp_Diat & + morz - graz_Det & - graz_Det_B & )+ feorgads & + feorgads_ca & + fecol - remife - expofe + impofe) 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 expoout_B = expoout_B + expo_B remiout_B = remiout_B + remi_B graz_Det_out_B = graz_Det_out_B + graz_Det_B nppout = nppout + npp morptout = morptout + morpt remiout = remiout + remi excrout = excrout + excr npp_Dout = npp_Dout + npp_D graz_Dout = graz_Dout + graz_D morp_Dout = morp_Dout + morp_D nfixout = nfixout + npp_D - no3upt_D expofeout = expofeout + expofe remifeout = remifeout + remife calproout = calproout + calpro calattout = calattout + calatt disslout = disslout + dissl expocaco3out = expocaco3out + expocaco3 npp_Cout = npp_Cout + npp_C graz_Cout = graz_Cout + graz_C morp_Cout = morp_Cout + morp_C morpt_Cout = morpt_Cout + morpt_C npp_Diatout = npp_Diatout + npp_Diat graz_Diatout = graz_Diatout + graz_Diat morp_Diatout = morp_Diatout + morp_Diat morpt_Diatout = morpt_Diatout + morpt_Diat 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 enddo ! This becomes the change in the tracer, not the tracer bioout(ibion) = bion - bioin(ibion) bioout(ibiop) = biop - bioin(ibiop) bioout(ibioz) = bioz - bioin(ibioz) bioout(ibiod) = biod - bioin(ibiod) bioout(ibiono3) = biono3 - bioin(ibiono3) bioout(ibiodiaz) = biodiaz - bioin(ibiodiaz) bioout(ibiodfe) = biodfe - bioin(ibiodfe) bioout(ibiodetrfe) = biodetrfe - bioin(ibiodetrfe) bioout(ibioc) = bioc - bioin(ibioc) bioout(ibiodiat) = biodiat - bioin(ibiodiat) bioout(ibiocaco3) = biocaco3 - bioin(ibiocaco3) bioout(ibiod_B) = biod_B - bioin(ibiod_B) return end