program driver c c This program runs a 1-D offline version of the Hood et al. (2000) c N-DON-P-T-Z-D ecosystem model. It has been specifically set up to be c driven with 1-D output (surface irradiance, mixed layer depth and c vertical velocities) from the McCreary et al. (2000) Arabian Sea 3-D c circulation model. The McCreary et al. model is a layer model, c the output from which has been converted to level (fixed vertical grid) c coordinates using the program l2l.f. c c This program incorporates subroutines for calculating the diffusive fluxes c of all biological quantities, downward sinking fluxes of detritus, and c vertical advective fluxes derived using vertical velocities from c McCreary's physical model. c c The boundaries are closed at the surface, and open at the bottom. c c To create an executable, compile with: c c 3) common_blocks.h (declares common variables) c 4) dimensions.h (declares array indexes/sizes) c c Last edited by R. R. Hood on August, 25 2000. c c Include common parameters and variables: c include 'dimensions.h' include 'common_blocks.h' c c Declare locally utilized real, integer and character variables c integer isteps,opt real wsink(1000),rkzn(1000),dout,endtime real badv(nbio),bflux(nbio),bsinkd,tmass,ttmass real ctime,cdays,dml,qi,sout,ajm1,djm1,bbc,wadv(1000) real bioi(nbio),hadvb(1000,nbio) logical first c c This logical if is used to specify euler step in first Runge-Kutta step c first = .true. c c Specify the domain, time/space resolution, no. of vertical points, c and output interval and run duration... c dz = 1.0 !vertical resolution (meters) delt1 = 1800.0 !time step for first step nz = 400 !number of points in the domain dout = 1.0 !output interval (days) sout = 1.0 !output interval to screen isteps = 500000 !max. no. of time steps c c Specify here whether or not to use the optional advective fields: c opt = 0 gives no horizontal advection c opt = 1 gives horizontal advection of all bio quantities c opt = 2 gives horizontal advection of DIN only c opt = 1 c c Specify here whether or not to use the initial conditions from c McCreary et al. (2000): c icond = 0 do not use the intitial conditions c icond = 1 use the initial conditions c icond = 1 c c Recall the order of the variables and integer designations c in the biological model... c c N = bio(k,in) where in = 1 c DON = bio(k,idon) where idon = 2 c P = bio(k,iph) where iph = 3 c T = bio(k,it) where it = 4 c Z = bio(k,ih) where ih = 5 c D = bio(k,id) where id = 6 c c c Sum mass in the system and print initial total c do j=1,nz tmass = bio(j,in)+bio(j,idon)+bio(j,iph) . +bio(j,it)+bio(j,ih)+bio(j,id) stmass = stmass + tmass*dz enddo print*,'total initial mass (mmoles/m^2) =',stmass c c c Set vertical diffusion profile. For now, its vertically uniform, c with the value set in biodat.f, but the diffusion coefficient can be made c to vary vertically if desired. c do k = 1,nz rkzn(k) = rnkz end do c c c Open the output files so you can write to it as you go: c open(unit=8,status='new',file='bio1.dat') !for vertical profiles open(unit=10,status='new',file='bio2.dat') !for vertical integrals c c c Open the forcing file (1-D output from the McCreary et al. model), c which includes the "required fields" time (cdays), surface radiation (qi), c depth of the mixed layer (dml), and the vertical velocity profile wadv(n) c open (unit=9,status='old',form='unformatted',file='levelcode.dat') c c If specified, open the optional field which specifies the horizontal c advection of biological quantities: c if ((opt.eq.1).or.(opt.eq.2)) then open (unit=11,status='old',form='unformatted', . file='bioadv_level.dat') endif c c If specified, open the optional field which specifies the initial c conditions: c c Set the initial values of the biological quantities if (icond.eq.0) then do k = 1,nz bio(k,in) = 10.0 bio(k,idon) = 0.1 bio(k,iph) = 0.1 bio(k,it) = 0.0 bio(k,ih) = 0.1 bio(k,id) = 0.1 enddo elseif(icond.eq.1) then open (unit=12,status='old',form='unformatted', . file='initbio_level.dat') read(12)(bio(n,in),n=1,nz+1) read(12)(bio(n,iph),n=1,nz+1) read(12)(bio(n,ih),n=1,nz+1) read(12)(bio(n,id),n=1,nz+1) do n=1,nz+1 bio(n,idon) = 0.15 !Init. don at mid-range value bio(n,it) = 0.0 !Don't initialize tricho enddo else print *,'Invalid value for icond, program stopped' stop endif c c**** start main loop over time (picking up forcing data as we go)**** c c Set delt1 = 1800 for the first step. After that use the time step c as specified by the forcing file c do 11 jmain=1,isteps istep = jmain cdaysl = cdays read(9,end=12) cdays,qi,dml read(9)(wadv(n),n=1,nz+1) do n=1,nz+1 wadv(n)=-wadv(n) enddo if ((opt.eq.1).or.(opt.eq.2)) then read(11)(hadvb(n,in),n=1,nz+1) read(11)(hadvb(n,iph),n=1,nz+1) read(11)(hadvb(n,ih),n=1,nz+1) read(11)(hadvb(n,id),n=1,nz+1) do n=1,nz+1 hadvb(n,idon) = hadvb(n,iph) !assume don grads similar to phyto hadvb(n,it) = 0 !assume tricho is horizont. uniform enddo endif if (istep.gt.1) delt1 = (cdays - cdaysl)*86400.0 qi = qi/0.4 !convert from PAR to total solar ndml = int(dml/dz) !calculate no. of cells in mixed layer c c First calculate changes in biological quantities due to biological c processes c call bio4nl(qi,first,istep) c c Now advect and diffuse phytoplankton c ajm1 = bio(1,iph) !set surface boundary condition bbc = bio(nz,iph) - (bio(nz-1,iph) - bio(nz,iph)) !set bot. bound. cond. if (bbc.lt.0.0) bbc = 0.0 call advectv (wadv,delt1,dz,nz,bbc,ajm1,bio(1,iph),badv(iph)) call diffusv (rkzn,delt1,dz,nz,bbc,bio(1,iph),bflux(iph)) if (opt.eq.1) call advecth (delt1,nz,bio(1,iph),hadvb(1,iph)) c c Now advect and diffuse heterotrophs c ajm1 = bio(1,ih) bbc = bio(nz,ih) - (bio(nz-1,ih) - bio(nz,ih)) if (bbc.lt.0.0) bbc = 0.0 call advectv (wadv,delt1,dz,nz,bbc,ajm1,bio(1,ih),badv(ih)) call diffusv (rkzn,delt1,dz,nz,bbc,bio(1,ih),bflux(ih)) if (opt.eq.1) call advecth (delt1,nz,bio(1,ih),hadvb(1,ih)) c c Now advect and diffuse Trichodesmium c ajm1 = bio(1,it) bbc = bio(nz,it) - (bio(nz-1,it) - bio(nz,it)) if (bbc.lt.0.0) bbc = 0.0 call advectv (wadv,delt1,dz,nz,bbc,ajm1,bio(1,it),badv(it)) call diffusv (rkzn,delt1,dz,nz,bbc,bio(1,it),bflux(it)) if (opt.eq.1) call advecth (delt1,nz,bio(1,it),hadvb(1,it)) c c Now advect and diffuse DON c ajm1 = bio(1,idon) bbc = bio(nz,idon) - (bio(nz-1,idon) - bio(nz,idon)) if (bbc.lt.0.0) bbc = 0.0 call advectv (wadv,delt1,dz,nz,bbc,ajm1,bio(1,idon),badv(idon)) call diffusv (rkzn,delt1,dz,nz,bbc,bio(1,idon),bflux(idon)) if (opt.eq.1) call advecth (delt1,nz,bio(1,idon),hadvb(1,idon)) c c Now advect and diffuse DIN (open bottom boundary) c ajm1 = bio(1,in) !surface boundary cond. for advection bbc = bnutz !bottom boundary cond. for adv. and diff. call advectv (wadv,delt1,dz,nz,bbc,ajm1,bio(1,in),badv(in)) call diffusv (rkzn,delt1,dz,nz,bbc,bio(1,in),bflux(in)) if ((opt.eq.1).or.(opt.eq.2)) call . advecth (delt1,nz,bio(1,in),hadvb(1,in)) c c Now advect, diffuse and sink detritus c ajm1 = bio(1,id) ! + (bio(1,id) - bio(2,id)) bbc = bio(nz,id) - (bio(nz-1,id) - bio(nz,id)) if (bbc.lt.0.0) bbc = 0.0 call advectv (wadv,delt1,dz,nz,bbc,ajm1,bio(1,id),badv(id)) call diffusv (rkzn,delt1,dz,nz,bbc,bio(1,id),bflux(id)) call sink (ws,delt1,dz,nz,bio(1,id),bsinkd) if (opt.eq.1) call advecth (delt1,nz,bio(1,id),hadvb(1,id)) c c Now uniformly mix everything down to dml c do inb = 1,nbio call mix1m(bio(1,inb),ndml) enddo c c Sum the bottom fluxes over time c stbflx = stbflx+bsinkd c Write to the output file if the time matches the output time c if(mod(cdays,dout).lt..00002) then do inb = 1,nbio do j = 1,nz bioi(inb) = bioi(inb) + bio(J,inb) enddo enddo write(10,101) cdays,bioi(in),bioi(idon), . bioi(iph),bioi(it),bioi(ih),bioi(id) ttmass = stmass - stbflx do k = 1,nz depth = real(k*dz) write(8,100) cdays,depth,bio(k,in),bio(k,idon), . bio(k,iph),bio(k,it),bio(k,ih),bio(k,id) enddo endif if(jmain.eq.1) then print*,'days ','dml ','DIN ','DON ','P ','T ','H ','D ' endif if(mod(cdays,sout).lt..00002) then ttmass = stmass + stbflx write(6,100) cdays,dml,bioi(in),bioi(idon),bioi(iph), . bioi(it),bioi(ih),bioi(id) endif do inb = 1,nbio bioi(inb) = 0.0 enddo 100 format(f7.2,1x,f7.2,1x,6(f9.3,x)) 101 format(f7.2,1x,1x,6(f9.3,x)) c c c 11 continue 12 continue c c Close open files c close (unit=8) close (unit=9) end c c subroutine diffusv(rkz,dt,dz,nz,bnz,a,bflux) real a(1),ar(799),br(799),cr(799),r(799),rkz(1),dt,dz,bnz, . bflux,dconst,anz,rkznzp1 C C This subroutine applies a Crank-Nicolson vertically variable diffusion C operation to the array A where the gradient is set to zero at the C upper boundary, and, depending upon the value of BNZ, applies either a C closed or open bottom boundary condition. By setting BNZ equal to the C value at NZ (just above) the bottom boundary can be closed (i.e., set C the gradient across the bottom to zero). Otherwise the bottom boundary C is open with the flux determined by the gradient between the C concentration below (BNZ) the bottom boundary and the dynamically C determined interior concentration just above the bottom boundary. C dconst = dt/(2.0*(dz**2)) nzm = nz - 1 c c do 1 j=2,nzm rkzm = (rkz(j-1) + rkz(j))/2.0 rkzp = (rkz(j+1) + rkz(j))/2.0 r(j) = dconst*rkzm*a(j-1) + (1.0-(dconst*(rkzp+rkzm)))*a(j) . + dconst*rkzp*a(j+1) 1 continue c rkzp = (rkz(1) + rkz(2))/2.0 r(1)=(1.0-(dconst*rkzp))*a(1) + dconst*rkzp*a(2) !For top boundary rkzm = (rkz(nz) + rkz(nz-1))/2.0 rkznzp1 = rkz(nz) + (rkz(nz) - rkz(nz-1)) !Extrapolate KZ to NZ+1 rkzp = (rkz(nz) + rkznzp1)/2.0 !Calc KZ at NZ+1/2 r(nz) = dconst*rkzm*a(nz-1) . + (1.0-(dconst*(rkzp+rkzm)))*a(nz) . + 2*dconst*rkzp*bnz !iFor bottom boundary C C Below sets coefficients for tridiagonal matrix which are different at C the top and bottom boundaries. C do 2 j=2,nzm rkzm = (rkz(j-1) + rkz(j))/2.0 rkzp = (rkz(j+1) + rkz(j))/2.0 ar(j) = -dconst*rkzm br(j) = 1.0 + (dconst*(rkzm + rkzp)) cr(j) = -dconst*rkzp 2 continue rkzp = (rkz(1) + rkz(2))/2.0 br(1) = 1.0 + (dconst*rkzp) cr(1) = -dconst*rkzp rkzm = (rkz(nz) + rkz(nz-1))/2.0 rkznzp1 = rkz(nz) + (rkz(nz) - rkz(nz-1)) !Extrapolate KZ to NZ+1 rkzp = (rkz(nz) + rkznzp1)/2.0 !Calc. KZ at NZ+1/2 ar(nz)= -dconst*rkzm br(nz)= 1.0 + (dconst*(rkzm+rkzp)) c anz = a(nz) call tridag(ar,br,cr,r,a,nz) C C Calculate the flux at the bottom: C bflux = dz*(dconst*rkzp)*((2*bnz)-a(nz)-anz) !FOR FLUX AT BOTTOM c return end C C C subroutine sink (w,dt,dz,nz,a,fbnet) real a(799),astar(799),w,dt,dz,wbarp,wbarm, . fbelow,fabove,e,fbstar,fbnet C C C This subroutine applies a modified second order flux-corrected C transport scheme (Smolarkiewics, 1983) to calculate the effect of detrital C sinking on the detritus array A. It is modified in that downward sinking C velocities, which are vertically invariant are assumed, which simplifies C the calculations. The top boundary is closed (no sinking input from above) C and the bottom boundary is open. (Note, this scheme can be greatly C simplified by dropping the correction for numerical diffusion. However, C without the correction, the scheme doesn't conserve as well). C C E = 10.0**(-15.0) C C C Calculate ASTAR using a simple upstream scheme at the top boundary: C fabove = 0.0 !No flux from above the top grid point fbelow = w*a(1)*dt/dz !Flux out of the the top grid point astar(1) = a(1) - (fbelow - fabove) C C Extrapolate one point below A(NZ) because you will need it to C calculate ASTAR(NZ+1) for the anti-diffusion velocity at the C bottom boundary: C a(nz+1) = a(nz) - (a(nz-1) - a(nz)) if(a(nz+1).lt.0.0) a(nz+1) = 0.0 C C Calculate ASTAR values using a standard upstream scheme through NZ+1: C do 1 j=2,nz+1 fabove = w*a(j-1)*dt/dz fbelow = w*a(j)*dt/dz astar(j) = a(j) - (fbelow - fabove) 1 continue fbstar = w*a(nz)*dt/dz C C Calculate anti-diffusion velocity for top boundary: C wbarp = ((w*dz - (dt*(w**2.0)))*(astar(2) - astar(1))) . /((astar(1) + astar(2) + e)*dz) C C Calculate correction fluxes using ASTAR values for the top boundary: C (note that the corrction flux can be in either direcion). C fabove = 0.0 fbelow = (((wbarp + abs(wbarp))*astar(1)) + . ((wbarp - abs(wbarp))*astar(2)))*(dt/(dz*2.)) c c fbstar = w*a(1)*(dt/dz) c print*,fbstar,fbelow C C And the new A value at the top boundary: C a(1) = astar(1) - (fbelow - fabove) C C Now do the rest of the column: C do 2 j=2,nz C C Calculate anti-diffusion velocities: C wbarp = ((w*dz - (dt*(w**2.0)))*(astar(j+1) - astar(j))) . /((astar(j) + astar(j+1) + e)*dz) wbarm = ((w*dz - (dt*(w**2.0)))*(astar(j) - astar(j-1))) . /((astar(j) + astar(j-1) + e)*dz) C C Calculate the correction fluxes using ASTAR values: C fabove = (((wbarm + abs(wbarm))*astar(j-1)) + . ((wbarm - abs(wbarm))*astar(j)))*(dt/(dz*2.)) fbelow = (((wbarp + abs(wbarp))*astar(j)) + . ((wbarp - abs(wbarp))*astar(j+1)))*(dt/(dz*2.)) C C Calculate the corrected A values: C a(j) = astar(j) - (fbelow - fabove) 2 continue C C Now calculate the flux out of the bottom: C wbarp = ((w*dz - (dt*(w**2.0)))*(astar(nz+1) . - astar(nz)))/((astar(nz) + astar(nz+1) + e)*dz) fbelow = (((wbarp + abs(wbarp))*astar(nz)) + . ((wbarp - abs(wbarp))*astar(nz+1)))*(dt/(dz*2.)) fbnet = (fbstar + fbelow)*dz c c c return end c c c subroutine advecth(dt,nz,a,b) real a(nz),b(nz),da,dt C C C This subroutine applies a simple, centered, finite-difference, gradient C form advection scheme to calculate the effect of horizontal advection on c the biological quantities using the product of the biomass gradient times c the horizontal advection provided from the McCreary et al. model. C C Calculate the horizontal advective changes in the biological quantities C over the entire domain: C do j=1,nz da = b(j)*dt a(j) = a(j) + da end do C C C return end C C C subroutine advectv(w,dt,dz,nz,anzp1,ajm1,a,bsink) real a(nz),da(1000),w(nz),dt,dz,nzp1,anzp1, . wnzp1,bsink,ajm1 C C C This subroutine applies a simple, centered, finite-difference, gradient C form advection scheme to calculate the effect of advection on the biological C quantities. The gradient form, as opposed to the flux form, is used here C to prevent accumulation or loss of mass convergence or divergence in 1-D. C C Set the vertical velocity just below the bottom boundary to C be the same as the vertical velocity at the bottom boundary C wnzp1 = w(nz) C C Calculate changes in the biological quantities due to advection or C sinking in the surface cell. For advection ajm1 is set to a(1) C (and w(1) should also be zero) which eliminates the first term in the C equation. C c ajm1 = a(1) da(1) = 0.5*((w(1)*(ajm1 - a(1))) . + (w(2)*(a(1) - a(2))))*dt/dz a(1) = a(1) + da(1) C C Calculate sinking or advective changes in the biological quantities C in the interior: C do j=2,nz-1 da(j) = 0.5*((w(j)*(a(j-1) - a(j))) . + (w(j+1)*(a(j) - a(j+1))))*dt/dz a(j) = a(j) + da(j) end do C C Calculate advective or sinking changes in the biological quantities in the C bottom cell. For both sinking and advection the bottom boundary condition C is set by extrapolating the gradient just above the bottom boundary to C that just below. C da(nz) = 0.5*((w(nz)* (a(nz-1) - a(nz))) . + (wnzp1*(a(nz) - anzp1)))*dt/dz a(nz) = a(nz) + da(nz) C C Calculate the net flux into the top and out of the bottom (note this C calculation only holds when you fold over both the bottom AND top C boundary conditions, and have vertically uniform sinking rate) C bsink = (w(1)*a(1) - wnzp1*a(nz))*dt C C C return end C C subroutine tridag(a,b,c,r,u,n) parameter (nmax=1000) real gam(nmax),a(n),b(n),c(n),r(n),u(n),bet C C Solves for a vector U of length N the tridiagonal linear set C given by the Crank-Nicolson scheme. A, B, C and R are input C vectors and are not modified. This subroutine is from Press et al., c numerical recipes, p. 40. C if(b(1).eq.0.)pause bet=b(1) u(1)=r(1)/bet do 1 j=2,n gam(j) = c(j-1)/bet bet=b(j)-a(j)*gam(j) if(bet.eq.0.)pause u(j)=(r(j)-a(j)*u(j-1))/bet 1 continue do 2 j=n-1,1,-1 u(j)=u(j)-gam(j+1)*u(j+1) 2 continue return end c c c subroutine mix1m(d,k) c C This subroutine completely mixes the array D down to level K. C real d(1),ds ds = 0. do 1 j=1,k ds = ds + d(j) 1 continue ds = ds/REAL(k) do 2 j=1,k d(j) = ds 2 continue return end C C block data biodat C include 'dimensions.h' include 'common_blocks.h' c-- in = din index c-- idon = don index c-- iph = phyto index c-- it = tricho index c-- ih = zoop index c-- id = detritus index data in,idon,iph,it,ih,id/1,2,3,4,5,6/ C C Note that all of the biological model parameters are in MKS C C Set the bio. model paramters: C DATA GEP/0.20/ !H growth efficiency on P . ,AEP/0.70/ !H assim. efficiency on P . ,GED/0.20/ !H growth efficiency on D . ,AED/0.70/ !H assim. efficiency on D . ,GEH/0.20/ !H growth efficiency on H . ,AEH/0.70/ !H assim. efficiency on H . ,GET/0.20/ !H growth efficiency on T . ,AET/0.70/ !H assim. efficiency on T . ,GEN/0.20/ !H growth efficiency on DON . ,AEN/1.0/ !H assim. efficiency on DON . ,CMAX/7.407E-5/ ! = 6.4/86400., max. Het. consumption rate . ,HKS/.80/ !H half sat. const. . ,PHIP/0.25/ !H preference for P . ,PHID/0.25/ !H preference for D . ,PHIH/0.25/ !H preference for Z . ,PHIT/0.0/ !H preference for T . ,PHIN/0.25/ !H preference for DON . ,UMAXP/3.727E-5/ != 3.22/86400., P max. growth rate . ,UMAXT/0/ !/1.967E-6/ != 0.17/86400., T max. growth rate . ,PKS/.5/ !P half sat. const. . ,TKS/.5/ !T half sat. const. . ,SP/1.157E-7/ != 0.01/86400., P Senescence rate . ,ST/2.89E-7/ != 0.025/86400., T Senescence rate . ,ST/3.13e-7/ != 0.027/86400., T Senescence rate . ,BETA/.25/ !Senesence fraction to DETR . ,ALPHA/.70/ !Partitioning of pripro. to T/P pools . ,RLTKP/40.0/ !Sat. const. for P WATTS/M^2 . ,RLTKT/80.0/ !Sat. const. for T WATTS/M^2 . ,GAMMA/0.75/ !Partitioning of excretion to DIN . ,RLTB/400.0/ !P photoinh. parameter C C Set light model paramters C data rs/0.6/ !Fraction of QI that is longwave . , wkr/1.6666/,wkb/0.03/ !Red and blue Ks for water C C Set sinking rate of detritus C data ws/1.16e-4/ !10m/day, c c Set the vertical diffusivity c data rnkz/1.0e-5/ !1 cm^2/sec c c Set bottom boundary DIN concentration c data bnutz/20/ !20 uM, after McCreary et al. (2000) c end subroutine bio4nl(qi,first,istep) c c c This subroutine initializes and time-steps a two-species N-DON-P-T-Z-D c ecosystem model, where DON is labile dissolved inorganic nitrogen and T c is the diazotrophic cyanobacterium, Trichodesmium. This model is derived c from bio2nl.f that was used in Hood et al. (2000), Deep-Sea Research, 2nd c Special Issue on the U.S. JGOFS Ocean Time Series Stations. c c This version was specifically adapted for making offline calculations at c fixed levels using output (surface radiation, mixed layer depth and vertical c velocities) adapted from the McCreary et al. (1996; 2000) 3-D, reduced c gravity (layer) circulation model of the Arabian Sea. And it includes a c simple 2-wavelength optical model (calclight.f) for calculating the c subsurface light field. c c Model parameter values are specified in biodat.f c c Last edited by R. Hood on July, 24 2000. c c Include common parameters and variables: c include 'dimensions.h' include 'common_blocks.h' c c Declare locally utilized real and integer variables, except c those specific to the subroutine derivs: c real yt(nbio),dyt(nbio+1),dym(nbio+1),qi real tempbio(nbio),dydtemp(nbio+1) integer istep c real hh,h6 logical first c c Split the light into long and shortwave components: c rior = abs(qi*rs) !Red fraction in watts/m^2 riob = abs(qi*(1-rs)) !Blue fraction in watts/m^2 c c *****Main loop over depth****** c do 9 k = 1,nz c c First calculate the depths of the box boundaries: c zl2 = k*dz !z of bottom of the box zl1 = zl2-dz !z of top of the box c c Calculate average light in the box under consideration: c (subroutine calculates new values for rior and riob for the c top of the next box when it is done) c call calclight(bio(k,iph),bio(k,it)) c c The following is a slightly modified version of RK4 from Numerical c Recipes where the first time step is Euler. c c Note that dydt also stores one additional term dydt(k,nbio+1) which is c the gross rate of N2-fixation. c c There is also a third array (dydti(x,nbio) which stores the individual c terms in each of the model equations. The value of x depends upon c the number of equations and the number of terms in each. For now, a c value of 50 will do. c c derivs(bio,dydt) is the subroutine that returns derivatives c at dydt. c do inb = 1,nbio tempbio(inb) = bio(k,inb) enddo hh=delt1*0.5 h6=delt1/6. if (first) then call derivs(tempbio,dydtemp) do 10 inb=1,nbio bio(k,inb)=bio(k,inb)+delt1*dydtemp(inb) !normal Euler step 10 continue c if(k.eq.1) print*,istep,bio(1,5),bio(1,6) else call derivs(tempbio,dydtemp) do 11 inb=1,nbio yt(inb)=bio(k,inb)+hh*dydtemp(inb) !bio mags. after half step 11 continue call derivs(yt,dyt) !gives derivs (dyt) at halfstep do 12 inb=1,nbio yt(inb)=bio(k,inb)+hh*dyt(inb) !mags. after half step w/dyt 12 continue call derivs(yt,dym) !gives deriv (dym) at half step do 13 inb=1,nbio yt(inb)=bio(k,inb)+delt1*dym(inb) !mags at full step w/dym dym(inb)=dyt(inb)+dym(inb) !sums two half step derivatives 13 continue call derivs(yt,dyt) !gives derivs at full step do 14 inb=1,nbio bio(k,inb)=bio(k,inb)+h6*(dydtemp(inb)+ . dyt(inb)+2.*dym(inb)) c c Print warnings if we get negative or excessively large bomasses in any c of the model compartments. c if (bio(k,inb) .lt. 0) then print *,'neg bio istep, k,inb',istep,k,inb,bio(k,inb) endif if (abs(bio(k,inb)) .gt. 50) then print *,'big! bio ',istep,k,inb,bio(k,inb) endif 14 continue endif c c PRINT*,'made it past DERIVS' c 9 continue if (first) first=.false. return end c c c Come here to calculate the derivatives: c subroutine derivs(y,dydtt) c c This subroutine calculates time change in each biological compartment. c c Note that arrays y(nbio) and dydtt(nbio) are dummy arrays which c are passed different values of bio and dydt depending on the step size c in question. c c Also note that the values of dydti (the individual rate terms) will be c specified, ultimately, by the last call to derives. c c Include common bio arrays. c include 'dimensions.h' !defines array sizes include 'common_blocks.h' !common variables and arrays and type c c declare variables specific to this subroutine: c real y(nbio),dydtt(nbio+1),thta,hp,hd,hh,ht,hn,up,ut c c Recall that for the current version of the model: c c N = bio(k,1) is y(in) c DON = bio(k,2) is y(idon) c P = bio(k,3) is y(iph) c T = bio(k,4) is y(it) c Z = bio(k,5) is y(ih) c D = bio(k,6) is y(id) c c Here are the grazing formulations for P, D, T and Z. c All hyperbolic saturation functions: c thta=(phip*y(iph)+phid*y(id)+phih*y(ih)+phit*y(it) . +phin*y(idon)+hks) C hp = phip/thta hd = phid/thta ht = phit/thta hh = phih/thta hn = phin/thta C C Here is the phytoplankton DIN-uptake equation: C up = umaxp * . ((1-exp(-rlt/rltkp))*exp(-rlt/rltb)) . * (y(in)/(y(in)+pks)) C C Here is the Trichodesmium DIN-uptake equation: C ut = umaxt * (1-exp(-rlt/rltkt)) . * (y(in)/(y(in)+tks)) C C Here is the Trichodesmium growth equation: C gt = umaxt * (1-exp(-rlt/rltkt)) C C C Here are the individual N terms: C dydti(1,in) = gamma*(aep-gep)*hp*cmax*y(ih)*y(iph) dydti(2,in) = gamma*(aed-ged)*hd*cmax*y(ih)*y(id) dydti(3,in) = gamma*(aet-get)*ht*cmax*y(ih)*y(it) dydti(4,in) = gamma*(aeh-geh)*hh*cmax*y(ih)*y(ih) dydti(5,in) = (aen-gen)*hn*cmax*y(ih)*y(idon) dydti(6,in) = - up * y(iph) dydti(7,in) = - ut * y(it) C C Here is the the net N rate: C dydtt(in) = dydti(1,in)+dydti(2,in)+dydti(3,in) . +dydti(4,in)+dydti(5,in)+dydti(6,in)+dydti(7,in) C C Here are the individual DON terms: C dydti(1,idon) = (1-gamma)*(aep-gep)*hp*cmax*y(ih)*y(iph) dydti(2,idon) = (1-gamma)*(aed-ged)*hd*cmax*y(ih)*y(id) dydti(3,idon) = (1-gamma)*(aet-get)*ht*cmax*y(ih)*y(it) dydti(4,idon) = (1-gamma)*(aeh-geh)*hh*cmax*y(ih)*y(ih) dydti(5,idon) = (0.0)*(aen-gen)*hn*cmax*y(ih)*y(idon) dydti(6,idon) = (1-alpha) * gt * y(it) dydti(7,idon) = (1-alpha) * up * y(iph) dydti(8,idon) = (1-beta) * sp * y(iph) dydti(9,idon) = (1-beta) * st * y(it) dydti(10,idon) = - y(ih) * hn * cmax * y(idon) C C Here is the net DON rate: C dydtt(idon) = dydti(1,idon)+dydti(2,idon)+dydti(3,idon) . +dydti(4,idon)+dydti(5,idon)+dydti(6,idon)+dydti(7,idon) . +dydti(8,idon)+dydti(9,idon)+dydti(10,idon) C C Here are the individual P terms: C dydti(1,iph) = alpha * up * y(iph) dydti(2,iph) = - sp * y(iph) dydti(3,iph) = - hp * y(ih) * cmax * y(iph) C C Here is the net P rate: C dydtt(iph) = dydti(1,iph) + dydti(2,iph) + dydti(3,iph) C C Here are the individual T terms: C dydti(1,it) = alpha * gt * y(it) dydti(2,it) = - st * y(it) dydti(3,it) = - ht * y(ih) * cmax * y(it) C C Here is the net T rate: C dydtt(it) = dydti(1,it) + dydti(2,it) + dydti(3,it) C C Here are the individual H terms: C dydti(1,ih) = gep * hp * y(ih) * cmax * y(iph) !growth on P dydti(2,ih) = get * ht * y(ih) * cmax * y(it) !growth on T dydti(3,ih) = geh * hh * y(ih) * cmax * y(ih) !growth on H dydti(4,ih) = ged * hd * y(ih) * cmax * y(id) !growth on D dydti(5,ih) = gen * hn * y(ih) * cmax * y(idon) !growth on DON dydti(6,ih) = - hh * y(ih) * cmax * y(ih) !self predation C C Here is the net H rate: C dydtt(ih) = dydti(1,ih)+dydti(2,ih)+dydti(3,ih) . +dydti(4,ih)+dydti(5,ih)+dydti(6,ih) C C Here are the individual D terms: C dydti(1,id) = (1 - aep) * hp * y(ih) * cmax * y(iph) dydti(2,id) = (1 - aet) * ht * y(ih) * cmax * y(it) dydti(3,id) = (1 - aeh) * hh * y(ih) * cmax * y(ih) dydti(4,id) = (1 - aen) * hn * y(ih) * cmax * y(idon) dydti(5,id) = (1 - aed) * hd * y(ih) * cmax * y(id) dydti(6,id) = beta * sp * y(iph) dydti(7,id) = beta * st * y(it) dydti(8,id) = - hd * y(ih) * cmax * y(id) C C Here is the net D rate: C dydtt(id) = dydti(1,id)+dydti(2,id)+dydti(3,id) . +dydti(4,id)+dydti(5,id)+dydti(6,id)+dydti(7,id) . +dydti(8,id) C C Calculate the Gross N2 fixation rate C dydtt(nbio+1) = (gt - ut) * y(it) return end C C C subroutine calclight(bioiph,bioit) c c Ths subroutine calculates the average light in the layer/level in c question using a dual-wavelength model which includes the effects c of light absorption by phytoplankton and Trichodesmium c c last modified by R. Hood, June 27, 2000 c c Include common bio arrays. c include 'dimensions.h' !defines array sizes include 'common_blocks.h' !common variables and arrays and type c c declare local variables that are not used elsewhere: c real tp,pkr,pkb,kr,kb,riravg,ribavg real bioiph,bioit c c First calculate the total phytoplankton conc (trichodesmium plus c phytoplankton) and the extinction coefficients due to phytoplankton c absorption of light: c tp = bioiph + bioit pkr = 0.0*(1.0/50.0)*12.0*(106.0/16.0)*tp !Set Kp red pkb = 0.014*(1.0/50.0)*12.0*(106.0/16.0)*tp !Set Kp blue c c pkr is for red and infrared light (set to zero here, assuming red c and infrared absorption by phytoplankton is small compared to water). c pkb is for blue light. Note that it is assumed here that Kc is the c same for the two species of phytoplankton. c c Conversion factors: c 0.014 is Kc (=0.014) (Kirk 1983, p. 273) c 1/50 is the wt/wt chlorophyll-a/carbon ratio (Parsons et al. 1984, p.48) c 106/16 is Redfield atomic carbon/nitrogen ratio (Redfield et al., The c Sea, V. 2, p. 27.) 12 is the molecular weight of carbon c kr = wkr + pkr !Calc. Kt (water plus phyto.) for red. kb = wkb + pkb !Calc. Kt (water plus phyto.) for blue. c c Now calculate the average light for the layer for the red and blue c spectral components. c riravg = - rior*(exp(-kr*zl2) - exp(-kr*zl1))/(kr*(zl2-zl1)) ribavg = - riob*(exp(-kb*zl2) - exp(-kb*zl1))/(kb*(zl2-zl1)) rlt = riravg + ribavg c c Reset long and shortwave radiations for the top of the next level c rior = rior*exp(-kr*(zl2-zl1)) riob = riob*exp(-kb*(zl2-zl1)) return end