ROMS
Loading...
Searching...
No Matches
nemuro.h
Go to the documentation of this file.
1 MODULE biology_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! Nemuro Lower Trophic Level Ecosystem Model. !
11! !
12! This routine computes the biological sources and sinks and adds !
13! then the global biological fields. Currently, the ecosystem has !
14! the following functional compartments: !
15! !
16! iSphy small phytoplankton biomass, nanophytoplankton !
17! iLphy large phytoplankton biomass, diatoms !
18! iSzoo small zooplankton biomass, microzooplankton (ciliates) !
19! iLzoo large zooplankton biomass, mesozooplankton (copepods) !
20! iPzoo predator zooplankton biomass (euphausiids, etc) !
21! iNO3_ nitrate concentration, NO3 !
22! iNH4_ ammonium concentration, NH4 !
23! iPON_ particulate organic nitrogen !
24! iDON_ dissolved organic nitrogen !
25! iSiOH silicate concentration, Si(OH)4 (silicic acid) !
26! iopal particulate organic silica !
27! !
28! Reference: !
29! !
30! Kishi, M. J., et all, 2007: Nemuro - a lower trophic level !
31! model for the North Pacific marine ecosystem, Ecological !
32! Modelling, 202, 12-25. !
33! !
34!=======================================================================
35!
36 implicit none
37!
38 PRIVATE
39 PUBLIC :: biology
40!
41 CONTAINS
42!
43!***********************************************************************
44 SUBROUTINE biology (ng,tile)
45!***********************************************************************
46!
47 USE mod_param
48 USE mod_forces
49 USE mod_grid
50 USE mod_ncparam
51 USE mod_ocean
52 USE mod_stepping
53!
54! Imported variable declarations.
55!
56 integer, intent(in) :: ng, tile
57!
58! Local variable declarations.
59!
60 character (len=*), parameter :: MyFile = &
61 & __FILE__
62!
63#include "tile.h"
64!
65! Set header file name.
66!
67#ifdef DISTRIBUTE
68 IF (lbiofile(inlm)) THEN
69#else
70 IF (lbiofile(inlm).and.(tile.eq.0)) THEN
71#endif
72 lbiofile(inlm)=.false.
73 bioname(inlm)=myfile
74 END IF
75!
76#ifdef PROFILE
77 CALL wclock_on (ng, inlm, 15, __line__, myfile)
78#endif
79 CALL nemuro_tile (ng, tile, &
80 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
81 & imins, imaxs, jmins, jmaxs, &
82 & nstp(ng), nnew(ng), &
83#ifdef MASKING
84 & grid(ng) % rmask, &
85#endif
86 & grid(ng) % Hz, &
87 & grid(ng) % z_r, &
88 & grid(ng) % z_w, &
89 & forces(ng) % srflx, &
90 & ocean(ng) % t)
91
92#ifdef PROFILE
93 CALL wclock_off (ng, inlm, 15, __line__, myfile)
94#endif
95!
96 RETURN
97 END SUBROUTINE biology
98!
99!-----------------------------------------------------------------------
100 SUBROUTINE nemuro_tile (ng, tile, &
101 & LBi, UBi, LBj, UBj, UBk, UBt, &
102 & IminS, ImaxS, JminS, JmaxS, &
103 & nstp, nnew, &
104#ifdef MASKING
105 & rmask, &
106#endif
107 & Hz, z_r, z_w, &
108 & srflx, &
109 & t)
110!-----------------------------------------------------------------------
111!
112 USE mod_param
113 USE mod_biology
114 USE mod_ncparam
115 USE mod_scalars
116!
117! Imported variable declarations.
118!
119 integer, intent(in) :: ng, tile
120 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
121 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
122 integer, intent(in) :: nstp, nnew
123
124#ifdef ASSUMED_SHAPE
125# ifdef MASKING
126 real(r8), intent(in) :: rmask(LBi:,LBj:)
127# endif
128 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
129 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
130 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
131 real(r8), intent(in) :: srflx(LBi:,LBj:)
132 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
133#else
134# ifdef MASKING
135 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
136# endif
137 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
138 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
139 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
140 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
141 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
142#endif
143!
144! Local variable declarations.
145!
146 integer, parameter :: Nsink = 2
147
148 integer :: Iter, ibio, indx, isink, itime, itrc, iTrcMax
149 integer :: i, j, k, ks
150
151 integer, dimension(Nsink) :: idsink
152
153 real(r8), parameter :: MinVal = 1.0e-6_r8
154
155 real(r8) :: AttL, AttS, IrrL, IrrS, KappaL, KappaS
156 real(r8) :: dtdays, dz
157 real(r8) :: GppAPS, GppAPL, GppNPS, GppNPL, GppPS, GppPL
158 real(r8) :: GraPL2ZL, GraPL2ZP, GraPS2ZL, GraPS2ZS
159 real(r8) :: GraZL2ZP, GraZS2ZL, GraZS2ZP
160 real(r8) :: EgeZL, EgeZP, EgeZS
161 real(r8) :: ExcPL, ExcPS, ExcZL, ExcZP, ExcZS
162 real(r8) :: MorPL, MorPS
163 real(r8) :: ResPL, ResPS
164 real(r8) :: RnewL, RnewS
165 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5, cff6, cff7
166 real(r8) :: fac, fac1, fac2, fac3, fac4, fac5, fac6, fac7
167 real(r8) :: cffL, cffR, cu, dltL, dltR
168
169 real(r8), dimension(Nsink) :: Wbio
170
171 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
172
173 real(r8), dimension(IminS:ImaxS) :: PARsur
174
175 real(r8), dimension(NT(ng),2) :: BioTrc
176
177 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
178 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
179
180 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
181
182 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
183 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
184 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
185 real(r8), dimension(IminS:ImaxS,N(ng)) :: LightL
186 real(r8), dimension(IminS:ImaxS,N(ng)) :: LightS
187 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
188 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
189 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
190 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
191 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
192
193#include "set_bounds.h"
194!
195!-----------------------------------------------------------------------
196! Add biological Source/Sink terms.
197!-----------------------------------------------------------------------
198!
199! Avoid computing source/sink terms if no biological iterations.
200!
201 IF (bioiter(ng).le.0) RETURN
202!
203! Set time-stepping size (days) according to the number of iterations.
204!
205 dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
206!
207! Set vertical sinking indentification vector.
208!
209 idsink(1)=ipon_ ! particulate organic nitrogen
210 idsink(2)=iopal ! particulate organic silica
211!
212! Set vertical sinking velocity vector in the same order as the
213! identification vector, IDSINK.
214!
215 wbio(1)=setvpon(ng) ! particulate organic nitrogen
216 wbio(2)=setvopal(ng) ! particulate organic silica
217!
218! Compute inverse thickness to avoid repeated divisions.
219!
220 j_loop : DO j=jstr,jend
221 DO k=1,n(ng)
222 DO i=istr,iend
223 hz_inv(i,k)=1.0_r8/hz(i,j,k)
224 END DO
225 END DO
226 DO k=1,n(ng)-1
227 DO i=istr,iend
228 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
229 END DO
230 END DO
231 DO k=2,n(ng)-1
232 DO i=istr,iend
233 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
234 END DO
235 END DO
236!
237! Extract biological variables from tracer arrays, place them into
238! scratch arrays, and restrict their values to be positive definite.
239! At input, all tracers (index nnew) from predictor step have
240! transport units (m Tunits) since we do not have yet the new
241! values for zeta and Hz. These are known after the 2D barotropic
242! time-stepping.
243!
244 DO itrc=1,nbt
245 ibio=idbio(itrc)
246 DO k=1,n(ng)
247 DO i=istr,iend
248 bio_old(i,k,ibio)=max(0.0_r8,t(i,j,k,nstp,ibio))
249 bio(i,k,ibio)=bio_old(i,k,ibio)
250 END DO
251 END DO
252 END DO
253!
254! Extract potential temperature and salinity.
255!
256 DO k=1,n(ng)
257 DO i=istr,iend
258 bio(i,k,itemp)=t(i,j,k,nstp,itemp)
259 END DO
260 END DO
261!
262! Calculate surface Photosynthetically Available Radiation (PAR). The
263! net shortwave radiation is scaled back to Watts/m2 and multiplied by
264! the fraction that is photosynthetically available, PARfrac.
265!
266 DO i=istr,iend
267 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
268 END DO
269!
270!=======================================================================
271! Start internal iterations to achieve convergence of the nonlinear
272! backward-implicit solution.
273!=======================================================================
274!
275! During the iterative procedure a series of fractional time steps are
276! performed in a chained mode (splitting by different biological
277! conversion processes) in sequence of the main food chain. In all
278! stages the concentration of the component being consumed is treated
279! in a fully implicit manner, so the algorithm guarantees non-negative
280! values, no matter how strong the concentration of active consuming
281! component (Phytoplankton or Zooplankton). The overall algorithm,
282! as well as any stage of it, is formulated in conservative form
283! (except explicit sinking) in sense that the sum of concentration of
284! all components is conserved.
285!
286! In the implicit algorithm, we have for example (N: nutrient,
287! P: phytoplankton),
288!
289! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
290! {Michaelis-Menten}
291! below, we set
292! The N in the numerator of
293! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
294! as N(new)
295!
296! so the time-stepping of the equations becomes:
297!
298! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
299! consuming, divide by (1 + cff)
300! and
301!
302! P(new) = P(old) + cff * N(new) (2) when adding a source term,
303! growing, add (cff * source)
304!
305! Notice that if you substitute (1) in (2), you will get:
306!
307! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
308!
309! If you add (1) and (3), you get
310!
311! N(new) + P(new) = N(old) + P(old)
312!
313! implying conservation regardless how "cff" is computed. Therefore,
314! this scheme is unconditionally stable regardless of the conversion
315! rate. It does not generate negative values since the constituent
316! to be consumed is always treated implicitly. It is also biased
317! toward damping oscillations.
318!
319! The iterative loop below is to iterate toward an universal Backward-
320! Euler treatment of all terms. So if there are oscillations in the
321! system, they are only physical oscillations. These iterations,
322! however, do not improve the accuaracy of the solution.
323!
324 iter_loop: DO iter=1,bioiter(ng)
325!
326! Compute light attenuation as function of depth.
327!
328 cff1=1.0/vmaxs(ng)
329 cff2=1.0/vmaxl(ng)
330 DO i=istr,iend
331 atts=parsur(i)
332 attl=parsur(i)
333 IF (parsur(i).gt.0.0_r8) THEN ! day time
334 DO k=n(ng),1,-1
335!
336! Attenuate the light to the center of the grid cell using the
337! Platt et al. (1980) photoinhibition formulation. Here, AttSW is
338! the light attenuation due to seawater and AttPS and AttPL is the
339! attenuation due to Small and Large Phytoplankton (self-shading
340! coefficient).
341!
342 dz=0.5_r8*(z_w(i,j,k)-z_w(i,j,k-1))
343 kappas=attsw(ng)+ &
344 & attps(ng)*(bio(i,k,isphy)+bio(i,k,ilphy))
345 kappal=attsw(ng)+ &
346 & attpl(ng)*(bio(i,k,isphy)+bio(i,k,ilphy))
347 irrs=exp(-kappas*dz)
348 irrl=exp(-kappal*dz)
349 atts=atts*irrs
350 attl=attl*irrl
351 lights(i,k)=(1.0_r8-exp(-alphaps(ng)*atts*cff1))* &
352 & exp(-betaps(ng)*atts*cff1)
353 lightl(i,k)=(1.0_r8-exp(-alphapl(ng)*attl*cff2))* &
354 & exp(-betapl(ng)*attl*cff2)
355!
356! Attenuate the light to the bottom of the grid cell.
357!
358 atts=atts*irrs
359 attl=attl*irrl
360 END DO
361 ELSE ! night time
362 DO k=1,n(ng)
363 lights(i,k)=0.0_r8
364 lightl(i,k)=0.0_r8
365 END DO
366 END IF
367 END DO
368!
369!-----------------------------------------------------------------------
370! Phytoplankton primary productivity.
371!-----------------------------------------------------------------------
372!
373! Gross primary production of Small Phytoplankton consisting of
374! nutrient uptake (NO3 and NH4) terms, temperature-dependend term,
375! and light limitation term. The Michaelis-Menten curve is used to
376! describe the change in uptake rate as a function of nutrient
377! concentration.
378!
379 cff=dtdays*vmaxs(ng)
380 DO k=1,n(ng)
381 DO i=istr,iend
382!
383! Small Phytoplankton gross primary productivity, GppPS.
384!
385 cff1=cff*exp(kgpps(ng)*bio(i,k,itemp))*lights(i,k)* &
386 & bio(i,k,isphy)
387 cff2=cff1*exp(-pusais(ng)*bio(i,k,inh4_))/ &
388 & (kno3s(ng)+bio(i,k,ino3_))
389 cff3=cff1/(knh4s(ng)+bio(i,k,inh4_))
390 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff2)
391 bio(i,k,inh4_)=bio(i,k,inh4_)/(1.0_r8+cff3)
392 gppnps=bio(i,k,ino3_)*cff2
393 gppaps=bio(i,k,inh4_)*cff3
394 gppps=gppnps+gppaps
395 bio(i,k,isphy)=bio(i,k,isphy)+gppps
396!
397! Small Phytoplankton respiration rate, ResPS, assumed to be
398! proportional to biomass. Use ratio of NO3 uptake to total update
399! (NO3 + NH4) to compute respiration contributions.
400!
401 rnews=gppnps/max(minval,gppps)
402 cff4=dtdays*resps0(ng)*exp(kresps(ng)*bio(i,k,itemp))
403 bio(i,k,isphy)=bio(i,k,isphy)/(1.0_r8+cff4)
404 resps=bio(i,k,isphy)*cff4
405 bio(i,k,ino3_)=bio(i,k,ino3_)+resps*rnews
406 bio(i,k,inh4_)=bio(i,k,inh4_)+resps*(1.0_r8-rnews)
407!
408! Small Phytoplankton extracellular excrection rate, ExcPS, assumed to
409! be proportional to production.
410!
411 excps=gppps*gammas(ng)
412 bio(i,k,isphy)=bio(i,k,isphy)-excps
413 bio(i,k,idon_)=bio(i,k,idon_)+excps
414 END DO
415 END DO
416!
417! Gross primary production of Large Phytoplankton consisting of
418! nutrient uptake (NO3, NH4, Silicate) terms, temperature-dependend
419! term, and light limitation term. Notice that there is a silicate
420! limiting term (which complicates the implicit algorithm). Again,
421! the Michaelis-Menten curve is used to describe the change in
422! uptake rate as a function of nutrient concentration.
423!
424 cff=dtdays*vmaxl(ng)
425 fac1=1.0/rsin(ng)
426 fac2=dtdays*respl0(ng)
427 DO k=1,n(ng)
428 DO i=istr,iend
429!
430! Large Phytoplankton gross primary productivity, GppPL. Notice that
431! the primary productivity is limited by previous time-step silicate
432! concentration.
433!
434 cff1=cff*exp(kgppl(ng)*bio(i,k,itemp))*lightl(i,k)* &
435 & bio(i,k,ilphy)
436 cff2=exp(-pusail(ng)*bio(i,k,inh4_))/ &
437 & (kno3l(ng)+bio(i,k,ino3_))
438 cff3=1.0_r8/(knh4l(ng)+bio(i,k,inh4_))
439 cff4=cff2*bio(i,k,ino3_)
440 cff5=cff3*bio(i,k,inh4_)
441 cff6=bio(i,k,isioh)/(ksil(ng)+bio(i,k,isioh))
442 cff7=cff6/max(minval,cff4+cff5)
443 cff4=cff1*cff2*min(1.0_r8,cff7) ! Si limitation on N03
444 cff5=cff1*cff3*min(1.0_r8,cff7) ! Si limitation on NH4
445 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff4)
446 bio(i,k,inh4_)=bio(i,k,inh4_)/(1.0_r8+cff5)
447 gppnpl=bio(i,k,ino3_)*cff4
448 gppapl=bio(i,k,inh4_)*cff5
449 gpppl=gppnpl+gppapl
450 bio(i,k,ilphy)=bio(i,k,ilphy)+gpppl
451 bio(i,k,isioh)=bio(i,k,isioh)-gpppl*rsin(ng)
452!
453! Large Phytoplankton respiration rate, ResPL, assumed to be
454! proportional to biomass. Use ratio of NO3 uptake to total update
455! (NO3 + NH4) to compute respiration contributions. Use Si:N ratio to
456! compute SiOH4 contribution.
457!
458 rnewl=gppnpl/max(minval,gpppl)
459 cff7=fac2*exp(krespl(ng)*bio(i,k,itemp))
460 bio(i,k,ilphy)=bio(i,k,ilphy)/(1.0_r8+cff7)
461 respl=bio(i,k,ilphy)*cff7
462 bio(i,k,ino3_)=bio(i,k,ino3_)+respl*rnewl
463 bio(i,k,inh4_)=bio(i,k,inh4_)+respl*(1.0_r8-rnewl)
464 bio(i,k,isioh)=bio(i,k,isioh)+respl*rsin(ng)
465!
466! Large Phytoplankton extracellular excrection rate, ExcPL, assumed to
467! be proportional to production.
468!
469 excpl=gpppl*gammal(ng)
470 bio(i,k,ilphy)=bio(i,k,ilphy)-excpl
471 bio(i,k,idon_)=bio(i,k,idon_)+excpl
472 bio(i,k,isioh)=bio(i,k,isioh)+excpl*rsin(ng)
473 END DO
474 END DO
475!
476!-----------------------------------------------------------------------
477! Phytoplankton mortality to particulate organic nitrogen.
478!-----------------------------------------------------------------------
479!
480 fac1=dtdays*morps0(ng)
481 fac2=dtdays*morpl0(ng)
482 DO k=1,n(ng)
483 DO i=istr,iend
484 cff1=fac1*bio(i,k,isphy)*exp(kmorps(ng)*bio(i,k,itemp))
485 cff2=fac2*bio(i,k,ilphy)*exp(kmorpl(ng)*bio(i,k,itemp))
486 bio(i,k,isphy)=bio(i,k,isphy)/(1.0_r8+cff1)
487 bio(i,k,ilphy)=bio(i,k,ilphy)/(1.0_r8+cff2)
488 morps=bio(i,k,isphy)*cff1
489 morpl=bio(i,k,ilphy)*cff2
490 bio(i,k,ipon_)=bio(i,k,ipon_)+morps+morpl
491 bio(i,k,iopal)=bio(i,k,iopal)+morpl*rsin(ng)
492 END DO
493 END DO
494!
495!-----------------------------------------------------------------------
496! Zooplankton grazing, egestion and excretion.
497!-----------------------------------------------------------------------
498
499#if defined IVLEV_EXPLICIT
500!
501! The rate of grazing by the zooplankton is modeled using an Ivlev
502! equation with a feeding threshold using an explicit (non-conserving)
503! algorithm to the original formulation. Notice that term is forced to
504! be positive using a MAX function.
505!
506#elif defined HOLLING_GRAZING
507!
508! The rate of grazing by the zooplankton is modeled using a Holling-
509! type s-shaped curve. It is known to be numerically more stable and
510! allows an implicit discretization.
511!
512! P(new) = P(old) - dt*mu*[P(old)/(Kp + P(old)^2)]*P(new)*Z(old)
513!
514! The implicit grazing term is then:
515!
516! P(new) = P(old) / (1 + G)
517!
518! were the grazing rate, G, is:
519!
520! G = dt * mu * [P(old)/(Kp + P(old)^2)] * Z(old)
521!
522#else
523# define IVLEV_IMPLICIT
524!
525! The rate of grazing by the zooplankton is modeled using an Ivlev
526! equation with a feeding threshold. An implicit algorithm is
527! achieved by multiplying grazing term by a unity factor, alpha:
528!
529! alpha = 1 = P(new)/(P(old)+deltaP)
530!
531! where
532!
533! deltaP = P(new) - P(old) = - dt*mu*[1-EXP(lambda*P(old))]
534!
535! The factor alpha can be approximated using Taylor series to:
536!
537! alpha = [P(new) / P(old)] * [1 - deltaP/P(old)]
538!
539! The discretized grazing term is then:
540!
541! P(new) = P(old) - dt*mu*[1-EXP(lambda*P)] * alpha * Z(old)
542!
543! Which can be approximated with an implicit algorithm to:
544!
545! P(new) = P(old) / (1 + G)
546!
547! were the grazing rate, G, is:
548!
549! G = [1 + P(old)/(dt*mu*(1 - EXP(lambda*P(old))))] * Z(old)
550!
551#endif
552 fac1=dtdays*grmaxsps(ng)
553 fac2=dtdays*grmaxlps(ng)
554 fac3=dtdays*grmaxlpl(ng)
555 fac4=dtdays*grmaxlzs(ng)
556 fac5=dtdays*grmaxppl(ng)
557 fac6=dtdays*grmaxpzs(ng)
558 fac7=dtdays*grmaxpzl(ng)
559 DO k=1,n(ng)
560 DO i=istr,iend
561!
562! Temperature-dependent term (Q10).
563!
564 cff1=exp(kgras(ng)*bio(i,k,itemp))
565 cff2=exp(kgral(ng)*bio(i,k,itemp))
566 cff3=exp(kgrap(ng)*bio(i,k,itemp))
567!
568! Small Zooplankton grazing on Small Phytoplankton, GraPS2ZS.
569!
570#if defined IVLEV_EXPLICIT
571 cff4=1.0_r8-exp(lams(ng)*(ps2zsstar(ng)-bio(i,k,isphy)))
572 graps2zs=fac1*cff1*max(0.0_r8,cff4)*bio(i,k,iszoo)
573 bio(i,k,isphy)=bio(i,k,isphy)-graps2zs
574 bio(i,k,iszoo)=bio(i,k,iszoo)+graps2zs
575#else
576# ifdef HOLLING_GRAZING
577 cff4=1.0_r8/(kps2zs(ng)+bio(i,k,isphy)*bio(i,k,isphy))
578 cff=fac1*cff1*cff4*bio(i,k,iszoo)*bio(i,k,isphy)
579# elif defined IVLEV_IMPLICIT
580 cff4=1.0_r8-exp(lams(ng)*(ps2zsstar(ng)-bio(i,k,isphy)))
581 cff5=1.0_r8/(fac1*cff4)
582 cff=(1.0_r8+bio(i,k,isphy)*cff5)*cff1*bio(i,k,iszoo)
583# endif
584 bio(i,k,isphy)=bio(i,k,isphy)/(1.0_r8+cff)
585 graps2zs=cff*bio(i,k,isphy)
586 bio(i,k,iszoo)=bio(i,k,iszoo)+graps2zs
587#endif
588!
589! Large Zooplankton grazing on Small Phytoplankton, GraPS2ZL.
590!
591#if defined IVLEV_EXPLICIT
592 cff4=1.0_r8-exp(laml(ng)*(ps2zlstar(ng)-bio(i,k,isphy)))
593 graps2zl=fac2*cff2*max(0.0_r8,cff4)*bio(i,k,ilzoo)
594 bio(i,k,isphy)=bio(i,k,isphy)-graps2zl
595 bio(i,k,ilzoo)=bio(i,k,ilzoo)+graps2zl
596#else
597# ifdef HOLLING_GRAZING
598 cff4=1.0_r8/(kps2zl(ng)+bio(i,k,isphy)*bio(i,k,isphy))
599 cff=fac2*cff2*cff4*bio(i,k,ilzoo)*bio(i,k,isphy)
600# elif defined IVLEV_IMPLICIT
601 cff4=1.0_r8-exp(laml(ng)*(ps2zlstar(ng)-bio(i,k,isphy)))
602 cff5=1.0_r8/(fac2*cff4)
603 cff=(1.0_r8+bio(i,k,isphy)*cff5)*cff2*bio(i,k,ilzoo)
604# endif
605 bio(i,k,isphy)=bio(i,k,isphy)/(1.0_r8+cff)
606 graps2zl=cff*bio(i,k,isphy)
607 bio(i,k,ilzoo)=bio(i,k,ilzoo)+graps2zl
608#endif
609!
610! Large Zooplankton grazing on Large Phytoplankton, GraPL2ZL.
611!
612#if defined IVLEV_EXPLICIT
613 cff4=1.0_r8-exp(laml(ng)*(pl2zlstar(ng)-bio(i,k,ilphy)))
614 grapl2zl=fac3*cff2*max(0.0_r8,cff4)*bio(i,k,ilzoo)
615 bio(i,k,ilphy)=bio(i,k,ilphy)-grapl2zl
616 bio(i,k,ilzoo)=bio(i,k,ilzoo)+grapl2zl
617#else
618# ifdef HOLLING_GRAZING
619 cff4=1.0_r8/(kpl2zl(ng)+bio(i,k,ilphy)*bio(i,k,ilphy))
620 cff=fac3*cff2*cff4*bio(i,k,ilzoo)*bio(i,k,ilphy)
621# elif defined IVLEV_IMPLICIT
622 cff4=1.0_r8-exp(laml(ng)*(pl2zlstar(ng)-bio(i,k,ilphy)))
623 cff5=1.0_r8/(fac3*cff4)
624 cff=(1.0_r8+bio(i,k,ilphy)*cff5)*cff2*bio(i,k,ilzoo)
625# endif
626 bio(i,k,ilphy)=bio(i,k,ilphy)/(1.0_r8+cff)
627 grapl2zl=cff*bio(i,k,ilphy)
628 bio(i,k,ilzoo)=bio(i,k,ilzoo)+grapl2zl
629#endif
630!
631! Large Zooplankton grazing on Small Zooplankton, GraZS2ZL.
632!
633#if defined IVLEV_EXPLICIT
634 cff4=1.0_r8-exp(laml(ng)*(zs2zlstar(ng)-bio(i,k,iszoo)))
635 grazs2zl=fac4*cff2*max(0.0_r8,cff4)*bio(i,k,ilzoo)
636 bio(i,k,iszoo)=bio(i,k,iszoo)-grazs2zl
637 bio(i,k,ilzoo)=bio(i,k,ilzoo)+grazs2zl
638#else
639# ifdef HOLLING_GRAZING
640 cff4=1.0_r8/(kzs2zl(ng)+bio(i,k,iszoo)*bio(i,k,iszoo))
641 cff=fac4*cff2*cff4*bio(i,k,ilzoo)*bio(i,k,iszoo)
642# elif defined IVLEV_IMPLICIT
643 cff4=1.0_r8-exp(laml(ng)*(zs2zlstar(ng)-bio(i,k,iszoo)))
644 cff5=1.0_r8/(fac4*cff4)
645 cff=(1.0_r8+bio(i,k,isphy)*cff5)*cff2*bio(i,k,ilzoo)
646# endif
647 bio(i,k,iszoo)=bio(i,k,iszoo)/(1.0_r8+cff)
648 grazs2zl=cff*bio(i,k,iszoo)
649 bio(i,k,ilzoo)=bio(i,k,ilzoo)+grazs2zl
650#endif
651!
652! Predactor Zooplankton grazing on Large Phytoplankton, GraPL2ZP.
653!
654#if defined IVLEV_EXPLICIT
655 cff4=1.0_r8-exp(lamp(ng)*(pl2zpstar(ng)-bio(i,k,ilphy)))
656 cff5=exp(-pusaipl(ng)*(bio(i,k,ilzoo)+bio(i,k,iszoo)))
657 grapl2zp=fac5*cff3*cff5*max(0.0_r8,cff4)*bio(i,k,ipzoo)
658 bio(i,k,ilphy)=bio(i,k,ilphy)-grapl2zp
659 bio(i,k,ipzoo)=bio(i,k,ipzoo)+grapl2zp
660#else
661# ifdef HOLLING_GRAZING
662 cff4=1.0_r8/(kpl2zp(ng)+bio(i,k,ilphy)*bio(i,k,ilphy))
663 cff5=exp(-pusaipl(ng)*(bio(i,k,ilzoo)+bio(i,k,iszoo)))
664 cff=fac5*cff3*cff4*cff5*bio(i,k,ipzoo)*bio(i,k,ilphy)
665# elif defined IVLEV_IMPLICIT
666 cff4=1.0_r8-exp(lamp(ng)*(pl2zpstar(ng)-bio(i,k,ilphy)))
667 cff5=exp(-pusaipl(ng)*(bio(i,k,ilzoo)+bio(i,k,iszoo)))
668 cff6=1.0_r8/(fac5*cff4)
669 cff=(1.0_r8+bio(i,k,ilphy)*cff6)*cff3*cff5*bio(i,k,ipzoo)
670# endif
671 bio(i,k,ilphy)=bio(i,k,ilphy)/(1.0_r8+cff)
672 grapl2zp=cff*bio(i,k,ilphy)
673 bio(i,k,ipzoo)=bio(i,k,ipzoo)+grapl2zp
674#endif
675!
676! Predactory Zooplankton grazing on Small Zooplankton, GraZS2ZP.
677!
678#if defined IVLEV_EXPLICIT
679 cff4=1.0_r8-exp(lamp(ng)*(zs2zpstar(ng)-bio(i,k,iszoo)))
680 cff5=exp(-pusaizs(ng)*bio(i,k,ilzoo))
681 grazs2zp=fac6*cff3*cff5*max(0.0_r8,cff4)*bio(i,k,ipzoo)
682 bio(i,k,iszoo)=bio(i,k,iszoo)-grazs2zp
683 bio(i,k,ipzoo)=bio(i,k,ipzoo)+grazs2zp
684#else
685# ifdef HOLLING_GRAZING
686 cff4=1.0_r8/(kzs2zp(ng)+bio(i,k,iszoo)*bio(i,k,iszoo))
687 cff5=exp(-pusaizs(ng)*bio(i,k,ilzoo))
688 cff=fac6*cff3*cff4*cff5*bio(i,k,ipzoo)*bio(i,k,iszoo)
689# elif defined IVLEV_IMPLICIT
690 cff4=1.0_r8-exp(lamp(ng)*(zs2zpstar(ng)-bio(i,k,iszoo)))
691 cff5=exp(-pusaizs(ng)*bio(i,k,ilzoo))
692 cff6=1.0_r8/(fac6*cff4)
693 cff=(1.0_r8+bio(i,k,iszoo)*cff6)*cff3*cff5*bio(i,k,ipzoo)
694# endif
695 bio(i,k,iszoo)=bio(i,k,iszoo)/(1.0_r8+cff)
696 grazs2zp=cff*bio(i,k,iszoo)
697 bio(i,k,ipzoo)=bio(i,k,ipzoo)+grazs2zp
698#endif
699!
700! Predactory Zooplankton grazing on Large Zooplankton, GraZL2ZP.
701!
702#if defined IVLEV_EXPLICIT
703 cff4=1.0_r8-exp(lamp(ng)*(zl2zpstar(ng)-bio(i,k,ilzoo)))
704 grazl2zp=fac7*cff3*max(0.0_r8,cff4)*bio(i,k,ipzoo)
705 bio(i,k,ilzoo)=bio(i,k,ilzoo)-grazl2zp
706 bio(i,k,ipzoo)=bio(i,k,ipzoo)+grazl2zp
707#else
708# ifdef HOLLING_GRAZING
709 cff4=1.0_r8/(kzl2zp(ng)+bio(i,k,ilzoo)*bio(i,k,ilzoo))
710 cff=fac7*cff3*cff4*bio(i,k,ipzoo)*bio(i,k,ilzoo)
711# elif defined IVLEV_IMPLICIT
712 cff4=1.0_r8-exp(lamp(ng)*(zl2zpstar(ng)-bio(i,k,ilzoo)))
713 cff5=1.0_r8/(fac7*cff4)
714 cff=(1.0_r8+bio(i,k,ilzoo)*cff5)*cff3*bio(i,k,ipzoo)
715# endif
716 bio(i,k,ilzoo)=bio(i,k,ilzoo)/(1.0_r8+cff)
717 grazl2zp=cff*bio(i,k,ilzoo)
718 bio(i,k,ipzoo)=bio(i,k,ipzoo)+grazl2zp
719#endif
720!
721! Zooplankton egestion to Particulate Organic Nitrogen (PON) and
722! Particulate Organic Silica (opal).
723!
724 egezs=(1.0_r8-alphazs(ng))* &
725 & graps2zs
726 egezl=(1.0_r8-alphazl(ng))* &
727 & (graps2zl+grapl2zl+grazs2zl)
728 egezp=(1.0_r8-alphazp(ng))* &
729 & (grapl2zp+grazs2zp+grazl2zp)
730 bio(i,k,iszoo)=bio(i,k,iszoo)-egezs
731 bio(i,k,ilzoo)=bio(i,k,ilzoo)-egezl
732 bio(i,k,ipzoo)=bio(i,k,ipzoo)-egezp
733 bio(i,k,ipon_)=bio(i,k,ipon_)+egezs+egezl+egezp
734 bio(i,k,iopal)=bio(i,k,iopal)+(grapl2zl+grapl2zp)*rsin(ng)
735!
736! Zooplankton excretion to NH4.
737!
738 exczs=(alphazs(ng)-betazs(ng))* &
739 & graps2zs
740 exczl=(alphazl(ng)-betazl(ng))* &
741 & (graps2zl+grapl2zl+grazs2zl)
742 exczp=(alphazp(ng)-betazp(ng))* &
743 & (grapl2zp+grazs2zp+grazl2zp)
744 bio(i,k,iszoo)=bio(i,k,iszoo)-exczs
745 bio(i,k,ilzoo)=bio(i,k,ilzoo)-exczl
746 bio(i,k,ipzoo)=bio(i,k,ipzoo)-exczp
747 bio(i,k,inh4_)=bio(i,k,inh4_)+exczs+exczl+exczp
748 END DO
749 END DO
750!
751!-----------------------------------------------------------------------
752! Zooplankton motality to particulate organic nitrogen.
753!-----------------------------------------------------------------------
754!
755 fac1=dtdays*morzs0(ng)
756 fac2=dtdays*morzl0(ng)
757 fac3=dtdays*morzp0(ng)
758 DO k=1,n(ng)
759 DO i=istr,iend
760 cff1=fac1*bio(i,k,iszoo)*exp(kmorzs(ng)*bio(i,k,itemp))
761 cff2=fac2*bio(i,k,ilzoo)*exp(kmorzl(ng)*bio(i,k,itemp))
762 cff3=fac3*bio(i,k,ipzoo)*exp(kmorzp(ng)*bio(i,k,itemp))
763 bio(i,k,iszoo)=bio(i,k,iszoo)/(1.0_r8+cff1)
764 bio(i,k,ilzoo)=bio(i,k,ilzoo)/(1.0_r8+cff2)
765 bio(i,k,ipzoo)=bio(i,k,ipzoo)/(1.0_r8+cff3)
766 bio(i,k,ipon_)=bio(i,k,ipon_)+ &
767 & bio(i,k,iszoo)*cff1+ &
768 & bio(i,k,ilzoo)*cff2+ &
769 & bio(i,k,ipzoo)*cff3
770 END DO
771 END DO
772!
773!-----------------------------------------------------------------------
774! Nutrient decomposition.
775!-----------------------------------------------------------------------
776!
777 fac1=dtdays*nit0(ng)
778 fac2=dtdays*vp2n0(ng)
779 fac3=dtdays*vp2d0(ng)
780 fac4=dtdays*vd2n0(ng)
781 fac5=dtdays*vo2s0(ng)
782 DO k=1,n(ng)
783 DO i=istr,iend
784!
785! Nitrification: NH4 to NO3.
786!
787 cff1=fac1*exp(knit(ng)*bio(i,k,itemp))
788 bio(i,k,inh4_)=bio(i,k,inh4_)/(1.0_r8+cff1)
789 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
790 & bio(i,k,inh4_)*cff1
791!
792! Decomposition: PON to NH4.
793!
794 cff2=fac2*exp(kp2n(ng)*bio(i,k,itemp))
795 bio(i,k,ipon_)=bio(i,k,ipon_)/(1.0_r8+cff2)
796 bio(i,k,inh4_)=bio(i,k,inh4_)+ &
797 & bio(i,k,ipon_)*cff2
798!
799! Decomposition: PON to DON.
800!
801 cff3=fac3*exp(kp2d(ng)*bio(i,k,itemp))
802 bio(i,k,ipon_)=bio(i,k,ipon_)/(1.0_r8+cff3)
803 bio(i,k,idon_)=bio(i,k,idon_)+ &
804 & bio(i,k,ipon_)*cff3
805!
806! Decomposition: DON to NH4.
807!
808 cff4=fac4*exp(kd2n(ng)*bio(i,k,itemp))
809 bio(i,k,idon_)=bio(i,k,idon_)/(1.0_r8+cff4)
810 bio(i,k,inh4_)=bio(i,k,inh4_)+ &
811 & bio(i,k,idon_)*cff4
812!
813! Decomposition: Opal to SiOH4.
814!
815 cff5=fac5*exp(ko2s(ng)*bio(i,k,itemp))
816 bio(i,k,iopal)=bio(i,k,iopal)/(1.0_r8+cff5)
817 bio(i,k,isioh)=bio(i,k,isioh)+ &
818 & bio(i,k,iopal)*cff5
819 END DO
820 END DO
821!
822!-----------------------------------------------------------------------
823! Vertical sinking terms: PON and Opal.
824!-----------------------------------------------------------------------
825!
826! Reconstruct vertical profile of selected biological constituents
827! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
828! grid box. Then, compute semi-Lagrangian flux due to sinking.
829!
830 sink_loop: DO isink=1,nsink
831 ibio=idsink(isink)
832!
833! Copy concentration of biological particulates into scratch array
834! "qc" (q-central, restrict it to be positive) which is hereafter
835! interpreted as a set of grid-box averaged values for biogeochemical
836! constituent concentration.
837!
838 DO k=1,n(ng)
839 DO i=istr,iend
840 qc(i,k)=bio(i,k,ibio)
841 END DO
842 END DO
843!
844 DO k=n(ng)-1,1,-1
845 DO i=istr,iend
846 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
847 END DO
848 END DO
849 DO k=2,n(ng)-1
850 DO i=istr,iend
851 dltr=hz(i,j,k)*fc(i,k)
852 dltl=hz(i,j,k)*fc(i,k-1)
853 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
854 cffr=cff*fc(i,k)
855 cffl=cff*fc(i,k-1)
856!
857! Apply PPM monotonicity constraint to prevent oscillations within the
858! grid box.
859!
860 IF ((dltr*dltl).le.0.0_r8) THEN
861 dltr=0.0_r8
862 dltl=0.0_r8
863 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
864 dltr=cffl
865 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
866 dltl=cffr
867 END IF
868!
869! Compute right and left side values (bR,bL) of parabolic segments
870! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
871!
872! NOTE: Although each parabolic segment is monotonic within its grid
873! box, monotonicity of the whole profile is not guaranteed,
874! because bL(k+1)-bR(k) may still have different sign than
875! qc(i,k+1)-qc(i,k). This possibility is excluded,
876! after bL and bR are reconciled using WENO procedure.
877!
878 cff=(dltr-dltl)*hz_inv3(i,k)
879 dltr=dltr-cff*hz(i,j,k+1)
880 dltl=dltl+cff*hz(i,j,k-1)
881 br(i,k)=qc(i,k)+dltr
882 bl(i,k)=qc(i,k)-dltl
883 wr(i,k)=(2.0_r8*dltr-dltl)**2
884 wl(i,k)=(dltr-2.0_r8*dltl)**2
885 END DO
886 END DO
887 cff=1.0e-14_r8
888 DO k=2,n(ng)-2
889 DO i=istr,iend
890 dltl=max(cff,wl(i,k ))
891 dltr=max(cff,wr(i,k+1))
892 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
893 bl(i,k+1)=br(i,k)
894 END DO
895 END DO
896 DO i=istr,iend
897 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
898#if defined LINEAR_CONTINUATION
899 bl(i,n(ng))=br(i,n(ng)-1)
900 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
901#elif defined NEUMANN
902 bl(i,n(ng))=br(i,n(ng)-1)
903 br(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
904#else
905 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
906 bl(i,n(ng))=qc(i,n(ng)) ! conditions
907 br(i,n(ng)-1)=qc(i,n(ng))
908#endif
909#if defined LINEAR_CONTINUATION
910 br(i,1)=bl(i,2)
911 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
912#elif defined NEUMANN
913 br(i,1)=bl(i,2)
914 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
915#else
916 bl(i,2)=qc(i,1) ! bottom grid boxes are
917 br(i,1)=qc(i,1) ! re-assumed to be
918 bl(i,1)=qc(i,1) ! piecewise constant.
919#endif
920 END DO
921!
922! Apply monotonicity constraint again, since the reconciled interfacial
923! values may cause a non-monotonic behavior of the parabolic segments
924! inside the grid box.
925!
926 DO k=1,n(ng)
927 DO i=istr,iend
928 dltr=br(i,k)-qc(i,k)
929 dltl=qc(i,k)-bl(i,k)
930 cffr=2.0_r8*dltr
931 cffl=2.0_r8*dltl
932 IF ((dltr*dltl).lt.0.0_r8) THEN
933 dltr=0.0_r8
934 dltl=0.0_r8
935 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
936 dltr=cffl
937 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
938 dltl=cffr
939 END IF
940 br(i,k)=qc(i,k)+dltr
941 bl(i,k)=qc(i,k)-dltl
942 END DO
943 END DO
944!
945! After this moment reconstruction is considered complete. The next
946! stage is to compute vertical advective fluxes, FC. It is expected
947! that sinking may occurs relatively fast, the algorithm is designed
948! to be free of CFL criterion, which is achieved by allowing
949! integration bounds for semi-Lagrangian advective flux to use as
950! many grid boxes in upstream direction as necessary.
951!
952! In the two code segments below, WL is the z-coordinate of the
953! departure point for grid box interface z_w with the same indices;
954! FC is the finite volume flux; ksource(:,k) is index of vertical
955! grid box which contains the departure point (restricted by N(ng)).
956! During the search: also add in content of whole grid boxes
957! participating in FC.
958!
959 cff=dtdays*abs(wbio(isink))
960 DO k=1,n(ng)
961 DO i=istr,iend
962 fc(i,k-1)=0.0_r8
963 wl(i,k)=z_w(i,j,k-1)+cff
964 wr(i,k)=hz(i,j,k)*qc(i,k)
965 ksource(i,k)=k
966 END DO
967 END DO
968 DO k=1,n(ng)
969 DO ks=k,n(ng)-1
970 DO i=istr,iend
971 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
972 ksource(i,k)=ks+1
973 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
974 END IF
975 END DO
976 END DO
977 END DO
978!
979! Finalize computation of flux: add fractional part.
980!
981 DO k=1,n(ng)
982 DO i=istr,iend
983 ks=ksource(i,k)
984 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
985 fc(i,k-1)=fc(i,k-1)+ &
986 & hz(i,j,ks)*cu* &
987 & (bl(i,ks)+ &
988 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
989 & (1.5_r8-cu)* &
990 & (br(i,ks)+bl(i,ks)- &
991 & 2.0_r8*qc(i,ks))))
992 END DO
993 END DO
994 DO k=1,n(ng)
995 DO i=istr,iend
996 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
997 END DO
998 END DO
999
1000#ifdef BIO_SEDIMENT
1001!
1002! Particulate fluxes reaching the seafloor are remineralized and returned
1003! to the dissolved nutrient pool. Without this conversion, particulate
1004! material falls out of the system. This is a temporary fix to restore
1005! total nutrient conservation. This may require a time delay
1006! remineralization in the future.
1007!
1008! HGA: The original Nemuro model has a restoring upwelling rate (UPW).
1009! The code below is an interpretation in terms of the semi-
1010! Lagrangian algorithm. What is the correct nutrient path from
1011! the benthos to the water column? NH4 to NO3?
1012!
1013 IF (ibio.eq.ipon_) THEN
1014 DO i=istr,iend
1015 cff1=fc(i,0)*hz_inv(i,1)
1016 bio(i,1,ino3_)=bio(i,1,ino3_)+cff1
1017 END DO
1018 ELSE IF (ibio.eq.iopal) THEN
1019 DO i=istr,iend
1020 cff1=fc(i,0)*hz_inv(i,1)
1021 bio(i,1,isioh)=bio(i,1,isioh)+cff1
1022 END DO
1023 END IF
1024#endif
1025
1026 END DO sink_loop
1027 END DO iter_loop
1028!
1029!-----------------------------------------------------------------------
1030! Update global tracer variables: Add increment due to BGC processes
1031! to tracer array in time index "nnew". Index "nnew" is solution after
1032! advection and mixing and has transport units (m Tunits) hence the
1033! increment is multiplied by Hz. Notice that we need to subtract
1034! original values "Bio_old" at the top of the routine to just account
1035! for the concentractions affected by BGC processes. This also takes
1036! into account any constraints (non-negative concentrations, carbon
1037! concentration range) specified before entering BGC kernel. If "Bio"
1038! were unchanged by BGC processes, the increment would be exactly
1039! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
1040! bounded >=0 so that we can preserve total inventory of nutrients
1041! when advection causes tracer concentration to go negative.
1042!-----------------------------------------------------------------------
1043!
1044 DO itrc=1,nbt
1045 ibio=idbio(itrc)
1046 DO k=1,n(ng)
1047 DO i=istr,iend
1048 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
1049 t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*hz(i,j,k)
1050 END DO
1051 END DO
1052 END DO
1053
1054 END DO j_loop
1055!
1056 RETURN
1057 END SUBROUTINE nemuro_tile
1058
1059 END MODULE biology_mod
subroutine, public biology(ng, tile)
Definition ecosim.h:57
subroutine nemuro_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nnew, rmask, hz, z_r, z_w, srflx, t)
Definition nemuro.h:110
real(r8), dimension(:), allocatable pl2zlstar
Definition nemuro_mod.h:259
integer isioh
Definition nemuro_mod.h:187
real(r8), dimension(:), allocatable kps2zs
Definition nemuro_mod.h:237
real(r8), dimension(:), allocatable parfrac
Definition fennel_mod.h:139
real(r8), dimension(:), allocatable kmorzl
Definition nemuro_mod.h:224
real(r8), dimension(:), allocatable kpl2zl
Definition nemuro_mod.h:235
integer ipzoo
Definition nemuro_mod.h:182
real(r8), dimension(:), allocatable ko2s
Definition nemuro_mod.h:232
real(r8), dimension(:), allocatable pusaipl
Definition nemuro_mod.h:256
real(r8), dimension(:), allocatable knit
Definition nemuro_mod.h:229
real(r8), dimension(:), allocatable pusais
Definition nemuro_mod.h:257
real(r8), dimension(:), allocatable kzs2zl
Definition nemuro_mod.h:243
real(r8), dimension(:), allocatable vp2n0
Definition nemuro_mod.h:273
integer ipon_
Definition nemuro_mod.h:185
real(r8), dimension(:), allocatable alphaps
Definition nemuro_mod.h:195
integer ilphy
Definition nemuro_mod.h:178
real(r8), dimension(:), allocatable vmaxs
Definition nemuro_mod.h:270
real(r8), dimension(:), allocatable vp2d0
Definition nemuro_mod.h:272
real(r8), dimension(:), allocatable ps2zsstar
Definition nemuro_mod.h:262
integer, dimension(:), allocatable bioiter
Definition ecosim_mod.h:343
real(r8), dimension(:), allocatable zs2zpstar
Definition nemuro_mod.h:276
real(r8), dimension(:), allocatable attsw
Definition fennel_mod.h:125
real(r8), dimension(:), allocatable knh4l
Definition nemuro_mod.h:227
real(r8), dimension(:), allocatable kgpps
Definition nemuro_mod.h:218
real(r8), dimension(:), allocatable grmaxppl
Definition nemuro_mod.h:212
integer ino3_
Definition ecosim_mod.h:277
integer iszoo
Definition nemuro_mod.h:181
real(r8), dimension(:), allocatable lams
Definition nemuro_mod.h:247
real(r8), dimension(:), allocatable kzs2zp
Definition nemuro_mod.h:244
real(r8), dimension(:), allocatable kps2zl
Definition nemuro_mod.h:236
real(r8), dimension(:), allocatable vd2n0
Definition nemuro_mod.h:268
real(r8), dimension(:), allocatable setvpon
Definition nemuro_mod.h:267
real(r8), dimension(:), allocatable kmorzs
Definition nemuro_mod.h:226
real(r8), dimension(:), allocatable vmaxl
Definition nemuro_mod.h:269
real(r8), dimension(:), allocatable morpl0
Definition nemuro_mod.h:248
real(r8), dimension(:), allocatable kmorps
Definition nemuro_mod.h:223
real(r8), dimension(:), allocatable rsin
Definition nemuro_mod.h:265
real(r8), dimension(:), allocatable grmaxpzs
Definition nemuro_mod.h:214
real(r8), dimension(:), allocatable vo2s0
Definition nemuro_mod.h:271
real(r8), dimension(:), allocatable pusail
Definition nemuro_mod.h:255
real(r8), dimension(:), allocatable kgppl
Definition nemuro_mod.h:217
real(r8), dimension(:), allocatable kresps
Definition nemuro_mod.h:240
real(r8), dimension(:), allocatable grmaxlzs
Definition nemuro_mod.h:211
real(r8), dimension(:), allocatable krespl
Definition nemuro_mod.h:239
real(r8), dimension(:), allocatable kgras
Definition nemuro_mod.h:221
real(r8), dimension(:), allocatable alphazl
Definition nemuro_mod.h:196
real(r8), dimension(:), allocatable kgral
Definition nemuro_mod.h:219
real(r8), dimension(:), allocatable betapl
Definition nemuro_mod.h:202
real(r8), dimension(:), allocatable resps0
Definition nemuro_mod.h:264
real(r8), dimension(:), allocatable attpl
Definition nemuro_mod.h:199
real(r8), dimension(:), allocatable kp2n
Definition nemuro_mod.h:234
integer isphy
Definition nemuro_mod.h:179
real(r8), dimension(:), allocatable attps
Definition nemuro_mod.h:200
real(r8), dimension(:), allocatable morzl0
Definition nemuro_mod.h:250
real(r8), dimension(:), allocatable betazp
Definition nemuro_mod.h:206
real(r8), dimension(:), allocatable morzs0
Definition nemuro_mod.h:252
real(r8), dimension(:), allocatable gammal
Definition nemuro_mod.h:207
real(r8), dimension(:), allocatable lamp
Definition nemuro_mod.h:246
integer inh4_
Definition ecosim_mod.h:278
integer, dimension(:), allocatable idbio
Definition ecosim_mod.h:256
real(r8), dimension(:), allocatable alphapl
Definition nemuro_mod.h:194
real(r8), dimension(:), allocatable kmorzp
Definition nemuro_mod.h:225
real(r8), dimension(:), allocatable kpl2zp
Definition nemuro_mod.h:238
real(r8), dimension(:), allocatable ps2zlstar
Definition nemuro_mod.h:261
real(r8), dimension(:), allocatable grmaxlps
Definition nemuro_mod.h:210
real(r8), dimension(:), allocatable kd2n
Definition nemuro_mod.h:216
real(r8), dimension(:), allocatable zl2zpstar
Definition nemuro_mod.h:274
real(r8), dimension(:), allocatable nit0
Definition nemuro_mod.h:253
real(r8), dimension(:), allocatable kno3l
Definition nemuro_mod.h:230
real(r8), dimension(:), allocatable alphazs
Definition nemuro_mod.h:198
real(r8), dimension(:), allocatable kzl2zp
Definition nemuro_mod.h:242
real(r8), dimension(:), allocatable betazl
Definition nemuro_mod.h:205
real(r8), dimension(:), allocatable gammas
Definition nemuro_mod.h:208
real(r8), dimension(:), allocatable setvopal
Definition nemuro_mod.h:266
real(r8), dimension(:), allocatable grmaxlpl
Definition nemuro_mod.h:209
real(r8), dimension(:), allocatable respl0
Definition nemuro_mod.h:263
real(r8), dimension(:), allocatable kno3s
Definition nemuro_mod.h:231
real(r8), dimension(:), allocatable morzp0
Definition nemuro_mod.h:251
real(r8), dimension(:), allocatable kmorpl
Definition nemuro_mod.h:222
real(r8), dimension(:), allocatable pl2zpstar
Definition nemuro_mod.h:260
real(r8), dimension(:), allocatable pusaizs
Definition nemuro_mod.h:258
integer iopal
Definition nemuro_mod.h:188
real(r8), dimension(:), allocatable alphazp
Definition nemuro_mod.h:197
real(r8), dimension(:), allocatable grmaxsps
Definition nemuro_mod.h:215
real(r8), dimension(:), allocatable kgrap
Definition nemuro_mod.h:220
real(r8), dimension(:), allocatable betazs
Definition nemuro_mod.h:204
real(r8), dimension(:), allocatable knh4s
Definition nemuro_mod.h:228
real(r8), dimension(:), allocatable ksil
Definition nemuro_mod.h:241
integer idon_
Definition nemuro_mod.h:186
real(r8), dimension(:), allocatable grmaxpzl
Definition nemuro_mod.h:213
real(r8), dimension(:), allocatable laml
Definition nemuro_mod.h:245
real(r8), dimension(:), allocatable morps0
Definition nemuro_mod.h:249
real(r8), dimension(:), allocatable kp2d
Definition nemuro_mod.h:233
real(r8), dimension(:), allocatable betaps
Definition nemuro_mod.h:203
real(r8), dimension(:), allocatable zs2zlstar
Definition nemuro_mod.h:275
integer ilzoo
Definition nemuro_mod.h:180
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
integer, parameter inlm
Definition mod_param.F:662
integer, dimension(:), allocatable n
Definition mod_param.F:479
integer nbt
Definition mod_param.F:509
integer, dimension(:), allocatable nt
Definition mod_param.F:489
real(dp), dimension(:), allocatable dt
real(dp) cp
real(dp), parameter sec2day
integer itemp
real(dp) rho0
integer, dimension(:), allocatable nnew
integer, dimension(:), allocatable nstp
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3