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