subroutine priprod c c======================================================================= c === c This routine computes the primary productivity (both no3pr and === c nh4pr) for the current row. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,k,l FLOAT * fac1,fac2,fac3,fac4,fac6, * fac7,fac9, #if defined bio_12_A & defined bioFasham * fac10,fac12, #endif * irrad,rl,aqph,phic,sqfscl FLOAT * irrad0(imt),irrad1(imt),irrad0s(imt,lm),irrad1s(imt,lm), * aph(lm) cjad need to move this stuff to input card or to blkdat.F FLOAT * chlcint,chlcslop c15aug parameter(chlcint=1.312,chlcslop=-0.357) parameter(chlcint=0.1335,chlcslop=0.006979/c100) c c======================================================================= c Begin executable code. c======================================================================= c c c----------------------------------------------------------------------- c Code common to all models: initialization of conversion factors, etc c----------------------------------------------------------------------- c sqfscl=c1e6/(rho0*cp) cjad - note: code for c:chl kinetics is only included for the c exponential model. Needs to be added to the others. c c----------------------------------------------------------------------- c Exponential model c----------------------------------------------------------------------- c if(ipmod.eq.ipmexp) then do 10 k=1,km do 10 i=1,imt 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)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac3=tb(i,k,ino3) fac4=tb(i,k,ichl) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif if (k.eq.1) irrad0(i)=srf(i)*parfrac*ht2qnta/sqfscl* * c1e6/avogadro fac1=(atth2o+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) if(c2chl.gt.c0) then rl=photorm*(c1-fac1)*fac2*c2chl else if((fac7+fac9).gt.c0) then rl=photorm*(c1-fac1)*fac2*fac4/(fac7+fac9)* * c2n/c12 else rl=c0 endif endif #if defined bioDuse & defined dusDroop if((fac7+fac9).gt.c0) then no3pr(i,k)=rl*fac7* * (1-nquota*fac4/(fac7+fac9)) nh4pr(i,k)=rl*fac9* * (1-nquota*fac4/(fac7+fac9)) else no3pr(i,k)=rl*fac7 nh4pr(i,k)=rl*fac9 endif no3pr(i,k)=max(c0,no3pr(i,k)) nh4pr(i,k)=max(c0,nh4pr(i,k)) #else no3pr(i,k)=(fac7+fac9)*rl*fac3*exp(-no3inh*fac6)/ * (hsno3+fac3) nh4pr(i,k)=(fac7+fac9)*rl*fac6/(hsnh4+fac6) #endif if(c2chl.lt.c0) then if(irrad.gt.c1) then c15Aug chlcinf(i,k)=chlcint+chlcslop*log10(irrad) chlcinf(i,k)=chlcint+chlcslop*tdepth(i,k,jrs) else cjad this line designed to prevent acclimation at night chlcinf(i,k) = c0 c if((fac7+fac9).gt.c0) then c chlcinf(i,k)=fac4/(fac7+fac9)*c2n c else c chlcinf(i,k)=c0 c endif endif endif irrad0(i)=irrad1(i) 10 continue c c----------------------------------------------------------------------- c Hyperbolic tangent model c----------------------------------------------------------------------- c elseif(ipmod.eq.ipmtnh) then do 20 k=1,km do 20 i=1,imt 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)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac3=tb(i,k,ino3) fac4=tb(i,k,ichl) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif if (k.eq.1) irrad0(i)=srf(i)*parfrac*ht2qnta/sqfscl* * c1e6/avogadro fac1=(atth2o+attphy*fac4)*dzqz(i,k,0) irrad1(i)=irrad0(i)*exp(-fac1) irrad=(irrad0(i)-irrad1(i))/fac1 rl=photorm*tanh(photor0*irrad/photorm)*c2chl #if defined bioDuse & defined dusDroop if((fac7+fac9).gt.c0) then no3pr(i,k)=rl*fac7* * (1-nquota*fac4/(fac7+fac9)) nh4pr(i,k)=rl*fac9* * (1-nquota*fac4/(fac7+fac9)) else no3pr(i,k)=rl*fac7 nh4pr(i,k)=rl*fac9 endif no3pr(i,k)=max(c0,no3pr(i,k)) nh4pr(i,k)=max(c0,nh4pr(i,k)) #else no3pr(i,k)=(fac7+fac9)*rl*fac3*exp(-no3inh*fac6)/ * (hsno3+fac3) nh4pr(i,k)=(fac7+fac9)*rl*fac6/(hsnh4+fac6) #endif irrad0(i)=irrad1(i) 20 continue c c elseif(ipmod.eq.ipmkm) then c c above is for Kiefer and Mitchell model not used yet c c c----------------------------------------------------------------------- c Smith hyperbolic model c----------------------------------------------------------------------- c elseif(ipmod.eq.ipmsmi) then do 40 k=1,km do 40 i=1,imt 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)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac3=tb(i,k,ino3) fac4=tb(i,k,ichl) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif if (k.eq.1) irrad0(i)=srf(i)*parfrac*ht2qnta/sqfscl* * c1e6/avogadro fac1=(atth2o+attphy*fac4)*dzqz(i,k,0) irrad1(i)=irrad0(i)*exp(-fac1) irrad=(irrad0(i)-irrad1(i))/fac1 rl=photorm*photor0*irrad/sqrt(photorm**c2+photor0**c2* * irrad**2)*c2chl #if defined bioDuse & defined dusDroop if((fac7+fac9).gt.c0) then no3pr(i,k)=rl*fac7* * (1-nquota*fac4/(fac7+fac9)) nh4pr(i,k)=rl*fac9* * (1-nquota*fac4/(fac7+fac9)) else no3pr(i,k)=rl*fac7 nh4pr(i,k)=rl*fac9 endif no3pr(i,k)=max(c0,no3pr(i,k)) nh4pr(i,k)=max(c0,nh4pr(i,k)) #else no3pr(i,k)=(fac7+fac9)*rl*fac3*exp(-no3inh*fac6)/ * (hsno3+fac3) nh4pr(i,k)=(fac7+fac9)*rl*fac6/(hsnh4+fac6) #endif irrad0(i)=irrad1(i) 40 continue c c----------------------------------------------------------------------- c Bidigare et al. (1987) Spectral, phi-based model c----------------------------------------------------------------------- c elseif(ipmod.eq.ipmbid) then do 50 k=1,km do 50 i=1,imt 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)) fac7=max(c0,tb(i,k,iqn3)) fac9=max(c0,tb(i,k,iqn4)) else fac3=tb(i,k,ino3) fac4=tb(i,k,ichl) fac6=tb(i,k,inh4) fac7=tb(i,k,iqn3) fac9=tb(i,k,iqn4) endif aqph=0.0 do 54 l=1,lm cjad currently using a semi-empirical model below cjad should change to allow for real data, when available if (k.eq.1) irrad0s(i,l)=srf(i)*wsj(l)*parfrac* * ht2qnta/sqfscl/ * avogadro fac1=(atth2os(l)+chil(l)*(fac4**elam(l)))*c1em2* * dzqz(i,k,0) irrad1s(i,l)=irrad0s(i,l)*exp(-fac1) irrad=(irrad0s(i,l)-irrad1s(i,l))/fac1 aph(l)=(achla(l)+cha2chb*achlb(l)+cha2chc*achlc(l)+ * cha2psc*apsc(l)+cha2ppc*appc(l)) cjad actual concentrations of accessory pigments can be used above cjad if available - otherwise use assumed ratio aqph=aqph+aph(l)*irrad 54 continue aqph=aqph*c12000 phic=phim*((photorm/phim)/(photorm/phim+aqph)) rl=phic*aqph*c2chl #if defined bioDuse & defined dusDroop if((fac7+fac9).gt.c0) then no3pr(i,k)=rl*fac7* * (1-nquota*fac4/(fac7+fac9)) nh4pr(i,k)=rl*fac9* * (1-nquota*fac4/(fac7+fac9)) else no3pr(i,k)=rl*fac7 nh4pr(i,k)=rl*fac9 endif no3pr(i,k)=max(c0,no3pr(i,k)) nh4pr(i,k)=max(c0,nh4pr(i,k)) #else no3pr(i,k)=(fac7+fac9)*rl*fac3*exp(-no3inh*fac6)/ * (hsno3+fac3) nh4pr(i,k)=(fac7+fac9)*rl*fac6/(hsnh4+fac6) #endif do 56 l=1,lm irrad0s(i,l)=irrad1s(i,l) 56 continue 50 continue c c additional productivity models can be inserted at this point... c else call exitus('PRIPROD') endif #if defined bio_12_A & defined bioFasham c c calculate bacterial uptake and/or additional fluxes c do 60 k=1,km do 60 i=1,imt if(biopos.eq.2) then fac6=max(c0,tb(i,k,inh4)) fac10=max(c0,tb(i,k,ibac)) fac12=max(c0,tb(i,k,idon)) else fac6=tb(i,k,inh4) fac10=tb(i,k,ibac) fac12=tb(i,k,idon) endif fac1=min(fac6,bacfrac*fac12) fac2=bacgrr*fac10/(bacgrk+fac1+fac12) bgrdon(i,k)=fac2*fac12 bgrnh4(i,k)=fac2*fac1 60 continue c c calculate zooplankton grazing rates c do 70 k=1,km do 70 i=1,imt if(biopos.eq.2) then fac5=max(c0,tb(i,k,izoo)) fac7=max(c0,tb(i,k,iqn3)) fac8=max(c0,tb(i,k,idet)) fac9=max(c0,tb(i,k,iqn4)) fac10=max(c0,tb(i,k,ibac)) else fac5=tb(i,k,izoo) fac7=tb(i,k,iqn3) fac8=tb(i,k,idet) fac9=tb(i,k,iqn4) fac10=tb(i,k,ibac) endif fac79i=civlevP*(fac7+fac9) fac8i=civlevD*fac8 fac10i=civlevB*fac10 if((fac79i).gt.0) then zgrphy(i,k)=fac5*grazrm * *((c1-exp(-fac79i))*(exp(-fac10i) * *exp(-fac8i)+(c1-exp(-fac10i)) * *exp(-fac8i)*fac79i/(fac79i+fac10i) * +(c1-exp(-fac8i))*exp(-fac10i) * *fac79i/(fac79i+fac8i)+(c1-exp(-fac81)) * *(c1-exp(-fac10i))*fac79i * /(fac79i+fac8i+fac10i))) else zgrphy(i,k)=c0 endif if(fac10i.gt.0) then zgrbac(i,k)=fac5*grazrm * *((c1-exp(-fac10i))*(exp(-fac79i) * *exp(-fac8i)+(c1-exp(-fac79i)) * *exp(-fac8i)*fac10i/(fac10i+fac79i) * +(c1-exp(-fac8i))*exp(-fac79i) * *fac10i/(fac10i+fac8i)+(c1-exp(-fac8i)) * *(c1-exp(-fac79i))*fac10i * /(fac10i+fac8i+fac79i))) else zgrbac(i,k)=c0 endif if(fac8i.gt.0) then zgrdet(i,k)=fac5*grazrm * *((c1-exp(-fac8i))*(exp(-fac79i) * *exp(-fac10i)+(c1-exp(-fac79i)) * *exp(-fac10i)*fac8i/(fac8i+fac79i) * +(c1-exp(-fac10i))*exp(-fac79i) * *fac8i/(fac8i+fac10i)+(c1-exp(-fac10i)) * *(c1-exp(-fac79i))*fac8i * /(fac8i+fac10i+fac79i))) else zgrdet(i,k)=c0 endif 70 continue #endif return end