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 Dusenberry Biological Model. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #ifdef rmdenbar # include #endif #include #include #include #include #if defined codunlim | defined codlim | defined cod_ing # include #endif c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,j,k,m FLOAT #ifdef dusDroop * fac3, #else cjad * chlcrtio,chlcrate, * chlcrtio,photosrc, #endif * fac4,fac4km1,fac5,fac6,fac7,fac7km1,fac8,fac8km1,fac9, * fac9km1, * wsnkphy2,WsnkTz,zgrfac #if ((defined codunlim | defined codlim) & defined codvadvect) | defined cod_ing * ,fx #endif #if ( defined codunlim | defined codlim ) & ! defined cod_ing * ,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 #if defined cod_ing FLOAT * cupw(imt,km),cutw(imt,km),cvts(imt,km),cvtn(imt,km) * ,cfuw(imt,km),cfvn(imt,km),cfvs(imt,km) * ,Ucx(imt,km),Vcy(imt,km),Wcz(imt,km),cw(imt,kmp1) * ,fy,fc,fxa,fxb # ifdef zoo_ing * ,zupw(imt,km),zutw(imt,km),zvts(imt,km),zvtn(imt,km) # endif #endif c c move below to blkdat.F at some point c parameter(chlcrate=3.89e-5) c c======================================================================= c Begin executable code. c======================================================================= c c----------------------------------------------------------------------- c Nitrate source/sink term. c----------------------------------------------------------------------- c if(m.eq.ino3) then #ifdef dusDroop c c ---------------------------------------------------------------- c Droop equations for nitrogen uptake/growth c ---------------------------------------------------------------- c do 10 k=1,km do 10 i=2,imtm1 if(biopos.eq.2) then fac3=max(c0,tb(i,k,ino3)) fac4=max(c0,tb(i,k,ichl)) fac6=max(c0,tb(i,k,inh4)) else fac3=tb(i,k,ino3) fac4=tb(i,k,ichl) fac6=tb(i,k,inh4) endif Tsrc(i,k)=-no3urm*fac4*fac3*exp(-no3inh*fac6)/(hsno3+fac3)+ * remnnh4*fac6 10 continue # else c c ---------------------------------------------------------------- c Constant (?) chl:N ratio (default) c ---------------------------------------------------------------- c do 10 k=1,km do 10 i=2,imtm1 if(biopos.eq.2) then fac6=max(c0,tb(i,k,inh4)) else fac6=tb(i,k,inh4) endif Tsrc(i,k)=-no3pr(i,k)+remnnh4*fac6 10 continue #endif c c----------------------------------------------------------------------- c Phytoplankton nitrate source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.iqn3) then #ifdef dusDroop c c ---------------------------------------------------------------- c Droop equations for nitrogen uptake/growth c ---------------------------------------------------------------- c do 20 k=1,km do 20 i=2,imtm1 if(biopos.eq.2) then fac3=max(c0,tb(i,k,ino3)) fac4=max(c0,tb(i,k,ichl)) fac5=max(c0,tb(i,k,izoo)) fac6=max(c0,tb(i,k,inh4)) fac7=max(c0,tb(i,k,iqn3)) if(k.gt.1) fac7km1=max(c0,tb(i,k-1,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac3=tb(i,k,ino3) fac4=tb(i,k,ichl) fac5=tb(i,k,izoo) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) if(k.gt.1) fac7km1=tb(i,k-1,iqn3) fac9=tb(i,k,iqn4) endif if(wsnkphy.lt.c0) then if(fac4.gt.c0) then wsnkphy2=c2*(chl2rad*fac4**chl2rade)*grav*cm2m* * (rhophy-rhos(i,k) # ifdef rmdenbar * -rhobar(i,j,k) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkphy2*(fac7km1-fac7)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkphy2*fac7*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkphy2*fac7km1*c2*dz2rqz(i,km,0) endif if((fac7+fac9).gt.0) then zgrfac=fac5*grazrm*fac7/(fac7+fac9) * *(c1-exp(-civlev*(fac7+fac9))) else zgrfac=c0 endif Tsrc(i,k)=no3urm*fac4*fac3*exp(-no3inh*fac6)/(hsno3+fac3) * -zgrfac+WsnkTz * -fac7*(phylr1+(fac7+fac9)*phylr2) 20 continue # else c c ---------------------------------------------------------------- c Constant (?) Chl:N ratio (default) c ---------------------------------------------------------------- c do 20 k=1,km do 20 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) if(k.gt.1) fac7km1=max(c0,tb(i,k-1,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac4=tb(i,k,ichl) fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) if(k.gt.1) fac7km1=tb(i,k-1,iqn3) fac9=tb(i,k,iqn4) endif if(wsnkphy.lt.c0) then if(fac4.gt.c0) then wsnkphy2=c2*(chl2rad*fac4**chl2rade)*grav*cm2m* * (rhophy-rhos(i,k) # ifdef rmdenbar * -rhobar(i,j,k) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkphy2*(fac7km1-fac7)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkphy2*fac7*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkphy2*fac7km1*c2*dz2rqz(i,km,0) endif if((fac7+fac9).gt.0) then zgrfac=fac5*grazrm*fac7/(fac7+fac9) * *(c1-exp(-civlev*(fac7+fac9))) else zgrfac=c0 endif Tsrc(i,k)=no3pr(i,k) * -zgrfac+WsnkTz * -fac7*(phylr1+(fac7+fac9)*phylr2) 20 continue #endif c c----------------------------------------------------------------------- c Zooplankton (heterotroph) source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.izoo) then #ifdef dusDroop c c ---------------------------------------------------------------- c Droop equations for nitrogen uptake/growth c ---------------------------------------------------------------- c do 30 k=1,km do 30 i=2,imtm1 if(biopos.eq.2) then fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif zgrphy(i,k)=fac5*grazrm*(c1-exp(-civlev*(fac7+fac9))) Tsrc(i,k)=(c1-zooexcn-zooexcd)*zgrphy(i,k) * -fac5*(zoolr1+zoolr2*fac5) 30 continue # else c c ---------------------------------------------------------------- c Constant (?) Chl:N ratio (default) c ---------------------------------------------------------------- c do 30 k=1,km do 30 i=2,imtm1 if(biopos.eq.2) then fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif zgrphy(i,k)=fac5*grazrm*(c1-exp(-civlev*(fac7+fac9))) Tsrc(i,k)=(c1-zooexcn-zooexcd)*zgrphy(i,k) * -fac5*(zoolr1+zoolr2*fac5) 30 continue #endif c c----------------------------------------------------------------------- c Ammonia source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.inh4) then #ifdef dusDroop c c ---------------------------------------------------------------- c Droop equations for nitrogen uptake/growth c ---------------------------------------------------------------- c do 40 k=1,km do 40 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) fac5=max(c0,tb(i,k,izoo)) fac6=max(c0,tb(i,k,inh4)) fac7=max(c0,tb(i,k,iqn3)) fac8=max(c0,tb(i,k,idet)) fac9=max(c0,tb(i,k,iqn4)) else fac4=tb(i,k,ichl) fac5=tb(i,k,izoo) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) fac8=tb(i,k,idet) fac9=tb(i,k,iqn4) endif zgrfac=fac5*grazrm*(c1-exp(-civlev*(fac7+fac9))) Tsrc(i,k)=-nh4urm*fac4*fac6/(hsnh4+fac6) * +zooexcn*zgrfac * +fac5*((c1-zoolf1)*zoolr1+(c1-zoolf2)*zoolr2*fac5) * -remnnh4*fac6+remndet*fac8 40 continue # else c c ---------------------------------------------------------------- c Constant (?) Chl:N ratio (default) c ---------------------------------------------------------------- c do 40 k=1,km do 40 i=2,imtm1 if(biopos.eq.2) then fac5=max(c0,tb(i,k,izoo)) fac6=max(c0,tb(i,k,inh4)) fac7=max(c0,tb(i,k,iqn3)) fac8=max(c0,tb(i,k,idet)) fac9=max(c0,tb(i,k,iqn4)) else fac5=tb(i,k,izoo) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) fac8=tb(i,k,idet) fac9=tb(i,k,iqn4) endif zgrfac=fac5*grazrm*(c1-exp(-civlev*(fac7+fac9))) Tsrc(i,k)=-nh4pr(i,k) * +zooexcn*zgrfac * +fac5*((c1-zoolf1)*zoolr1+(c1-zoolf2)*zoolr2*fac5) * -remnnh4*fac6+remndet*fac8 40 continue #endif c c----------------------------------------------------------------------- c Detritus source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.idet) then #ifdef dusDroop c c ---------------------------------------------------------------- c Droop equations for nitrogen uptake/growth c ---------------------------------------------------------------- c do 50 k=1,km do 50 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac8=max(c0,tb(i,k,idet)) if(k.gt.1) fac8km1=max(c0,tb(i,k-1,idet)) fac9=max(c0,tb(i,k,iqn4)) else fac4=tb(i,k,ichl) fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac8=tb(i,k,idet) if(k.gt.1) fac8km1=tb(i,k-1,idet) fac9=tb(i,k,iqn4) endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkdet*(fac8km1-fac8)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkdet*fac8*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkdet*fac8km1*c2*dz2rqz(i,km,0) endif zgrfac=fac5*grazrm*(c1-exp(-civlev*(fac7+fac9))) Tsrc(i,k)=-remndet*fac8+zooexcd*zgrfac * +fac5*(zoolf1*zoolr1+zoolf2*zoolr2*fac5)+WsnkTz * +(fac7+fac9)*(phylr1+(fac7+fac9)*phylr2) 50 continue # else c c ---------------------------------------------------------------- c Constant (?) Chl:N ratio (default) c ---------------------------------------------------------------- c do 50 k=1,km do 50 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac8=max(c0,tb(i,k,idet)) if(k.gt.1) fac8km1=max(c0,tb(i,k-1,idet)) fac9=max(c0,tb(i,k,iqn4)) else fac4=tb(i,k,ichl) fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac8=tb(i,k,idet) if(k.gt.1) fac8km1=tb(i,k-1,idet) fac9=tb(i,k,iqn4) endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkdet*(fac8km1-fac8)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkdet*fac8*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkdet*fac8km1*c2*dz2rqz(i,km,0) endif zgrfac=fac5*grazrm*(c1-exp(-civlev*(fac7+fac9))) Tsrc(i,k)=-remndet*fac8+zooexcd*zgrfac * +fac5*(zoolf1*zoolr1+zoolf2*zoolr2*fac5)+WsnkTz * +(fac7+fac9)*(phylr1+(fac7+fac9)*phylr2) 50 continue #endif c c----------------------------------------------------------------------- c Phytoplankton chlorophyll source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.ichl) then #ifdef dusDroop c c ---------------------------------------------------------------- c Droop equations for nitrogen uptake/growth c ---------------------------------------------------------------- c do 60 k=1,km do 60 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) if(k.gt.1) fac4km1=max(c0,tb(i,k-1,ichl)) fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac4=tb(i,k,ichl) if(k.gt.1) fac4km1=tb(i,k-1,ichl) fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif if(wsnkphy.lt.c0) then if(fac4.gt.c0) then wsnkphy2=c2*(chl2rad*fac4**chl2rade)*grav*cm2m* * (rhophy-rhos(i,k) # ifdef rmdenbar * -rhobar(i,j,k) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkphy2*(fac4km1-fac4)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkphy2*fac4*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkphy2*fac4km1*c2*dz2rqz(i,km,0) endif if((fac7+fac9).gt.0) then zgrfac=fac5*grazrm*fac4/(fac7+fac9) * *(c1-exp(-civlev*(fac7+fac9))) else zgrfac=c0 endif if((fac7+fac9).gt.0) then Tsrc(i,k)=(no3pr(i,k)+nh4pr(i,k))*fac4/(fac7+fac9) * -zgrfac+WsnkTz * -fac4*(phylr1+(fac7+fac9)*phylr2) else c chl with no biomass attached?? guard against this... Tsrc(i,k)=WsnkTz-fac4*(phylr1+(fac7+fac9)*phylr2) endif 60 continue # else c c ---------------------------------------------------------------- c Constant (?) Chl:N ratio (default) c ---------------------------------------------------------------- c do 60 k=1,km do 60 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) if(k.gt.1) fac4km1=max(c0,tb(i,k-1,ichl)) fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac4=tb(i,k,ichl) if(k.gt.1) fac4km1=tb(i,k-1,ichl) fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif if(wsnkphy.lt.c0) then if(fac4.gt.c0) then wsnkphy2=c2*(chl2rad*fac4**chl2rade)*grav*cm2m* * (rhophy-rhos(i,k) # ifdef rmdenbar * -rhobar(i,j,k) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkphy2*(fac4km1-fac4)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkphy2*fac4*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkphy2*fac4km1*c2*dz2rqz(i,km,0) endif if((fac7+fac9).gt.0) then if((c2chl.lt.c0).and.(chlcinf(i,k).gt.c0)) then chlcrtio=fac4/(fac7+fac9)*c2n photosrc=chlcrate*fac4*(c1 * -(chlcrtio/chlcinf(i,k))) else photosrc=c0 endif zgrfac=fac5*grazrm*fac4/(fac7+fac9) * *(c1-exp(-civlev*(fac7+fac9))) Tsrc(i,k)=(no3pr(i,k)+nh4pr(i,k))*fac4/(fac7+fac9) * -zgrfac+WsnkTz * -fac4*(phylr1+(fac7+fac9)*phylr2) * +photosrc else zgrfac=c0 Tsrc(i,k)=WsnkTz-fac4*(phylr1+(fac7+fac9)*phylr2) endif 60 continue #endif c c----------------------------------------------------------------------- c Phytoplankton ammonia source/sink term. c----------------------------------------------------------------------- c elseif(m.eq.iqn4) then #ifdef dusDroop c c ---------------------------------------------------------------- c Droop equations for nitrogen uptake/growth c ---------------------------------------------------------------- c do 70 k=1,km do 70 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) fac5=max(c0,tb(i,k,izoo)) fac6=max(c0,tb(i,k,inh4)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) if(k.gt.1) fac9km1=max(c0,tb(i,k-1,iqn4)) else fac4=tb(i,k,ichl) fac5=tb(i,k,izoo) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) if(k.gt.1) fac9km1=tb(i,k-1,iqn4) endif if(wsnkphy.lt.c0) then if(fac4.gt.c0) then wsnkphy2=c2*(chl2rad*fac4**chl2rade)*grav*cm2m* * (rhophy-rhos(i,k) # ifdef rmdenbar * -rhobar(i,j,k) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkphy2*(fac9km1-fac9)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkphy2*fac9*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkphy2*fac9km1*c2*dz2rqz(i,km,0) endif if((fac7+fac9).gt.0) then zgrfac=fac5*grazrm*fac9/(fac7+fac9) * *(c1-exp(-civlev*(fac7+fac9))) else zgrfac=c0 endif Tsrc(i,k)=nh4urm*fac4*fac6/(hsnh4+fac6)- * zgrfac+WsnkTz * -fac9*(phylr1+(fac7+fac9)*phylr2) 70 continue endif # else c c ---------------------------------------------------------------- c Constant (?) Chl:N ratio (default) c ---------------------------------------------------------------- c do 70 k=1,km do 70 i=2,imtm1 if(biopos.eq.2) then fac4=max(c0,tb(i,k,ichl)) fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) if(k.gt.1) fac9km1=max(c0,tb(i,k-1,iqn4)) else fac4=tb(i,k,ichl) fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) if(k.gt.1) fac9km1=tb(i,k-1,iqn4) endif if(wsnkphy.lt.c0) then if(fac4.gt.c0) then wsnkphy2=c2*(chl2rad*fac4**chl2rade)*grav*cm2m* * (rhophy-rhos(i,k) # ifdef rmdenbar * -rhobar(i,j,k) # endif * )/(c9*visch2o) else wsnkphy2=c0 endif else wsnkphy2=wsnkphy endif if((k.gt.1).and.(k.lt.km)) then WsnkTz=wsnkphy2*(fac9km1-fac9)* * c2*dz2rqz(i,k,0) elseif(k.eq.1) then if(k.eq.km) then WsnkTz = c0 else WsnkTz=-wsnkphy2*fac9*c2*dz2rqz(i,1,0) endif else WsnkTz=wsnkphy2*fac9km1*c2*dz2rqz(i,km,0) endif if((fac7+fac9).gt.0) then zgrfac=fac5*grazrm*fac9/(fac7+fac9) * *(c1-exp(-civlev*(fac7+fac9))) else zgrfac=c0 endif Tsrc(i,k)=nh4pr(i,k) * -zgrfac+WsnkTz * -fac9*(phylr1+(fac7+fac9)*phylr2) 70 continue endif #endif c #if ( defined codunlim | defined codlim ) & ! defined cod_ing c----------------------------------------------------------------------- c Cod source/sink term. c----------------------------------------------------------------------- c if (m.eq.icod) then c c Temperature Response Velocities. c do 80 k = 1, km do 80 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) 80 continue c c Depth Response Velocities. c do 90 k = 1, km do 90 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) 90 continue c c Prey Response Velocities. c do 100 k = 1, km do 100 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 100 continue c c Total Response Fluxes. c do 110 k = 1, km do 110 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)) 110 continue c c Cod Response flux. c do 120 k = 1, km do 120 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) 120 continue c c Hardwire cod to the floor. c # ifndef codvadvect do 130 k = 1, km do 130 i = 1, imt WTz(i,k) = c0 # ifndef codvmix Tzz(i,k) = c0 # endif 130 continue # else fx = cdwmax*c2 do 140 i = 1, imt do 130 k = 2, kmm1 if (biopos.ne.2) then WTz(i,k) = WTz(i,k) + & fx*(t(i,k-1,icod)-t(i,k,icod))*dz2rqz(i,k,0) else WTz(i,k) = WTz(i,k) + & fx*(max(c0,t(i,k-1,icod))- & max(c0,t(i,k ,icod)))*dz2rqz(i,k,0) end if # ifndef codvmix Tzz(i,k) = c0 # endif 130 continue if (biopos.ne.2) then 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) else WTz(i,1) = WTz(i,1) - & fx*max(c0,t(i,1,icod))*dz2rqz(i,k,0) WTz(i,km) = WTz(i,km) + & fx*max(c0,t(i,kmm1,icod))*dz2rqz(i,k,0) end if # ifndef codvmix Tzz(i,1) = c0 Tzz(i,km) = c0 # endif 140 continue # endif end if c #elif defined cod_ing c----------------------------------------------------------------------- c Cod source/sink term. c----------------------------------------------------------------------- c if (m.eq.icod) then c c Temperature Response Velocities. c c c -- cod speed at u-points c do 80 k = 1, km do 80 i = 2, imt fx= p25*(t(i-1,k,1)+t(i,k,1)+tm(i,k,1)+tm(i-1,k,1)) # ifndef cod_square fx = cdspd*(fx-cdtmin)/cdtmax # else fx=sign(fx-cdtmin,cdspd*((fx-cdtmin)/cdtmax)**2) # endif if(abs(fx).gt.cdspd)fx=sign(fx,cdspd) cutw(i-1,k) = -fx fx = p25*(t(i-1,k,1)+t(i,k,1)+tp(i,k,1)+tp(i-1,k,1)) # ifndef cod_square fx = cdspd*(fx-cdtmin)/cdtmax # else fx=sign(fx-cdtmin,cdspd*((fx-cdtmin)/cdtmax)**2) # endif if(abs(fx).gt.cdspd)fx=sign(fx,cdspd) cupw(i-1,k) = -fx # ifdef zoo_ing CDKZ Minimum tolerable prey value CDKP Optimal prey value for Cod fx= p25*(t(i-1,k,izoo)+t(i,k,izoo)+ $ tm(i,k,izoo)+tm(i-1,k,izoo)) if(fx.lt.cdkz) then fx=cdspd elseif(fx.gt.cdkp) then fx=c1em4 else fx = cdspd*(fx-cdkz)/(cdkp-cdkz)+c1em4 endif zutw(i-1,k) = fx fx = p25*(t(i-1,k,izoo)+t(i,k,izoo)+ $ tp(i,k,izoo)+tp(i-1,k,izoo)) fx = cdspd*(fx-cdtmin)/cdtmax if(fx.lt.cdkz) then fx=cdspd elseif(fx.gt.cdkp) then fx=c1em4 else fx = cdspd*(fx-cdkz)/(cdkp-cdkz)+c1em4 endif zupw(i-1,k) = fx # endif 80 continue c c -- cod direction u-points c do 90 k = 1, km do 90 i = 2, imt fx= -t(i-1,k,1)+t(i,k,1)+tm(i,k,1)-tm(i-1,k,1) fy= t(i-1,k,1)+t(i,k,1)-tm(i,k,1)-tm(i-1,k,1) fc= sqrt(fx*fx+fy*fy) if(fc.ne.0) then fx=fx/fc fy=fy/fc else fx=1.0 fy=0.0 endif cvts(i-1,k)=cutw(i-1,k)*fy cutw(i-1,k)=cutw(i-1,k)*fx # ifdef zoo_ing fx= -t(i-1,k,izoo)+t(i,k,izoo)+ $ tm(i,k,izoo)-tm(i-1,k,izoo) fy= t(i-1,k,izoo)+t(i,k,izoo)- $ tm(i,k,izoo)-tm(i-1,k,izoo) fc= sqrt(fx*fx+fy*fy) if(fc.ne.0) then fx=fx/fc fy=fy/fc else fx=1.0 fy=0.0 endif zvts(i-1,k)=zutw(i-1,k)*fy zutw(i-1,k)=zutw(i-1,k)*fx c c -- combine enviromental influences c cvts(i-1,k)=(c1-cdkt)*cvts(i-1,k)+cdkt*zvts(i-1,k) cutw(i-1,k)=(c1-cdkt)*cutw(i-1,k)+cdkt*zutw(i-1,k) # endif fx= -t(i-1,k,1)+t(i,k,1)+tp(i,k,1)-tp(i-1,k,1) fy= -t(i-1,k,1)-t(i,k,1)+tp(i,k,1)+tp(i-1,k,1) fc= sqrt(fx*fx+fy*fy) if(fc.ne.0) then fx=fx/fc fy=fy/fc else fx=1.0 fy=0.0 endif cvtn(i-1,k)=cupw(i-1,k)*fy cupw(i-1,k)=cupw(i-1,k)*fx # ifdef zoo_ing fx= -t(i-1,k,izoo)+t(i,k,izoo)+tp(i,k,izoo)-tp(i-1,k,izoo) fy= -t(i-1,k,izoo)-t(i,k,izoo)+tp(i,k,izoo)+tp(i-1,k,izoo) fc= sqrt(fx*fx+fy*fy) if(fc.ne.0) then fx=fx/fc fy=fy/fc else fx=1.0 fy=0.0 endif zvtn(i-1,k)=zupw(i-1,k)*fy zupw(i-1,k)=zupw(i-1,k)*fx c c -- combine enviromental influences c cvtn(i-1,k)=(c1-cdkt)*cvtn(i-1,k)+cdkt*zvtn(i-1,k) cupw(i-1,k)=(c1-cdkt)*cupw(i-1,k)+cdkt*zupw(i-1,k) # endif 90 continue c c----------------------------------------------------------------------- c Compute the advective coefficients CFUW at west face of T grid box c CFVN at the north face of T grid box and CFVS at the south face c DO NOT BOOTSTRAP (yet) c----------------------------------------------------------------------- c fxa=cstr(j)*dytr(j) fxb=fxa*cs(j) do 110 k=1,km do 110 i=2,imt cfuw(i,k)=(cupw(i-1,k)*dyu(j )*dzvqz(i-1,k,0)+ * cutw(i-1,k)*dyu(j-1)*dzvqz(i-1,k,1))*fxa cfvn(i,k)=(cvtn(i ,k)*dxuq(i ,k)*dzvqz(i ,k,0)+ * cvtn(i-1,k)*dxuq(i-1,k)*dzvqz(i-1,k,0))*fxb*dxt4rq(i,k) cfvs(i,k)=(cvts(i ,k)*dxuq(i ,k)*dzvqz(i ,k,1)+ * cvts(i-1,k)*dxuq(i-1,k)*dzvqz(i-1,k,1))*fxb*dxt4rq(i,k) 110 continue c c----------------------------------------------------------------------- c Compute "omega" vertical velocity in T columns. c----------------------------------------------------------------------- c c Set "omega" vertical velocity at the surface to zero. c do 120 i=2,imtm1 cw(i,1)=c0 120 continue c c Compute omega due to horizontal velocity contributions c do 140 k=1,km-1 do 140 i=2,imtm1 fx=(cupw(i-1,k)+cupw(i,k)+cutw(i-1,k)+cutw(i,k)) * *cstr(j)*dxt4rq(i,k)*p5* * ( (vdepth(i ,k,jrn)- * vdepth(i-1,k,jrn) )+ * (vdepth(i ,k,jrs) - * vdepth(i-1,k,jrs) ) ) + * (cvtn(i-1,k)+cvtn(i,k)+cvts(i-1,k)+cvts(i,k)) * *dyt4r(j)*p5* * ( (vdepth(i ,k,jrn)- * vdepth(i ,k,jrs)) + * (vdepth(i-1,k,jrn)- * vdepth(i-1,k,jrs)) ) fxa=(cupw(i-1,k+1)+cupw(i,k+1)+cutw(i-1,k+1)+cutw(i,k+1)) * *cstr(j)*dxt4rq(i,k+1)*p5* * ( (vdepth(i ,k+1,jrn)- * vdepth(i-1,k+1,jrn) )+ * (vdepth(i ,k+1,jrs) - * vdepth(i-1,k+1,jrs) ) ) + * (cvtn(i-1,k+1)+cvtn(i,k+1)+cvts(i-1,k+1)+cvts(i,k+1)) * *dyt4r(j)*p5* * ( (vdepth(i ,k+1,jrn)- * vdepth(i ,k+1,jrs)) + * (vdepth(i-1,k+1,jrn)- * vdepth(i-1,k+1,jrs)) ) cw(i,k+1)=p5*(fx+fxa) 140 continue c c Set bottom omega velocity to zero c k=km-1 do 150 i=2,imtm1 cw(i,km+1)= * (cupw(i-1,k+1)+cupw(i,k+1)+cutw(i-1,k+1)+cutw(i,k+1)) * *cstr(j)*dxt4rq(i,k+1)*p5* * ( (vdepth(i ,k+1,jrn)- * vdepth(i-1,k+1,jrn) )+ * (vdepth(i ,k+1,jrs) - * vdepth(i-1,k+1,jrs) ) ) + * (cvtn(i-1,k+1)+cvtn(i,k+1)+cvts(i-1,k+1)+cvts(i,k+1)) * *dyt4r(j)*p5* * ( (vdepth(i ,k+1,jrn)- * vdepth(i ,k+1,jrs)) + * (vdepth(i-1,k+1,jrn)- * vdepth(i-1,k+1,jrs)) ) 150 continue c c----------------------------------------------------------------------- c Compute total advection of tracers. c----------------------------------------------------------------------- c c Compute zonal advection of tracer. c do 200 k=1,km do 200 i=2,imtm1 Ucx(i,k)=(cfuw(i ,k)*(t(i ,k,m)+t(i-1,k,m))- * cfuw(i+1,k)*(t(i+1,k,m)+t(i ,k,m))) * *dxt4rq(i,k)/dzqz(i,k,0) 200 continue c c Compute meridional advection of tracer. c do 210 k=1,km do 210 i=2,imtm1 Vcy(i,k)=( cfvs(i,k)*(t (i,k,m)+tm(i,k,m)) * -cfvn (i,k)*(tp(i,k,m)+t (i,k,m))) * /dzqz(i,k,0) 210 continue c c Compute vertical advection of tracer. c do 220 k=2,km do 220 i=2,imtm1 tempb(i,k)=cw(i,k)*(t(i,k-1,m)+t(i,k,m)) 220 continue do 230 i=2,imtm1 tempb(i, 1)=cw(i, 1)*t(i, 1,m) tempb(i,kmp1)=cw(i,kmp1)*t(i,km,m) 230 continue c do 240 k=1,km do 240 i=2,imtm1 Wcz(i,k)=(tempb(i,k+1)-tempb(i,k))*dz2rqz(i,k,0) 240 continue c c Cod Response flux. c do 250 k = 1, km do 250 i = 2, imtm1 Tsrc(i,k) = Ucx(i,k)+Vcy(i,k)+Wcz(i,k) 250 continue c c Hardwire cod to the floor. c # ifndef codvadvect do 1300 k = 1, km do 1300 i = 1, imt WTz(i,k) = c0 # ifndef codvmix Tzz(i,k) = c0 # endif 1300 continue # else fx = cdwmax*c2 do 1400 i = 1, imt do 1300 k = 2, kmm1 if (biopos.ne.2) then WTz(i,k) = WTz(i,k) + & fx*(t(i,k-1,icod)-t(i,k,icod))*dz2rqz(i,k,0) else WTz(i,k) = WTz(i,k) + & fx*(max(c0,t(i,k-1,icod))- & max(c0,t(i,k ,icod)))*dz2rqz(i,k,0) end if # ifndef codvmix Tzz(i,k) = c0 # endif 1300 continue if (biopos.ne.2) then 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) else WTz(i,1) = WTz(i,1) - & fx*max(c0,t(i,1,icod))*dz2rqz(i,k,0) WTz(i,km) = WTz(i,km) + & fx*max(c0,t(i,kmm1,icod))*dz2rqz(i,k,0) end if # ifndef codvmix Tzz(i,1) = c0 Tzz(i,km) = c0 # endif 1400 continue # endif end if c #endif return end