! source file: /nfs/tape_cache/smomw258/latin_hypercube_nobio/template/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) 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 wmp0 = wmp0*1.e2 wdz0 = wdz0*1.e2 alpha = alpha/daylen alpha_D = alpha_D/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 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] !uniform plastic rising rates (for now) wmp(k) = wmp0/dzt(k) ![s-1] wdz(k) = (wdz0+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 do k=km,1,-1 if (zw(k) .gt. 24000) kmfe = k-1 enddo ! read in the dissolved Fe conc from O_dissolved_fe.nc fe_dissolved(:,:,:,:) = 0. ib(:) = 1 ic(:) = imtm2 ic(2) = jmtm2 ic(3) = kmfe ic(4) = 12 fname = new_file_name ("O_fe_dissolved.nc") inquire (file=trim(fname), exist=exists) if (exists) then c1e9 = 1000000000 call openfile (trim(fname), iou) call getvara ('O_dissolved_fe', iou, ic(1)*ic(2)*ic(3)*ic(4) &, ib, ic, tmpijkm, c1e9, c0) fe_dissolved(2:imtm1,2:jmtm1,1:kmfe,:) = tmpijkm(1:imtm2 &, 1:jmtm2,1:kmfe,:) do m=1,12 do k=1,kmfe do j=1,jmt fe_dissolved(1,j,k,m) = fe_dissolved(imtm1,j,k,m) fe_dissolved(imt,j,k,m) = fe_dissolved(2,j,k,m) enddo do i=1,imt fe_dissolved(i,1,k,m) = fe_dissolved(i,2,k,m) fe_dissolved(i,jmt,k,m) = fe_dissolved(i,jmtm1,k,m) enddo enddo enddo else print*,"Warning => Cannot find", trim(fname) endif !--------------------------------------------------------------------- ! 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 &, agg_con, poo_con,expompout_a,impomp_a &, expompout_p,impomp_p,agg_uptout,zoop_uptout &, agg_relout, zoop_relout &, wwdz, remizout, expozout, impoz &, 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 &, felimit_C &, impocaco3, wwc, expocaco3out, disslout &, calproout &, calattout,dissk1 &, avej_out, avej_D_out, gmax_out, no3P_out &, po4P_out, po4_D_out &, 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:4) = N,P,Z,D [mmol m-3] ! bioin(5) = nitrate [mmol m-3] ! bioin(6) = diazotrophs [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 ! felimit = Fe limitation parameter ! felmit_D = Fe limitation parameter for diazotrophs !======================================================================= 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 agg_upt, agg_rel,agg_uptout,zoop_uptout, agg_relout real agg_con, poo_con, biomp_a, mar_snow, zoop_relout real biomp_p, mpaflag, mppflag, ing_MP, g_MP, graz_MP real biomp, mpflag, impomp_a,expomp_a,expompout_a real impomp_p, expomp_p, expompout_p, zoop_upt, zoop_rel real wwdz, biodz, expoz, expozout real impoz, remiz, remizout, dzflag real felimit, felimit_D 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 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 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) ! 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)) & ) & - kc_c*max(trcmin,bioin(ibiocaco3)) & )*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*alpha/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)) & ) & + 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*felimit_C gd_C = jmax_C*dayfrac u1 = max(gl*alpha_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)) & ) & + kc_c*max(trcmin,bioin(ibiocaco3)) & )*dzt) jmax_D = max(0.,abio*(bct - 2.6)*felimit_D)*jdiar gd_D = max(1.e-14,jmax_D*dayfrac) u1 = max(gl*alpha/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)) & ) & + 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 zprefP = (1.-(zprefMP+zprefDiaz))*0.25 zprefDet = (1.-(zprefMP+zprefDiaz))*0.25 zprefC = (1.-(zprefMP+zprefDiaz))*0.25 zprefZ = 1.-(zprefMP+zprefDiaz+zprefC+zprefDet & +zprefP+zprefDiat) 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 END IF nupt = nupt0*bct bioout(:) = 0.0 bion = bioin(ibion) biop = bioin(ibiop) bioz = bioin(ibioz) biod = bioin(ibiod) bioc = bioin(ibioc) biono3 = bioin(ibiono3) biodiaz = bioin(ibiodiaz) nuct = nuct0*bct npp_Cout = 0.0 graz_Cout = 0.0 morp_Cout = 0.0 morpt_Cout = 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 biomp = bioin(ibiomp) biomp_a = bioin(ibiompa) biomp_p = bioin(ibiompp) expompout_a = 0.0 expompout_p = 0.0 agg_uptout = 0.0 zoop_uptout = 0.0 agg_relout = 0.0 zoop_relout = 0.0 mar_snow = 0.0 biodz = bioin(ibiodz) remizout = 0.0 expozout = 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 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 ! 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))) ! 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))) 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 & +biodz & ) & +zprefZ*max(trcmin,bioz) + kzoo & +zprefDiaz*max(trcmin,biodiaz) & + zprefC*max(trcmin,bioc) & + zprefDet*max(trcmin,biod_B) & + zprefMP*max(trcmin,biomp)*poo_con !KK test ing_P = zprefP/thetaZ ing_Det = zprefDet/thetaZ ing_Z = zprefZ/thetaZ ing_D = zprefDiaz/thetaZ ing_C = zprefC/thetaZ ! in N units ing_MP = zprefMP/thetaZ npp = u_P*biop npp_C = u_C*bioc 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 & +biodz & ) graz_Det = g_Det*bioz g_C = gmax*ing_C*bioc graz_C = g_C*bioz morp_C = nuc*bioc morpt_C = nuct*bioc 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 g_MP = 0.!gmax*ing_MP*biomp*poo_con graz_MP = 0.!g_MP*bioz zoop_upt = 0.!(gmax*ing_MP*bioz)/poo_con*biomp zoop_rel = 0.!nud*bct*biomp_p expomp_p = 0.!wwdz*biomp_p expomp_a = 0.!wwd*biomp_a agg_rel = 0.!nud*bct*biomp_a remiz = nud*bct*biodz expoz = wwdz*biodz ! 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) mpflag = 0.5 + sign(0.5,biomp - trcmin) mppflag = 0.5 + sign(0.5,biomp_p - trcmin) mpaflag = 0.5 + sign(0.5,biomp_a - trcmin) dzflag = 0.5 + sign(0.5,biodz - 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) 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 expomp_a = expomp_a*mpaflag expomp_p = expomp_p*mppflag zoop_upt = zoop_upt*zflag*mpflag zoop_rel = zoop_rel*mppflag agg_rel = agg_rel*mpaflag remiz = remiz*dzflag expoz = expoz*dzflag npp = npp*nflag*no3flag*pflag npp_D = npp_D*nflag*diazflag graz_D = graz_D*diazflag*zflag morp_D = morp_D*diazflag no3upt_D = no3upt_D*no3flag*diazflag 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 graz_C = graz_C*cflag*zflag morp_C = morp_C*cflag morpt_C = morpt_C*cflag avej = avej*pflag avej_D = avej_D*diazflag gmax = gmax*zflag no3P = no3P*no3flag po4P = po4P*nflag po4_D = po4_D*nflag !plastic taken up in phytoplankton aggregates ! assume 3% aggregation rate !!marine snow quantity in bio units- total S/S of MS mar_snow = ((1.-gamma1)* & (1.-rpoo)* & (graz + graz_Det & + (graz_Z & - bapr*graz_Z & ) & + graz_D & + (graz_C & - bapr*graz_C & ) & ) & + morp_D & + morp & + (morp_C & - bapr*morp_C & ) & + (morz & - bapr*morz & )) ! & + remi_B!*agg_con ! max potential uptake limited by availability agg_upt = 0.!min(mar_snow*agg_con*biomp/(kp+biomp),biomp)*mpflag ! agg_upt = (mar_snow*biomp/(kp+biomp))*mpflag ! Excretion is the difference between ! the zooplankton assimilation and growth efficiencies excr = (gamma1-geZ)*(graz+graz_Z+graz_Det & + graz_D & + graz_C & + graz_Det_B & ) ! nutrients equation bion = bion + tsb*redptn*(remi + excr - npp + morpt & + remiz & - npp_D & - npp_C + morpt_C & ) ! phytoplankton equation biop = biop + tsb*(npp - morp - graz - morpt) ! coccolithophores equation bioc = bioc + tsb*(npp_C - morp_C - graz_C & - morpt_C) ! zooplankton equation bioz = bioz + tsb*(geZ*(graz + graz_Det + graz_Z & + graz_D & + graz_C & + graz_Det_B & ) & - morz - graz_Z) ! detritus equations biod = biod + tsb*((1.-gamma1)* & (1.-rpoo)* & (graz + graz_Det & + (graz_Z & - bapr*graz_Z & ) & + graz_D & + (graz_C & - bapr*graz_C & ) & ) & + morp_D & + (morp_C & - bapr*morp_C & ) & + (morz & - bapr*morz & ) & + morp - remi - graz_Det - expo + impo) & + remi_B biodz = biodz + tsb*((1.-gamma1)* & rpoo* & (graz + graz_Det & + (graz_Z & - bapr*graz_Z & ) & + graz_D & + (graz_C & - bapr*graz_C & ) & ) & - remiz - expoz + impoz) biomp = biomp +tsb*(zoop_rel - zoop_upt + agg_rel-agg_upt) biomp_a = 0.!biomp_a + tsb*(impomp_a-expomp_a-agg_rel+agg_upt) biomp_p = 0.!biomp_p + tsb*(zoop_upt-zoop_rel+impomp_p-expomp_p) 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 + 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 & + remiz & - npp_C + morpt_C & ) ! diazotroph equation biodiaz = biodiaz + tsb*(npp_D - morp_D - graz_D) expoout = expoout + expo expompout_a = expompout_a + expomp_a expompout_p = expompout_p + expomp_p agg_uptout = agg_uptout + agg_upt zoop_uptout = zoop_uptout + zoop_upt agg_relout = agg_relout + agg_rel zoop_relout = zoop_relout + zoop_rel 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 expozout = expozout + expoz remizout = remizout + remiz 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 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 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(ibiompa) = max(0.,biomp_a) - bioin(ibiompa) ! bioout(ibiompp) = max(0.,biomp_p) - bioin(ibiompp) ! bioout(ibiomp) = max(0.,biomp) - bioin(ibiomp) bioout(ibiompa) = biomp_a - bioin(ibiompa) bioout(ibiompp) = biomp_p - bioin(ibiompp) bioout(ibiomp) = biomp - bioin(ibiomp) bioout(ibiodz) = biodz - bioin(ibiodz) bioout(ibiono3) = biono3 - bioin(ibiono3) bioout(ibiodiaz) = biodiaz - bioin(ibiodiaz) bioout(ibioc) = bioc - bioin(ibioc) bioout(ibiocaco3) = biocaco3 - bioin(ibiocaco3) bioout(ibiod_B) = biod_B - bioin(ibiod_B) return end