ROMS
Loading...
Searching...
No Matches
tl_npzd_Franks.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! !
12! This routine computes the biological sources and sinks and adds !
13! then the global biological fields. !
14! !
15! Reference: !
16! !
17! Franks et al, 1986: Behavior of simple plankton model with !
18! food-level acclimation by herbivores, Marine Biology, 91, !
19! 121-129. !
20! !
21!=======================================================================
22!
23 implicit none
24!
25 PRIVATE
26 PUBLIC :: tl_biology
27!
28 CONTAINS
29!
30!***********************************************************************
31 SUBROUTINE tl_biology (ng,tile)
32!***********************************************************************
33!
34 USE mod_param
35 USE mod_grid
36 USE mod_ncparam
37 USE mod_ocean
38 USE mod_stepping
39!
40! Imported variable declarations.
41!
42 integer, intent(in) :: ng, tile
43!
44! Local variable declarations.
45!
46 character (len=*), parameter :: myfile = &
47 & __FILE__
48!
49#include "tile.h"
50!
51! Set header file name.
52!
53#ifdef DISTRIBUTE
54 IF (lbiofile(itlm)) THEN
55#else
56 IF (lbiofile(itlm).and.(tile.eq.0)) THEN
57#endif
58 lbiofile(itlm)=.false.
59 bioname(itlm)=myfile
60 END IF
61!
62#ifdef PROFILE
63 CALL wclock_on (ng, itlm, 15, __line__, myfile)
64#endif
65 CALL tl_npzd_franks_tile (ng, tile, &
66 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
67 & imins, imaxs, jmins, jmaxs, &
68 & nstp(ng), nnew(ng), &
69#ifdef MASKING
70 & grid(ng) % rmask, &
71#endif
72 & grid(ng) % Hz, &
73 & grid(ng) % tl_Hz, &
74 & grid(ng) % z_r, &
75 & grid(ng) % tl_z_r, &
76 & grid(ng) % z_w, &
77 & grid(ng) % tl_z_w, &
78 & ocean(ng) % t, &
79 & ocean(ng) % tl_t)
80
81#ifdef PROFILE
82 CALL wclock_off (ng, itlm, 15, __line__, myfile)
83#endif
84!
85 RETURN
86 END SUBROUTINE tl_biology
87!
88!-----------------------------------------------------------------------
89 SUBROUTINE tl_npzd_franks_tile (ng, tile, &
90 & LBi, UBi, LBj, UBj, UBk, UBt, &
91 & IminS, ImaxS, JminS, JmaxS, &
92 & nstp, nnew, &
93#ifdef MASKING
94 & rmask, &
95#endif
96 & Hz, tl_Hz, &
97 & z_r, tl_z_r, &
98 & z_w, tl_z_w, &
99 & t, tl_t)
100!-----------------------------------------------------------------------
101!
102 USE mod_param
103 USE mod_biology
104 USE mod_ncparam
105 USE mod_scalars
106!
107! Imported variable declarations.
108!
109 integer, intent(in) :: ng, tile
110 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
111 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
112 integer, intent(in) :: nstp, nnew
113
114#ifdef ASSUMED_SHAPE
115# ifdef MASKING
116 real(r8), intent(in) :: rmask(LBi:,LBj:)
117# endif
118 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
119 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
120 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
121 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
122
123 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
124 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
125 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
126 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
127#else
128# ifdef MASKING
129 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
130# endif
131 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
132 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
133 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
134 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
135
136 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,UBk)
137 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,UBk)
138 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:UBk)
139 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,UBk,3,UBt)
140#endif
141!
142! Local variable declarations.
143!
144 integer, parameter :: Nsink = 1
145
146 integer :: Iter, i, ibio, isink, itrc, itrmx, j, k, ks
147 integer :: Iteradj
148
149 integer, dimension(Nsink) :: idsink
150
151 real(r8), parameter :: eps = 1.0e-16_r8
152
153 real(r8) :: cff, cff1, cff2, cff3, dtdays
154 real(r8) :: tl_cff, tl_cff1
155 real(r8) :: cffL, cffR, cu, dltL, dltR
156 real(r8) :: tl_cffL, tl_cffR, tl_cu, tl_dltL, tl_dltR
157
158 real(r8), dimension(Nsink) :: Wbio
159 real(r8), dimension(Nsink) :: tl_Wbio
160
161 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
162
163 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
164 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio1
165 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
166
167 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: tl_Bio
168 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: tl_Bio_old
169
170 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
171 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FC
172
173 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
174 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
175 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
176 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
177 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
178 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
179 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL1
180 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
181 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR1
182 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
183
184 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv
185 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv2
186 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hz_inv3
187 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_WL
188 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_WR
189 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bL
190 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bR
191 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_qc
192
193#include "set_bounds.h"
194!
195!-----------------------------------------------------------------------
196! Add biological Source/Sink terms.
197!-----------------------------------------------------------------------
198!
199! Avoid computing source/sink terms if no biological iterations.
200!
201 IF (bioiter(ng).le.0) RETURN
202!
203! Set time-stepping according to the number of iterations.
204!
205 dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
206!
207! Set vertical sinking indentification vector.
208!
209 idsink(1)=isdet ! Small detritus
210!
211! Set vertical sinking velocity vector in the same order as the
212! identification vector, IDSINK.
213!
214 wbio(1)=wdet(ng) ! Small detritus
215 tl_wbio(1)=tl_wdet(ng) ! Small detritus
216!
217 j_loop : DO j=jstr,jend
218!
219! Compute inverse thickness to avoid repeated divisions.
220!
221 DO k=1,n(ng)
222 DO i=istr,iend
223 hz_inv(i,k)=1.0_r8/hz(i,j,k)
224 tl_hz_inv(i,k)=-hz_inv(i,k)*hz_inv(i,k)*tl_hz(i,j,k)
225 END DO
226 END DO
227 DO k=1,n(ng)-1
228 DO i=istr,iend
229 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
230 tl_hz_inv2(i,k)=-hz_inv2(i,k)*hz_inv2(i,k)* &
231 & (tl_hz(i,j,k)+tl_hz(i,j,k+1))
232 END DO
233 END DO
234 DO k=2,n(ng)-1
235 DO i=istr,iend
236 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
237 tl_hz_inv3(i,k)=-hz_inv3(i,k)*hz_inv3(i,k)* &
238 & (tl_hz(i,j,k-1)+tl_hz(i,j,k)+ &
239 & tl_hz(i,j,k+1))
240 END DO
241 END DO
242!
243! Clear tl_Bio and Bio arrays.
244!
245 DO itrc=1,nbt
246 ibio=idbio(itrc)
247 DO k=1,n(ng)
248 DO i=istr,iend
249 bio(i,k,ibio)=0.0_r8
250 bio1(i,k,ibio)=0.0_r8
251 tl_bio(i,k,ibio)=0.0_r8
252 END DO
253 END DO
254 END DO
255!
256! Extract biological variables from tracer arrays, place them into
257! scratch arrays, and restrict their values to be positive definite.
258! At input, all tracers (index nnew) from predictor step have
259! transport units (m Tunits) since we do not have yet the new
260! values for zeta and Hz. These are known after the 2D barotropic
261! time-stepping. In this routine, this is not a problem because
262! we only use index nstp in the right-hand-side equations.
263!
264 DO itrc=1,nbt
265 ibio=idbio(itrc)
266 DO k=1,n(ng)
267 DO i=istr,iend
268 bio_old(i,k,ibio)=t(i,j,k,nstp,ibio)
269 tl_bio_old(i,k,ibio)=tl_t(i,j,k,nstp,ibio)
270 END DO
271 END DO
272 END DO
273!
274! Determine Correction for negativity.
275!
276 DO k=1,n(ng)
277 DO i=istr,iend
278 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
279 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
280 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
281 & max(0.0_r8,eps-bio_old(i,k,isdet))
282 tl_cff1=-(0.5_r8-sign(0.5_r8,bio_old(i,k,ino3_)-eps))* &
283 & tl_bio_old(i,k,ino3_)- &
284 & (0.5_r8-sign(0.5_r8,bio_old(i,k,iphyt)-eps))* &
285 & tl_bio_old(i,k,iphyt)- &
286 & (0.5_r8-sign(0.5_r8,bio_old(i,k,izoop)-eps))* &
287 & tl_bio_old(i,k,izoop)- &
288 & (0.5_r8-sign(0.5_r8,bio_old(i,k,isdet)-eps))* &
289 & tl_bio_old(i,k,isdet)
290!
291! If correction needed, determine the largest pool to debit.
292!
293 IF (cff1.gt.0.0) THEN
294 itrmx=idbio(1)
295 cff=t(i,j,k,nstp,itrmx)
296 DO ibio=idbio(2),idbio(nbt)
297 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
298 itrmx=ibio
299 cff=t(i,j,k,nstp,ibio)
300 END IF
301 END DO
302!
303! Update new values.
304!
305 DO itrc=1,nbt
306 ibio=idbio(itrc)
307 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
308 & cff1* &
309 & (sign(0.5_r8, real(itrmx-ibio,r8)**2)+ &
310 & sign(0.5_r8,-real(itrmx-ibio,r8)**2))
311 tl_bio(i,k,ibio)=(0.5_r8- &
312 & sign(0.5_r8,eps-bio_old(i,k,ibio)))* &
313 & tl_bio_old(i,k,ibio)- &
314 & tl_cff1* &
315 & (sign(0.5_r8, real(itrmx-ibio,r8)**2)+ &
316 & sign(0.5_r8,-real(itrmx-ibio,r8)**2))
317 END DO
318 ELSE
319 DO itrc=1,nbt
320 ibio=idbio(itrc)
321 bio(i,k,ibio)=bio_old(i,k,ibio)
322 tl_bio(i,k,ibio)=tl_bio_old(i,k,ibio)
323 END DO
324 END IF
325 END DO
326 END DO
327!
328!=======================================================================
329! Start internal iterations to achieve convergence of the nonlinear
330! backward-implicit solution.
331!=======================================================================
332!
333! During the iterative procedure a series of fractional time steps are
334! performed in a chained mode (splitting by different biological
335! conversion processes) in sequence of the main food chain. In all
336! stages the concentration of the component being consumed is treated
337! in a fully implicit manner, so the algorithm guarantees non-negative
338! values, no matter how strong the concentration of active consuming
339! component (Phytoplankton or Zooplankton). The overall algorithm,
340! as well as any stage of it, is formulated in conservative form
341! (except explicit sinking) in sense that the sum of concentration of
342! all components is conserved.
343!
344! In the implicit algorithm, we have for example (N: nutrient,
345! P: phytoplankton),
346!
347! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
348! {Michaelis-Menten}
349! below, we set
350! The N in the numerator of
351! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
352! as N(new)
353!
354! so the time-stepping of the equations becomes:
355!
356! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
357! consuming, divide by (1 + cff)
358! and
359!
360! P(new) = P(old) + cff * N(new) (2) when adding a source term,
361! growing, add (cff * source)
362!
363! Notice that if you substitute (1) in (2), you will get:
364!
365! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
366!
367! If you add (1) and (3), you get
368!
369! N(new) + P(new) = N(old) + P(old)
370!
371! implying conservation regardless how "cff" is computed. Therefore,
372! this scheme is unconditionally stable regardless of the conversion
373! rate. It does not generate negative values since the constituent
374! to be consumed is always treated implicitly. It is also biased
375! toward damping oscillations.
376!
377! The iterative loop below is to iterate toward an universal Backward-
378! Euler treatment of all terms. So if there are oscillations in the
379! system, they are only physical oscillations. These iterations,
380! however, do not improve the accuaracy of the solution.
381!
382 iter_loop: DO iter=1,bioiter(ng)
383!
384! Compute appropriate basic state arrays I.
385!
386! Determine Correction for negativity.
387!
388 DO k=1,n(ng)
389 DO i=istr,iend
390 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
391 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
392 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
393 & max(0.0_r8,eps-bio_old(i,k,isdet))
394!
395! If correction needed, determine the largest pool to debit.
396!
397 IF (cff1.gt.0.0) THEN
398 itrmx=idbio(1)
399 cff=t(i,j,k,nstp,itrmx)
400 DO ibio=idbio(2),idbio(nbt)
401 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
402 itrmx=ibio
403 cff=t(i,j,k,nstp,ibio)
404 END IF
405 END DO
406!
407! Update new values.
408!
409 DO itrc=1,nbt
410 ibio=idbio(itrc)
411 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
412 & cff1*(sign(0.5_r8, &
413 & real(itrmx-ibio,r8)**2)+ &
414 & sign(0.5_r8, &
415 & -real(itrmx-ibio,r8)**2))
416 END DO
417 ELSE
418 DO itrc=1,nbt
419 ibio=idbio(itrc)
420 bio(i,k,ibio)=bio_old(i,k,ibio)
421 END DO
422 END IF
423 END DO
424 END DO
425!
426!=======================================================================
427! Start internal iterations to achieve convergence of the nonlinear
428! backward-implicit solution.
429!=======================================================================
430!
431 DO iteradj=1,iter
432!
433! Nutrient uptake by phytoplankton.
434!
435 cff1=dtdays*vm_no3(ng)
436 DO k=1,n(ng)
437 DO i=istr,iend
438 cff=bio(i,k,iphyt)* &
439 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
440 & (k_no3(ng)+bio(i,k,ino3_))
441 bio1(i,k,ino3_)=bio(i,k,ino3_)
442 bio(i,k,ino3_)=bio(i,k,ino3_)/ &
443 & (1.0_r8+cff)
444 bio1(i,k,iphyt)=bio(i,k,iphyt)
445 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
446 & bio(i,k,ino3_)*cff
447 END DO
448 END DO
449!
450 IF (iteradj.ne.iter) THEN
451!
452! Phytoplankton grazing by Zooplankton and mortality to Detritus
453! (rate: PhyMR).
454!
455 cff1=dtdays*zoogr(ng)
456 cff2=dtdays*phymr(ng)
457 cff3=k_phy(ng)*k_phy(ng)
458 DO k=1,n(ng)
459 DO i=istr,iend
460 cff=bio(i,k,izoop)*bio(i,k,iphyt)*cff1/ &
461 & (cff3+bio(i,k,iphyt)*bio(i,k,iphyt))
462 bio1(i,k,iphyt)=bio(i,k,iphyt)
463 bio(i,k,iphyt)=bio(i,k,iphyt)/ &
464 & (1.0_r8+cff+cff2)
465 bio1(i,k,izoop)=bio(i,k,izoop)
466 bio(i,k,izoop)=bio(i,k,izoop)+ &
467 & bio(i,k,iphyt)*cff*(1.0_r8-zooga(ng))
468 bio(i,k,isdet)=bio(i,k,isdet)+ &
469 & bio(i,k,iphyt)* &
470 & (cff2+cff*(zooga(ng)-zooec(ng)))
471 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
472 & bio(i,k,iphyt)*cff*zooec(ng)
473 END DO
474 END DO
475!
476! Zooplankton excretion to nutrients and mortality to Detritus.
477!
478 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
479 cff2=dtdays*zoomr(ng)
480 cff3=dtdays*zoomd(ng)
481 DO k=1,n(ng)
482 DO i=istr,iend
483 bio(i,k,izoop)=bio(i,k,izoop)*cff1
484 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
485 & bio(i,k,izoop)*cff2
486 bio(i,k,isdet)=bio(i,k,isdet)+ &
487 & bio(i,k,izoop)*cff3
488 END DO
489 END DO
490!
491! Detritus breakdown to nutrients.
492!
493 cff1=dtdays*detrr(ng)
494 cff2=1.0_r8/(1.0_r8+cff1)
495 DO k=1,n(ng)
496 DO i=istr,iend
497 bio(i,k,isdet)=bio(i,k,isdet)*cff2
498 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
499 & bio(i,k,isdet)*cff1
500 END DO
501 END DO
502!
503!-----------------------------------------------------------------------
504! Vertical sinking terms.
505!-----------------------------------------------------------------------
506!
507! Reconstruct vertical profile of selected biological constituents
508! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
509! grid box. Then, compute semi-Lagrangian flux due to sinking.
510!
511 DO isink=1,nsink
512 ibio=idsink(isink)
513!
514! Copy concentration of biological particulates into scratch array
515! "qc" (q-central, restrict it to be positive) which is hereafter
516! interpreted as a set of grid-box averaged values for biogeochemical
517! constituent concentration.
518!
519 DO k=1,n(ng)
520 DO i=istr,iend
521 qc(i,k)=bio(i,k,ibio)
522 END DO
523 END DO
524!
525 DO k=n(ng)-1,1,-1
526 DO i=istr,iend
527 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
528 END DO
529 END DO
530 DO k=2,n(ng)-1
531 DO i=istr,iend
532 dltr=hz(i,j,k)*fc(i,k)
533 dltl=hz(i,j,k)*fc(i,k-1)
534 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
535 cffr=cff*fc(i,k)
536 cffl=cff*fc(i,k-1)
537!
538! Apply PPM monotonicity constraint to prevent oscillations within the
539! grid box.
540!
541 IF ((dltr*dltl).le.0.0_r8) THEN
542 dltr=0.0_r8
543 dltl=0.0_r8
544 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
545 dltr=cffl
546 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
547 dltl=cffr
548 END IF
549!
550! Compute right and left side values (bR,bL) of parabolic segments
551! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
552!
553! NOTE: Although each parabolic segment is monotonic within its grid
554! box, monotonicity of the whole profile is not guaranteed,
555! because bL(k+1)-bR(k) may still have different sign than
556! qc(i,k+1)-qc(i,k). This possibility is excluded,
557! after bL and bR are reconciled using WENO procedure.
558!
559 cff=(dltr-dltl)*hz_inv3(i,k)
560 dltr=dltr-cff*hz(i,j,k+1)
561 dltl=dltl+cff*hz(i,j,k-1)
562 br(i,k)=qc(i,k)+dltr
563 bl(i,k)=qc(i,k)-dltl
564 wr(i,k)=(2.0_r8*dltr-dltl)**2
565 wl(i,k)=(dltr-2.0_r8*dltl)**2
566 END DO
567 END DO
568 cff=1.0e-14_r8
569 DO k=2,n(ng)-2
570 DO i=istr,iend
571 dltl=max(cff,wl(i,k ))
572 dltr=max(cff,wr(i,k+1))
573 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
574 bl(i,k+1)=br(i,k)
575 END DO
576 END DO
577 DO i=istr,iend
578 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
579#if defined LINEAR_CONTINUATION
580 bl(i,n(ng))=br(i,n(ng)-1)
581 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
582#elif defined NEUMANN
583 bl(i,n(ng))=br(i,n(ng)-1)
584 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
585#else
586 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
587 bl(i,n(ng))=qc(i,n(ng)) ! conditions
588 br(i,n(ng)-1)=qc(i,n(ng))
589#endif
590#if defined LINEAR_CONTINUATION
591 br(i,1)=bl(i,2)
592 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
593#elif defined NEUMANN
594 br(i,1)=bl(i,2)
595 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
596#else
597 bl(i,2)=qc(i,1) ! bottom grid boxes are
598 br(i,1)=qc(i,1) ! re-assumed to be
599 bl(i,1)=qc(i,1) ! piecewise constant.
600#endif
601 END DO
602!
603! Apply monotonicity constraint again, since the reconciled interfacial
604! values may cause a non-monotonic behavior of the parabolic segments
605! inside the grid box.
606!
607 DO k=1,n(ng)
608 DO i=istr,iend
609 dltr=br(i,k)-qc(i,k)
610 dltl=qc(i,k)-bl(i,k)
611 cffr=2.0_r8*dltr
612 cffl=2.0_r8*dltl
613 IF ((dltr*dltl).lt.0.0_r8) THEN
614 dltr=0.0_r8
615 dltl=0.0_r8
616 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
617 dltr=cffl
618 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
619 dltl=cffr
620 END IF
621 br(i,k)=qc(i,k)+dltr
622 bl(i,k)=qc(i,k)-dltl
623 END DO
624 END DO
625!
626! After this moment reconstruction is considered complete. The next
627! stage is to compute vertical advective fluxes, FC. It is expected
628! that sinking may occurs relatively fast, the algorithm is designed
629! to be free of CFL criterion, which is achieved by allowing
630! integration bounds for semi-Lagrangian advective flux to use as
631! many grid boxes in upstream direction as necessary.
632!
633! In the two code segments below, WL is the z-coordinate of the
634! departure point for grid box interface z_w with the same indices;
635! FC is the finite volume flux; ksource(:,k) is index of vertical
636! grid box which contains the departure point (restricted by N(ng)).
637! During the search: also add in content of whole grid boxes
638! participating in FC.
639!
640 cff=dtdays*abs(wbio(isink))
641 DO k=1,n(ng)
642 DO i=istr,iend
643 fc(i,k-1)=0.0_r8
644 wl(i,k)=z_w(i,j,k-1)+cff
645 wr(i,k)=hz(i,j,k)*qc(i,k)
646 ksource(i,k)=k
647 END DO
648 END DO
649 DO k=1,n(ng)
650 DO ks=k,n(ng)-1
651 DO i=istr,iend
652 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
653 ksource(i,k)=ks+1
654 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
655 END IF
656 END DO
657 END DO
658 END DO
659!
660! Finalize computation of flux: add fractional part.
661!
662 DO k=1,n(ng)
663 DO i=istr,iend
664 ks=ksource(i,k)
665 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
666 fc(i,k-1)=fc(i,k-1)+ &
667 & hz(i,j,ks)*cu* &
668 & (bl(i,ks)+ &
669 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
670 & (1.5_r8-cu)* &
671 & (br(i,ks)+bl(i,ks)- &
672 & 2.0_r8*qc(i,ks))))
673 END DO
674 END DO
675 DO k=1,n(ng)
676 DO i=istr,iend
677 bio(i,k,ibio)=qc(i,k)+ &
678 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
679 END DO
680 END DO
681 END DO
682 END IF
683 END DO
684!
685! End of compute basic state arrays I.
686!
687! Nutrient uptake by phytoplankton.
688!
689 cff1=dtdays*vm_no3(ng)
690 DO k=1,n(ng)
691 DO i=istr,iend
692 cff=bio1(i,k,iphyt)* &
693 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
694 & (k_no3(ng)+bio1(i,k,ino3_))
695 tl_cff=(tl_bio(i,k,iphyt)* &
696 & cff1*exp(k_ext(ng)*z_r(i,j,k))- &
697 & tl_bio(i,k,ino3_)*cff)/ &
698 & (k_no3(ng)+bio1(i,k,ino3_))+ &
699 & k_ext(ng)*tl_z_r(i,j,k)*cff
700!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)/ &
701!^ & (1.0_r8+cff)
702!^
703 tl_bio(i,k,ino3_)=(tl_bio(i,k,ino3_)- &
704 & tl_cff*bio(i,k,ino3_))/ &
705 & (1.0_r8+cff)
706!^ Bio(i,k,iPhyt)=Bio(i,k,iPhyt)+ &
707!^ & Bio(i,k,iNO3_)*cff
708!^
709 tl_bio(i,k,iphyt)=tl_bio(i,k,iphyt)+ &
710 & tl_bio(i,k,ino3_)*cff+ &
711 & bio(i,k,ino3_)*tl_cff
712 END DO
713 END DO
714!
715! Compute appropriate basic state arrays II.
716!
717! Determine Correction for negativity.
718!
719 DO k=1,n(ng)
720 DO i=istr,iend
721 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
722 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
723 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
724 & max(0.0_r8,eps-bio_old(i,k,isdet))
725!
726! If correction needed, determine the largest pool to debit.
727!
728 IF (cff1.gt.0.0) THEN
729 itrmx=idbio(1)
730 cff=t(i,j,k,nstp,itrmx)
731 DO ibio=idbio(2),idbio(nbt)
732 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
733 itrmx=ibio
734 cff=t(i,j,k,nstp,ibio)
735 END IF
736 END DO
737!
738! Update new values.
739!
740 DO itrc=1,nbt
741 ibio=idbio(itrc)
742 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
743 & cff1*(sign(0.5_r8, &
744 & real(itrmx-ibio,r8)**2)+ &
745 & sign(0.5_r8, &
746 & -real(itrmx-ibio,r8)**2))
747 END DO
748 ELSE
749 DO itrc=1,nbt
750 ibio=idbio(itrc)
751 bio(i,k,ibio)=bio_old(i,k,ibio)
752 END DO
753 END IF
754 END DO
755 END DO
756!
757!=======================================================================
758! Start internal iterations to achieve convergence of the nonlinear
759! backward-implicit solution.
760!=======================================================================
761!
762 DO iteradj=1,iter
763!
764! Nutrient uptake by phytoplankton.
765!
766 cff1=dtdays*vm_no3(ng)
767 DO k=1,n(ng)
768 DO i=istr,iend
769 cff=bio(i,k,iphyt)* &
770 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
771 & (k_no3(ng)+bio(i,k,ino3_))
772 bio1(i,k,ino3_)=bio(i,k,ino3_)
773 bio(i,k,ino3_)=bio(i,k,ino3_)/ &
774 & (1.0_r8+cff)
775 bio1(i,k,iphyt)=bio(i,k,iphyt)
776 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
777 & bio(i,k,ino3_)*cff
778 END DO
779 END DO
780!
781! Phytoplankton grazing by Zooplankton and mortality to Detritus
782! (rate: PhyMR).
783!
784 cff1=dtdays*zoogr(ng)
785 cff2=dtdays*phymr(ng)
786 cff3=k_phy(ng)*k_phy(ng)
787 DO k=1,n(ng)
788 DO i=istr,iend
789 cff=bio(i,k,izoop)*bio(i,k,iphyt)*cff1/ &
790 & (cff3+bio(i,k,iphyt)*bio(i,k,iphyt))
791 bio1(i,k,iphyt)=bio(i,k,iphyt)
792 bio(i,k,iphyt)=bio(i,k,iphyt)/ &
793 & (1.0_r8+cff+cff2)
794 bio1(i,k,izoop)=bio(i,k,izoop)
795 bio(i,k,izoop)=bio(i,k,izoop)+ &
796 & bio(i,k,iphyt)*cff*(1.0_r8-zooga(ng))
797 bio(i,k,isdet)=bio(i,k,isdet)+ &
798 & bio(i,k,iphyt)* &
799 & (cff2+cff*(zooga(ng)-zooec(ng)))
800 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
801 & bio(i,k,iphyt)*cff*zooec(ng)
802 END DO
803 END DO
804!
805 IF (iteradj.ne.iter) THEN
806!
807! Zooplankton excretion to nutrients and mortality to Detritus.
808!
809 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
810 cff2=dtdays*zoomr(ng)
811 cff3=dtdays*zoomd(ng)
812 DO k=1,n(ng)
813 DO i=istr,iend
814 bio(i,k,izoop)=bio(i,k,izoop)*cff1
815 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
816 & bio(i,k,izoop)*cff2
817 bio(i,k,isdet)=bio(i,k,isdet)+ &
818 & bio(i,k,izoop)*cff3
819 END DO
820 END DO
821!
822! Detritus breakdown to nutrients.
823!
824 cff1=dtdays*detrr(ng)
825 cff2=1.0_r8/(1.0_r8+cff1)
826 DO k=1,n(ng)
827 DO i=istr,iend
828 bio(i,k,isdet)=bio(i,k,isdet)*cff2
829 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
830 & bio(i,k,isdet)*cff1
831 END DO
832 END DO
833!
834!-----------------------------------------------------------------------
835! Vertical sinking terms.
836!-----------------------------------------------------------------------
837!
838! Reconstruct vertical profile of selected biological constituents
839! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
840! grid box. Then, compute semi-Lagrangian flux due to sinking.
841!
842 DO isink=1,nsink
843 ibio=idsink(isink)
844!
845! Copy concentration of biological particulates into scratch array
846! "qc" (q-central, restrict it to be positive) which is hereafter
847! interpreted as a set of grid-box averaged values for biogeochemical
848! constituent concentration.
849!
850 DO k=1,n(ng)
851 DO i=istr,iend
852 qc(i,k)=bio(i,k,ibio)
853 END DO
854 END DO
855!
856 DO k=n(ng)-1,1,-1
857 DO i=istr,iend
858 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
859 END DO
860 END DO
861 DO k=2,n(ng)-1
862 DO i=istr,iend
863 dltr=hz(i,j,k)*fc(i,k)
864 dltl=hz(i,j,k)*fc(i,k-1)
865 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
866 cffr=cff*fc(i,k)
867 cffl=cff*fc(i,k-1)
868!
869! Apply PPM monotonicity constraint to prevent oscillations within the
870! grid box.
871!
872 IF ((dltr*dltl).le.0.0_r8) THEN
873 dltr=0.0_r8
874 dltl=0.0_r8
875 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
876 dltr=cffl
877 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
878 dltl=cffr
879 END IF
880!
881! Compute right and left side values (bR,bL) of parabolic segments
882! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
883!
884! NOTE: Although each parabolic segment is monotonic within its grid
885! box, monotonicity of the whole profile is not guaranteed,
886! because bL(k+1)-bR(k) may still have different sign than
887! qc(i,k+1)-qc(i,k). This possibility is excluded,
888! after bL and bR are reconciled using WENO procedure.
889!
890 cff=(dltr-dltl)*hz_inv3(i,k)
891 dltr=dltr-cff*hz(i,j,k+1)
892 dltl=dltl+cff*hz(i,j,k-1)
893 br(i,k)=qc(i,k)+dltr
894 bl(i,k)=qc(i,k)-dltl
895 wr(i,k)=(2.0_r8*dltr-dltl)**2
896 wl(i,k)=(dltr-2.0_r8*dltl)**2
897 END DO
898 END DO
899 cff=1.0e-14_r8
900 DO k=2,n(ng)-2
901 DO i=istr,iend
902 dltl=max(cff,wl(i,k ))
903 dltr=max(cff,wr(i,k+1))
904 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
905 bl(i,k+1)=br(i,k)
906 END DO
907 END DO
908 DO i=istr,iend
909 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
910#if defined LINEAR_CONTINUATION
911 bl(i,n(ng))=br(i,n(ng)-1)
912 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
913#elif defined NEUMANN
914 bl(i,n(ng))=br(i,n(ng)-1)
915 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
916#else
917 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
918 bl(i,n(ng))=qc(i,n(ng)) ! conditions
919 br(i,n(ng)-1)=qc(i,n(ng))
920#endif
921#if defined LINEAR_CONTINUATION
922 br(i,1)=bl(i,2)
923 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
924#elif defined NEUMANN
925 br(i,1)=bl(i,2)
926 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
927#else
928 bl(i,2)=qc(i,1) ! bottom grid boxes are
929 br(i,1)=qc(i,1) ! re-assumed to be
930 bl(i,1)=qc(i,1) ! piecewise constant.
931#endif
932 END DO
933!
934! Apply monotonicity constraint again, since the reconciled interfacial
935! values may cause a non-monotonic behavior of the parabolic segments
936! inside the grid box.
937!
938 DO k=1,n(ng)
939 DO i=istr,iend
940 dltr=br(i,k)-qc(i,k)
941 dltl=qc(i,k)-bl(i,k)
942 cffr=2.0_r8*dltr
943 cffl=2.0_r8*dltl
944 IF ((dltr*dltl).lt.0.0_r8) THEN
945 dltr=0.0_r8
946 dltl=0.0_r8
947 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
948 dltr=cffl
949 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
950 dltl=cffr
951 END IF
952 br(i,k)=qc(i,k)+dltr
953 bl(i,k)=qc(i,k)-dltl
954 END DO
955 END DO
956!
957! After this moment reconstruction is considered complete. The next
958! stage is to compute vertical advective fluxes, FC. It is expected
959! that sinking may occurs relatively fast, the algorithm is designed
960! to be free of CFL criterion, which is achieved by allowing
961! integration bounds for semi-Lagrangian advective flux to use as
962! many grid boxes in upstream direction as necessary.
963!
964! In the two code segments below, WL is the z-coordinate of the
965! departure point for grid box interface z_w with the same indices;
966! FC is the finite volume flux; ksource(:,k) is index of vertical
967! grid box which contains the departure point (restricted by N(ng)).
968! During the search: also add in content of whole grid boxes
969! participating in FC.
970!
971 cff=dtdays*abs(wbio(isink))
972 DO k=1,n(ng)
973 DO i=istr,iend
974 fc(i,k-1)=0.0_r8
975 wl(i,k)=z_w(i,j,k-1)+cff
976 wr(i,k)=hz(i,j,k)*qc(i,k)
977 ksource(i,k)=k
978 END DO
979 END DO
980 DO k=1,n(ng)
981 DO ks=k,n(ng)-1
982 DO i=istr,iend
983 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
984 ksource(i,k)=ks+1
985 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
986 END IF
987 END DO
988 END DO
989 END DO
990!
991! Finalize computation of flux: add fractional part.
992!
993 DO k=1,n(ng)
994 DO i=istr,iend
995 ks=ksource(i,k)
996 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
997 fc(i,k-1)=fc(i,k-1)+ &
998 & hz(i,j,ks)*cu* &
999 & (bl(i,ks)+ &
1000 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1001 & (1.5_r8-cu)* &
1002 & (br(i,ks)+bl(i,ks)- &
1003 & 2.0_r8*qc(i,ks))))
1004 END DO
1005 END DO
1006 DO k=1,n(ng)
1007 DO i=istr,iend
1008 bio(i,k,ibio)=qc(i,k)+ &
1009 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1010 END DO
1011 END DO
1012 END DO
1013 END IF
1014 END DO
1015!
1016! End of compute basic state arrays II.
1017!
1018! Phytoplankton grazing by Zooplankton and mortality to Detritus
1019! (rate: PhyMR).
1020!
1021 cff1=dtdays*zoogr(ng)
1022 cff2=dtdays*phymr(ng)
1023 cff3=k_phy(ng)*k_phy(ng)
1024 DO k=1,n(ng)
1025 DO i=istr,iend
1026 cff=bio1(i,k,izoop)*bio1(i,k,iphyt)*cff1/ &
1027 & (cff3+bio1(i,k,iphyt)*bio1(i,k,iphyt))
1028 tl_cff=((tl_bio(i,k,izoop)*bio1(i,k,iphyt)+ &
1029 & bio1(i,k,izoop)*tl_bio(i,k,iphyt))*cff1- &
1030 & 2.0_r8*bio1(i,k,iphyt)*tl_bio(i,k,iphyt)*cff)/ &
1031 & (cff3+bio1(i,k,iphyt)*bio1(i,k,iphyt))
1032!^ Bio(i,k,iPhyt)=Bio(i,k,iPhyt)/ &
1033!^ & (1.0_r8+cff+cff2)
1034!^
1035 tl_bio(i,k,iphyt)=(tl_bio(i,k,iphyt)- &
1036 & tl_cff*bio(i,k,iphyt))/ &
1037 & (1.0_r8+cff+cff2)
1038!^ Bio(i,k,iZoop)=Bio(i,k,iZoop)+ &
1039!^ & Bio(i,k,iPhyt)*cff*(1.0_r8-ZooGA(ng))
1040!^
1041 tl_bio(i,k,izoop)=tl_bio(i,k,izoop)+ &
1042 & tl_bio(i,k,iphyt)* &
1043 & cff*(1.0_r8-zooga(ng))+ &
1044 & bio(i,k,iphyt)* &
1045 & tl_cff*(1.0_r8-zooga(ng))
1046!^ Bio(i,k,iSDet)=Bio(i,k,iSDet)+ &
1047!^ & Bio(i,k,iPhyt)* &
1048!^ & (cff2+cff*(ZooGA(ng)-ZooEC(ng)))
1049!^
1050 tl_bio(i,k,isdet)=tl_bio(i,k,isdet)+ &
1051 & tl_bio(i,k,iphyt)* &
1052 & (cff2+cff*(zooga(ng)-zooec(ng)))+ &
1053 & bio(i,k,iphyt)* &
1054 & tl_cff*(zooga(ng)-zooec(ng))
1055!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
1056!^ & Bio(i,k,iPhyt)*cff*ZooEC(ng)
1057!^
1058 tl_bio(i,k,ino3_)=tl_bio(i,k,ino3_)+ &
1059 & tl_bio(i,k,iphyt)*cff*zooec(ng)+ &
1060 & bio(i,k,iphyt)*tl_cff*zooec(ng)
1061 END DO
1062 END DO
1063!
1064! Zooplankton excretion to nutrients and mortality to Detritus.
1065!
1066 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
1067 cff2=dtdays*zoomr(ng)
1068 cff3=dtdays*zoomd(ng)
1069 DO k=1,n(ng)
1070 DO i=istr,iend
1071!^ Bio(i,k,iZoop)=Bio(i,k,iZoop)*cff1
1072!^
1073 tl_bio(i,k,izoop)=tl_bio(i,k,izoop)*cff1
1074!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
1075!^ & Bio(i,k,iZoop)*cff2
1076!^
1077 tl_bio(i,k,ino3_)=tl_bio(i,k,ino3_)+ &
1078 & tl_bio(i,k,izoop)*cff2
1079!^ Bio(i,k,iSDet)=Bio(i,k,iSDet)+ &
1080!^ & Bio(i,k,iZoop)*cff3
1081!^
1082 tl_bio(i,k,isdet)=tl_bio(i,k,isdet)+ &
1083 & tl_bio(i,k,izoop)*cff3
1084 END DO
1085 END DO
1086!
1087! Detritus breakdown to nutrients.
1088!
1089 cff1=dtdays*detrr(ng)
1090 cff2=1.0_r8/(1.0_r8+cff1)
1091 DO k=1,n(ng)
1092 DO i=istr,iend
1093!^ Bio(i,k,iSDet)=Bio(i,k,iSDet)*cff2
1094!^
1095 tl_bio(i,k,isdet)=tl_bio(i,k,isdet)*cff2
1096!^ Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+ &
1097!^ & Bio(i,k,iSDet)*cff1
1098!^
1099 tl_bio(i,k,ino3_)=tl_bio(i,k,ino3_)+ &
1100 & tl_bio(i,k,isdet)*cff1
1101 END DO
1102 END DO
1103!
1104! Compute appropriate basic state arrays III.
1105!
1106! Determine Correction for negativity.
1107!
1108 DO k=1,n(ng)
1109 DO i=istr,iend
1110 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
1111 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
1112 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
1113 & max(0.0_r8,eps-bio_old(i,k,isdet))
1114!
1115! If correction needed, determine the largest pool to debit.
1116!
1117 IF (cff1.gt.0.0) THEN
1118 itrmx=idbio(1)
1119 cff=t(i,j,k,nstp,itrmx)
1120 DO ibio=idbio(2),idbio(nbt)
1121 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
1122 itrmx=ibio
1123 cff=t(i,j,k,nstp,ibio)
1124 END IF
1125 END DO
1126!
1127! Update new values.
1128!
1129 DO itrc=1,nbt
1130 ibio=idbio(itrc)
1131 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
1132 & cff1*(sign(0.5_r8, &
1133 & real(itrmx-ibio,r8)**2)+ &
1134 & sign(0.5_r8, &
1135 & -real(itrmx-ibio,r8)**2))
1136 END DO
1137 ELSE
1138 DO itrc=1,nbt
1139 ibio=idbio(itrc)
1140 bio(i,k,ibio)=bio_old(i,k,ibio)
1141 END DO
1142 END IF
1143 END DO
1144 END DO
1145!
1146!=======================================================================
1147! Start internal iterations to achieve convergence of the nonlinear
1148! backward-implicit solution.
1149!=======================================================================
1150!
1151 DO iteradj=1,iter
1152!
1153! Nutrient uptake by phytoplankton.
1154!
1155 cff1=dtdays*vm_no3(ng)
1156 DO k=1,n(ng)
1157 DO i=istr,iend
1158 cff=bio(i,k,iphyt)* &
1159 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
1160 & (k_no3(ng)+bio(i,k,ino3_))
1161 bio1(i,k,ino3_)=bio(i,k,ino3_)
1162 bio(i,k,ino3_)=bio(i,k,ino3_)/ &
1163 & (1.0_r8+cff)
1164 bio1(i,k,iphyt)=bio(i,k,iphyt)
1165 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
1166 & bio(i,k,ino3_)*cff
1167 END DO
1168 END DO
1169!
1170! Phytoplankton grazing by Zooplankton and mortality to Detritus
1171! (rate: PhyMR).
1172!
1173 cff1=dtdays*zoogr(ng)
1174 cff2=dtdays*phymr(ng)
1175 cff3=k_phy(ng)*k_phy(ng)
1176 DO k=1,n(ng)
1177 DO i=istr,iend
1178 cff=bio(i,k,izoop)*bio(i,k,iphyt)*cff1/ &
1179 & (cff3+bio(i,k,iphyt)*bio(i,k,iphyt))
1180 bio1(i,k,iphyt)=bio(i,k,iphyt)
1181 bio(i,k,iphyt)=bio(i,k,iphyt)/ &
1182 & (1.0_r8+cff+cff2)
1183 bio1(i,k,izoop)=bio(i,k,izoop)
1184 bio(i,k,izoop)=bio(i,k,izoop)+ &
1185 & bio(i,k,iphyt)*cff*(1.0_r8-zooga(ng))
1186 bio(i,k,isdet)=bio(i,k,isdet)+ &
1187 & bio(i,k,iphyt)* &
1188 & (cff2+cff*(zooga(ng)-zooec(ng)))
1189 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1190 & bio(i,k,iphyt)*cff*zooec(ng)
1191 END DO
1192 END DO
1193!
1194! Zooplankton excretion to nutrients and mortality to Detritus.
1195!
1196 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
1197 cff2=dtdays*zoomr(ng)
1198 cff3=dtdays*zoomd(ng)
1199 DO k=1,n(ng)
1200 DO i=istr,iend
1201 bio(i,k,izoop)=bio(i,k,izoop)*cff1
1202 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1203 & bio(i,k,izoop)*cff2
1204 bio(i,k,isdet)=bio(i,k,isdet)+ &
1205 & bio(i,k,izoop)*cff3
1206 END DO
1207 END DO
1208!
1209! Detritus breakdown to nutrients.
1210!
1211 cff1=dtdays*detrr(ng)
1212 cff2=1.0_r8/(1.0_r8+cff1)
1213 DO k=1,n(ng)
1214 DO i=istr,iend
1215 bio(i,k,isdet)=bio(i,k,isdet)*cff2
1216 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1217 & bio(i,k,isdet)*cff1
1218 END DO
1219 END DO
1220!
1221 IF (iteradj.ne.iter) THEN
1222!
1223!-----------------------------------------------------------------------
1224! Vertical sinking terms.
1225!-----------------------------------------------------------------------
1226!
1227! Reconstruct vertical profile of selected biological constituents
1228! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
1229! grid box. Then, compute semi-Lagrangian flux due to sinking.
1230!
1231 DO isink=1,nsink
1232 ibio=idsink(isink)
1233!
1234! Copy concentration of biological particulates into scratch array
1235! "qc" (q-central, restrict it to be positive) which is hereafter
1236! interpreted as a set of grid-box averaged values for biogeochemical
1237! constituent concentration.
1238!
1239 DO k=1,n(ng)
1240 DO i=istr,iend
1241 qc(i,k)=bio(i,k,ibio)
1242 END DO
1243 END DO
1244!
1245 DO k=n(ng)-1,1,-1
1246 DO i=istr,iend
1247 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1248 END DO
1249 END DO
1250 DO k=2,n(ng)-1
1251 DO i=istr,iend
1252 dltr=hz(i,j,k)*fc(i,k)
1253 dltl=hz(i,j,k)*fc(i,k-1)
1254 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1255 cffr=cff*fc(i,k)
1256 cffl=cff*fc(i,k-1)
1257!
1258! Apply PPM monotonicity constraint to prevent oscillations within the
1259! grid box.
1260!
1261 IF ((dltr*dltl).le.0.0_r8) THEN
1262 dltr=0.0_r8
1263 dltl=0.0_r8
1264 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1265 dltr=cffl
1266 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1267 dltl=cffr
1268 END IF
1269!
1270! Compute right and left side values (bR,bL) of parabolic segments
1271! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1272!
1273! NOTE: Although each parabolic segment is monotonic within its grid
1274! box, monotonicity of the whole profile is not guaranteed,
1275! because bL(k+1)-bR(k) may still have different sign than
1276! qc(i,k+1)-qc(i,k). This possibility is excluded,
1277! after bL and bR are reconciled using WENO procedure.
1278!
1279 cff=(dltr-dltl)*hz_inv3(i,k)
1280 dltr=dltr-cff*hz(i,j,k+1)
1281 dltl=dltl+cff*hz(i,j,k-1)
1282 br(i,k)=qc(i,k)+dltr
1283 bl(i,k)=qc(i,k)-dltl
1284 wr(i,k)=(2.0_r8*dltr-dltl)**2
1285 wl(i,k)=(dltr-2.0_r8*dltl)**2
1286 END DO
1287 END DO
1288 cff=1.0e-14_r8
1289 DO k=2,n(ng)-2
1290 DO i=istr,iend
1291 dltl=max(cff,wl(i,k ))
1292 dltr=max(cff,wr(i,k+1))
1293 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1294 bl(i,k+1)=br(i,k)
1295 END DO
1296 END DO
1297 DO i=istr,iend
1298 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1299#if defined LINEAR_CONTINUATION
1300 bl(i,n(ng))=br(i,n(ng)-1)
1301 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1302#elif defined NEUMANN
1303 bl(i,n(ng))=br(i,n(ng)-1)
1304 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1305#else
1306 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1307 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1308 br(i,n(ng)-1)=qc(i,n(ng))
1309#endif
1310#if defined LINEAR_CONTINUATION
1311 br(i,1)=bl(i,2)
1312 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1313#elif defined NEUMANN
1314 br(i,1)=bl(i,2)
1315 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1316#else
1317 bl(i,2)=qc(i,1) ! bottom grid boxes are
1318 br(i,1)=qc(i,1) ! re-assumed to be
1319 bl(i,1)=qc(i,1) ! piecewise constant.
1320#endif
1321 END DO
1322!
1323! Apply monotonicity constraint again, since the reconciled interfacial
1324! values may cause a non-monotonic behavior of the parabolic segments
1325! inside the grid box.
1326!
1327 DO k=1,n(ng)
1328 DO i=istr,iend
1329 dltr=br(i,k)-qc(i,k)
1330 dltl=qc(i,k)-bl(i,k)
1331 cffr=2.0_r8*dltr
1332 cffl=2.0_r8*dltl
1333 IF ((dltr*dltl).lt.0.0_r8) THEN
1334 dltr=0.0_r8
1335 dltl=0.0_r8
1336 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1337 dltr=cffl
1338 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1339 dltl=cffr
1340 END IF
1341 br(i,k)=qc(i,k)+dltr
1342 bl(i,k)=qc(i,k)-dltl
1343 END DO
1344 END DO
1345!
1346! After this moment reconstruction is considered complete. The next
1347! stage is to compute vertical advective fluxes, FC. It is expected
1348! that sinking may occurs relatively fast, the algorithm is designed
1349! to be free of CFL criterion, which is achieved by allowing
1350! integration bounds for semi-Lagrangian advective flux to use as
1351! many grid boxes in upstream direction as necessary.
1352!
1353! In the two code segments below, WL is the z-coordinate of the
1354! departure point for grid box interface z_w with the same indices;
1355! FC is the finite volume flux; ksource(:,k) is index of vertical
1356! grid box which contains the departure point (restricted by N(ng)).
1357! During the search: also add in content of whole grid boxes
1358! participating in FC.
1359!
1360 cff=dtdays*abs(wbio(isink))
1361 DO k=1,n(ng)
1362 DO i=istr,iend
1363 fc(i,k-1)=0.0_r8
1364 wl(i,k)=z_w(i,j,k-1)+cff
1365 wr(i,k)=hz(i,j,k)*qc(i,k)
1366 ksource(i,k)=k
1367 END DO
1368 END DO
1369 DO k=1,n(ng)
1370 DO ks=k,n(ng)-1
1371 DO i=istr,iend
1372 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1373 ksource(i,k)=ks+1
1374 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1375 END IF
1376 END DO
1377 END DO
1378 END DO
1379!
1380! Finalize computation of flux: add fractional part.
1381!
1382 DO k=1,n(ng)
1383 DO i=istr,iend
1384 ks=ksource(i,k)
1385 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1386 fc(i,k-1)=fc(i,k-1)+ &
1387 & hz(i,j,ks)*cu* &
1388 & (bl(i,ks)+ &
1389 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1390 & (1.5_r8-cu)* &
1391 & (br(i,ks)+bl(i,ks)- &
1392 & 2.0_r8*qc(i,ks))))
1393 END DO
1394 END DO
1395 DO k=1,n(ng)
1396 DO i=istr,iend
1397 bio(i,k,ibio)=qc(i,k)+ &
1398 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1399 END DO
1400 END DO
1401 END DO
1402 END IF
1403 END DO
1404!
1405! End of compute basic state arrays III.
1406!
1407!-----------------------------------------------------------------------
1408! Vertical sinking terms.
1409!-----------------------------------------------------------------------
1410!
1411! Reconstruct vertical profile of selected biological constituents
1412! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
1413! grid box. Then, compute semi-Lagrangian flux due to sinking.
1414!
1415 sink_loop: DO isink=1,nsink
1416 ibio=idsink(isink)
1417!
1418! Copy concentration of biological particulates into scratch array
1419! "qc" (q-central, restrict it to be positive) which is hereafter
1420! interpreted as a set of grid-box averaged values for biogeochemical
1421! constituent concentration.
1422!
1423 DO k=1,n(ng)
1424 DO i=istr,iend
1425 qc(i,k)=bio(i,k,ibio)
1426 tl_qc(i,k)=tl_bio(i,k,ibio)
1427 END DO
1428 END DO
1429!
1430 DO k=n(ng)-1,1,-1
1431 DO i=istr,iend
1432 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1433 tl_fc(i,k)=(tl_qc(i,k+1)-tl_qc(i,k))*hz_inv2(i,k)+ &
1434 & (qc(i,k+1)-qc(i,k))*tl_hz_inv2(i,k)
1435 END DO
1436 END DO
1437 DO k=2,n(ng)-1
1438 DO i=istr,iend
1439 dltr=hz(i,j,k)*fc(i,k)
1440 tl_dltr=tl_hz(i,j,k)*fc(i,k)+hz(i,j,k)*tl_fc(i,k)
1441 dltl=hz(i,j,k)*fc(i,k-1)
1442 tl_dltl=tl_hz(i,j,k)*fc(i,k-1)+hz(i,j,k)*tl_fc(i,k-1)
1443 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1444 tl_cff=tl_hz(i,j,k-1)+2.0_r8*tl_hz(i,j,k)+tl_hz(i,j,k+1)
1445 cffr=cff*fc(i,k)
1446 tl_cffr=tl_cff*fc(i,k)+cff*tl_fc(i,k)
1447 cffl=cff*fc(i,k-1)
1448 tl_cffl=tl_cff*fc(i,k-1)+cff*tl_fc(i,k-1)
1449!
1450! Apply PPM monotonicity constraint to prevent oscillations within the
1451! grid box.
1452!
1453 IF ((dltr*dltl).le.0.0_r8) THEN
1454 dltr=0.0_r8
1455 tl_dltr=0.0_r8
1456 dltl=0.0_r8
1457 tl_dltl=0.0_r8
1458 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1459 dltr=cffl
1460 tl_dltr=tl_cffl
1461 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1462 dltl=cffr
1463 tl_dltl=tl_cffr
1464 END IF
1465!
1466! Compute right and left side values (bR,bL) of parabolic segments
1467! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1468!
1469! NOTE: Although each parabolic segment is monotonic within its grid
1470! box, monotonicity of the whole profile is not guaranteed,
1471! because bL(k+1)-bR(k) may still have different sign than
1472! qc(i,k+1)-qc(i,k). This possibility is excluded,
1473! after bL and bR are reconciled using WENO procedure.
1474!
1475 cff=(dltr-dltl)*hz_inv3(i,k)
1476 tl_cff=(tl_dltr-tl_dltl)*hz_inv3(i,k)+ &
1477 & (dltr-dltl)*tl_hz_inv3(i,k)
1478 dltr=dltr-cff*hz(i,j,k+1)
1479 tl_dltr=tl_dltr-tl_cff*hz(i,j,k+1)-cff*tl_hz(i,j,k+1)
1480 dltl=dltl+cff*hz(i,j,k-1)
1481 tl_dltl=tl_dltl+tl_cff*hz(i,j,k-1)+cff*tl_hz(i,j,k-1)
1482 br(i,k)=qc(i,k)+dltr
1483 tl_br(i,k)=tl_qc(i,k)+tl_dltr
1484 bl(i,k)=qc(i,k)-dltl
1485 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
1486 wr(i,k)=(2.0_r8*dltr-dltl)**2
1487 tl_wr(i,k)=2.0_r8*(2.0_r8*dltr-dltl)* &
1488 & (2.0_r8*tl_dltr-tl_dltl)
1489 wl(i,k)=(dltr-2.0_r8*dltl)**2
1490 tl_wl(i,k)=2.0_r8*(dltr-2.0_r8*dltl)* &
1491 & (tl_dltr-2.0_r8*tl_dltl)
1492 END DO
1493 END DO
1494 cff=1.0e-14_r8
1495 DO k=2,n(ng)-2
1496 DO i=istr,iend
1497 dltl=max(cff,wl(i,k ))
1498 tl_dltl=(0.5_r8-sign(0.5_r8,cff-wl(i,k )))* &
1499 & tl_wl(i,k )
1500 dltr=max(cff,wr(i,k+1))
1501 tl_dltr=(0.5_r8-sign(0.5_r8,cff-wr(i,k+1)))* &
1502 & tl_wr(i,k+1)
1503 br1(i,k)=br(i,k)
1504 bl1(i,k+1)=bl(i,k+1)
1505 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1506 tl_br(i,k)=(tl_dltr*br1(i,k )+dltr*tl_br(i,k )+ &
1507 & tl_dltl*bl1(i,k+1)+dltl*tl_bl(i,k+1))/ &
1508 & (dltr+dltl)- &
1509 & (tl_dltr+tl_dltl)*br(i,k)/(dltr+dltl)
1510 bl(i,k+1)=br(i,k)
1511 tl_bl(i,k+1)=tl_br(i,k)
1512 END DO
1513 END DO
1514 DO i=istr,iend
1515 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1516 tl_fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1517#if defined LINEAR_CONTINUATION
1518 bl(i,n(ng))=br(i,n(ng)-1)
1519 tl_bl(i,n(ng))=tl_br(i,n(ng)-1)
1520 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1521 tl_br(i,n(ng))=2.0_r8*tl_qc(i,n(ng))-tl_bl(i,n(ng))
1522#elif defined NEUMANN
1523 bl(i,n(ng))=br(i,n(ng)-1)
1524 tl_bl(i,n(ng))=tl_br(i,n(ng)-1)
1525 br(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1526 tl_br(i,n(ng))=1.5_r8*tl_qc(i,n(ng))-0.5_r8*tl_bl(i,n(ng))
1527#else
1528 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1529 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1530 br(i,n(ng)-1)=qc(i,n(ng))
1531 tl_br(i,n(ng))=tl_qc(i,n(ng)) ! default strictly monotonic
1532 tl_bl(i,n(ng))=tl_qc(i,n(ng)) ! conditions
1533 tl_br(i,n(ng)-1)=tl_qc(i,n(ng))
1534#endif
1535#if defined LINEAR_CONTINUATION
1536 br(i,1)=bl(i,2)
1537 tl_br(i,1)=tl_bl(i,2)
1538 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1539 tl_bl(i,1)=2.0_r8*tl_qc(i,1)-tl_br(i,1)
1540#elif defined NEUMANN
1541 br(i,1)=bl(i,2)
1542 tl_br(i,1)=tl_bl(i,2)
1543 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1544 tl_bl(i,1)=1.5_r8*tl_qc(i,1)-0.5_r8*tl_br(i,1)
1545#else
1546 bl(i,2)=qc(i,1) ! bottom grid boxes are
1547 br(i,1)=qc(i,1) ! re-assumed to be
1548 bl(i,1)=qc(i,1) ! piecewise constant.
1549 tl_bl(i,2)=tl_qc(i,1) ! bottom grid boxes are
1550 tl_br(i,1)=tl_qc(i,1) ! re-assumed to be
1551 tl_bl(i,1)=tl_qc(i,1) ! piecewise constant.
1552#endif
1553 END DO
1554!
1555! Apply monotonicity constraint again, since the reconciled interfacial
1556! values may cause a non-monotonic behavior of the parabolic segments
1557! inside the grid box.
1558!
1559 DO k=1,n(ng)
1560 DO i=istr,iend
1561 dltr=br(i,k)-qc(i,k)
1562 tl_dltr=tl_br(i,k)-tl_qc(i,k)
1563 dltl=qc(i,k)-bl(i,k)
1564 tl_dltl=tl_qc(i,k)-tl_bl(i,k)
1565 cffr=2.0_r8*dltr
1566 tl_cffr=2.0_r8*tl_dltr
1567 cffl=2.0_r8*dltl
1568 tl_cffl=2.0_r8*tl_dltl
1569 IF ((dltr*dltl).lt.0.0_r8) THEN
1570 dltr=0.0_r8
1571 tl_dltr=0.0_r8
1572 dltl=0.0_r8
1573 tl_dltl=0.0_r8
1574 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1575 dltr=cffl
1576 tl_dltr=tl_cffl
1577 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1578 dltl=cffr
1579 tl_dltl=tl_cffr
1580 END IF
1581 br(i,k)=qc(i,k)+dltr
1582 tl_br(i,k)=tl_qc(i,k)+tl_dltr
1583 bl(i,k)=qc(i,k)-dltl
1584 tl_bl(i,k)=tl_qc(i,k)-tl_dltl
1585 END DO
1586 END DO
1587!
1588! After this moment reconstruction is considered complete. The next
1589! stage is to compute vertical advective fluxes, FC. It is expected
1590! that sinking may occurs relatively fast, the algorithm is designed
1591! to be free of CFL criterion, which is achieved by allowing
1592! integration bounds for semi-Lagrangian advective flux to use as
1593! many grid boxes in upstream direction as necessary.
1594!
1595! In the two code segments below, WL is the z-coordinate of the
1596! departure point for grid box interface z_w with the same indices;
1597! FC is the finite volume flux; ksource(:,k) is index of vertical
1598! grid box which contains the departure point (restricted by N(ng)).
1599! During the search: also add in content of whole grid boxes
1600! participating in FC.
1601!
1602 cff=dtdays*abs(wbio(isink))
1603 tl_cff=dtdays*sign(1.0_r8,wbio(isink))*tl_wbio(isink)
1604 DO k=1,n(ng)
1605 DO i=istr,iend
1606 fc(i,k-1)=0.0_r8
1607 tl_fc(i,k-1)=0.0_r8
1608 wl(i,k)=z_w(i,j,k-1)+cff
1609 tl_wl(i,k)=tl_z_w(i,j,k-1)+tl_cff
1610 wr(i,k)=hz(i,j,k)*qc(i,k)
1611 tl_wr(i,k)=tl_hz(i,j,k)*qc(i,k)+hz(i,j,k)*tl_qc(i,k)
1612 ksource(i,k)=k
1613 END DO
1614 END DO
1615 DO k=1,n(ng)
1616 DO ks=k,n(ng)-1
1617 DO i=istr,iend
1618 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1619 ksource(i,k)=ks+1
1620 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1621 tl_fc(i,k-1)=tl_fc(i,k-1)+tl_wr(i,ks)
1622 END IF
1623 END DO
1624 END DO
1625 END DO
1626!
1627! Finalize computation of flux: add fractional part.
1628!
1629 DO k=1,n(ng)
1630 DO i=istr,iend
1631 ks=ksource(i,k)
1632 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1633 tl_cu=(0.5_r8+sign(0.5_r8, &
1634 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
1635 & hz_inv(i,ks))))* &
1636 & ((tl_wl(i,k)-tl_z_w(i,j,ks-1))*hz_inv(i,ks)+ &
1637 & (wl(i,k)-z_w(i,j,ks-1))*tl_hz_inv(i,ks))
1638 fc(i,k-1)=fc(i,k-1)+ &
1639 & hz(i,j,ks)*cu* &
1640 & (bl(i,ks)+ &
1641 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1642 & (1.5_r8-cu)* &
1643 & (br(i,ks)+bl(i,ks)- &
1644 & 2.0_r8*qc(i,ks))))
1645 tl_fc(i,k-1)=tl_fc(i,k-1)+ &
1646 & (tl_hz(i,j,ks)*cu+hz(i,j,ks)*tl_cu)* &
1647 & (bl(i,ks)+ &
1648 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1649 & (1.5_r8-cu)* &
1650 & (br(i,ks)+bl(i,ks)- &
1651 & 2.0_r8*qc(i,ks))))+ &
1652 & hz(i,j,ks)*cu* &
1653 & (tl_bl(i,ks)+ &
1654 & tl_cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1655 & (1.5_r8-cu)* &
1656 & (br(i,ks)+bl(i,ks)- &
1657 & 2.0_r8*qc(i,ks)))+ &
1658 & cu*(0.5_r8*(tl_br(i,ks)-tl_bl(i,ks))+ &
1659 & tl_cu* &
1660 & (br(i,ks)+bl(i,ks)-2.0_r8*qc(i,ks))- &
1661 & (1.5_r8-cu)* &
1662 & (tl_br(i,ks)+tl_bl(i,ks)- &
1663 & 2.0_r8*tl_qc(i,ks))))
1664 END DO
1665 END DO
1666 DO k=1,n(ng)
1667 DO i=istr,iend
1668 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1669 tl_bio(i,k,ibio)=tl_qc(i,k)+ &
1670 & (tl_fc(i,k)-tl_fc(i,k-1))*hz_inv(i,k)+ &
1671 & (fc(i,k)-fc(i,k-1))*tl_hz_inv(i,k)
1672 END DO
1673 END DO
1674
1675 END DO sink_loop
1676 END DO iter_loop
1677!
1678!-----------------------------------------------------------------------
1679! Update global tracer variables: Add increment due to BGC processes
1680! to tracer array in time index "nnew". Index "nnew" is solution after
1681! advection and mixing and has transport units (m Tunits) hence the
1682! increment is multiplied by Hz. Notice that we need to subtract
1683! original values "Bio_old" at the top of the routine to just account
1684! for the concentractions affected by BGC processes. This also takes
1685! into account any constraints (non-negative concentrations, carbon
1686! concentration range) specified before entering BGC kernel. If "Bio"
1687! were unchanged by BGC processes, the increment would be exactly
1688! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
1689! bounded >=0 so that we can preserve total inventory of nutrients
1690! when advection causes tracer concentration to go negative.
1691!-----------------------------------------------------------------------
1692!
1693 DO itrc=1,nbt
1694 ibio=idbio(itrc)
1695 DO k=1,n(ng)
1696 DO i=istr,iend
1697 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
1698 tl_cff=tl_bio(i,k,ibio)-tl_bio_old(i,k,ibio)
1699!^ t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*Hz(i,j,k)
1700!^
1701 tl_t(i,j,k,nnew,ibio)=tl_t(i,j,k,nnew,ibio)+ &
1702 & tl_cff*hz(i,j,k)+cff*tl_hz(i,j,k)
1703 END DO
1704 END DO
1705 END DO
1706
1707 END DO j_loop
1708!
1709 RETURN
1710 END SUBROUTINE tl_npzd_franks_tile
1711
1712 END MODULE tl_biology_mod
mod_biology
Definition ecosim_mod.h:1
mod_biology::zooga
real(r8), dimension(:), allocatable zooga
Definition npzd_Franks_mod.h:67
mod_biology::tl_wdet
real(r8), dimension(:), allocatable tl_wdet
Definition npzd_Franks_mod.h:59
mod_biology::zoomr
real(r8), dimension(:), allocatable zoomr
Definition fennel_mod.h:162
mod_biology::k_phy
real(r8), dimension(:), allocatable k_phy
Definition fennel_mod.h:135
mod_biology::detrr
real(r8), dimension(:), allocatable detrr
Definition npzd_Franks_mod.h:51
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::zoomd
real(r8), dimension(:), allocatable zoomd
Definition npzd_Franks_mod.h:66
mod_biology::bioiter
integer, dimension(:), allocatable bioiter
Definition ecosim_mod.h:343
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::k_no3
real(r8), dimension(:), allocatable k_no3
Definition fennel_mod.h:133
mod_biology::phymr
real(r8), dimension(:), allocatable phymr
Definition fennel_mod.h:145
mod_biology::iphyt
integer iphyt
Definition fennel_mod.h:81
mod_biology::zooec
real(r8), dimension(:), allocatable zooec
Definition npzd_Franks_mod.h:68
mod_biology::idbio
integer, dimension(:), allocatable idbio
Definition ecosim_mod.h:256
mod_biology::vm_no3
real(r8), dimension(:), allocatable vm_no3
Definition npzd_Franks_mod.h:56
mod_biology::izoop
integer izoop
Definition fennel_mod.h:82
mod_biology::k_ext
real(r8), dimension(:), allocatable k_ext
Definition npzd_Franks_mod.h:52
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::sec2day
real(dp), parameter sec2day
Definition mod_scalars.F:817
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_franks_tile
subroutine tl_npzd_franks_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nnew, rmask, hz, tl_hz, z_r, tl_z_r, z_w, tl_z_w, t, tl_t)
Definition tl_npzd_Franks.h:100
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