subroutine biosource(j,m) c c======================================================================= c === c This routine computes the tracers source/sink term Tsrc at row J === c and for tracer M. Notice that the biological tracers are forced === c to be non-negative. If the switch BIOPOS=1, they are forced to === c be non-negative in routine STEP. === c === c This particular version computes the biological source terms for === c the Fasham et al. Biological Model. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #include #include #if defined codunlim | defined codlim # include #endif c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,j,k,m FLOAT * fac1,fac2,fac3,fac4,fac5,fac6,fac7,fac8,fac9,remnfac, * WsnkTz,fac0a,fac1a,pref1,pref2,pref3,deptop,depbot,rl #if (defined codunlim | defined codlim) & defined codvadvect * ,fx #endif FLOAT * irrad0(imt),irrad1(imt) #if defined codunlim | defined codlim * ,cuw(imt,km),cvs(imt,km),cvn(imt,km) * ,cupw(imt,km),cvps(imt,km),cvpn(imt,km) * ,cutw(imt,km),cvts(imt,km),cvtn(imt,km) * ,cuzw(imt,km),cvzs(imt,km),cvzn(imt,km) FLOAT * cdrspv #endif c parameter(remnfac=-0.858) c c======================================================================= c Begin executable code. c======================================================================= c c----------------------------------------------------------------------- c Nitrate source/sink term. Also compute recycled production. c----------------------------------------------------------------------- c if(m.eq.ino3) then do 10 i=2,imtm1 do 10 k=1,km fac6=max(c0,tb(i,k,inh4)) if(tdepth(i,k,jrs).le.sbldep) then fac3=max(c0,tb(i,k,ino3)) fac4=max(c0,tb(i,k,iphy)) if (k.eq.1) irrad0(i)=srf(i)*parfrac fac1=(atth2o+attphy*fac4)*dzqz(i,k,0) fac2=grthpa*(grthpb**(grthpc*tb(i,k,1))) irrad1(i)=irrad0(i)*exp(-fac1) c c Option 1: integrate I(z) over box depth (= IRRAD), and use in c equation for RL. Not integrated over timestep. c c irrad=(irrad0(i)-irrad1(i))/fac1 c rl=(fac2*phygrth*irrad)/ c * (sqrt(fac2*fac2+phygrth*phygrth*irrad*irrad)) c c Option 2: use I(z) in the equation for rl(z), and integrate c RL(z) over box depth. Not integrated over timestep. c fac0a=irrad0(i)+ * sqrt(irrad0(i)*irrad0(i)+fac2*fac2/(phygrth*phygrth)) fac1a=irrad1(i)+ * sqrt(irrad1(i)*irrad1(i)+fac2*fac2/(phygrth*phygrth)) rl=(fac2/fac1)*log(fac0a/fac1a) no3pr(i,k)=fac4*rl*fac3*exp(-no3inh*fac6)/(hsno3+fac3) nh4pr(i,k)=fac4*rl*fac6/(hsnh4+fac6) Tsrc(i,k)=-no3pr(i,k) irrad0(i)=irrad1(i) else Tsrc(i,k)=rmnnh4*fac6 endif 10 continue c c----------------------------------------------------------------------- c Phytoplankton source/sink term. Also compute zooplankton grazing. c----------------------------------------------------------------------- c elseif(m.eq.iphy) then do 20 i=2,imtm1 do 20 k=1,km fac4=max(c0,tb(i,k,iphy)) if(tdepth(i,k,jrs).le.sbldep) then fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,ipon)) fac9=max(c0,tb(i,k,ibac)) fac1=prefphy*fac4+prefbac*fac9+prefpon*fac7 pref1=prefphy*fac4/fac1 pref2=prefbac*fac9/fac1 pref3=prefpon*fac7/fac1 fac2=zgrazr*fac5/(zgrazk+pref1*fac4+pref2*fac9+pref3*fac7) zgrphy(i,k)=fac2*pref1*fac4 zgrbac(i,k)=fac2*pref2*fac9 zgrpon(i,k)=fac2*pref3*fac7 Tsrc(i,k)=(c1-phyexud)*(no3pr(i,k)+nh4pr(i,k)) * -zgrphy(i,k)-fac4*phylr1 else Tsrc(i,k)=-rmnbio*fac4 endif 20 continue c c----------------------------------------------------------------------- c Zooplankton source/sink term. Also compute export flux. c----------------------------------------------------------------------- c elseif(m.eq.izoo) then do 30 i=2,imtm1 do 30 k=1,km fac5=max(c0,tb(i,k,izoo)) if(tdepth(i,k,jrs).le.sbldep) then Tsrc(i,k)=zasseff*(zgrphy(i,k)+zgrbac(i,k)+zgrpon(i,k)) * -(zexcrr+zmortr)*fac5 else Tsrc(i,k)=-rmnbio*fac5 endif 30 continue c c Compute fast-sinking particle flux at base of tracer boxes in c 10^-5 moles nitrogen m-2 second-1 (fpflux), and remineralization c rate in micromoles nitrogen liter-1 second-1 (fpremn). c Remineralize the residual flux in the bottom box. c FPFLUX(imt,1:kmp1) is defined on the wT grid (i.e. box interfaces). c do 60 i=2,imtm1 fpflux(i,1)=c0 deptop=c0 do 50 k=1,km if(tdepth(i,k,jrs).le.sbldep) then fac5=max(c0,tb(i,k,izoo)) fpflux(i,k+1)=fpflux(i,k)+ * zmortex*zmortr*fac5*dzqz(i,k,jrs) fpremn(i,k)=c0 deptop=tdepth(i,k,jrs)+p5*dzqz(i,k,jrs) else depbot=tdepth(i,k,jrs)+p5*dzqz(i,k,jrs) fpflux(i,k+1)=fpflux(i,k)*((depbot/deptop)**remnfac) fpremn(i,k)=(fpflux(i,k)-fpflux(i,k+1))/dzqz(i,k,jrs) deptop=depbot endif 50 continue fpremn(i,km)=fpremn(i,km)+fpflux(i,kmp1)/dzqz(i,km,jrs) fpflux(i,kmp1)=c0 60 continue c c----------------------------------------------------------------------- c Ammonia source/sink term. Also compute bacterial uptake. c----------------------------------------------------------------------- c elseif(m.eq.inh4) then do 70 i=2,imtm1 do 70 k=1,km fac5=max(c0,tb(i,k,izoo)) fac6=max(c0,tb(i,k,inh4)) fac8=max(c0,tb(i,k,idon)) fac9=max(c0,tb(i,k,ibac)) if(tdepth(i,k,jrs).le.sbldep) then fac1=min(fac6,bacfrac*fac8) fac2=bacgrr*fac9/(bacgrk+fac1+fac8) bgrdon(i,k)=fac2*fac8 bgrnh4(i,k)=fac2*fac1 Tsrc(i,k)=-nh4pr(i,k)-bgrnh4(i,k) * +(zexcnh4*zexcrr+(c1-zmortex)*zmortr)*fac5 * +bexcrr*fac9 else fac4=max(c0,tb(i,k,iphy)) fac7=max(c0,tb(i,k,ipon)) Tsrc(i,k)=-rmnnh4*fac6+rmnbio*(fac4+fac5+fac8+fac9) * +fac7*(rmnpon1+fac7*rmnpon2)+fpremn(i,k) endif 70 continue c c----------------------------------------------------------------------- c Particulate organic nitrogen source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.ipon) then do 80 i=2,imtm1 do 80 k=1,km fac7=max(c0,tb(i,k,ipon)) #ifndef barotropic if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkpon*(tb(i,k-1,ipon)-tb(i,k,ipon))* * c2*dz2rqz(i,k,0) elseif (k.eq.1) then WsnkTz=-wsnkpon*tb(i,1,ipon)*c2*dz2rqz(i,1,0) else WsnkTz=wsnkpon*tb(i,kmm1,ipon)*c2*dz2rqz(i,km,0) endif # else WsnkTz = c0 #endif if(tdepth(i,k,jrs).le.sbldep) then fac4=max(c0,tb(i,k,iphy)) Tsrc(i,k)=(c1-zasseff)*(zgrphy(i,k)+zgrbac(i,k)+zgrpon(i,k)) * +phylr1*fac4-zgrpon(i,k)-brknpon*fac7+WsnkTz else Tsrc(i,k)=-fac7*(rmnpon1+fac7*rmnpon2)+WsnkTz endif 80 continue c c----------------------------------------------------------------------- c Dissolved organic nitrogen source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.idon) then do 90 i=2,imtm1 do 90 k=1,km if(tdepth(i,k,jrs).le.sbldep) then fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,ipon)) Tsrc(i,k)=phyexud*(no3pr(i,k)+nh4pr(i,k)) * +(c1-zexcnh4)*zexcrr*fac5+brknpon*fac7-bgrdon(i,k) else fac8=max(c0,tb(i,k,idon)) Tsrc(i,k)=-rmnbio*fac8 endif 90 continue c c----------------------------------------------------------------------- c Bacteria source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.ibac) then do 100 i=2,imtm1 do 100 k=1,km fac9=max(c0,tb(i,k,ibac)) if(tdepth(i,k,jrs).le.sbldep) then Tsrc(i,k)=bgrnh4(i,k)+bgrdon(i,k)-zgrbac(i,k)-bexcrr*fac9 else Tsrc(i,k)=-rmnbio*fac9 endif 100 continue endif c #if defined codunlim | defined codlim c----------------------------------------------------------------------- c Cod source/sink term. c----------------------------------------------------------------------- c if (m.eq.icod) then c c Temperature Response Velocities. c do 110 k = 1, km do 110 i = 2, imt cutw(i,k) = cdrspv (t(i-1,k,1),t(i,k,1),cdtmin,cdtmax, & cstr(j)*dxur(i-1),cdkt) if (j.eq.2) then cvts(i,k) = cdrspv (tm(i,k,1),t(i,k,1),cdtmin,cdtmax, & dyur(j-1),cdkt) end if cvtn(i,k) = cdrspv (t(i,k,1),tp(i,k,1),cdtmin,cdtmax, & dyur(j),cdkt) 110 continue c c Depth Response Velocities. c do 120 k = 1, km do 120 i = 2, imt cuzw(i,k) = cdrspv (hd(i-1,j),hd(i,j),cdzmin,cdzmax, & cstr(j)*dxur(i-1),cdkz) if (j.eq.2) then cvzs(i,k) = cdrspv (hd(i,j-1),hd(i,j),cdzmin,cdzmax, & dyur(j-1),cdkz) end if cvzn(i,k) = cdrspv (hd(i,j),hd(i,j+1),cdzmin,cdzmax, & dyur(j),cdkz) 120 continue c c Prey Response Velocities. c do 130 k = 1, km do 130 i = 2, imt cupw(i,k) = (t(i,k,izoo)-t(i-1,k,izoo))* & cstr(j)*dxur(i-1)*cdkp if (j.eq.2) then cvps(i,k) = (t(i,k,izoo)-tm(i,k,izoo))*dyur(j-1)*cdkp end if cvpn(i,k) = (tp(i,k,izoo)-t(i,k,izoo))*dyur(j)*cdkp # ifdef codlim cupw(i,k) = tanh(cupw(i,k)) if (j.eq.2) cvps(i,k) = tanh(cvps(i,k)) cvpn(i,k) = tanh(cvpn(i,k)) # endif 130 continue c c Total Response Fluxes. c do 140 k = 1, km do 140 i = 2, imt cuw(i,k) = cupw(i,k) + cutw(i,k) + cuzw(i,k) # ifdef codlim cuw(i,k) = tanh(cuw(i,k))*cdspd # endif cuw(i,k) = cuw(i,k)*(t(i,k,icod)+t(i-1,k,icod)) c if (j.gt.2) then cvs(i,k) = cvn(i,k) else cvs(i,k) = cvps(i,k) + cvts(i,k) + cvzs(i,k) # ifdef codlim cvs(i,k) = tanh(cvs(i,k))*cdspd # endif cvs(i,k) = cvs(i,k)*(t(i,k,icod)+tm(i,k,icod)) end if c cvn(i,k) = cvpn(i,k) + cvtn(i,k) + cvzn(i,k) # ifdef codlim cvn(i,k) = tanh(cvn(i,k))*cdspd # endif cvn(i,k) = cvn(i,k)*(t(i,k,icod)+tp(i,k,icod)) 140 continue c c Cod Response flux. c do 150 k = 1, km do 150 i = 2, imtm1 Tsrc(i,k) = (cuw(i-1,k)-cuw(i,k))*cstr(j)*dxt2r(i) + & (cvs(i,k)-cvn(i,k))*dyt2r(j) 150 continue c c Hardwire cod to the floor. c # ifndef codvadvect do 160 k = 1, km do 160 i = 1, imt WTz(i,k) = c0 # ifndef codvmix Tzz(i,k) = c0 # endif 160 continue # else fx = cdwmax*c2 do 170 i = 1, imt do 160 k = 2, kmm1 WTz(i,k) = WTz(i,k) + & fx*(t(i,k-1,icod)-t(i,k,icod))*dz2rqz(i,k,0) # ifndef codvmix Tzz(i,k) = c0 # endif 160 continue WTz(i,1) = WTz(i,1) - fx*t(i,1,icod)*dz2rqz(i,k,0) WTz(i,km) = WTz(i,km) + fx*t(i,kmm1,icod)*dz2rqz(i,k,0) # ifndef codvmix Tzz(i,1) = c0 Tzz(i,km) = c0 # endif 170 continue # endif end if c #endif return end