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