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