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 McGillicuddy 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,irrad,rl,WsnkTz,dep #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 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 k=1,km do 10 i=2,imtm1 fac3=max(c0,tb(i,k,ino3)) fac4=max(c0,tb(i,k,iphy)) fac6=max(c0,tb(i,k,inh4)) if (k.eq.1) irrad0(i)=srf(i)*parfrac fac1=(atth2o1+attphy*fac4)*dzqz(i,k,0) irrad1(i)=irrad0(i)*exp(-fac1) irrad=(irrad0(i)-irrad1(i))/fac1 fac1=exp(-photor0 *irrad/photorm) fac2=exp(-photoinh*irrad/photorm) rl=photorm*(c1-fac1)*fac2 q1(i,k)=fac3*exp(-no3inh*fac6)/(hsno3+fac3) q2(i,k)=fac6/(hsnh4+fac6) no3pr(i,k)=fac4*rl*q1(i,k) nh4pr(i,k)=fac4*rl*q2(i,k) Tsrc(i,k)=-no3pr(i,k) irrad0(i)=irrad1(i) 10 continue c c----------------------------------------------------------------------- c Phytoplankton (autotroph) source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.iphy) then do 40 i=2,imtm1 c c f-ratio dependent sinking rate. c if(sinkfac.ne.c0) then fave(i)=c0 k=0 20 continue k=k+1 dep=tdepth(i,k,jrs)+p5*dzqz(i,k,jrs) if(dep.lt.sbldep) then if((q1(i,k)+q2(i,k)).gt.c0) then fave(i)=fave(i)+dzqz(i,k,jrs)* * (q1(i,k)/(q1(i,k)+q2(i,k))) endif goto 20 else if((q1(i,k)+q2(i,k)).gt.c0) then fave(i)=fave(i)+(sbldep-dep)* * (q1(i,k)/(q1(i,k)+q2(i,k))) endif endif fave(i)=fave(i)/sbldep endif wsnkphy=sinkphy+fave(i)*sinkfac c do 30 k=1,km fac4=max(c0,tb(i,k,iphy)) fac5=max(c0,tb(i,k,izoo)) #ifndef barotropic if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkphy*(tb(i,k-1,iphy)-tb(i,k,iphy))* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then WsnkTz=-wsnkphy*tb(i,1,iphy)*c2*dz2rqz(i,1,0) else WsnkTz=wsnkphy*tb(i,kmm1,iphy)*c2*dz2rqz(i,km,0) endif # else WsnkTz = c0 #endif zgrphy(i,k)=fac5*grazrm*(c1-exp(-civlev*fac4)) Tsrc(i,k)=no3pr(i,k)+nh4pr(i,k)-zgrphy(i,k)+WsnkTz 30 continue 40 continue c c----------------------------------------------------------------------- c Zooplankton (heterotroph) source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.izoo) then do 50 k=1,km do 50 i=2,imtm1 fac5=max(c0,tb(i,k,izoo)) Tsrc(i,k)=(c1-zooexc1)*zgrphy(i,k) * -fac5*(zoolr1+zoolr2*fac5) 50 continue c c----------------------------------------------------------------------- c Ammonia source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.inh4) then do 60 k=1,km do 60 i=2,imtm1 fac5=max(c0,tb(i,k,izoo)) Tsrc(i,k)=-nh4pr(i,k)+zooexc1*zgrphy(i,k) * +fac5*((c1-zooexp1)*zoolr1+(c1-zooexp2)*zoolr2*fac5) c c Keep track of export term. c export(i,k)=fac5*(zooexp1*zoolr1+zooexp2*zoolr2*fac5) 60 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 70 k = 1, km do 70 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) 70 continue c c Depth Response Velocities. c do 80 k = 1, km do 80 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) 80 continue c c Prey Response Velocities. c do 90 k = 1, km do 90 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 90 continue c c Total Response Fluxes. c do 100 k = 1, km do 100 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)) 100 continue c c Cod Response flux. c do 110 k = 1, km do 110 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) 110 continue c c Hardwire cod to the floor. c # ifndef codvadvect do 120 k = 1, km do 120 i = 1, imt WTz(i,k) = c0 # ifndef codvmix Tzz(i,k) = c0 # endif 120 continue # else fx = cdwmax*c2 do 130 i = 1, imt do 120 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 120 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 130 continue # endif end if c #endif return end