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