ROMS
Loading...
Searching...
No Matches
tl_npzd_iron.h
Go to the documentation of this file.
1 MODULE tl_biology_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! Nutrient-Phytoplankton-Zooplankton-Detritus Model, !
11! including Iron Limitation on Phytoplankton Growth. !
12! !
13! This routine computes the biological sources and sinks and adds !
14! then the global biological fields. !
15! !
16! Reference: !
17! !
18! Fiechter, J., A.M. Moore, C.A. Edwards, K.W. Bruland, !
19! E. Di Lorenzo, C.V.W. Lewis, T.M. Powell, E. Curchitser !
20! and K. Hedstrom, 2009: Modeling iron limitation of primary !
21! production in the coastal Gulf of Alaska, Deep Sea Res. II !
22! !
23!=======================================================================
24!
25 implicit none
26!
27 PRIVATE
28 PUBLIC :: tl_biology
29!
30 CONTAINS
31!
32!***********************************************************************
33 SUBROUTINE tl_biology (ng,tile)
34!***********************************************************************
35!
36 USE mod_param
37 USE mod_forces
38 USE mod_grid
39 USE mod_ncparam
40 USE mod_ocean
41 USE mod_stepping
42!
43! Imported variable declarations.
44!
45 integer, intent(in) :: ng, tile
46!
47! Local variable declarations.
48!
49 character (len=*), parameter :: MyFile = &
50 & __FILE__
51!
52#include "tile.h"
53!
54! Set header file name.
55!
56#ifdef DISTRIBUTE
57 IF (lbiofile(itlm)) THEN
58#else
59 IF (lbiofile(itlm).and.(tile.eq.0)) THEN
60#endif
61 lbiofile(itlm)=.false.
62 bioname(itlm)=myfile
63 END IF
64!
65#ifdef PROFILE
66 CALL wclock_on (ng, itlm, 15, __line__, myfile)
67#endif
68 CALL tl_npzd_iron_tile (ng, tile, &
69 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
70 & imins, imaxs, jmins, jmaxs, &
71 & nstp(ng), nnew(ng), &
72#ifdef MASKING
73 & grid(ng) % rmask, &
74#endif
75#if defined IRON_LIMIT && defined IRON_RELAX
76 & grid(ng) % h, &
77#endif
78 & grid(ng) % Hz, &
79 & grid(ng) % tl_Hz, &
80 & grid(ng) % z_r, &
81 & grid(ng) % tl_z_r, &
82 & grid(ng) % z_w, &
83 & grid(ng) % tl_z_w, &
84 & forces(ng) % srflx, &
85 & forces(ng) % tl_srflx, &
86 & ocean(ng) % t, &
87 & ocean(ng) % tl_t)
88
89#ifdef PROFILE
90 CALL wclock_off (ng, itlm, 15, __line__, myfile)
91#endif
92!
93 RETURN
94 END SUBROUTINE tl_biology
95!
96!-----------------------------------------------------------------------
97 SUBROUTINE tl_npzd_iron_tile (ng, tile, &
98 & LBi, UBi, LBj, UBj, UBk, UBt, &
99 & IminS, ImaxS, JminS, JmaxS, &
100 & nstp, nnew, &
101#ifdef MASKING
102 & rmask, &
103#endif
104#if defined IRON_LIMIT && defined IRON_RELAX
105 & h, &
106#endif
107 & Hz, tl_Hz, &
108 & z_r, tl_z_r, &
109 & z_w, tl_z_w, &
110 & srflx, tl_srflx, &
111 & t, tl_t)
112!-----------------------------------------------------------------------
113!
114 USE mod_param
115 USE mod_biology
116 USE mod_ncparam
117 USE mod_scalars
118!
119! Imported variable declarations.
120!
121 integer, intent(in) :: ng, tile
122 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
123 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
124 integer, intent(in) :: nstp, nnew
125
126#ifdef ASSUMED_SHAPE
127# ifdef MASKING
128 real(r8), intent(in) :: rmask(LBi:,LBj:)
129# endif
130# if defined IRON_LIMIT && defined IRON_RELAX
131 real(r8), intent(in) :: h(LBi:,LBj:)
132# endif
133 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
134 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
135 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
136 real(r8), intent(in) :: srflx(LBi:,LBj:)
137 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
138
139 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
140 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
141 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
142 real(r8), intent(in) :: tl_srflx(LBi:,LBj:)
143 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
144#else
145# ifdef MASKING
146 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
147# endif
148# if defined IRON_LIMIT && defined IRON_RELAX
149 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
150# endif
151 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
152 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
153 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
154 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
155 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
156
157 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,UBk)
158 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,UBk)
159 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:UBk)
160 real(r8), intent(in) :: tl_srflx(LBi:UBi,LBj:UBj)
161 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,UBk,3,UBt)
162#endif
163!
164! Local variable declarations.
165!
166 integer, parameter :: Nsink = 2
167
168 integer :: Iter, i, ibio, isink, itime, itrc, iTrcMax, j, k, ks
169 integer :: Iteradj
170
171 integer, dimension(Nsink) :: idsink
172
173 real(r8), parameter :: MinVal = 1.0e-6_r8
174
175 real(r8) :: Att, ExpAtt, Itop, PAR
176 real(r8) :: tl_Att, tl_ExpAtt, tl_Itop, tl_PAR
177 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5, cff6, dtdays
178 real(r8) :: tl_cff, tl_cff1, tl_cff4, tl_cff5, tl_cff6
179 real(r8) :: cffL, cffR, cu, dltL, dltR
180 real(r8) :: tl_cffL, tl_cffR, tl_cu, tl_dltL, tl_dltR
181 real(r8) :: fac, fac1, fac2
182 real(r8) :: tl_fac, tl_fac1, tl_fac2
183#ifdef IRON_LIMIT
184 real(r8) :: Nlimit, FNlim
185 real(r8) :: tl_Nlimit, tl_FNlim
186 real(r8) :: FNratio, FCratio, FCratioE, Flimit
187 real(r8) :: tl_FNratio, tl_FCratio, tl_FCratioE, tl_Flimit
188 real(r8) :: FeC2FeN, FeN2FeC
189# ifdef IRON_RELAX
190 real(r8) :: FeNudgCoef
191# endif
192#endif
193 real(r8), dimension(Nsink) :: Wbio
194 real(r8), dimension(Nsink) :: tl_Wbio
195
196 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
197
198 real(r8), dimension(IminS:ImaxS) :: PARsur
199 real(r8), dimension(IminS:ImaxS) :: tl_PARsur
200
201 real(r8), dimension(NT(ng),2) :: BioTrc
202 real(r8), dimension(NT(ng),2) :: BioTrc1
203 real(r8), dimension(NT(ng),2) :: tl_BioTrc
204 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
205 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio1
206 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio2
207 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
208
209 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: tl_Bio
210 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: tl_Bio_old
211
212 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
213 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FC
214
215 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
216 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
217 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
218 real(r8), dimension(IminS:ImaxS,N(ng)) :: Light
219 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
220 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
221 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
222 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL1
223 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
224 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR1
225 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
226
227 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv
228 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv2
229 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv3
230 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Light
231 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_WL
232 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_WR
233 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bL
234 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bR
235 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_qc
236
237#include "set_bounds.h"
238!
239!-----------------------------------------------------------------------
240! Add biological Source/Sink terms.
241!-----------------------------------------------------------------------
242!
243! Avoid computing source/sink terms if no biological iterations.
244!
245 IF (bioiter(ng).le.0) RETURN
246!
247! Set time-stepping size (days) according to the number of iterations.
248!
249 dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
250
251#if defined IRON_LIMIT && defined IRON_RELAX
252!
253! Set nudging coefficient for dissolved iron over the shelf.
254!
255 fenudgcoef=dt(ng)/(fenudgtime(ng)*86400.0_r8)
256#endif
257#ifdef IRON_LIMIT
258!
259! Set Fe:N and Fe:C conversion ratio and its inverse.
260!
261 fen2fec=(16.0_r8/106.0_r8)*1.0e3_r8
262 fec2fen=(106.0_r8/16.0_r8)*1.0e-3_r8
263#endif
264!
265! Set vertical sinking indentification vector.
266!
267 idsink(1)=iphyt ! Phytoplankton
268 idsink(2)=isdet ! Small detritus
269!
270! Set vertical sinking velocity vector in the same order as the
271! identification vector, IDSINK.
272!
273 wbio(1)=wphy(ng) ! Phytoplankton
274 tl_wbio(1)=tl_wphy(ng) ! Phytoplankton
275 wbio(2)=wdet(ng) ! Small detritus
276 tl_wbio(2)=tl_wdet(ng) ! Small detritus
277!
278 j_loop : DO j=jstr,jend
279!
280! Compute inverse thickness to avoid repeated divisions.
281!
282 DO k=1,n(ng)
283 DO i=istr,iend
284 hz_inv(i,k)=1.0_r8/hz(i,j,k)
285 tl_hz_inv(i,k)=-hz_inv(i,k)*hz_inv(i,k)*tl_hz(i,j,k)
286 END DO
287 END DO
288 DO k=1,n(ng)-1
289 DO i=istr,iend
290 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
291 tl_hz_inv2(i,k)=-hz_inv2(i,k)*hz_inv2(i,k)* &
292 & (tl_hz(i,j,k)+tl_hz(i,j,k+1))
293 END DO
294 END DO
295 DO k=2,n(ng)-1
296 DO i=istr,iend
297 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
298 tl_hz_inv3(i,k)=-hz_inv3(i,k)*hz_inv3(i,k)* &
299 & (tl_hz(i,j,k-1)+tl_hz(i,j,k)+ &
300 & tl_hz(i,j,k+1))
301 END DO
302 END DO
303!
304! Clear tl_Bio and Bio arrays.
305!
306 DO itrc=1,nbt
307 ibio=idbio(itrc)
308 DO k=1,n(ng)
309 DO i=istr,iend
310 bio(i,k,ibio)=0.0_r8
311 bio1(i,k,ibio)=0.0_r8
312 bio2(i,k,ibio)=0.0_r8
313 tl_bio(i,k,ibio)=0.0_r8
314 END DO
315 END DO
316 END DO
317!
318! Restrict biological tracer to be positive definite. If a negative
319! concentration is detected, nitrogen is drawn from the most abundant
320! pool to supplement the negative pools to a lower limit of MinVal
321! which is set to 1E-6 above.
322!
323 DO k=1,n(ng)
324 DO i=istr,iend
325!
326! At input, all tracers (index nnew) from predictor step have
327! transport units (m Tunits) since we do not have yet the new
328! values for zeta and Hz. These are known after the 2D barotropic
329! time-stepping.
330!
331! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
332! tracer times depth. However the basic state (nstp and nnew
333! indices) that is read from the forward file is in units of
334! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
335! use t instead of t*Hz_inv.
336!
337 DO itrc=1,nbt
338 ibio=idbio(itrc)
339!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
340!^
341 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
342 tl_biotrc(ibio,nstp)=tl_t(i,j,k,nstp,ibio)
343!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
344!^
345 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
346 tl_biotrc(ibio,nnew)=tl_t(i,j,k,nnew,ibio)* &
347 & hz_inv(i,k)+ &
348 & t(i,j,k,nnew,ibio)*hz(i,j,k)* &
349 & tl_hz_inv(i,k)
350 END DO
351!
352! Impose positive definite concentrations.
353!
354 cff2=0.0_r8
355 DO itime=1,2
356 cff1=0.0_r8
357 tl_cff1=0.0_r8
358 itrcmax=idbio(1)
359#ifdef IRON_LIMIT
360 DO itrc=1,nbt-2
361#else
362 DO itrc=1,nbt
363#endif
364 ibio=idbio(itrc)
365 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
366 tl_cff1=tl_cff1- &
367 & (0.5_r8-sign(0.5_r8, &
368 & biotrc(ibio,itime)-minval))* &
369 & tl_biotrc(ibio,itime)
370 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
371 itrcmax=ibio
372 END IF
373 biotrc1(ibio,itime)=biotrc(ibio,itime)
374 biotrc(ibio,itime)=max(minval,biotrc1(ibio,itime))
375 tl_biotrc(ibio,itime)=(0.5_r8- &
376 & sign(0.5_r8, &
377 & minval- &
378 & biotrc1(ibio,itime)))* &
379 & tl_biotrc(ibio,itime)
380 END DO
381 IF (biotrc(itrcmax,itime).gt.cff1) THEN
382 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
383 tl_biotrc(itrcmax,itime)=tl_biotrc(itrcmax,itime)- &
384 & tl_cff1
385 END IF
386#ifdef IRON_LIMIT
387 DO itrc=nbt-1,nbt
388 ibio=idbio(itrc)
389 biotrc1(ibio,itime)=biotrc(ibio,itime)
390 biotrc(ibio,itime)=max(minval,biotrc1(ibio,itime))
391 tl_biotrc(ibio,itime)=(0.5_r8- &
392 & sign(0.5_r8, &
393 & minval- &
394 & biotrc1(ibio,itime)))* &
395 & tl_biotrc(ibio,itime)
396 END DO
397#endif
398 END DO
399!
400! Load biological tracers into local arrays.
401!
402 DO itrc=1,nbt
403 ibio=idbio(itrc)
404 bio_old(i,k,ibio)=biotrc(ibio,nstp)
405 tl_bio_old(i,k,ibio)=tl_biotrc(ibio,nstp)
406 bio(i,k,ibio)=biotrc(ibio,nstp)
407 tl_bio(i,k,ibio)=tl_biotrc(ibio,nstp)
408 END DO
409
410#if defined IRON_LIMIT && defined IRON_RELAX
411!
412! Relax dissolved iron at coast (h <= FeHim) to a constant value
413! (FeMax) over a time scale (FeNudgTime; days) to simulate sources
414! at the shelf.
415!
416 IF (h(i,j).le.fehmin(ng)) THEN
417!^ Bio(i,k,iFdis)=Bio(i,k,iFdis)+ &
418!^ & FeNudgCoef*(FeMax(ng)-Bio(i,k,iFdis))
419!^
420 tl_bio(i,k,ifdis)=tl_bio(i,k,ifdis)- &
421 & fenudgcoef*tl_bio(i,k,ifdis)
422 END IF
423#endif
424 END DO
425 END DO
426!
427! Calculate surface Photosynthetically Available Radiation (PAR). The
428! net shortwave radiation is scaled back to Watts/m2 and multiplied by
429! the fraction that is photosynthetically available, PARfrac.
430!
431 DO i=istr,iend
432#ifdef CONST_PAR
433!
434! Specify constant surface irradiance a la Powell and Spitz.
435!
436 parsur(i)=158.075_r8
437 tl_parsur(i)=0.0_r8
438#else
439 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
440 tl_parsur(i)=(tl_parfrac(ng)*srflx(i,j)+ &
441 & parfrac(ng)*tl_srflx(i,j))*rho0*cp
442#endif
443 END DO
444!
445!=======================================================================
446! Start internal iterations to achieve convergence of the nonlinear
447! backward-implicit solution.
448!=======================================================================
449!
450! During the iterative procedure a series of fractional time steps are
451! performed in a chained mode (splitting by different biological
452! conversion processes) in sequence of the main food chain. In all
453! stages the concentration of the component being consumed is treated
454! in a fully implicit manner, so the algorithm guarantees non-negative
455! values, no matter how strong the concentration of active consuming
456! component (Phytoplankton or Zooplankton). The overall algorithm,
457! as well as any stage of it, is formulated in conservative form
458! (except explicit sinking) in sense that the sum of concentration of
459! all components is conserved.
460!
461! In the implicit algorithm, we have for example (N: nutrient,
462! P: phytoplankton),
463!
464! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
465! {Michaelis-Menten}
466! below, we set
467! The N in the numerator of
468! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
469! as N(new)
470!
471! so the time-stepping of the equations becomes:
472!
473! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
474! consuming, divide by (1 + cff)
475! and
476!
477! P(new) = P(old) + cff * N(new) (2) when adding a source term,
478! growing, add (cff * source)
479!
480! Notice that if you substitute (1) in (2), you will get:
481!
482! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
483!
484! If you add (1) and (3), you get
485!
486! N(new) + P(new) = N(old) + P(old)
487!
488! implying conservation regardless how "cff" is computed. Therefore,
489! this scheme is unconditionally stable regardless of the conversion
490! rate. It does not generate negative values since the constituent
491! to be consumed is always treated implicitly. It is also biased
492! toward damping oscillations.
493!
494! The iterative loop below is to iterate toward an universal Backward-
495! Euler treatment of all terms. So if there are oscillations in the
496! system, they are only physical oscillations. These iterations,
497! however, do not improve the accuaracy of the solution.
498!
499 iter_loop: DO iter=1,bioiter(ng)
500!
501! Compute appropriate basic state arrays I.
502!
503 DO k=1,n(ng)
504 DO i=istr,iend
505!
506! At input, all tracers (index nnew) from predictor step have
507! transport units (m Tunits) since we do not have yet the new
508! values for zeta and Hz. These are known after the 2D barotropic
509! time-stepping.
510!
511! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
512! tracer times depth. However the basic state (nstp and nnew
513! indices) that is read from the forward file is in units of
514! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
515! use t instead of t*Hz_inv.
516!
517 DO itrc=1,nbt
518 ibio=idbio(itrc)
519!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
520!^
521 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
522!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
523!^
524 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
525 END DO
526!
527! Impose positive definite concentrations.
528!
529 cff2=0.0_r8
530 DO itime=1,2
531 cff1=0.0_r8
532 itrcmax=idbio(1)
533#ifdef IRON_LIMIT
534 DO itrc=1,nbt-2
535#else
536 DO itrc=1,nbt
537#endif
538 ibio=idbio(itrc)
539 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
540 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
541 itrcmax=ibio
542 END IF
543 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
544 END DO
545 IF (biotrc(itrcmax,itime).gt.cff1) THEN
546 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
547 END IF
548#ifdef IRON_LIMIT
549 DO itrc=nbt-1,nbt
550 ibio=idbio(itrc)
551 biotrc1(ibio,itime)=biotrc(ibio,itime)
552 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
553 END DO
554#endif
555 END DO
556!
557! Load biological tracers into local arrays.
558!
559 DO itrc=1,nbt
560 ibio=idbio(itrc)
561 bio_old(i,k,ibio)=biotrc(ibio,nstp)
562 bio(i,k,ibio)=biotrc(ibio,nstp)
563 END DO
564
565#if defined IRON_LIMIT && defined IRON_RELAX
566!
567! Relax dissolved iron at coast (h <= FeHim) to a constant value
568! (FeMax) over a time scale (FeNudgTime; days) to simulate sources
569! at the shelf.
570!
571 IF (h(i,j).le.fehmin(ng)) THEN
572 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
573 & fenudgcoef*(femax(ng)-bio(i,k,ifdis))
574 END IF
575#endif
576 END DO
577 END DO
578!
579! Calculate surface Photosynthetically Available Radiation (PAR). The
580! net shortwave radiation is scaled back to Watts/m2 and multiplied by
581! the fraction that is photosynthetically available, PARfrac.
582!
583 DO i=istr,iend
584#ifdef CONST_PAR
585!
586! Specify constant surface irradiance a la Powell and Spitz.
587!
588 parsur(i)=158.075_r8
589#else
590 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
591#endif
592 END DO
593!
594!=======================================================================
595! Start internal iterations to achieve convergence of the nonlinear
596! backward-implicit solution.
597!=======================================================================
598!
599 DO iteradj=1,iter
600!
601! Compute light attenuation as function of depth.
602!
603 DO i=istr,iend
604 par=parsur(i)
605 IF (parsur(i).gt.0.0_r8) THEN ! day time
606 DO k=n(ng),1,-1
607!
608! Compute average light attenuation for each grid cell. Here, AttSW is
609! the light attenuation due to seawater and AttPhy is the attenuation
610! due to phytoplankton (self-shading coefficient).
611!
612 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
613 & (z_w(i,j,k)-z_w(i,j,k-1))
614 expatt=exp(-att)
615 itop=par
616 par=itop*(1.0_r8-expatt)/att ! average at cell center
617 light(i,k)=par
618!
619! Light attenuation at the bottom of the grid cell. It is the starting
620! PAR value for the next (deeper) vertical grid cell.
621!
622 par=itop*expatt
623 END DO
624 ELSE ! night time
625 DO k=1,n(ng)
626 light(i,k)=0.0_r8
627 END DO
628 END IF
629 END DO
630!
631! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
632! The Michaelis-Menten curve is used to describe the change in uptake
633! rate as a function of nitrate concentration. Here, PhyIS is the
634! initial slope of the P-I curve and K_NO3 is the half saturation of
635! phytoplankton nitrate uptake.
636#ifdef IRON_LIMIT
637!
638! Growth reduction factors due to iron limitation:
639!
640! FNratio current Fe:N ratio [umol-Fe/mmol-N]
641! FCratio current Fe:C ratio [umol-Fe/mol-C]
642! (umol-Fe/mmol-N)*(16 M-N/106 M-C)*(1E3 mmol-C/mol-C)
643! FCratioE empirical Fe:C ratio
644! Flimit Phytoplankton growth reduction factor due to Fe
645! limitation based on Fe:C ratio
646!
647#endif
648!
649 cff1=dtdays*vm_no3(ng)*phyis(ng)
650 cff2=vm_no3(ng)*vm_no3(ng)
651 cff3=phyis(ng)*phyis(ng)
652 DO k=1,n(ng)
653 DO i=istr,iend
654#ifdef IRON_LIMIT
655!
656! Calculate growth reduction factor due to iron limitation.
657!
658 fnratio=bio(i,k,ifphy)/max(minval,bio(i,k,iphyt))
659 fcratio=fnratio*fen2fec
660 fcratioe=b_fe(ng)*bio(i,k,ifdis)**a_fe(ng)
661 flimit=fcratio*fcratio/ &
662 & (fcratio*fcratio+k_fec(ng)*k_fec(ng))
663
664 nlimit=1.0_r8/(k_no3(ng)+bio(i,k,ino3_))
665 fnlim=min(1.0_r8,flimit/(bio(i,k,ino3_)*nlimit))
666#endif
667 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
668 cff=bio(i,k,iphyt)* &
669#ifdef IRON_LIMIT
670 & cff1*cff4*light(i,k)*fnlim*nlimit
671#else
672 & cff1*cff4*light(i,k)/ &
673 & (k_no3(ng)+bio(i,k,ino3_))
674#endif
675 bio1(i,k,ino3_)=bio(i,k,ino3_)
676 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
677 bio1(i,k,iphyt)=bio(i,k,iphyt)
678 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
679 & bio(i,k,ino3_)*cff
680
681#ifdef IRON_LIMIT
682!
683! Iron uptake proportional to growth.
684!
685 fac=cff*bio(i,k,ino3_)*fnratio/ &
686 & max(minval,bio(i,k,ifdis))
687 bio1(i,k,ifdis)=bio(i,k,ifdis)
688 bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+fac)
689 bio2(i,k,ifdis)=bio(i,k,ifdis)
690 bio1(i,k,ifphy)=bio(i,k,ifphy)
691 bio(i,k,ifphy)=bio(i,k,ifphy)+ &
692 & bio(i,k,ifdis)*fac
693 bio2(i,k,ifphy)=bio(i,k,ifphy)
694!
695! Iron uptake to reach appropriate Fe:C ratio.
696!
697 cff5=dtdays*(fcratioe-fcratio)/t_fe(ng)
698 cff6=bio(i,k,iphyt)*cff5*fec2fen
699 IF (cff6.ge.0.0_r8) THEN
700 cff=cff6/max(minval,bio(i,k,ifdis))
701 bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+cff)
702 bio(i,k,ifphy)=bio(i,k,ifphy)+ &
703 & bio(i,k,ifdis)*cff
704 ELSE
705 cff=-cff6/max(minval,bio(i,k,ifphy))
706 bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
707 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
708 & bio(i,k,ifphy)*cff
709 END IF
710#endif
711 END DO
712 END DO
713!
714 IF (iteradj.ne.iter) THEN
715!
716! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
717! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
718! pool as function of "sloppy feeding" and metabolic processes
719! (ZooEEN and ZooEED fractions).
720#ifdef IRON_LIMIT
721! The lost of phytoplankton to the dissolve iron pool is scale by the
722! remineralization rate (FeRR).
723#endif
724!
725 cff1=dtdays*zoogr(ng)
726 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
727 DO k=1,n(ng)
728 DO i=istr,iend
729 cff=bio(i,k,izoop)* &
730 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
731 & bio(i,k,iphyt)
732 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
733 bio(i,k,izoop)=bio(i,k,izoop)+ &
734 & bio(i,k,iphyt)*cff2*cff
735 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
736 & bio(i,k,iphyt)*zooeen(ng)*cff
737 bio(i,k,isdet)=bio(i,k,isdet)+ &
738 & bio(i,k,iphyt)*zooeed(ng)*cff
739#ifdef IRON_LIMIT
740 bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
741 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
742 & bio(i,k,ifphy)*cff*ferr(ng)
743#endif
744 END DO
745 END DO
746!
747! Phytoplankton mortality to nutrients (PhyMRNro rate), detritus
748! (PhyMRD rate), and if applicable dissolved iron (FeRR rate).
749!
750 cff3=dtdays*phymrd(ng)
751 cff2=dtdays*phymrn(ng)
752 cff1=1.0_r8/(1.0_r8+cff2+cff3)
753 DO k=1,n(ng)
754 DO i=istr,iend
755 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
756 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
757 & bio(i,k,iphyt)*cff2
758 bio(i,k,isdet)=bio(i,k,isdet)+ &
759 & bio(i,k,iphyt)*cff3
760#ifdef IRON_LIMIT
761 bio(i,k,ifphy)=bio(i,k,ifphy)*cff1
762 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
763 & bio(i,k,ifphy)*(cff2+cff3)*ferr(ng)
764#endif
765 END DO
766 END DO
767!
768! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
769! (ZooMRD rate).
770!
771 cff3=dtdays*zoomrd(ng)
772 cff2=dtdays*zoomrn(ng)
773 cff1=1.0_r8/(1.0_r8+cff2+cff3)
774 DO k=1,n(ng)
775 DO i=istr,iend
776 bio(i,k,izoop)=bio(i,k,izoop)*cff1
777 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
778 & bio(i,k,izoop)*cff2
779 bio(i,k,isdet)=bio(i,k,isdet)+ &
780 & bio(i,k,izoop)*cff3
781 END DO
782 END DO
783!
784! Detritus breakdown to nutrients: remineralization (DetRR rate).
785!
786 cff2=dtdays*detrr(ng)
787 cff1=1.0_r8/(1.0_r8+cff2)
788 DO k=1,n(ng)
789 DO i=istr,iend
790 bio(i,k,isdet)=bio(i,k,isdet)*cff1
791 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
792 & bio(i,k,isdet)*cff2
793 END DO
794 END DO
795!
796!-----------------------------------------------------------------------
797! Vertical sinking terms: Phytoplankton and Detritus
798!-----------------------------------------------------------------------
799!
800! Reconstruct vertical profile of selected biological constituents
801! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
802! grid box. Then, compute semi-Lagrangian flux due to sinking.
803!
804 DO isink=1,nsink
805 ibio=idsink(isink)
806!
807! Copy concentration of biological particulates into scratch array
808! "qc" (q-central, restrict it to be positive) which is hereafter
809! interpreted as a set of grid-box averaged values for biogeochemical
810! constituent concentration.
811!
812 DO k=1,n(ng)
813 DO i=istr,iend
814 qc(i,k)=bio(i,k,ibio)
815 END DO
816 END DO
817!
818 DO k=n(ng)-1,1,-1
819 DO i=istr,iend
820 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
821 END DO
822 END DO
823 DO k=2,n(ng)-1
824 DO i=istr,iend
825 dltr=hz(i,j,k)*fc(i,k)
826 dltl=hz(i,j,k)*fc(i,k-1)
827 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
828 cffr=cff*fc(i,k)
829 cffl=cff*fc(i,k-1)
830!
831! Apply PPM monotonicity constraint to prevent oscillations within the
832! grid box.
833!
834 IF ((dltr*dltl).le.0.0_r8) THEN
835 dltr=0.0_r8
836 dltl=0.0_r8
837 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
838 dltr=cffl
839 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
840 dltl=cffr
841 END IF
842!
843! Compute right and left side values (bR,bL) of parabolic segments
844! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
845!
846! NOTE: Although each parabolic segment is monotonic within its grid
847! box, monotonicity of the whole profile is not guaranteed,
848! because bL(k+1)-bR(k) may still have different sign than
849! qc(i,k+1)-qc(i,k). This possibility is excluded,
850! after bL and bR are reconciled using WENO procedure.
851!
852 cff=(dltr-dltl)*hz_inv3(i,k)
853 dltr=dltr-cff*hz(i,j,k+1)
854 dltl=dltl+cff*hz(i,j,k-1)
855 br(i,k)=qc(i,k)+dltr
856 bl(i,k)=qc(i,k)-dltl
857 wr(i,k)=(2.0_r8*dltr-dltl)**2
858 wl(i,k)=(dltr-2.0_r8*dltl)**2
859 END DO
860 END DO
861 cff=1.0e-14_r8
862 DO k=2,n(ng)-2
863 DO i=istr,iend
864 dltl=max(cff,wl(i,k ))
865 dltr=max(cff,wr(i,k+1))
866 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
867 bl(i,k+1)=br(i,k)
868 END DO
869 END DO
870 DO i=istr,iend
871 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
872#if defined LINEAR_CONTINUATION
873 bl(i,n(ng))=br(i,n(ng)-1)
874 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
875#elif defined NEUMANN
876 bl(i,n(ng))=br(i,n(ng)-1)
877 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
878#else
879 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
880 bl(i,n(ng))=qc(i,n(ng)) ! conditions
881 br(i,n(ng)-1)=qc(i,n(ng))
882#endif
883#if defined LINEAR_CONTINUATION
884 br(i,1)=bl(i,2)
885 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
886#elif defined NEUMANN
887 br(i,1)=bl(i,2)
888 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
889#else
890 bl(i,2)=qc(i,1) ! bottom grid boxes are
891 br(i,1)=qc(i,1) ! re-assumed to be
892 bl(i,1)=qc(i,1) ! piecewise constant.
893#endif
894 END DO
895!
896! Apply monotonicity constraint again, since the reconciled interfacial
897! values may cause a non-monotonic behavior of the parabolic segments
898! inside the grid box.
899!
900 DO k=1,n(ng)
901 DO i=istr,iend
902 dltr=br(i,k)-qc(i,k)
903 dltl=qc(i,k)-bl(i,k)
904 cffr=2.0_r8*dltr
905 cffl=2.0_r8*dltl
906 IF ((dltr*dltl).lt.0.0_r8) THEN
907 dltr=0.0_r8
908 dltl=0.0_r8
909 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
910 dltr=cffl
911 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
912 dltl=cffr
913 END IF
914 br(i,k)=qc(i,k)+dltr
915 bl(i,k)=qc(i,k)-dltl
916 END DO
917 END DO
918!
919! After this moment reconstruction is considered complete. The next
920! stage is to compute vertical advective fluxes, FC. It is expected
921! that sinking may occurs relatively fast, the algorithm is designed
922! to be free of CFL criterion, which is achieved by allowing
923! integration bounds for semi-Lagrangian advective flux to use as
924! many grid boxes in upstream direction as necessary.
925!
926! In the two code segments below, WL is the z-coordinate of the
927! departure point for grid box interface z_w with the same indices;
928! FC is the finite volume flux; ksource(:,k) is index of vertical
929! grid box which contains the departure point (restricted by N(ng)).
930! During the search: also add in content of whole grid boxes
931! participating in FC.
932!
933 cff=dtdays*abs(wbio(isink))
934 DO k=1,n(ng)
935 DO i=istr,iend
936 fc(i,k-1)=0.0_r8
937 wl(i,k)=z_w(i,j,k-1)+cff
938 wr(i,k)=hz(i,j,k)*qc(i,k)
939 ksource(i,k)=k
940 END DO
941 END DO
942 DO k=1,n(ng)
943 DO ks=k,n(ng)-1
944 DO i=istr,iend
945 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
946 ksource(i,k)=ks+1
947 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
948 END IF
949 END DO
950 END DO
951 END DO
952!
953! Finalize computation of flux: add fractional part.
954!
955 DO k=1,n(ng)
956 DO i=istr,iend
957 ks=ksource(i,k)
958 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
959 fc(i,k-1)=fc(i,k-1)+ &
960 & hz(i,j,ks)*cu* &
961 & (bl(i,ks)+ &
962 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
963 & (1.5_r8-cu)* &
964 & (br(i,ks)+bl(i,ks)- &
965 & 2.0_r8*qc(i,ks))))
966 END DO
967 END DO
968 DO k=1,n(ng)
969 DO i=istr,iend
970 bio(i,k,ibio)=qc(i,k)+ &
971 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
972 END DO
973 END DO
974 END DO
975 END IF
976 END DO
977!
978! End of compute basic state arrays I.
979!
980! Compute light attenuation as function of depth.
981!
982 DO i=istr,iend
983 par=parsur(i)
984 tl_par=tl_parsur(i)
985 IF (parsur(i).gt.0.0_r8) THEN ! day time
986 DO k=n(ng),1,-1
987!
988! Compute average light attenuation for each grid cell. Here, AttSW is
989! the light attenuation due to seawater and AttPhy is the attenuation
990! due to phytoplankton (self-shading coefficient).
991!
992 att=(attsw(ng)+attphy(ng)*bio1(i,k,iphyt))* &
993 & (z_w(i,j,k)-z_w(i,j,k-1))
994 tl_att=attphy(ng)*tl_bio(i,k,iphyt)* &
995 & (z_w(i,j,k)-z_w(i,j,k-1))+ &
996 & (attsw(ng)+attphy(ng)*bio1(i,k,iphyt))* &
997 & (tl_z_w(i,j,k)-tl_z_w(i,j,k-1))
998 expatt=exp(-att)
999 tl_expatt=-expatt*tl_att
1000 itop=par
1001 tl_itop=tl_par
1002 par=itop*(1.0_r8-expatt)/att ! average at cell center
1003 tl_par=(-tl_att*par+tl_itop*(1.0_r8-expatt)- &
1004 & itop*tl_expatt)/att
1005!^ Light(i,k)=PAR
1006!^
1007 tl_light(i,k)=tl_par
1008!
1009! Light attenuation at the bottom of the grid cell. It is the starting
1010! PAR value for the next (deeper) vertical grid cell.
1011!
1012 par=itop*expatt
1013 tl_par=tl_itop*expatt+itop*tl_expatt
1014 END DO
1015 ELSE ! night time
1016 DO k=1,n(ng)
1017!^ Light(i,k)=0.0_r8
1018!^
1019 tl_light(i,k)=0.0_r8
1020 END DO
1021 END IF
1022 END DO
1023!
1024! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
1025! The Michaelis-Menten curve is used to describe the change in uptake
1026! rate as a function of nitrate concentration. Here, PhyIS is the
1027! initial slope of the P-I curve and K_NO3 is the half saturation of
1028! phytoplankton nitrate uptake.
1029#ifdef IRON_LIMIT
1030!
1031! Growth reduction factors due to iron limitation:
1032!
1033! FNratio current Fe:N ratio [umol-Fe/mmol-N]
1034! FCratio current Fe:C ratio [umol-Fe/mol-C]
1035! (umol-Fe/mmol-N)*(16 M-N/106 M-C)*(1E3 mmol-C/mol-C)
1036! FCratioE empirical Fe:C ratio
1037! Flimit Phytoplankton growth reduction factor due to Fe
1038! limitation based on Fe:C ratio
1039!
1040#endif
1041!
1042 cff1=dtdays*vm_no3(ng)*phyis(ng)
1043 cff2=vm_no3(ng)*vm_no3(ng)
1044 cff3=phyis(ng)*phyis(ng)
1045 DO k=1,n(ng)
1046 DO i=istr,iend
1047#ifdef IRON_LIMIT
1048!
1049! Calculate growth reduction factor due to iron limitation.
1050!
1051!^ FNratio=Bio(i,k,iFphy)/MAX(MinVal,Bio(i,k,iPhyt))
1052!^
1053 fac1=max(minval,bio1(i,k,iphyt))
1054 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio1(i,k,iphyt)))* &
1055 & tl_bio(i,k,iphyt)
1056 fnratio=bio1(i,k,ifphy)/fac1
1057 tl_fnratio=(tl_bio(i,k,ifphy)-tl_fac1*fnratio)/fac1
1058 fcratio=fnratio*fen2fec
1059 tl_fcratio=tl_fnratio*fen2fec
1060 fcratioe=b_fe(ng)*bio1(i,k,ifdis)**a_fe(ng)
1061 tl_fcratioe=a_fe(ng)*b_fe(ng)* &
1062 & bio1(i,k,ifdis)**(a_fe(ng)-1.0_r8)* &
1063 & tl_bio(i,k,ifdis)
1064 flimit=fcratio*fcratio/ &
1065 & (fcratio*fcratio+k_fec(ng)*k_fec(ng))
1066 tl_flimit=2.0_r8*(tl_fcratio*fcratio- &
1067 & tl_fcratio*fcratio*flimit)/ &
1068 & (fcratio*fcratio+k_fec(ng)*k_fec(ng))
1069 nlimit=1.0_r8/(k_no3(ng)+bio1(i,k,ino3_))
1070 tl_nlimit=-tl_bio(i,k,ino3_)*nlimit*nlimit
1071!^ FNlim=MIN(1.0_r8,Flimit/(Bio1(i,k,iNO3_)*Nlimit))
1072!^
1073 fac1=flimit/(bio1(i,k,ino3_)*nlimit)
1074 tl_fac1=tl_flimit/(bio1(i,k,ino3_)*nlimit)- &
1075 & (tl_bio(i,k,ino3_)*nlimit+ &
1076 & bio1(i,k,ino3_)*tl_nlimit)*fac1/ &
1077 & (bio1(i,k,ino3_)*nlimit)
1078 fnlim=min(1.0_r8,fac1)
1079 tl_fnlim=(0.5_r8+sign(0.5_r8,1.0_r8-fac1))*tl_fac1
1080#endif
1081 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1082 tl_cff4=-cff3*tl_light(i,k)*light(i,k)*cff4*cff4*cff4
1083 cff=bio1(i,k,iphyt)* &
1084#ifdef IRON_LIMIT
1085 & cff1*cff4*light(i,k)*fnlim*nlimit
1086#else
1087 & cff1*cff4*light(i,k)/ &
1088 & (k_no3(ng)+bio1(i,k,ino3_))
1089#endif
1090#ifdef IRON_LIMIT
1091 tl_cff=tl_bio(i,k,iphyt)* &
1092 & cff1*cff4*light(i,k)*fnlim*nlimit+ &
1093 & bio1(i,k,iphyt)*cff1*cff4* &
1094 & (tl_light(i,k)*fnlim*nlimit+ &
1095 & light(i,k)*tl_fnlim*nlimit+ &
1096 & light(i,k)*fnlim*tl_nlimit)+ &
1097 & bio1(i,k,iphyt)*cff1*tl_cff4* &
1098 & light(i,k)*fnlim*nlimit
1099#else
1100 tl_cff=(tl_bio(i,k,iphyt)*cff1*cff4*light(i,k)+ &
1101 & bio1(i,k,iphyt)*cff1* &
1102 & (tl_cff4*light(i,k)+cff4*tl_light(i,k))- &
1103 & tl_bio(i,k,ino3_)*cff)/ &
1104 & (k_no3(ng)+bio1(i,k,ino3_))
1105#endif
1106!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)/(1.0_r8+cff)
1107!^
1108 tl_bio(i,k,ino3_)=(tl_bio(i,k,ino3_)- &
1109 & tl_cff*bio(i,k,ino3_))/ &
1110 & (1.0_r8+cff)
1111!^ Bio(i,k,iPhyt)=Bio(i,k,iPhyt)+ &
1112!^ & Bio(i,k,iNO3_)*cff
1113!^
1114 tl_bio(i,k,iphyt)=tl_bio(i,k,iphyt)+ &
1115 & tl_bio(i,k,ino3_)*cff+ &
1116 & bio(i,k,ino3_)*tl_cff
1117
1118#ifdef IRON_LIMIT
1119!
1120! Iron uptake proportional to growth.
1121!
1122!^ fac=cff*Bio(i,k,iNO3_)*FNratio/MAX(MinVal,Bio1(i,k,iFdis))
1123!^
1124 fac1=max(minval,bio1(i,k,ifdis))
1125 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio1(i,k,ifdis)))* &
1126 & tl_bio(i,k,ifdis)
1127 fac2=1.0_r8/fac1
1128 tl_fac2=-fac2*fac2*tl_fac1
1129 fac=cff*bio(i,k,ino3_)*fnratio*fac2
1130 tl_fac=fnratio*fac2*(tl_cff*bio(i,k,ino3_)+ &
1131 & cff*tl_bio(i,k,ino3_))+ &
1132 & cff*bio(i,k,ino3_)*(tl_fnratio*fac2+ &
1133 & fnratio*tl_fac2)
1134!^ Bio(i,k,iFdis)=Bio(i,k,iFdis)/(1.0_r8+fac)
1135!^
1136 tl_bio(i,k,ifdis)=(tl_bio(i,k,ifdis)- &
1137 & tl_fac*bio2(i,k,ifdis))/ &
1138 & (1.0_r8+fac)
1139!^ Bio(i,k,iFphy)=Bio(i,k,iFphy)+ &
1140!^ & Bio(i,k,iFdis)*fac
1141!^
1142 tl_bio(i,k,ifphy)=tl_bio(i,k,ifphy)+ &
1143 & tl_bio(i,k,ifdis)*fac+ &
1144 & bio2(i,k,ifdis)*tl_fac
1145!
1146! Iron uptake to reach appropriate Fe:C ratio.
1147!
1148 cff5=dtdays*(fcratioe-fcratio)/t_fe(ng)
1149 tl_cff5=dtdays*(tl_fcratioe-tl_fcratio)/t_fe(ng)
1150 cff6=bio(i,k,iphyt)*cff5*fec2fen
1151 tl_cff6=(tl_bio(i,k,iphyt)*cff5+ &
1152 & bio(i,k,iphyt)*tl_cff5)*fec2fen
1153 IF (cff6.ge.0.0_r8) THEN
1154!^ cff=cff6/MAX(MinVal,Bio2(i,k,iFdis))
1155!^
1156 fac1=max(minval,bio2(i,k,ifdis))
1157 tl_fac1=(0.5_r8-sign(0.5_r8,minval-bio2(i,k,ifdis)))* &
1158 & tl_bio(i,k,ifdis)
1159 cff=cff6/fac1
1160 tl_cff=(tl_cff6-tl_fac1*cff)/fac1
1161!^ Bio(i,k,iFdis)=Bio(i,k,iFdis)/(1.0_r8+cff)
1162!^
1163 tl_bio(i,k,ifdis)=(tl_bio(i,k,ifdis)- &
1164 & tl_cff*bio(i,k,ifdis))/ &
1165 & (1.0_r8+cff)
1166!^ Bio(i,k,iFphy)=Bio(i,k,iFphy)+ &
1167!^ & Bio(i,k,iFdis)*cff
1168!^
1169 tl_bio(i,k,ifphy)=tl_bio(i,k,ifphy)+ &
1170 & tl_bio(i,k,ifdis)*cff+ &
1171 & bio(i,k,ifdis)*tl_cff
1172 ELSE
1173!^ cff=-cff6/MAX(MinVal,Bio2(i,k,iFphy))
1174!^
1175 fac1=-max(minval,bio2(i,k,ifphy))
1176 tl_fac1=-(0.5_r8-sign(0.5_r8,minval-bio2(i,k,ifphy)))* &
1177 & tl_bio(i,k,ifphy)
1178 cff=cff6/fac1
1179 tl_cff=(tl_cff6-tl_fac1*cff)/fac1
1180!^ Bio(i,k,iFphy)=Bio(i,k,iFphy)/(1.0_r8+cff)
1181!^
1182 tl_bio(i,k,ifphy)=(tl_bio(i,k,ifphy)- &
1183 & tl_cff*bio(i,k,ifphy))/ &
1184 & (1.0_r8+cff)
1185!^ Bio(i,k,iFdis)=Bio(i,k,iFdis)+ &
1186!^ & Bio(i,k,iFphy)*cff
1187!^
1188 tl_bio(i,k,ifdis)=tl_bio(i,k,ifdis)+ &
1189 & tl_bio(i,k,ifphy)*cff+ &
1190 & bio(i,k,ifphy)*tl_cff
1191 END IF
1192#endif
1193 END DO
1194 END DO
1195!
1196! Compute appropriate basic state arrays II.
1197!
1198 DO k=1,n(ng)
1199 DO i=istr,iend
1200!
1201! At input, all tracers (index nnew) from predictor step have
1202! transport units (m Tunits) since we do not have yet the new
1203! values for zeta and Hz. These are known after the 2D barotropic
1204! time-stepping.
1205!
1206! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
1207! tracer times depth. However the basic state (nstp and nnew
1208! indices) that is read from the forward file is in units of
1209! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
1210! use t instead of t*Hz_inv.
1211!
1212 DO itrc=1,nbt
1213 ibio=idbio(itrc)
1214!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1215!^
1216 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1217!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
1218!^
1219 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
1220 END DO
1221!
1222! Impose positive definite concentrations.
1223!
1224 cff2=0.0_r8
1225 DO itime=1,2
1226 cff1=0.0_r8
1227 itrcmax=idbio(1)
1228#ifdef IRON_LIMIT
1229 DO itrc=1,nbt-2
1230#else
1231 DO itrc=1,nbt
1232#endif
1233 ibio=idbio(itrc)
1234 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
1235 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
1236 itrcmax=ibio
1237 END IF
1238 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1239 END DO
1240 IF (biotrc(itrcmax,itime).gt.cff1) THEN
1241 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
1242 END IF
1243#ifdef IRON_LIMIT
1244 DO itrc=nbt-1,nbt
1245 ibio=idbio(itrc)
1246 biotrc1(ibio,itime)=biotrc(ibio,itime)
1247 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1248 END DO
1249#endif
1250 END DO
1251!
1252! Load biological tracers into local arrays.
1253!
1254 DO itrc=1,nbt
1255 ibio=idbio(itrc)
1256 bio_old(i,k,ibio)=biotrc(ibio,nstp)
1257 bio(i,k,ibio)=biotrc(ibio,nstp)
1258 END DO
1259
1260#if defined IRON_LIMIT && defined IRON_RELAX
1261!
1262! Relax dissolved iron at coast (h <= FeHim) to a constant value
1263! (FeMax) over a time scale (FeNudgTime; days) to simulate sources
1264! at the shelf.
1265!
1266 IF (h(i,j).le.fehmin(ng)) THEN
1267 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
1268 & fenudgcoef*(femax(ng)-bio(i,k,ifdis))
1269 END IF
1270#endif
1271 END DO
1272 END DO
1273!
1274! Calculate surface Photosynthetically Available Radiation (PAR). The
1275! net shortwave radiation is scaled back to Watts/m2 and multiplied by
1276! the fraction that is photosynthetically available, PARfrac.
1277!
1278 DO i=istr,iend
1279#ifdef CONST_PAR
1280!
1281! Specify constant surface irradiance a la Powell and Spitz.
1282!
1283 parsur(i)=158.075_r8
1284#else
1285 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
1286#endif
1287 END DO
1288!
1289!=======================================================================
1290! Start internal iterations to achieve convergence of the nonlinear
1291! backward-implicit solution.
1292!=======================================================================
1293!
1294 DO iteradj=1,iter
1295!
1296! Compute light attenuation as function of depth.
1297!
1298 DO i=istr,iend
1299 par=parsur(i)
1300 IF (parsur(i).gt.0.0_r8) THEN ! day time
1301 DO k=n(ng),1,-1
1302!
1303! Compute average light attenuation for each grid cell. Here, AttSW is
1304! the light attenuation due to seawater and AttPhy is the attenuation
1305! due to phytoplankton (self-shading coefficient).
1306!
1307 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
1308 & (z_w(i,j,k)-z_w(i,j,k-1))
1309 expatt=exp(-att)
1310 itop=par
1311 par=itop*(1.0_r8-expatt)/att ! average at cell center
1312 light(i,k)=par
1313!
1314! Light attenuation at the bottom of the grid cell. It is the starting
1315! PAR value for the next (deeper) vertical grid cell.
1316!
1317 par=itop*expatt
1318 END DO
1319 ELSE ! night time
1320 DO k=1,n(ng)
1321 light(i,k)=0.0_r8
1322 END DO
1323 END IF
1324 END DO
1325!
1326! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
1327! The Michaelis-Menten curve is used to describe the change in uptake
1328! rate as a function of nitrate concentration. Here, PhyIS is the
1329! initial slope of the P-I curve and K_NO3 is the half saturation of
1330! phytoplankton nitrate uptake.
1331#ifdef IRON_LIMIT
1332!
1333! Growth reduction factors due to iron limitation:
1334!
1335! FNratio current Fe:N ratio [umol-Fe/mmol-N]
1336! FCratio current Fe:C ratio [umol-Fe/mol-C]
1337! (umol-Fe/mmol-N)*(16 M-N/106 M-C)*(1E3 mmol-C/mol-C)
1338! FCratioE empirical Fe:C ratio
1339! Flimit Phytoplankton growth reduction factor due to Fe
1340! limitation based on Fe:C ratio
1341!
1342#endif
1343!
1344 cff1=dtdays*vm_no3(ng)*phyis(ng)
1345 cff2=vm_no3(ng)*vm_no3(ng)
1346 cff3=phyis(ng)*phyis(ng)
1347 DO k=1,n(ng)
1348 DO i=istr,iend
1349#ifdef IRON_LIMIT
1350!
1351! Calculate growth reduction factor due to iron limitation.
1352!
1353 fnratio=bio(i,k,ifphy)/max(minval,bio(i,k,iphyt))
1354 fcratio=fnratio*fen2fec
1355 fcratioe=b_fe(ng)*bio(i,k,ifdis)**a_fe(ng)
1356 flimit=fcratio*fcratio/ &
1357 & (fcratio*fcratio+k_fec(ng)*k_fec(ng))
1358
1359 nlimit=1.0_r8/(k_no3(ng)+bio(i,k,ino3_))
1360 fnlim=min(1.0_r8,flimit/(bio(i,k,ino3_)*nlimit))
1361#endif
1362 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1363 cff=bio(i,k,iphyt)* &
1364#ifdef IRON_LIMIT
1365 & cff1*cff4*light(i,k)*fnlim*nlimit
1366#else
1367 & cff1*cff4*light(i,k)/ &
1368 & (k_no3(ng)+bio(i,k,ino3_))
1369#endif
1370 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
1371 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
1372 & bio(i,k,ino3_)*cff
1373
1374#ifdef IRON_LIMIT
1375!
1376! Iron uptake proportional to growth.
1377!
1378 fac=cff*bio(i,k,ino3_)*fnratio/ &
1379 & max(minval,bio(i,k,ifdis))
1380 bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+fac)
1381 bio(i,k,ifphy)=bio(i,k,ifphy)+ &
1382 & bio(i,k,ifdis)*fac
1383!
1384! Iron uptake to reach appropriate Fe:C ratio.
1385!
1386 cff5=dtdays*(fcratioe-fcratio)/t_fe(ng)
1387 cff6=bio(i,k,iphyt)*cff5*fec2fen
1388 IF (cff6.ge.0.0_r8) THEN
1389 cff=cff6/max(minval,bio(i,k,ifdis))
1390 bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+cff)
1391 bio(i,k,ifphy)=bio(i,k,ifphy)+ &
1392 & bio(i,k,ifdis)*cff
1393 ELSE
1394 cff=-cff6/max(minval,bio(i,k,ifphy))
1395 bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
1396 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
1397 & bio(i,k,ifphy)*cff
1398 END IF
1399#endif
1400 END DO
1401 END DO
1402!
1403! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
1404! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
1405! pool as function of "sloppy feeding" and metabolic processes
1406! (ZooEEN and ZooEED fractions).
1407#ifdef IRON_LIMIT
1408! The lost of phytoplankton to the dissolve iron pool is scale by the
1409! remineralization rate (FeRR).
1410#endif
1411!
1412 cff1=dtdays*zoogr(ng)
1413 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
1414 DO k=1,n(ng)
1415 DO i=istr,iend
1416 cff=bio(i,k,izoop)* &
1417 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
1418 & bio(i,k,iphyt)
1419 bio1(i,k,iphyt)=bio(i,k,iphyt)
1420 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
1421 bio1(i,k,izoop)=bio(i,k,izoop)
1422 bio(i,k,izoop)=bio(i,k,izoop)+ &
1423 & bio(i,k,iphyt)*cff2*cff
1424 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1425 & bio(i,k,iphyt)*zooeen(ng)*cff
1426 bio(i,k,isdet)=bio(i,k,isdet)+ &
1427 & bio(i,k,iphyt)*zooeed(ng)*cff
1428#ifdef IRON_LIMIT
1429 bio1(i,k,ifphy)=bio(i,k,ifphy)
1430 bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
1431 bio2(i,k,ifphy)=bio(i,k,ifphy)
1432 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
1433 & bio(i,k,ifphy)*cff*ferr(ng)
1434#endif
1435 END DO
1436 END DO
1437!
1438! Phytoplankton mortality to nutrients (PhyMRNro rate), detritus
1439! (PhyMRD rate), and if applicable dissolved iron (FeRR rate).
1440!
1441 cff3=dtdays*phymrd(ng)
1442 cff2=dtdays*phymrn(ng)
1443 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1444 DO k=1,n(ng)
1445 DO i=istr,iend
1446 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
1447 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1448 & bio(i,k,iphyt)*cff2
1449 bio(i,k,isdet)=bio(i,k,isdet)+ &
1450 & bio(i,k,iphyt)*cff3
1451#ifdef IRON_LIMIT
1452 bio(i,k,ifphy)=bio(i,k,ifphy)*cff1
1453 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
1454 & bio(i,k,ifphy)*(cff2+cff3)*ferr(ng)
1455#endif
1456 END DO
1457 END DO
1458!
1459 IF (iteradj.ne.iter) THEN
1460!
1461! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
1462! (ZooMRD rate).
1463!
1464 cff3=dtdays*zoomrd(ng)
1465 cff2=dtdays*zoomrn(ng)
1466 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1467 DO k=1,n(ng)
1468 DO i=istr,iend
1469 bio(i,k,izoop)=bio(i,k,izoop)*cff1
1470 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1471 & bio(i,k,izoop)*cff2
1472 bio(i,k,isdet)=bio(i,k,isdet)+ &
1473 & bio(i,k,izoop)*cff3
1474 END DO
1475 END DO
1476!
1477! Detritus breakdown to nutrients: remineralization (DetRR rate).
1478!
1479 cff2=dtdays*detrr(ng)
1480 cff1=1.0_r8/(1.0_r8+cff2)
1481 DO k=1,n(ng)
1482 DO i=istr,iend
1483 bio(i,k,isdet)=bio(i,k,isdet)*cff1
1484 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1485 & bio(i,k,isdet)*cff2
1486 END DO
1487 END DO
1488!
1489!-----------------------------------------------------------------------
1490! Vertical sinking terms: Phytoplankton and Detritus
1491!-----------------------------------------------------------------------
1492!
1493! Reconstruct vertical profile of selected biological constituents
1494! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
1495! grid box. Then, compute semi-Lagrangian flux due to sinking.
1496!
1497 DO isink=1,nsink
1498 ibio=idsink(isink)
1499!
1500! Copy concentration of biological particulates into scratch array
1501! "qc" (q-central, restrict it to be positive) which is hereafter
1502! interpreted as a set of grid-box averaged values for biogeochemical
1503! constituent concentration.
1504!
1505 DO k=1,n(ng)
1506 DO i=istr,iend
1507 qc(i,k)=bio(i,k,ibio)
1508 END DO
1509 END DO
1510!
1511 DO k=n(ng)-1,1,-1
1512 DO i=istr,iend
1513 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1514 END DO
1515 END DO
1516 DO k=2,n(ng)-1
1517 DO i=istr,iend
1518 dltr=hz(i,j,k)*fc(i,k)
1519 dltl=hz(i,j,k)*fc(i,k-1)
1520 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1521 cffr=cff*fc(i,k)
1522 cffl=cff*fc(i,k-1)
1523!
1524! Apply PPM monotonicity constraint to prevent oscillations within the
1525! grid box.
1526!
1527 IF ((dltr*dltl).le.0.0_r8) THEN
1528 dltr=0.0_r8
1529 dltl=0.0_r8
1530 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1531 dltr=cffl
1532 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1533 dltl=cffr
1534 END IF
1535!
1536! Compute right and left side values (bR,bL) of parabolic segments
1537! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1538!
1539! NOTE: Although each parabolic segment is monotonic within its grid
1540! box, monotonicity of the whole profile is not guaranteed,
1541! because bL(k+1)-bR(k) may still have different sign than
1542! qc(i,k+1)-qc(i,k). This possibility is excluded,
1543! after bL and bR are reconciled using WENO procedure.
1544!
1545 cff=(dltr-dltl)*hz_inv3(i,k)
1546 dltr=dltr-cff*hz(i,j,k+1)
1547 dltl=dltl+cff*hz(i,j,k-1)
1548 br(i,k)=qc(i,k)+dltr
1549 bl(i,k)=qc(i,k)-dltl
1550 wr(i,k)=(2.0_r8*dltr-dltl)**2
1551 wl(i,k)=(dltr-2.0_r8*dltl)**2
1552 END DO
1553 END DO
1554 cff=1.0e-14_r8
1555 DO k=2,n(ng)-2
1556 DO i=istr,iend
1557 dltl=max(cff,wl(i,k ))
1558 dltr=max(cff,wr(i,k+1))
1559 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1560 bl(i,k+1)=br(i,k)
1561 END DO
1562 END DO
1563 DO i=istr,iend
1564 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1565#if defined LINEAR_CONTINUATION
1566 bl(i,n(ng))=br(i,n(ng)-1)
1567 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1568#elif defined NEUMANN
1569 bl(i,n(ng))=br(i,n(ng)-1)
1570 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1571#else
1572 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1573 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1574 br(i,n(ng)-1)=qc(i,n(ng))
1575#endif
1576#if defined LINEAR_CONTINUATION
1577 br(i,1)=bl(i,2)
1578 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1579#elif defined NEUMANN
1580 br(i,1)=bl(i,2)
1581 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1582#else
1583 bl(i,2)=qc(i,1) ! bottom grid boxes are
1584 br(i,1)=qc(i,1) ! re-assumed to be
1585 bl(i,1)=qc(i,1) ! piecewise constant.
1586#endif
1587 END DO
1588!
1589! Apply monotonicity constraint again, since the reconciled interfacial
1590! values may cause a non-monotonic behavior of the parabolic segments
1591! inside the grid box.
1592!
1593 DO k=1,n(ng)
1594 DO i=istr,iend
1595 dltr=br(i,k)-qc(i,k)
1596 dltl=qc(i,k)-bl(i,k)
1597 cffr=2.0_r8*dltr
1598 cffl=2.0_r8*dltl
1599 IF ((dltr*dltl).lt.0.0_r8) THEN
1600 dltr=0.0_r8
1601 dltl=0.0_r8
1602 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1603 dltr=cffl
1604 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1605 dltl=cffr
1606 END IF
1607 br(i,k)=qc(i,k)+dltr
1608 bl(i,k)=qc(i,k)-dltl
1609 END DO
1610 END DO
1611!
1612! After this moment reconstruction is considered complete. The next
1613! stage is to compute vertical advective fluxes, FC. It is expected
1614! that sinking may occurs relatively fast, the algorithm is designed
1615! to be free of CFL criterion, which is achieved by allowing
1616! integration bounds for semi-Lagrangian advective flux to use as
1617! many grid boxes in upstream direction as necessary.
1618!
1619! In the two code segments below, WL is the z-coordinate of the
1620! departure point for grid box interface z_w with the same indices;
1621! FC is the finite volume flux; ksource(:,k) is index of vertical
1622! grid box which contains the departure point (restricted by N(ng)).
1623! During the search: also add in content of whole grid boxes
1624! participating in FC.
1625!
1626 cff=dtdays*abs(wbio(isink))
1627 DO k=1,n(ng)
1628 DO i=istr,iend
1629 fc(i,k-1)=0.0_r8
1630 wl(i,k)=z_w(i,j,k-1)+cff
1631 wr(i,k)=hz(i,j,k)*qc(i,k)
1632 ksource(i,k)=k
1633 END DO
1634 END DO
1635 DO k=1,n(ng)
1636 DO ks=k,n(ng)-1
1637 DO i=istr,iend
1638 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1639 ksource(i,k)=ks+1
1640 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1641 END IF
1642 END DO
1643 END DO
1644 END DO
1645!
1646! Finalize computation of flux: add fractional part.
1647!
1648 DO k=1,n(ng)
1649 DO i=istr,iend
1650 ks=ksource(i,k)
1651 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1652 fc(i,k-1)=fc(i,k-1)+ &
1653 & hz(i,j,ks)*cu* &
1654 & (bl(i,ks)+ &
1655 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1656 & (1.5_r8-cu)* &
1657 & (br(i,ks)+bl(i,ks)- &
1658 & 2.0_r8*qc(i,ks))))
1659 END DO
1660 END DO
1661 DO k=1,n(ng)
1662 DO i=istr,iend
1663 bio(i,k,ibio)=qc(i,k)+ &
1664 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1665 END DO
1666 END DO
1667 END DO
1668 END IF
1669 END DO
1670!
1671! End of compute basic state arrays II.
1672!
1673! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
1674! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
1675! pool as function of "sloppy feeding" and metabolic processes
1676! (ZooEEN and ZooEED fractions).
1677#ifdef IRON_LIMIT
1678! The lost of phytoplankton to the dissolve iron pool is scale by the
1679! remineralization rate (FeRR).
1680#endif
1681!
1682 cff1=dtdays*zoogr(ng)
1683 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
1684 DO k=1,n(ng)
1685 DO i=istr,iend
1686 cff=bio1(i,k,izoop)* &
1687 & cff1*(1.0_r8-exp(-ivlev(ng)*bio1(i,k,iphyt)))/ &
1688 & bio1(i,k,iphyt)
1689 tl_cff=(tl_bio(i,k,izoop)* &
1690 & cff1*(1.0_r8-exp(-ivlev(ng)*bio1(i,k,iphyt)))+ &
1691 & bio1(i,k,izoop)*ivlev(ng)*tl_bio(i,k,iphyt)*cff1* &
1692 & exp(-ivlev(ng)*bio1(i,k,iphyt))- &
1693 & tl_bio(i,k,iphyt)*cff)/ &
1694 & bio1(i,k,iphyt)
1695!^ Bio(i,k,iPhyt)=Bio(i,k,iPhyt)/(1.0_r8+cff)
1696!^
1697 tl_bio(i,k,iphyt)=(tl_bio(i,k,iphyt)- &
1698 & tl_cff*bio(i,k,iphyt))/ &
1699 & (1.0_r8+cff)
1700!^ Bio(i,k,iZoop)=Bio(i,k,iZoop)+ &
1701!^ & Bio(i,k,iPhyt)*cff2*cff
1702!^
1703 tl_bio(i,k,izoop)=tl_bio(i,k,izoop)+ &
1704 & cff2*(tl_bio(i,k,iphyt)*cff+ &
1705 & bio(i,k,iphyt)*tl_cff)
1706!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
1707!^ & Bio(i,k,iPhyt)*ZooEEN(ng)*cff
1708!^
1709 tl_bio(i,k,ino3_)=tl_bio(i,k,ino3_)+ &
1710 & zooeen(ng)*(tl_bio(i,k,iphyt)*cff+ &
1711 & bio(i,k,iphyt)*tl_cff)
1712!^ Bio(i,k,iSDet)=Bio(i,k,iSDet)+ &
1713!^ & Bio(i,k,iPhyt)*ZooEED(ng)*cff
1714!^
1715 tl_bio(i,k,isdet)=tl_bio(i,k,isdet)+ &
1716 & zooeed(ng)*(tl_bio(i,k,iphyt)*cff+ &
1717 & bio(i,k,iphyt)*tl_cff)
1718
1719#ifdef IRON_LIMIT
1720!^ Bio(i,k,iFphy)=Bio(i,k,iFphy)/(1.0_r8+cff)
1721!^
1722 tl_bio(i,k,ifphy)=(tl_bio(i,k,ifphy)- &
1723 & tl_cff*bio2(i,k,ifphy))/ &
1724 & (1.0_r8+cff)
1725!^ Bio(i,k,iFdis)=Bio(i,k,iFdis)+ &
1726!^ & Bio(i,k,iFphy)*cff*FeRR(ng)
1727!^
1728 tl_bio(i,k,ifdis)=tl_bio(i,k,ifdis)+ &
1729 & (tl_bio(i,k,ifphy)*cff+ &
1730 & bio2(i,k,ifphy)*tl_cff)*ferr(ng)
1731#endif
1732 END DO
1733 END DO
1734!
1735! Phytoplankton mortality to nutrients (PhyMRNro rate), detritus
1736! (PhyMRD rate), and if applicable dissolved iron (FeRR rate).
1737!
1738 cff3=dtdays*phymrd(ng)
1739 cff2=dtdays*phymrn(ng)
1740 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1741 DO k=1,n(ng)
1742 DO i=istr,iend
1743!^ Bio(i,k,iPhyt)=Bio(i,k,iPhyt)*cff1
1744!^
1745 tl_bio(i,k,iphyt)=tl_bio(i,k,iphyt)*cff1
1746!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
1747!^ & Bio(i,k,iPhyt)*cff2
1748!^
1749 tl_bio(i,k,ino3_)=tl_bio(i,k,ino3_)+ &
1750 & tl_bio(i,k,iphyt)*cff2
1751!^ Bio(i,k,iSDet)=Bio(i,k,iSDet)+ &
1752!^ & Bio(i,k,iPhyt)*cff3
1753!^
1754 tl_bio(i,k,isdet)=tl_bio(i,k,isdet)+ &
1755 & tl_bio(i,k,iphyt)*cff3
1756
1757#ifdef IRON_LIMIT
1758!^ Bio(i,k,iFphy)=Bio(i,k,iFphy)*cff1
1759!^
1760 tl_bio(i,k,ifphy)=tl_bio(i,k,ifphy)*cff1
1761!^ Bio(i,k,iFdis)=Bio(i,k,iFdis)+ &
1762!^ & Bio(i,k,iFphy)*(cff2+cff3)*FeRR(ng)
1763!^
1764 tl_bio(i,k,ifdis)=tl_bio(i,k,ifdis)+ &
1765 & tl_bio(i,k,ifphy)*(cff2+cff3)*ferr(ng)
1766#endif
1767 END DO
1768 END DO
1769!
1770! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
1771! (ZooMRD rate).
1772!
1773 cff3=dtdays*zoomrd(ng)
1774 cff2=dtdays*zoomrn(ng)
1775 cff1=1.0_r8/(1.0_r8+cff2+cff3)
1776 DO k=1,n(ng)
1777 DO i=istr,iend
1778!^ Bio(i,k,iZoop)=Bio(i,k,iZoop)*cff1
1779!^
1780 tl_bio(i,k,izoop)=tl_bio(i,k,izoop)*cff1
1781!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
1782!^ & Bio(i,k,iZoop)*cff2
1783!^
1784 tl_bio(i,k,ino3_)=tl_bio(i,k,ino3_)+ &
1785 & tl_bio(i,k,izoop)*cff2
1786!^ Bio(i,k,iSDet)=Bio(i,k,iSDet)+ &
1787!^ & Bio(i,k,iZoop)*cff3
1788!^
1789 tl_bio(i,k,isdet)=tl_bio(i,k,isdet)+ &
1790 & tl_bio(i,k,izoop)*cff3
1791 END DO
1792 END DO
1793!
1794! Detritus breakdown to nutrients: remineralization (DetRR rate).
1795!
1796 cff2=dtdays*detrr(ng)
1797 cff1=1.0_r8/(1.0_r8+cff2)
1798 DO k=1,n(ng)
1799 DO i=istr,iend
1800!^ Bio(i,k,iSDet)=Bio(i,k,iSDet)*cff1
1801!^
1802 tl_bio(i,k,isdet)=tl_bio(i,k,isdet)*cff1
1803!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
1804!^ & Bio(i,k,iSDet)*cff2
1805!^
1806 tl_bio(i,k,ino3_)=tl_bio(i,k,ino3_)+ &
1807 & tl_bio(i,k,isdet)*cff2
1808 END DO
1809 END DO
1810!
1811! Compute appropriate basic state arrays III.
1812!
1813 DO k=1,n(ng)
1814 DO i=istr,iend
1815!
1816! At input, all tracers (index nnew) from predictor step have
1817! transport units (m Tunits) since we do not have yet the new
1818! values for zeta and Hz. These are known after the 2D barotropic
1819! time-stepping.
1820!
1821! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
1822! tracer times depth. However the basic state (nstp and nnew
1823! indices) that is read from the forward file is in units of
1824! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
1825! use t instead of t*Hz_inv.
1826!
1827 DO itrc=1,nbt
1828 ibio=idbio(itrc)
1829!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1830!^
1831 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1832!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
1833!^
1834 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
1835 END DO
1836!
1837! Impose positive definite concentrations.
1838!
1839 cff2=0.0_r8
1840 DO itime=1,2
1841 cff1=0.0_r8
1842 itrcmax=idbio(1)
1843#ifdef IRON_LIMIT
1844 DO itrc=1,nbt-2
1845#else
1846 DO itrc=1,nbt
1847#endif
1848 ibio=idbio(itrc)
1849 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
1850 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
1851 itrcmax=ibio
1852 END IF
1853 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1854 END DO
1855 IF (biotrc(itrcmax,itime).gt.cff1) THEN
1856 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
1857 END IF
1858#ifdef IRON_LIMIT
1859 DO itrc=nbt-1,nbt
1860 ibio=idbio(itrc)
1861 biotrc1(ibio,itime)=biotrc(ibio,itime)
1862 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
1863 END DO
1864#endif
1865 END DO
1866!
1867! Load biological tracers into local arrays.
1868!
1869 DO itrc=1,nbt
1870 ibio=idbio(itrc)
1871 bio_old(i,k,ibio)=biotrc(ibio,nstp)
1872 bio(i,k,ibio)=biotrc(ibio,nstp)
1873 END DO
1874
1875#if defined IRON_LIMIT && defined IRON_RELAX
1876!
1877! Relax dissolved iron at coast (h <= FeHim) to a constant value
1878! (FeMax) over a time scale (FeNudgTime; days) to simulate sources
1879! at the shelf.
1880!
1881 IF (h(i,j).le.fehmin(ng)) THEN
1882 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
1883 & fenudgcoef*(femax(ng)-bio(i,k,ifdis))
1884 END IF
1885#endif
1886 END DO
1887 END DO
1888!
1889! Calculate surface Photosynthetically Available Radiation (PAR). The
1890! net shortwave radiation is scaled back to Watts/m2 and multiplied by
1891! the fraction that is photosynthetically available, PARfrac.
1892!
1893 DO i=istr,iend
1894#ifdef CONST_PAR
1895!
1896! Specify constant surface irradiance a la Powell and Spitz.
1897!
1898 parsur(i)=158.075_r8
1899#else
1900 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
1901#endif
1902 END DO
1903!
1904!=======================================================================
1905! Start internal iterations to achieve convergence of the nonlinear
1906! backward-implicit solution.
1907!=======================================================================
1908!
1909 DO iteradj=1,iter
1910!
1911! Compute light attenuation as function of depth.
1912!
1913 DO i=istr,iend
1914 par=parsur(i)
1915 IF (parsur(i).gt.0.0_r8) THEN ! day time
1916 DO k=n(ng),1,-1
1917!
1918! Compute average light attenuation for each grid cell. Here, AttSW is
1919! the light attenuation due to seawater and AttPhy is the attenuation
1920! due to phytoplankton (self-shading coefficient).
1921!
1922 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
1923 & (z_w(i,j,k)-z_w(i,j,k-1))
1924 expatt=exp(-att)
1925 itop=par
1926 par=itop*(1.0_r8-expatt)/att ! average at cell center
1927 light(i,k)=par
1928!
1929! Light attenuation at the bottom of the grid cell. It is the starting
1930! PAR value for the next (deeper) vertical grid cell.
1931!
1932 par=itop*expatt
1933 END DO
1934 ELSE ! night time
1935 DO k=1,n(ng)
1936 light(i,k)=0.0_r8
1937 END DO
1938 END IF
1939 END DO
1940!
1941! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
1942! The Michaelis-Menten curve is used to describe the change in uptake
1943! rate as a function of nitrate concentration. Here, PhyIS is the
1944! initial slope of the P-I curve and K_NO3 is the half saturation of
1945! phytoplankton nitrate uptake.
1946#ifdef IRON_LIMIT
1947!
1948! Growth reduction factors due to iron limitation:
1949!
1950! FNratio current Fe:N ratio [umol-Fe/mmol-N]
1951! FCratio current Fe:C ratio [umol-Fe/mol-C]
1952! (umol-Fe/mmol-N)*(16 M-N/106 M-C)*(1E3 mmol-C/mol-C)
1953! FCratioE empirical Fe:C ratio
1954! Flimit Phytoplankton growth reduction factor due to Fe
1955! limitation based on Fe:C ratio
1956!
1957#endif
1958!
1959 cff1=dtdays*vm_no3(ng)*phyis(ng)
1960 cff2=vm_no3(ng)*vm_no3(ng)
1961 cff3=phyis(ng)*phyis(ng)
1962 DO k=1,n(ng)
1963 DO i=istr,iend
1964#ifdef IRON_LIMIT
1965 fnratio=bio(i,k,ifphy)/max(minval,bio(i,k,iphyt))
1966 fcratio=fnratio*fen2fec
1967 fcratioe=b_fe(ng)*bio(i,k,ifdis)**a_fe(ng)
1968 flimit=fcratio*fcratio/ &
1969 & (fcratio*fcratio+k_fec(ng)*k_fec(ng))
1970
1971 nlimit=1.0_r8/(k_no3(ng)+bio(i,k,ino3_))
1972 fnlim=min(1.0_r8,flimit/(bio(i,k,ino3_)*nlimit))
1973#endif
1974 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
1975 cff=bio(i,k,iphyt)* &
1976#ifdef IRON_LIMIT
1977 & cff1*cff4*light(i,k)*fnlim*nlimit
1978#else
1979 & cff1*cff4*light(i,k)/ &
1980 & (k_no3(ng)+bio(i,k,ino3_))
1981#endif
1982 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
1983 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
1984 & bio(i,k,ino3_)*cff
1985
1986#ifdef IRON_LIMIT
1987!
1988! Iron uptake proportional to growth.
1989!
1990 fac=cff*bio(i,k,ino3_)*fnratio/ &
1991 & max(minval,bio(i,k,ifdis))
1992 bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+fac)
1993 bio(i,k,ifphy)=bio(i,k,ifphy)+ &
1994 & bio(i,k,ifdis)*fac
1995!
1996! Iron uptake to reach appropriate Fe:C ratio.
1997!
1998 cff5=dtdays*(fcratioe-fcratio)/t_fe(ng)
1999 cff6=bio(i,k,iphyt)*cff5*fec2fen
2000 IF (cff6.ge.0.0_r8) THEN
2001 cff=cff6/max(minval,bio(i,k,ifdis))
2002 bio(i,k,ifdis)=bio(i,k,ifdis)/(1.0_r8+cff)
2003 bio(i,k,ifphy)=bio(i,k,ifphy)+ &
2004 & bio(i,k,ifdis)*cff
2005 ELSE
2006 cff=-cff6/max(minval,bio(i,k,ifphy))
2007 bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
2008 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
2009 & bio(i,k,ifphy)*cff
2010 END IF
2011#endif
2012 END DO
2013 END DO
2014!
2015! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
2016! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
2017! pool as function of "sloppy feeding" and metabolic processes
2018! (ZooEEN and ZooEED fractions).
2019#ifdef IRON_LIMIT
2020! The lost of phytoplankton to the dissolve iron pool is scale by the
2021! remineralization rate (FeRR).
2022#endif
2023!
2024 cff1=dtdays*zoogr(ng)
2025 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
2026 DO k=1,n(ng)
2027 DO i=istr,iend
2028 cff=bio(i,k,izoop)* &
2029 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
2030 & bio(i,k,iphyt)
2031 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
2032 bio(i,k,izoop)=bio(i,k,izoop)+ &
2033 & bio(i,k,iphyt)*cff2*cff
2034 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2035 & bio(i,k,iphyt)*zooeen(ng)*cff
2036 bio(i,k,isdet)=bio(i,k,isdet)+ &
2037 & bio(i,k,iphyt)*zooeed(ng)*cff
2038#ifdef IRON_LIMIT
2039 bio(i,k,ifphy)=bio(i,k,ifphy)/(1.0_r8+cff)
2040 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
2041 & bio(i,k,ifphy)*cff*ferr(ng)
2042#endif
2043 END DO
2044 END DO
2045!
2046! Phytoplankton mortality to nutrients (PhyMRNro rate), detritus
2047! (PhyMRD rate), and if applicable dissolved iron (FeRR rate).
2048!
2049 cff3=dtdays*phymrd(ng)
2050 cff2=dtdays*phymrn(ng)
2051 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2052 DO k=1,n(ng)
2053 DO i=istr,iend
2054 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
2055 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2056 & bio(i,k,iphyt)*cff2
2057 bio(i,k,isdet)=bio(i,k,isdet)+ &
2058 & bio(i,k,iphyt)*cff3
2059#ifdef IRON_LIMIT
2060 bio(i,k,ifphy)=bio(i,k,ifphy)*cff1
2061 bio(i,k,ifdis)=bio(i,k,ifdis)+ &
2062 & bio(i,k,ifphy)*(cff2+cff3)*ferr(ng)
2063#endif
2064 END DO
2065 END DO
2066!
2067! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
2068! (ZooMRD rate).
2069!
2070 cff3=dtdays*zoomrd(ng)
2071 cff2=dtdays*zoomrn(ng)
2072 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2073 DO k=1,n(ng)
2074 DO i=istr,iend
2075 bio(i,k,izoop)=bio(i,k,izoop)*cff1
2076 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2077 & bio(i,k,izoop)*cff2
2078 bio(i,k,isdet)=bio(i,k,isdet)+ &
2079 & bio(i,k,izoop)*cff3
2080 END DO
2081 END DO
2082!
2083! Detritus breakdown to nutrients: remineralization (DetRR rate).
2084!
2085 cff2=dtdays*detrr(ng)
2086 cff1=1.0_r8/(1.0_r8+cff2)
2087 DO k=1,n(ng)
2088 DO i=istr,iend
2089 bio(i,k,isdet)=bio(i,k,isdet)*cff1
2090 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2091 & bio(i,k,isdet)*cff2
2092 END DO
2093 END DO
2094!
2095 IF (iteradj.ne.iter) THEN
2096!
2097!-----------------------------------------------------------------------
2098! Vertical sinking terms: Phytoplankton and Detritus
2099!-----------------------------------------------------------------------
2100!
2101! Reconstruct vertical profile of selected biological constituents
2102! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
2103! grid box. Then, compute semi-Lagrangian flux due to sinking.
2104!
2105 DO isink=1,nsink
2106 ibio=idsink(isink)
2107!
2108! Copy concentration of biological particulates into scratch array
2109! "qc" (q-central, restrict it to be positive) which is hereafter
2110! interpreted as a set of grid-box averaged values for biogeochemical
2111! constituent concentration.
2112!
2113 DO k=1,n(ng)
2114 DO i=istr,iend
2115 qc(i,k)=bio(i,k,ibio)
2116 END DO
2117 END DO
2118!
2119 DO k=n(ng)-1,1,-1
2120 DO i=istr,iend
2121 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2122 END DO
2123 END DO
2124 DO k=2,n(ng)-1
2125 DO i=istr,iend
2126 dltr=hz(i,j,k)*fc(i,k)
2127 dltl=hz(i,j,k)*fc(i,k-1)
2128 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2129 cffr=cff*fc(i,k)
2130 cffl=cff*fc(i,k-1)
2131!
2132! Apply PPM monotonicity constraint to prevent oscillations within the
2133! grid box.
2134!
2135 IF ((dltr*dltl).le.0.0_r8) THEN
2136 dltr=0.0_r8
2137 dltl=0.0_r8
2138 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2139 dltr=cffl
2140 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2141 dltl=cffr
2142 END IF
2143!
2144! Compute right and left side values (bR,bL) of parabolic segments
2145! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
2146!
2147! NOTE: Although each parabolic segment is monotonic within its grid
2148! box, monotonicity of the whole profile is not guaranteed,
2149! because bL(k+1)-bR(k) may still have different sign than
2150! qc(i,k+1)-qc(i,k). This possibility is excluded,
2151! after bL and bR are reconciled using WENO procedure.
2152!
2153 cff=(dltr-dltl)*hz_inv3(i,k)
2154 dltr=dltr-cff*hz(i,j,k+1)
2155 dltl=dltl+cff*hz(i,j,k-1)
2156 br(i,k)=qc(i,k)+dltr
2157 bl(i,k)=qc(i,k)-dltl
2158 wr(i,k)=(2.0_r8*dltr-dltl)**2
2159 wl(i,k)=(dltr-2.0_r8*dltl)**2
2160 END DO
2161 END DO
2162 cff=1.0e-14_r8
2163 DO k=2,n(ng)-2
2164 DO i=istr,iend
2165 dltl=max(cff,wl(i,k ))
2166 dltr=max(cff,wr(i,k+1))
2167 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2168 bl(i,k+1)=br(i,k)
2169 END DO
2170 END DO
2171 DO i=istr,iend
2172 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
2173#if defined LINEAR_CONTINUATION
2174 bl(i,n(ng))=br(i,n(ng)-1)
2175 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
2176#elif defined NEUMANN
2177 bl(i,n(ng))=br(i,n(ng)-1)
2178 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
2179#else
2180 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
2181 bl(i,n(ng))=qc(i,n(ng)) ! conditions
2182 br(i,n(ng)-1)=qc(i,n(ng))
2183#endif
2184#if defined LINEAR_CONTINUATION
2185 br(i,1)=bl(i,2)
2186 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2187#elif defined NEUMANN
2188 br(i,1)=bl(i,2)
2189 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2190#else
2191 bl(i,2)=qc(i,1) ! bottom grid boxes are
2192 br(i,1)=qc(i,1) ! re-assumed to be
2193 bl(i,1)=qc(i,1) ! piecewise constant.
2194#endif
2195 END DO
2196!
2197! Apply monotonicity constraint again, since the reconciled interfacial
2198! values may cause a non-monotonic behavior of the parabolic segments
2199! inside the grid box.
2200!
2201 DO k=1,n(ng)
2202 DO i=istr,iend
2203 dltr=br(i,k)-qc(i,k)
2204 dltl=qc(i,k)-bl(i,k)
2205 cffr=2.0_r8*dltr
2206 cffl=2.0_r8*dltl
2207 IF ((dltr*dltl).lt.0.0_r8) THEN
2208 dltr=0.0_r8
2209 dltl=0.0_r8
2210 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2211 dltr=cffl
2212 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2213 dltl=cffr
2214 END IF
2215 br(i,k)=qc(i,k)+dltr
2216 bl(i,k)=qc(i,k)-dltl
2217 END DO
2218 END DO
2219!
2220! After this moment reconstruction is considered complete. The next
2221! stage is to compute vertical advective fluxes, FC. It is expected
2222! that sinking may occurs relatively fast, the algorithm is designed
2223! to be free of CFL criterion, which is achieved by allowing
2224! integration bounds for semi-Lagrangian advective flux to use as
2225! many grid boxes in upstream direction as necessary.
2226!
2227! In the two code segments below, WL is the z-coordinate of the
2228! departure point for grid box interface z_w with the same indices;
2229! FC is the finite volume flux; ksource(:,k) is index of vertical
2230! grid box which contains the departure point (restricted by N(ng)).
2231! During the search: also add in content of whole grid boxes
2232! participating in FC.
2233!
2234 cff=dtdays*abs(wbio(isink))
2235 DO k=1,n(ng)
2236 DO i=istr,iend
2237 fc(i,k-1)=0.0_r8
2238 wl(i,k)=z_w(i,j,k-1)+cff
2239 wr(i,k)=hz(i,j,k)*qc(i,k)
2240 ksource(i,k)=k
2241 END DO
2242 END DO
2243 DO k=1,n(ng)
2244 DO ks=k,n(ng)-1
2245 DO i=istr,iend
2246 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2247 ksource(i,k)=ks+1
2248 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2249 END IF
2250 END DO
2251 END DO
2252 END DO
2253!
2254! Finalize computation of flux: add fractional part.
2255!
2256 DO k=1,n(ng)
2257 DO i=istr,iend
2258 ks=ksource(i,k)
2259 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2260 fc(i,k-1)=fc(i,k-1)+ &
2261 & hz(i,j,ks)*cu* &
2262 & (bl(i,ks)+ &
2263 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2264 & (1.5_r8-cu)* &
2265 & (br(i,ks)+bl(i,ks)- &
2266 & 2.0_r8*qc(i,ks))))
2267 END DO
2268 END DO
2269 DO k=1,n(ng)
2270 DO i=istr,iend
2271 bio(i,k,ibio)=qc(i,k)+ &
2272 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2273 END DO
2274 END DO
2275 END DO
2276 END IF
2277 END DO
2278!
2279! End of compute basic state arrays III.
2280!
2281!-----------------------------------------------------------------------
2282! Tangent linear vertical sinking terms.
2283!-----------------------------------------------------------------------
2284!
2285! Reconstruct vertical profile of selected biological constituents
2286! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
2287! grid box. Then, compute semi-Lagrangian flux due to sinking.
2288!
2289 sink_loop: DO isink=1,nsink
2290 ibio=idsink(isink)
2291!
2292! Copy concentration of biological particulates into scratch array
2293! "qc" (q-central, restrict it to be positive) which is hereafter
2294! interpreted as a set of grid-box averaged values for biogeochemical
2295! constituent concentration.
2296!
2297 DO k=1,n(ng)
2298 DO i=istr,iend
2299 qc(i,k)=bio(i,k,ibio)
2300 tl_qc(i,k)=tl_bio(i,k,ibio)
2301 END DO
2302 END DO
2303!
2304 DO k=n(ng)-1,1,-1
2305 DO i=istr,iend
2306 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2307 tl_fc(i,k)=(tl_qc(i,k+1)-tl_qc(i,k))*hz_inv2(i,k)+ &
2308 & (qc(i,k+1)-qc(i,k))*tl_hz_inv2(i,k)
2309 END DO
2310 END DO
2311 DO k=2,n(ng)-1
2312 DO i=istr,iend
2313 dltr=hz(i,j,k)*fc(i,k)
2314 tl_dltr=tl_hz(i,j,k)*fc(i,k)+hz(i,j,k)*tl_fc(i,k)
2315 dltl=hz(i,j,k)*fc(i,k-1)
2316 tl_dltl=tl_hz(i,j,k)*fc(i,k-1)+hz(i,j,k)*tl_fc(i,k-1)
2317 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2318 tl_cff=tl_hz(i,j,k-1)+2.0_r8*tl_hz(i,j,k)+tl_hz(i,j,k+1)
2319 cffr=cff*fc(i,k)
2320 tl_cffr=tl_cff*fc(i,k)+cff*tl_fc(i,k)
2321 cffl=cff*fc(i,k-1)
2322 tl_cffl=tl_cff*fc(i,k-1)+cff*tl_fc(i,k-1)
2323!
2324! Apply PPM monotonicity constraint to prevent oscillations within the
2325! grid box.
2326!
2327 IF ((dltr*dltl).le.0.0_r8) THEN
2328 dltr=0.0_r8
2329 tl_dltr=0.0_r8
2330 dltl=0.0_r8
2331 tl_dltl=0.0_r8
2332 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2333 dltr=cffl
2334 tl_dltr=tl_cffl
2335 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2336 dltl=cffr
2337 tl_dltl=tl_cffr
2338 END IF
2339!
2340! Compute right and left side values (bR,bL) of parabolic segments
2341! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
2342!
2343! NOTE: Although each parabolic segment is monotonic within its grid
2344! box, monotonicity of the whole profile is not guaranteed,
2345! because bL(k+1)-bR(k) may still have different sign than
2346! qc(i,k+1)-qc(i,k). This possibility is excluded,
2347! after bL and bR are reconciled using WENO procedure.
2348!
2349 cff=(dltr-dltl)*hz_inv3(i,k)
2350 tl_cff=(tl_dltr-tl_dltl)*hz_inv3(i,k)+ &
2351 & (dltr-dltl)*tl_hz_inv3(i,k)
2352 dltr=dltr-cff*hz(i,j,k+1)
2353 tl_dltr=tl_dltr-tl_cff*hz(i,j,k+1)-cff*tl_hz(i,j,k+1)
2354 dltl=dltl+cff*hz(i,j,k-1)
2355 tl_dltl=tl_dltl+tl_cff*hz(i,j,k-1)+cff*tl_hz(i,j,k-1)
2356 br(i,k)=qc(i,k)+dltr
2357 tl_br(i,k)=tl_qc(i,k)+tl_dltr
2358 bl(i,k)=qc(i,k)-dltl
2359 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
2360 wr(i,k)=(2.0_r8*dltr-dltl)**2
2361 tl_wr(i,k)=2.0_r8*(2.0_r8*dltr-dltl)* &
2362 & (2.0_r8*tl_dltr-tl_dltl)
2363 wl(i,k)=(dltr-2.0_r8*dltl)**2
2364 tl_wl(i,k)=2.0_r8*(dltr-2.0_r8*dltl)* &
2365 & (tl_dltr-2.0_r8*tl_dltl)
2366 END DO
2367 END DO
2368 cff=1.0e-14_r8
2369 DO k=2,n(ng)-2
2370 DO i=istr,iend
2371 dltl=max(cff,wl(i,k ))
2372 tl_dltl=(0.5_r8-sign(0.5_r8,cff-wl(i,k )))* &
2373 & tl_wl(i,k )
2374 dltr=max(cff,wr(i,k+1))
2375 tl_dltr=(0.5_r8-sign(0.5_r8,cff-wr(i,k+1)))* &
2376 & tl_wr(i,k+1)
2377 br1(i,k)=br(i,k)
2378 bl1(i,k+1)=bl(i,k+1)
2379 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2380 tl_br(i,k)=(tl_dltr*br1(i,k )+dltr*tl_br(i,k )+ &
2381 & tl_dltl*bl1(i,k+1)+dltl*tl_bl(i,k+1))/ &
2382 & (dltr+dltl)- &
2383 & (tl_dltr+tl_dltl)*br(i,k)/(dltr+dltl)
2384 bl(i,k+1)=br(i,k)
2385 tl_bl(i,k+1)=tl_br(i,k)
2386 END DO
2387 END DO
2388 DO i=istr,iend
2389 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
2390 tl_fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
2391#if defined LINEAR_CONTINUATION
2392 bl(i,n(ng))=br(i,n(ng)-1)
2393 tl_bl(i,n(ng))=tl_br(i,n(ng)-1)
2394 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
2395 tl_br(i,n(ng))=2.0_r8*tl_qc(i,n(ng))-tl_bl(i,n(ng))
2396#elif defined NEUMANN
2397 bl(i,n(ng))=br(i,n(ng)-1)
2398 tl_bl(i,n(ng))=tl_br(i,n(ng)-1)
2399 br(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
2400 tl_br(i,n(ng))=1.5_r8*tl_qc(i,n(ng))-0.5_r8*tl_bl(i,n(ng))
2401#else
2402 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
2403 bl(i,n(ng))=qc(i,n(ng)) ! conditions
2404 br(i,n(ng)-1)=qc(i,n(ng))
2405 tl_br(i,n(ng))=tl_qc(i,n(ng)) ! default strictly monotonic
2406 tl_bl(i,n(ng))=tl_qc(i,n(ng)) ! conditions
2407 tl_br(i,n(ng)-1)=tl_qc(i,n(ng))
2408#endif
2409#if defined LINEAR_CONTINUATION
2410 br(i,1)=bl(i,2)
2411 tl_br(i,1)=tl_bl(i,2)
2412 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2413 tl_bl(i,1)=2.0_r8*tl_qc(i,1)-tl_br(i,1)
2414#elif defined NEUMANN
2415 br(i,1)=bl(i,2)
2416 tl_br(i,1)=tl_bl(i,2)
2417 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2418 tl_bl(i,1)=1.5_r8*tl_qc(i,1)-0.5_r8*tl_br(i,1)
2419#else
2420 bl(i,2)=qc(i,1) ! bottom grid boxes are
2421 br(i,1)=qc(i,1) ! re-assumed to be
2422 bl(i,1)=qc(i,1) ! piecewise constant.
2423 tl_bl(i,2)=tl_qc(i,1) ! bottom grid boxes are
2424 tl_br(i,1)=tl_qc(i,1) ! re-assumed to be
2425 tl_bl(i,1)=tl_qc(i,1) ! piecewise constant.
2426#endif
2427 END DO
2428!
2429! Apply monotonicity constraint again, since the reconciled interfacial
2430! values may cause a non-monotonic behavior of the parabolic segments
2431! inside the grid box.
2432!
2433 DO k=1,n(ng)
2434 DO i=istr,iend
2435 dltr=br(i,k)-qc(i,k)
2436 tl_dltr=tl_br(i,k)-tl_qc(i,k)
2437 dltl=qc(i,k)-bl(i,k)
2438 tl_dltl=tl_qc(i,k)-tl_bl(i,k)
2439 cffr=2.0_r8*dltr
2440 tl_cffr=2.0_r8*tl_dltr
2441 cffl=2.0_r8*dltl
2442 tl_cffl=2.0_r8*tl_dltl
2443 IF ((dltr*dltl).lt.0.0_r8) THEN
2444 dltr=0.0_r8
2445 tl_dltr=0.0_r8
2446 dltl=0.0_r8
2447 tl_dltl=0.0_r8
2448 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2449 dltr=cffl
2450 tl_dltr=tl_cffl
2451 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2452 dltl=cffr
2453 tl_dltl=tl_cffr
2454 END IF
2455 br(i,k)=qc(i,k)+dltr
2456 tl_br(i,k)=tl_qc(i,k)+tl_dltr
2457 bl(i,k)=qc(i,k)-dltl
2458 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
2459 END DO
2460 END DO
2461!
2462! After this moment reconstruction is considered complete. The next
2463! stage is to compute vertical advective fluxes, FC. It is expected
2464! that sinking may occurs relatively fast, the algorithm is designed
2465! to be free of CFL criterion, which is achieved by allowing
2466! integration bounds for semi-Lagrangian advective flux to use as
2467! many grid boxes in upstream direction as necessary.
2468!
2469! In the two code segments below, WL is the z-coordinate of the
2470! departure point for grid box interface z_w with the same indices;
2471! FC is the finite volume flux; ksource(:,k) is index of vertical
2472! grid box which contains the departure point (restricted by N(ng)).
2473! During the search: also add in content of whole grid boxes
2474! participating in FC.
2475!
2476 cff=dtdays*abs(wbio(isink))
2477 tl_cff=dtdays*sign(1.0_r8,wbio(isink))*tl_wbio(isink)
2478 DO k=1,n(ng)
2479 DO i=istr,iend
2480 fc(i,k-1)=0.0_r8
2481 tl_fc(i,k-1)=0.0_r8
2482 wl(i,k)=z_w(i,j,k-1)+cff
2483 tl_wl(i,k)=tl_z_w(i,j,k-1)+tl_cff
2484 wr(i,k)=hz(i,j,k)*qc(i,k)
2485 tl_wr(i,k)=tl_hz(i,j,k)*qc(i,k)+hz(i,j,k)*tl_qc(i,k)
2486 ksource(i,k)=k
2487 END DO
2488 END DO
2489 DO k=1,n(ng)
2490 DO ks=k,n(ng)-1
2491 DO i=istr,iend
2492 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2493 ksource(i,k)=ks+1
2494 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2495 tl_fc(i,k-1)=tl_fc(i,k-1)+tl_wr(i,ks)
2496 END IF
2497 END DO
2498 END DO
2499 END DO
2500!
2501! Finalize computation of flux: add fractional part.
2502!
2503 DO k=1,n(ng)
2504 DO i=istr,iend
2505 ks=ksource(i,k)
2506 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2507 tl_cu=(0.5_r8+sign(0.5_r8, &
2508 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
2509 & hz_inv(i,ks))))* &
2510 & ((tl_wl(i,k)-tl_z_w(i,j,ks-1))*hz_inv(i,ks)+ &
2511 & (wl(i,k)-z_w(i,j,ks-1))*tl_hz_inv(i,ks))
2512 fc(i,k-1)=fc(i,k-1)+ &
2513 & hz(i,j,ks)*cu* &
2514 & (bl(i,ks)+ &
2515 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2516 & (1.5_r8-cu)* &
2517 & (br(i,ks)+bl(i,ks)- &
2518 & 2.0_r8*qc(i,ks))))
2519 tl_fc(i,k-1)=tl_fc(i,k-1)+ &
2520 & (tl_hz(i,j,ks)*cu+hz(i,j,ks)*tl_cu)* &
2521 & (bl(i,ks)+ &
2522 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2523 & (1.5_r8-cu)* &
2524 & (br(i,ks)+bl(i,ks)- &
2525 & 2.0_r8*qc(i,ks))))+ &
2526 & hz(i,j,ks)*cu* &
2527 & (tl_bl(i,ks)+ &
2528 & tl_cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2529 & (1.5_r8-cu)* &
2530 & (br(i,ks)+bl(i,ks)- &
2531 & 2.0_r8*qc(i,ks)))+ &
2532 & cu*(0.5_r8*(tl_br(i,ks)-tl_bl(i,ks))+ &
2533 & tl_cu* &
2534 & (br(i,ks)+bl(i,ks)-2.0_r8*qc(i,ks))- &
2535 & (1.5_r8-cu)* &
2536 & (tl_br(i,ks)+tl_bl(i,ks)- &
2537 & 2.0_r8*tl_qc(i,ks))))
2538 END DO
2539 END DO
2540 DO k=1,n(ng)
2541 DO i=istr,iend
2542 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2543 tl_bio(i,k,ibio)=tl_qc(i,k)+ &
2544 & (tl_fc(i,k)-tl_fc(i,k-1))*hz_inv(i,k)+ &
2545 & (fc(i,k)-fc(i,k-1))*tl_hz_inv(i,k)
2546 END DO
2547 END DO
2548
2549 END DO sink_loop
2550 END DO iter_loop
2551!
2552!-----------------------------------------------------------------------
2553! Update global tracer variables: Add increment due to BGC processes
2554! to tracer array in time index "nnew". Index "nnew" is solution after
2555! advection and mixing and has transport units (m Tunits) hence the
2556! increment is multiplied by Hz. Notice that we need to subtract
2557! original values "Bio_old" at the top of the routine to just account
2558! for the concentractions affected by BGC processes. This also takes
2559! into account any constraints (non-negative concentrations, carbon
2560! concentration range) specified before entering BGC kernel. If "Bio"
2561! were unchanged by BGC processes, the increment would be exactly
2562! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
2563! bounded >=0 so that we can preserve total inventory of nutrients
2564! when advection causes tracer concentration to go negative.
2565!-----------------------------------------------------------------------
2566!
2567 DO itrc=1,nbt
2568 ibio=idbio(itrc)
2569 DO k=1,n(ng)
2570 DO i=istr,iend
2571 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
2572 tl_cff=tl_bio(i,k,ibio)-tl_bio_old(i,k,ibio)
2573!^ t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*Hz(i,j,k)
2574!^
2575 tl_t(i,j,k,nnew,ibio)=tl_t(i,j,k,nnew,ibio)+ &
2576 & tl_cff*hz(i,j,k)+cff*tl_hz(i,j,k)
2577 END DO
2578 END DO
2579 END DO
2580
2581 END DO j_loop
2582!
2583 RETURN
2584 END SUBROUTINE tl_npzd_iron_tile
2585
2586 END MODULE tl_biology_mod
mod_biology
Definition ecosim_mod.h:1
mod_biology::tl_wdet
real(r8), dimension(:), allocatable tl_wdet
Definition npzd_Franks_mod.h:59
mod_biology::parfrac
real(r8), dimension(:), allocatable parfrac
Definition fennel_mod.h:139
mod_biology::phymrd
real(r8), dimension(:), allocatable phymrd
Definition npzd_iron_mod.h:102
mod_biology::zooeed
real(r8), dimension(:), allocatable zooeed
Definition npzd_iron_mod.h:119
mod_biology::ifphy
integer ifphy
Definition npzd_iron_mod.h:68
mod_biology::detrr
real(r8), dimension(:), allocatable detrr
Definition npzd_Franks_mod.h:51
mod_biology::ferr
real(r8), dimension(:), allocatable ferr
Definition npzd_iron_mod.h:86
mod_biology::wdet
real(r8), dimension(:), allocatable wdet
Definition npzd_Franks_mod.h:57
mod_biology::isdet
integer isdet
Definition npzd_Franks_mod.h:42
mod_biology::bioiter
integer, dimension(:), allocatable bioiter
Definition ecosim_mod.h:343
mod_biology::phyis
real(r8), dimension(:), allocatable phyis
Definition fennel_mod.h:143
mod_biology::attsw
real(r8), dimension(:), allocatable attsw
Definition fennel_mod.h:125
mod_biology::a_fe
real(r8), dimension(:), allocatable a_fe
Definition npzd_iron_mod.h:84
mod_biology::ino3_
integer ino3_
Definition ecosim_mod.h:277
mod_biology::zoogr
real(r8), dimension(:), allocatable zoogr
Definition fennel_mod.h:160
mod_biology::ifdis
integer ifdis
Definition npzd_iron_mod.h:69
mod_biology::k_no3
real(r8), dimension(:), allocatable k_no3
Definition fennel_mod.h:133
mod_biology::iphyt
integer iphyt
Definition fennel_mod.h:81
mod_biology::tl_parfrac
real(r8), dimension(:), allocatable tl_parfrac
Definition npzd_iron_mod.h:96
mod_biology::tl_wphy
real(r8), dimension(:), allocatable tl_wphy
Definition npzd_iron_mod.h:114
mod_biology::fenudgtime
real(r8), dimension(:), allocatable fenudgtime
Definition npzd_iron_mod.h:90
mod_biology::ivlev
real(r8), dimension(:), allocatable ivlev
Definition npzd_iron_mod.h:81
mod_biology::attphy
real(r8), dimension(:), allocatable attphy
Definition npzd_iron_mod.h:77
mod_biology::b_fe
real(r8), dimension(:), allocatable b_fe
Definition npzd_iron_mod.h:85
mod_biology::idbio
integer, dimension(:), allocatable idbio
Definition ecosim_mod.h:256
mod_biology::wphy
real(r8), dimension(:), allocatable wphy
Definition fennel_mod.h:154
mod_biology::vm_no3
real(r8), dimension(:), allocatable vm_no3
Definition npzd_Franks_mod.h:56
mod_biology::fehmin
real(r8), dimension(:), allocatable fehmin
Definition npzd_iron_mod.h:88
mod_biology::k_fec
real(r8), dimension(:), allocatable k_fec
Definition npzd_iron_mod.h:92
mod_biology::t_fe
real(r8), dimension(:), allocatable t_fe
Definition npzd_iron_mod.h:93
mod_biology::zoomrn
real(r8), dimension(:), allocatable zoomrn
Definition npzd_iron_mod.h:123
mod_biology::femax
real(r8), dimension(:), allocatable femax
Definition npzd_iron_mod.h:89
mod_biology::phymrn
real(r8), dimension(:), allocatable phymrn
Definition npzd_iron_mod.h:103
mod_biology::izoop
integer izoop
Definition fennel_mod.h:82
mod_biology::zoomrd
real(r8), dimension(:), allocatable zoomrd
Definition npzd_iron_mod.h:122
mod_biology::zooeen
real(r8), dimension(:), allocatable zooeen
Definition npzd_iron_mod.h:120
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ncparam
Definition mod_ncparam.F:2
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_param::nbt
integer nbt
Definition mod_param.F:509
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_param::nt
integer, dimension(:), allocatable nt
Definition mod_param.F:489
mod_scalars
Definition mod_scalars.F:2
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::cp
real(dp) cp
Definition mod_scalars.F:456
mod_scalars::sec2day
real(dp), parameter sec2day
Definition mod_scalars.F:817
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
tl_biology_mod
Definition tl_biology.F:26
tl_biology_mod::tl_npzd_iron_tile
subroutine tl_npzd_iron_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nnew, rmask, h, hz, tl_hz, z_r, tl_z_r, z_w, tl_z_w, srflx, tl_srflx, t, tl_t)
Definition tl_npzd_iron.h:112
tl_biology_mod::tl_biology
subroutine, public tl_biology(ng, tile)
Definition tl_npzd_Franks.h:32
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3