ROMS
Loading...
Searching...
No Matches
ad_npzd_Franks.h
Go to the documentation of this file.
1 MODULE ad_biology_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! Nutrient-Phytoplankton-Zooplankton-Detritus Model. !
11! !
12! This routine computes the biological sources and sinks and adds !
13! then the global biological fields. !
14! !
15! Reference: !
16! !
17! Franks et al, 1986: Behavior of simple plankton model with !
18! food-level acclimation by herbivores, Marine Biology, 91, !
19! 121-129. !
20! !
21!=======================================================================
22!
23 implicit none
24!
25 PRIVATE
26 PUBLIC :: ad_biology
27!
28 CONTAINS
29!
30!***********************************************************************
31 SUBROUTINE ad_biology (ng,tile)
32!***********************************************************************
33!
34 USE mod_param
35 USE mod_grid
36 USE mod_ncparam
37 USE mod_ocean
38 USE mod_stepping
39!
40! Imported variable declarations.
41!
42 integer, intent(in) :: ng, tile
43!
44! Local variable declarations.
45!
46 character (len=*), parameter :: myfile = &
47 & __FILE__
48!
49#include "tile.h"
50!
51! Set header file name.
52!
53#ifdef DISTRIBUTE
54 IF (lbiofile(iadm)) THEN
55#else
56 IF (lbiofile(iadm).and.(tile.eq.0)) THEN
57#endif
58 lbiofile(iadm)=.false.
59 bioname(iadm)=myfile
60 END IF
61!
62#ifdef PROFILE
63 CALL wclock_on (ng, iadm, 15, __line__, myfile)
64#endif
65 CALL ad_npzd_franks_tile (ng, tile, &
66 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
67 & imins, imaxs, jmins, jmaxs, &
68 & nstp(ng), nnew(ng), &
69#ifdef MASKING
70 & grid(ng) % rmask, &
71#endif
72 & grid(ng) % Hz, &
73 & grid(ng) % ad_Hz, &
74 & grid(ng) % z_r, &
75 & grid(ng) % ad_z_r, &
76 & grid(ng) % z_w, &
77 & grid(ng) % ad_z_w, &
78 & ocean(ng) % t, &
79 & ocean(ng) % ad_t)
80
81#ifdef PROFILE
82 CALL wclock_off (ng, iadm, 15, __line__, myfile)
83#endif
84 RETURN
85 END SUBROUTINE ad_biology
86!
87!-----------------------------------------------------------------------
88 SUBROUTINE ad_npzd_franks_tile (ng, tile, &
89 & LBi, UBi, LBj, UBj, UBk, UBt, &
90 & IminS, ImaxS, JminS, JmaxS, &
91 & nstp, nnew, &
92#ifdef MASKING
93 & rmask, &
94#endif
95 & Hz, ad_Hz, z_r, ad_z_r, &
96 & z_w, ad_z_w, t, ad_t)
97!-----------------------------------------------------------------------
98!
99 USE mod_param
100 USE mod_biology
101 USE mod_ncparam
102 USE mod_scalars
103!
104! Imported variable declarations.
105!
106 integer, intent(in) :: ng, tile
107 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
108 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
109 integer, intent(in) :: nstp, nnew
110
111#ifdef ASSUMED_SHAPE
112# ifdef MASKING
113 real(r8), intent(in) :: rmask(LBi:,LBj:)
114# endif
115 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
116 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
117 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
118 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
119
120 real(r8), intent(inout) :: ad_Hz(LBi:,LBj:,:)
121 real(r8), intent(inout) :: ad_z_r(LBi:,LBj:,:)
122 real(r8), intent(inout) :: ad_z_w(LBi:,LBj:,0:)
123 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
124#else
125# ifdef MASKING
126 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
127# endif
128 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
129 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
130 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
131 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
132
133 real(r8), intent(inout) :: ad_Hz(LBi:UBi,LBj:UBj,UBk)
134 real(r8), intent(inout) :: ad_z_r(LBi:UBi,LBj:UBj,UBk)
135 real(r8), intent(inout) :: ad_z_w(LBi:UBi,LBj:UBj,0:UBk)
136 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,UBk,3,UBt)
137#endif
138!
139! Local variable declarations.
140!
141 integer, parameter :: Nsink = 1
142
143 integer :: Iter, i, ibio, isink, itrc, itrmx, j, k, ks
144 integer :: Iteradj
145
146 integer, dimension(Nsink) :: idsink
147
148 real(r8), parameter :: eps = 1.0e-16_r8
149
150 real(r8) :: cff, cff1, cff2, cff3, dtdays
151 real(r8) :: ad_cff, ad_cff1
152 real(r8) :: adfac, adfac1, adfac2, adfac3
153 real(r8) :: cffL, cffR, cu, dltL, dltR
154 real(r8) :: ad_cffL, ad_cffR, ad_cu, ad_dltL, ad_dltR
155
156 real(r8), dimension(Nsink) :: Wbio
157 real(r8), dimension(Nsink) :: ad_Wbio
158
159 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
160
161 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
162 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio1
163 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
164
165 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: ad_Bio
166 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: ad_Bio_old
167
168 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
169 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_FC
170
171 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
172 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
173 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
174 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
175 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
176 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
177 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL1
178 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
179 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR1
180 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
181
182 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv
183 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv2
184 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv3
185 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_WL
186 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_WR
187 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_bL
188 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_bR
189 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_qc
190
191#include "set_bounds.h"
192!
193!-----------------------------------------------------------------------
194! Add biological Source/Sink terms.
195!-----------------------------------------------------------------------
196!
197! Avoid computing source/sink terms if no biological iterations.
198!
199 IF (bioiter(ng).le.0) RETURN
200!
201! Set time-stepping according to the number of iterations.
202!
203 dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
204!
205! Set vertical sinking indentification vector.
206!
207 idsink(1)=isdet ! Small detritus
208!
209! Set vertical sinking velocity vector in the same order as the
210! identification vector, IDSINK.
211!
212 wbio(1)=wdet(ng) ! Small detritus
213!
214 j_loop : DO j=jstr,jend
215!
216!-----------------------------------------------------------------------
217! Initialize adjoint private variables.
218!-----------------------------------------------------------------------
219!
220 ad_dltl=0.0_r8
221 ad_dltr=0.0_r8
222 ad_cu=0.0_r8
223 ad_cff=0.0_r8
224 ad_cff1=0.0_r8
225 ad_cffl=0.0_r8
226 ad_cffr=0.0_r8
227 adfac=0.0_r8
228 adfac1=0.0_r8
229 adfac2=0.0_r8
230 adfac3=0.0_r8
231 ad_wdet(ng)=0.0_r8
232 ad_wbio(1)=0.0_r8
233!
234 DO k=1,n(ng)
235 DO i=imins,imaxs
236 ad_hz_inv(i,k)=0.0_r8
237 ad_hz_inv2(i,k)=0.0_r8
238 ad_hz_inv3(i,k)=0.0_r8
239 ad_wl(i,k)=0.0_r8
240 ad_wr(i,k)=0.0_r8
241 ad_bl(i,k)=0.0_r8
242 ad_br(i,k)=0.0_r8
243 ad_qc(i,k)=0.0_r8
244 END DO
245 END DO
246 DO k=0,n(ng)
247 DO i=imins,imaxs
248 ad_fc(i,k)=0.0_r8
249 END DO
250 END DO
251!
252! Clear ad_Bio and Bio arrays.
253!
254 DO itrc=1,nbt
255 ibio=idbio(itrc)
256 DO k=1,n(ng)
257 DO i=istr,iend
258 bio(i,k,ibio)=0.0_r8
259 bio1(i,k,ibio)=0.0_r8
260 ad_bio(i,k,ibio)=0.0_r8
261 ad_bio_old(i,k,ibio)=0.0_r8
262 END DO
263 END DO
264 END DO
265!
266! Compute inverse thickness to avoid repeated divisions.
267!
268 DO k=1,n(ng)
269 DO i=istr,iend
270 hz_inv(i,k)=1.0_r8/hz(i,j,k)
271 END DO
272 END DO
273 DO k=1,n(ng)-1
274 DO i=istr,iend
275 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
276 END DO
277 END DO
278 DO k=2,n(ng)-1
279 DO i=istr,iend
280 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
281 END DO
282 END DO
283!
284! Compute required basic state variables Bio_old and Bio.
285!
286! Extract biological variables from tracer arrays, place them into
287! scratch arrays, and restrict their values to be positive definite.
288! At input, all tracers (index nnew) from predictor step have
289! transport units (m Tunits) since we do not have yet the new
290! values for zeta and Hz. These are known after the 2D barotropic
291! time-stepping. In this routine, this is not a problem because
292! we only use index nstp in the right-hand-side equations.
293!
294 DO itrc=1,nbt
295 ibio=idbio(itrc)
296 DO k=1,n(ng)
297 DO i=istr,iend
298 bio_old(i,k,ibio)=t(i,j,k,nstp,ibio)
299 END DO
300 END DO
301 END DO
302!
303! Determine Correction for negativity.
304!
305 DO k=1,n(ng)
306 DO i=istr,iend
307 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
308 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
309 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
310 & max(0.0_r8,eps-bio_old(i,k,isdet))
311!
312! If correction needed, determine the largest pool to debit.
313!
314 IF (cff1.gt.0.0) THEN
315 itrmx=idbio(1)
316 cff=t(i,j,k,nstp,itrmx)
317 DO ibio=idbio(2),idbio(nbt)
318 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
319 itrmx=ibio
320 cff=t(i,j,k,nstp,ibio)
321 END IF
322 END DO
323!
324! Update new values.
325!
326 DO itrc=1,nbt
327 ibio=idbio(itrc)
328 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
329 & cff1*(sign(0.5_r8, &
330 & real(itrmx-ibio,r8)**2)+ &
331 & sign(0.5_r8, &
332 & -real(itrmx-ibio,r8)**2))
333 END DO
334 ELSE
335 DO itrc=1,nbt
336 ibio=idbio(itrc)
337 bio(i,k,ibio)=bio_old(i,k,ibio)
338 END DO
339 END IF
340 END DO
341 END DO
342!
343!=======================================================================
344! Start internal iterations to achieve convergence of the nonlinear
345! backward-implicit solution.
346!=======================================================================
347!
348! During the iterative procedure a series of fractional time steps are
349! performed in a chained mode (splitting by different biological
350! conversion processes) in sequence of the main food chain. In all
351! stages the concentration of the component being consumed is treated
352! in a fully implicit manner, so the algorithm guarantees non-negative
353! values, no matter how strong the concentration of active consuming
354! component (Phytoplankton or Zooplankton). The overall algorithm,
355! as well as any stage of it, is formulated in conservative form
356! (except explicit sinking) in sense that the sum of concentration of
357! all components is conserved.
358!
359! In the implicit algorithm, we have for example (N: nutrient,
360! P: phytoplankton),
361!
362! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
363! {Michaelis-Menten}
364! below, we set
365! The N in the numerator of
366! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
367! as N(new)
368!
369! so the time-stepping of the equations becomes:
370!
371! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
372! consuming, divide by (1 + cff)
373! and
374!
375! P(new) = P(old) + cff * N(new) (2) when adding a source term,
376! growing, add (cff * source)
377!
378! Notice that if you substitute (1) in (2), you will get:
379!
380! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
381!
382! If you add (1) and (3), you get
383!
384! N(new) + P(new) = N(old) + P(old)
385!
386! implying conservation regardless how "cff" is computed. Therefore,
387! this scheme is unconditionally stable regardless of the conversion
388! rate. It does not generate negative values since the constituent
389! to be consumed is always treated implicitly. It is also biased
390! toward damping oscillations.
391!
392! The iterative loop below is to iterate toward an universal Backward-
393! Euler treatment of all terms. So if there are oscillations in the
394! system, they are only physical oscillations. These iterations,
395! however, do not improve the accuaracy of the solution.
396!
397 DO iter=1,bioiter(ng)
398!
399! Nutrient uptake by phytoplankton.
400!
401 cff1=dtdays*vm_no3(ng)
402 DO k=1,n(ng)
403 DO i=istr,iend
404 cff=bio(i,k,iphyt)* &
405 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
406 & (k_no3(ng)+bio(i,k,ino3_))
407 bio(i,k,ino3_)=bio(i,k,ino3_)/ &
408 & (1.0_r8+cff)
409 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
410 & bio(i,k,ino3_)*cff
411 END DO
412 END DO
413!
414! Phytoplankton grazing by Zooplankton and mortality to Detritus
415! (rate: PhyMR).
416!
417 cff1=dtdays*zoogr(ng)
418 cff2=dtdays*phymr(ng)
419 cff3=k_phy(ng)*k_phy(ng)
420 DO k=1,n(ng)
421 DO i=istr,iend
422 cff=bio(i,k,izoop)*bio(i,k,iphyt)*cff1/ &
423 & (cff3+bio(i,k,iphyt)*bio(i,k,iphyt))
424 bio(i,k,iphyt)=bio(i,k,iphyt)/ &
425 & (1.0_r8+cff+cff2)
426 bio(i,k,izoop)=bio(i,k,izoop)+ &
427 & bio(i,k,iphyt)*cff*(1.0_r8-zooga(ng))
428 bio(i,k,isdet)=bio(i,k,isdet)+ &
429 & bio(i,k,iphyt)* &
430 & (cff2+cff*(zooga(ng)-zooec(ng)))
431 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
432 & bio(i,k,iphyt)*cff*zooec(ng)
433 END DO
434 END DO
435!
436! Zooplankton excretion to nutrients and mortality to Detritus.
437!
438 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
439 cff2=dtdays*zoomr(ng)
440 cff3=dtdays*zoomd(ng)
441 DO k=1,n(ng)
442 DO i=istr,iend
443 bio(i,k,izoop)=bio(i,k,izoop)*cff1
444 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
445 & bio(i,k,izoop)*cff2
446 bio(i,k,isdet)=bio(i,k,isdet)+ &
447 & bio(i,k,izoop)*cff3
448 END DO
449 END DO
450!
451! Detritus breakdown to nutrients.
452!
453 cff1=dtdays*detrr(ng)
454 cff2=1.0_r8/(1.0_r8+cff1)
455 DO k=1,n(ng)
456 DO i=istr,iend
457 bio(i,k,isdet)=bio(i,k,isdet)*cff2
458 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
459 & bio(i,k,isdet)*cff1
460 END DO
461 END DO
462!
463!-----------------------------------------------------------------------
464! Vertical sinking terms.
465!-----------------------------------------------------------------------
466!
467! Reconstruct vertical profile of selected biological constituents
468! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
469! grid box. Then, compute semi-Lagrangian flux due to sinking.
470!
471 DO isink=1,nsink
472 ibio=idsink(isink)
473!
474! Copy concentration of biological particulates into scratch array
475! "qc" (q-central, restrict it to be positive) which is hereafter
476! interpreted as a set of grid-box averaged values for biogeochemical
477! constituent concentration.
478!
479 DO k=1,n(ng)
480 DO i=istr,iend
481 qc(i,k)=bio(i,k,ibio)
482 END DO
483 END DO
484!
485 DO k=n(ng)-1,1,-1
486 DO i=istr,iend
487 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
488 END DO
489 END DO
490 DO k=2,n(ng)-1
491 DO i=istr,iend
492 dltr=hz(i,j,k)*fc(i,k)
493 dltl=hz(i,j,k)*fc(i,k-1)
494 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
495 cffr=cff*fc(i,k)
496 cffl=cff*fc(i,k-1)
497!
498! Apply PPM monotonicity constraint to prevent oscillations within the
499! grid box.
500!
501 IF ((dltr*dltl).le.0.0_r8) THEN
502 dltr=0.0_r8
503 dltl=0.0_r8
504 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
505 dltr=cffl
506 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
507 dltl=cffr
508 END IF
509!
510! Compute right and left side values (bR,bL) of parabolic segments
511! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
512!
513! NOTE: Although each parabolic segment is monotonic within its grid
514! box, monotonicity of the whole profile is not guaranteed,
515! because bL(k+1)-bR(k) may still have different sign than
516! qc(i,k+1)-qc(i,k). This possibility is excluded,
517! after bL and bR are reconciled using WENO procedure.
518!
519 cff=(dltr-dltl)*hz_inv3(i,k)
520 dltr=dltr-cff*hz(i,j,k+1)
521 dltl=dltl+cff*hz(i,j,k-1)
522 br(i,k)=qc(i,k)+dltr
523 bl(i,k)=qc(i,k)-dltl
524 wr(i,k)=(2.0_r8*dltr-dltl)**2
525 wl(i,k)=(dltr-2.0_r8*dltl)**2
526 END DO
527 END DO
528 cff=1.0e-14_r8
529 DO k=2,n(ng)-2
530 DO i=istr,iend
531 dltl=max(cff,wl(i,k ))
532 dltr=max(cff,wr(i,k+1))
533 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
534 bl(i,k+1)=br(i,k)
535 END DO
536 END DO
537 DO i=istr,iend
538 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
539#if defined LINEAR_CONTINUATION
540 bl(i,n(ng))=br(i,n(ng)-1)
541 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
542#elif defined NEUMANN
543 bl(i,n(ng))=br(i,n(ng)-1)
544 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
545#else
546 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
547 bl(i,n(ng))=qc(i,n(ng)) ! conditions
548 br(i,n(ng)-1)=qc(i,n(ng))
549#endif
550#if defined LINEAR_CONTINUATION
551 br(i,1)=bl(i,2)
552 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
553#elif defined NEUMANN
554 br(i,1)=bl(i,2)
555 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
556#else
557 bl(i,2)=qc(i,1) ! bottom grid boxes are
558 br(i,1)=qc(i,1) ! re-assumed to be
559 bl(i,1)=qc(i,1) ! piecewise constant.
560#endif
561 END DO
562!
563! Apply monotonicity constraint again, since the reconciled interfacial
564! values may cause a non-monotonic behavior of the parabolic segments
565! inside the grid box.
566!
567 DO k=1,n(ng)
568 DO i=istr,iend
569 dltr=br(i,k)-qc(i,k)
570 dltl=qc(i,k)-bl(i,k)
571 cffr=2.0_r8*dltr
572 cffl=2.0_r8*dltl
573 IF ((dltr*dltl).lt.0.0_r8) THEN
574 dltr=0.0_r8
575 dltl=0.0_r8
576 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
577 dltr=cffl
578 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
579 dltl=cffr
580 END IF
581 br(i,k)=qc(i,k)+dltr
582 bl(i,k)=qc(i,k)-dltl
583 END DO
584 END DO
585!
586! After this moment reconstruction is considered complete. The next
587! stage is to compute vertical advective fluxes, FC. It is expected
588! that sinking may occurs relatively fast, the algorithm is designed
589! to be free of CFL criterion, which is achieved by allowing
590! integration bounds for semi-Lagrangian advective flux to use as
591! many grid boxes in upstream direction as necessary.
592!
593! In the two code segments below, WL is the z-coordinate of the
594! departure point for grid box interface z_w with the same indices;
595! FC is the finite volume flux; ksource(:,k) is index of vertical
596! grid box which contains the departure point (restricted by N(ng)).
597! During the search: also add in content of whole grid boxes
598! participating in FC.
599!
600 cff=dtdays*abs(wbio(isink))
601 DO k=1,n(ng)
602 DO i=istr,iend
603 fc(i,k-1)=0.0_r8
604 wl(i,k)=z_w(i,j,k-1)+cff
605 wr(i,k)=hz(i,j,k)*qc(i,k)
606 ksource(i,k)=k
607 END DO
608 END DO
609 DO k=1,n(ng)
610 DO ks=k,n(ng)-1
611 DO i=istr,iend
612 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
613 ksource(i,k)=ks+1
614 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
615 END IF
616 END DO
617 END DO
618 END DO
619!
620! Finalize computation of flux: add fractional part.
621!
622 DO k=1,n(ng)
623 DO i=istr,iend
624 ks=ksource(i,k)
625 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
626 fc(i,k-1)=fc(i,k-1)+ &
627 & hz(i,j,ks)*cu* &
628 & (bl(i,ks)+ &
629 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
630 & (1.5_r8-cu)* &
631 & (br(i,ks)+bl(i,ks)- &
632 & 2.0_r8*qc(i,ks))))
633 END DO
634 END DO
635 DO k=1,n(ng)
636 DO i=istr,iend
637 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
638 END DO
639 END DO
640
641 END DO
642 END DO
643!
644!-----------------------------------------------------------------------
645! Update global tracer variables: Add increment due to BGC processes
646! to tracer array in time index "nnew". Index "nnew" is solution after
647! advection and mixing and has transport units (m Tunits) hence the
648! increment is multiplied by Hz. Notice that we need to subtract
649! original values "Bio_old" at the top of the routine to just account
650! for the concentractions affected by BGC processes. This also takes
651! into account any constraints (non-negative concentrations, carbon
652! concentration range) specified before entering BGC kernel. If "Bio"
653! were unchanged by BGC processes, the increment would be exactly
654! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
655! bounded >=0 so that we can preserve total inventory of nutrients
656! when advection causes tracer concentration to go negative.
657!-----------------------------------------------------------------------
658!
659 DO itrc=1,nbt
660 ibio=idbio(itrc)
661 DO k=1,n(ng)
662 DO i=istr,iend
663 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
664!^ tl_t(i,j,k,nnew,ibio)=tl_t(i,j,k,nnew,ibio)+ &
665!^ & tl_cff*Hz(i,j,k)+cff*tl_Hz(i,j,k)
666!^
667 ad_hz(i,j,k)=ad_hz(i,j,k)+cff*ad_t(i,j,k,nnew,ibio)
668 ad_cff=ad_cff+hz(i,j,k)*ad_t(i,j,k,nnew,ibio)
669!^ tl_cff=tl_Bio(i,k,ibio)-tl_Bio_old(i,k,ibio)
670!^
671 ad_bio_old(i,k,ibio)=ad_bio_old(i,k,ibio)-ad_cff
672 ad_bio(i,k,ibio)=ad_bio(i,k,ibio)+ad_cff
673 ad_cff=0.0_r8
674 END DO
675 END DO
676 END DO
677!
678!=======================================================================
679! Start adjoint internal iterations to achieve convergence of the
680! nonlinear backward-implicit solution.
681!=======================================================================
682!
683! During the iterative procedure a series of fractional time steps are
684! performed in a chained mode (splitting by different biological
685! conversion processes) in sequence of the main food chain. In all
686! stages the concentration of the component being consumed is treated
687! in fully implicit manner, so the algorithm guarantees non-negative
688! values, no matter how strong s the concentration of active consuming
689! component (Phytoplankton or Zooplankton). The overall algorithm,
690! as well as any stage of it, is formulated in conservative form
691! (except explicit sinking) in sense that the sum of concentration of
692! all components is conserved.
693!
694 iter_loop: DO iter=bioiter(ng),1,-1
695!
696!-----------------------------------------------------------------------
697! Adjoint vertical sinking terms.
698!-----------------------------------------------------------------------
699!
700! Reconstruct vertical profile of selected biological constituents
701! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
702! grid box. Then, compute semi-Lagrangian flux due to sinking.
703!
704!
705! Compute appropriate basic state arrays III.
706!
707! Determine Correction for negativity.
708!
709 DO k=1,n(ng)
710 DO i=istr,iend
711 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
712 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
713 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
714 & max(0.0_r8,eps-bio_old(i,k,isdet))
715!
716! If correction needed, determine the largest pool to debit.
717!
718 IF (cff1.gt.0.0) THEN
719 itrmx=idbio(1)
720 cff=t(i,j,k,nstp,itrmx)
721 DO ibio=idbio(2),idbio(nbt)
722 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
723 itrmx=ibio
724 cff=t(i,j,k,nstp,ibio)
725 END IF
726 END DO
727!
728! Update new values.
729!
730 DO itrc=1,nbt
731 ibio=idbio(itrc)
732 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
733 & cff1*(sign(0.5_r8, &
734 & real(itrmx-ibio,r8)**2)+ &
735 & sign(0.5_r8, &
736 & -real(itrmx-ibio,r8)**2))
737 END DO
738 ELSE
739 DO itrc=1,nbt
740 ibio=idbio(itrc)
741 bio(i,k,ibio)=bio_old(i,k,ibio)
742 END DO
743 END IF
744 END DO
745 END DO
746!
747!=======================================================================
748! Start internal iterations to achieve convergence of the nonlinear
749! backward-implicit solution.
750!=======================================================================
751!
752 DO iteradj=1,iter
753!
754! Nutrient uptake by phytoplankton.
755!
756 cff1=dtdays*vm_no3(ng)
757 DO k=1,n(ng)
758 DO i=istr,iend
759 cff=bio(i,k,iphyt)* &
760 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
761 & (k_no3(ng)+bio(i,k,ino3_))
762 bio1(i,k,ino3_)=bio(i,k,ino3_)
763 bio(i,k,ino3_)=bio(i,k,ino3_)/ &
764 & (1.0_r8+cff)
765 bio1(i,k,iphyt)=bio(i,k,iphyt)
766 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
767 & bio(i,k,ino3_)*cff
768 END DO
769 END DO
770!
771! Phytoplankton grazing by Zooplankton and mortality to Detritus
772! (rate: PhyMR).
773!
774 cff1=dtdays*zoogr(ng)
775 cff2=dtdays*phymr(ng)
776 cff3=k_phy(ng)*k_phy(ng)
777 DO k=1,n(ng)
778 DO i=istr,iend
779 cff=bio(i,k,izoop)*bio(i,k,iphyt)*cff1/ &
780 & (cff3+bio(i,k,iphyt)*bio(i,k,iphyt))
781 bio1(i,k,iphyt)=bio(i,k,iphyt)
782 bio(i,k,iphyt)=bio(i,k,iphyt)/ &
783 & (1.0_r8+cff+cff2)
784 bio1(i,k,izoop)=bio(i,k,izoop)
785 bio(i,k,izoop)=bio(i,k,izoop)+ &
786 & bio(i,k,iphyt)*cff*(1.0_r8-zooga(ng))
787 bio(i,k,isdet)=bio(i,k,isdet)+ &
788 & bio(i,k,iphyt)* &
789 & (cff2+cff*(zooga(ng)-zooec(ng)))
790 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
791 & bio(i,k,iphyt)*cff*zooec(ng)
792 END DO
793 END DO
794!
795! Zooplankton excretion to nutrients and mortality to Detritus.
796!
797 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
798 cff2=dtdays*zoomr(ng)
799 cff3=dtdays*zoomd(ng)
800 DO k=1,n(ng)
801 DO i=istr,iend
802 bio(i,k,izoop)=bio(i,k,izoop)*cff1
803 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
804 & bio(i,k,izoop)*cff2
805 bio(i,k,isdet)=bio(i,k,isdet)+ &
806 & bio(i,k,izoop)*cff3
807 END DO
808 END DO
809!
810! Detritus breakdown to nutrients.
811!
812 cff1=dtdays*detrr(ng)
813 cff2=1.0_r8/(1.0_r8+cff1)
814 DO k=1,n(ng)
815 DO i=istr,iend
816 bio(i,k,isdet)=bio(i,k,isdet)*cff2
817 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
818 & bio(i,k,isdet)*cff1
819 END DO
820 END DO
821!
822 IF (iteradj.ne.iter) THEN
823!
824!-----------------------------------------------------------------------
825! Vertical sinking terms.
826!-----------------------------------------------------------------------
827!
828! Reconstruct vertical profile of selected biological constituents
829! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
830! grid box. Then, compute semi-Lagrangian flux due to sinking.
831!
832 DO isink=1,nsink
833 ibio=idsink(isink)
834!
835! Copy concentration of biological particulates into scratch array
836! "qc" (q-central, restrict it to be positive) which is hereafter
837! interpreted as a set of grid-box averaged values for biogeochemical
838! constituent concentration.
839!
840 DO k=1,n(ng)
841 DO i=istr,iend
842 qc(i,k)=bio(i,k,ibio)
843 END DO
844 END DO
845!
846 DO k=n(ng)-1,1,-1
847 DO i=istr,iend
848 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
849 END DO
850 END DO
851 DO k=2,n(ng)-1
852 DO i=istr,iend
853 dltr=hz(i,j,k)*fc(i,k)
854 dltl=hz(i,j,k)*fc(i,k-1)
855 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
856 cffr=cff*fc(i,k)
857 cffl=cff*fc(i,k-1)
858!
859! Apply PPM monotonicity constraint to prevent oscillations within the
860! grid box.
861!
862 IF ((dltr*dltl).le.0.0_r8) THEN
863 dltr=0.0_r8
864 dltl=0.0_r8
865 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
866 dltr=cffl
867 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
868 dltl=cffr
869 END IF
870!
871! Compute right and left side values (bR,bL) of parabolic segments
872! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
873!
874! NOTE: Although each parabolic segment is monotonic within its grid
875! box, monotonicity of the whole profile is not guaranteed,
876! because bL(k+1)-bR(k) may still have different sign than
877! qc(i,k+1)-qc(i,k). This possibility is excluded,
878! after bL and bR are reconciled using WENO procedure.
879!
880 cff=(dltr-dltl)*hz_inv3(i,k)
881 dltr=dltr-cff*hz(i,j,k+1)
882 dltl=dltl+cff*hz(i,j,k-1)
883 br(i,k)=qc(i,k)+dltr
884 bl(i,k)=qc(i,k)-dltl
885 wr(i,k)=(2.0_r8*dltr-dltl)**2
886 wl(i,k)=(dltr-2.0_r8*dltl)**2
887 END DO
888 END DO
889 cff=1.0e-14_r8
890 DO k=2,n(ng)-2
891 DO i=istr,iend
892 dltl=max(cff,wl(i,k ))
893 dltr=max(cff,wr(i,k+1))
894 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
895 bl(i,k+1)=br(i,k)
896 END DO
897 END DO
898 DO i=istr,iend
899 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
900#if defined LINEAR_CONTINUATION
901 bl(i,n(ng))=br(i,n(ng)-1)
902 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
903#elif defined NEUMANN
904 bl(i,n(ng))=br(i,n(ng)-1)
905 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
906#else
907 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
908 bl(i,n(ng))=qc(i,n(ng)) ! conditions
909 br(i,n(ng)-1)=qc(i,n(ng))
910#endif
911#if defined LINEAR_CONTINUATION
912 br(i,1)=bl(i,2)
913 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
914#elif defined NEUMANN
915 br(i,1)=bl(i,2)
916 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
917#else
918 bl(i,2)=qc(i,1) ! bottom grid boxes are
919 br(i,1)=qc(i,1) ! re-assumed to be
920 bl(i,1)=qc(i,1) ! piecewise constant.
921#endif
922 END DO
923!
924! Apply monotonicity constraint again, since the reconciled interfacial
925! values may cause a non-monotonic behavior of the parabolic segments
926! inside the grid box.
927!
928 DO k=1,n(ng)
929 DO i=istr,iend
930 dltr=br(i,k)-qc(i,k)
931 dltl=qc(i,k)-bl(i,k)
932 cffr=2.0_r8*dltr
933 cffl=2.0_r8*dltl
934 IF ((dltr*dltl).lt.0.0_r8) THEN
935 dltr=0.0_r8
936 dltl=0.0_r8
937 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
938 dltr=cffl
939 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
940 dltl=cffr
941 END IF
942 br(i,k)=qc(i,k)+dltr
943 bl(i,k)=qc(i,k)-dltl
944 END DO
945 END DO
946!
947! After this moment reconstruction is considered complete. The next
948! stage is to compute vertical advective fluxes, FC. It is expected
949! that sinking may occurs relatively fast, the algorithm is designed
950! to be free of CFL criterion, which is achieved by allowing
951! integration bounds for semi-Lagrangian advective flux to use as
952! many grid boxes in upstream direction as necessary.
953!
954! In the two code segments below, WL is the z-coordinate of the
955! departure point for grid box interface z_w with the same indices;
956! FC is the finite volume flux; ksource(:,k) is index of vertical
957! grid box which contains the departure point (restricted by N(ng)).
958! During the search: also add in content of whole grid boxes
959! participating in FC.
960!
961 cff=dtdays*abs(wbio(isink))
962 DO k=1,n(ng)
963 DO i=istr,iend
964 fc(i,k-1)=0.0_r8
965 wl(i,k)=z_w(i,j,k-1)+cff
966 wr(i,k)=hz(i,j,k)*qc(i,k)
967 ksource(i,k)=k
968 END DO
969 END DO
970 DO k=1,n(ng)
971 DO ks=k,n(ng)-1
972 DO i=istr,iend
973 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
974 ksource(i,k)=ks+1
975 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
976 END IF
977 END DO
978 END DO
979 END DO
980!
981! Finalize computation of flux: add fractional part.
982!
983 DO k=1,n(ng)
984 DO i=istr,iend
985 ks=ksource(i,k)
986 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
987 fc(i,k-1)=fc(i,k-1)+ &
988 & hz(i,j,ks)*cu* &
989 & (bl(i,ks)+ &
990 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
991 & (1.5_r8-cu)* &
992 & (br(i,ks)+bl(i,ks)- &
993 & 2.0_r8*qc(i,ks))))
994 END DO
995 END DO
996 DO k=1,n(ng)
997 DO i=istr,iend
998 bio(i,k,ibio)=qc(i,k)+ &
999 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1000 END DO
1001 END DO
1002 END DO
1003 END IF
1004 END DO
1005!
1006! End compute basic state arrays III.
1007!
1008!
1009 sink_loop: DO isink=1,nsink
1010 ibio=idsink(isink)
1011!
1012! Compute required flux arrays.
1013!
1014! Copy concentration of biological particulates into scratch array
1015! "qc" (q-central, restrict it to be positive) which is hereafter
1016! interpreted as a set of grid-box averaged values for biogeochemical
1017! constituent concentration.
1018!
1019 DO k=1,n(ng)
1020 DO i=istr,iend
1021 qc(i,k)=bio(i,k,ibio)
1022 END DO
1023 END DO
1024!
1025 DO k=n(ng)-1,1,-1
1026 DO i=istr,iend
1027 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1028 END DO
1029 END DO
1030 DO k=2,n(ng)-1
1031 DO i=istr,iend
1032 dltr=hz(i,j,k)*fc(i,k)
1033 dltl=hz(i,j,k)*fc(i,k-1)
1034 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1035 cffr=cff*fc(i,k)
1036 cffl=cff*fc(i,k-1)
1037!
1038! Apply PPM monotonicity constraint to prevent oscillations within the
1039! grid box.
1040!
1041 IF ((dltr*dltl).le.0.0_r8) THEN
1042 dltr=0.0_r8
1043 dltl=0.0_r8
1044 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1045 dltr=cffl
1046 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1047 dltl=cffr
1048 END IF
1049!
1050! Compute right and left side values (bR,bL) of parabolic segments
1051! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1052!
1053! NOTE: Although each parabolic segment is monotonic within its grid
1054! box, monotonicity of the whole profile is not guaranteed,
1055! because bL(k+1)-bR(k) may still have different sign than
1056! qc(i,k+1)-qc(i,k). This possibility is excluded,
1057! after bL and bR are reconciled using WENO procedure.
1058!
1059 cff=(dltr-dltl)*hz_inv3(i,k)
1060 dltr=dltr-cff*hz(i,j,k+1)
1061 dltl=dltl+cff*hz(i,j,k-1)
1062 br(i,k)=qc(i,k)+dltr
1063 bl(i,k)=qc(i,k)-dltl
1064 wr(i,k)=(2.0_r8*dltr-dltl)**2
1065 wl(i,k)=(dltr-2.0_r8*dltl)**2
1066 END DO
1067 END DO
1068 cff=1.0e-14_r8
1069 DO k=2,n(ng)-2
1070 DO i=istr,iend
1071 dltl=max(cff,wl(i,k ))
1072 dltr=max(cff,wr(i,k+1))
1073 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1074 bl(i,k+1)=br(i,k)
1075 END DO
1076 END DO
1077 DO i=istr,iend
1078 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1079#if defined LINEAR_CONTINUATION
1080 bl(i,n(ng))=br(i,n(ng)-1)
1081 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1082#elif defined NEUMANN
1083 bl(i,n(ng))=br(i,n(ng)-1)
1084 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1085#else
1086 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1087 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1088 br(i,n(ng)-1)=qc(i,n(ng))
1089#endif
1090#if defined LINEAR_CONTINUATION
1091 br(i,1)=bl(i,2)
1092 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1093#elif defined NEUMANN
1094 br(i,1)=bl(i,2)
1095 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1096#else
1097 bl(i,2)=qc(i,1) ! bottom grid boxes are
1098 br(i,1)=qc(i,1) ! re-assumed to be
1099 bl(i,1)=qc(i,1) ! piecewise constant.
1100#endif
1101 END DO
1102!
1103! Apply monotonicity constraint again, since the reconciled interfacial
1104! values may cause a non-monotonic behavior of the parabolic segments
1105! inside the grid box.
1106!
1107 DO k=1,n(ng)
1108 DO i=istr,iend
1109 dltr=br(i,k)-qc(i,k)
1110 dltl=qc(i,k)-bl(i,k)
1111 cffr=2.0_r8*dltr
1112 cffl=2.0_r8*dltl
1113 IF ((dltr*dltl).lt.0.0_r8) THEN
1114 dltr=0.0_r8
1115 dltl=0.0_r8
1116 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1117 dltr=cffl
1118 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1119 dltl=cffr
1120 END IF
1121 br(i,k)=qc(i,k)+dltr
1122 bl(i,k)=qc(i,k)-dltl
1123 END DO
1124 END DO
1125!
1126! After this moment reconstruction is considered complete. The next
1127! stage is to compute vertical advective fluxes, FC. It is expected
1128! that sinking may occurs relatively fast, the algorithm is designed
1129! to be free of CFL criterion, which is achieved by allowing
1130! integration bounds for semi-Lagrangian advective flux to use as
1131! many grid boxes in upstream direction as necessary.
1132!
1133! In the two code segments below, WL is the z-coordinate of the
1134! departure point for grid box interface z_w with the same indices;
1135! FC is the finite volume flux; ksource(:,k) is index of vertical
1136! grid box which contains the departure point (restricted by N(ng)).
1137! During the search: also add in content of whole grid boxes
1138! participating in FC.
1139!
1140 cff=dtdays*abs(wbio(isink))
1141 DO k=1,n(ng)
1142 DO i=istr,iend
1143 fc(i,k-1)=0.0_r8
1144 wl(i,k)=z_w(i,j,k-1)+cff
1145 wr(i,k)=hz(i,j,k)*qc(i,k)
1146 ksource(i,k)=k
1147 END DO
1148 END DO
1149 DO k=1,n(ng)
1150 DO ks=k,n(ng)-1
1151 DO i=istr,iend
1152 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1153 ksource(i,k)=ks+1
1154 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1155 END IF
1156 END DO
1157 END DO
1158 END DO
1159!
1160! Finalize computation of flux: add fractional part.
1161!
1162 DO k=1,n(ng)
1163 DO i=istr,iend
1164 ks=ksource(i,k)
1165 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1166 fc(i,k-1)=fc(i,k-1)+ &
1167 & hz(i,j,ks)*cu* &
1168 & (bl(i,ks)+ &
1169 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1170 & (1.5_r8-cu)* &
1171 & (br(i,ks)+bl(i,ks)- &
1172 & 2.0_r8*qc(i,ks))))
1173 END DO
1174 END DO
1175 DO k=1,n(ng)
1176 DO i=istr,iend
1177!^ tl_Bio(i,k,ibio)=tl_qc(i,k)+ &
1178!^ & (tl_FC(i,k)-tl_FC(i,k-1))*Hz_inv(i,k)+ &
1179!^ & (FC(i,k)-FC(i,k-1))*tl_Hz_inv(i,k)
1180!^
1181 adfac=hz_inv(i,k)*ad_bio(i,k,ibio)
1182 ad_fc(i,k-1)=ad_fc(i,k-1)-adfac
1183 ad_fc(i,k )=ad_fc(i,k )+adfac
1184 ad_hz_inv(i,k)=ad_hz_inv(i,k)+ &
1185 & (fc(i,k)-fc(i,k-1))*ad_bio(i,k,ibio)
1186 ad_qc(i,k)=ad_qc(i,k)+ad_bio(i,k,ibio)
1187 ad_bio(i,k,ibio)=0.0_r8
1188 END DO
1189 END DO
1190!
1191! Adjoint of final computation of flux: add fractional part.
1192!
1193 DO k=1,n(ng)
1194 DO i=istr,iend
1195 ks=ksource(i,k)
1196 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1197!^ tl_FC(i,k-1)=tl_FC(i,k-1)+ &
1198!^ & (tl_Hz(i,j,ks)*cu+Hz(i,j,ks)*tl_cu)* &
1199!^ & (bL(i,ks)+ &
1200!^ & cu*(0.5_r8*(bR(i,ks)-bL(i,ks))- &
1201!^ & (1.5_r8-cu)* &
1202!^ & (bR(i,ks)+bL(i,ks)- &
1203!^ & 2.0_r8*qc(i,ks))))+ &
1204!^ & Hz(i,j,ks)*cu* &
1205!^ & (tl_bL(i,ks)+ &
1206!^ & tl_cu*(0.5_r8*(bR(i,ks)-bL(i,ks))- &
1207!^ & (1.5_r8-cu)* &
1208!^ & (bR(i,ks)+bL(i,ks)- &
1209!^ & 2.0_r8*qc(i,ks)))+ &
1210!^ & cu*(0.5_r8*(tl_bR(i,ks)-tl_bL(i,ks))+ &
1211!^ & tl_cu* &
1212!^ & (bR(i,ks)+bL(i,ks)-2.0_r8*qc(i,ks))- &
1213!^ & (1.5_r8-cu)* &
1214!^ & (tl_bR(i,ks)+tl_bL(i,ks)- &
1215!^ & 2.0_r8*tl_qc(i,ks))))
1216!^
1217 adfac=(bl(i,ks)+ &
1218 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1219 & (1.5_r8-cu)* &
1220 & (br(i,ks)+bl(i,ks)- &
1221 & 2.0_r8*qc(i,ks))))*ad_fc(i,k-1)
1222 adfac1=hz(i,j,ks)*cu*ad_fc(i,k-1)
1223 adfac2=adfac1*cu
1224 adfac3=adfac2*(1.5_r8-cu)
1225 ad_hz(i,j,ks)=ad_hz(i,j,ks)+cu*adfac
1226 ad_cu=ad_cu+hz(i,j,ks)*adfac
1227 ad_bl(i,ks)=ad_bl(i,ks)+adfac1
1228 ad_cu=ad_cu+ &
1229 & adfac1*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1230 & (1.5_r8-cu)* &
1231 & (br(i,ks)+bl(i,ks)- &
1232 & 2.0_r8*qc(i,ks)))+ &
1233 & adfac2*(br(i,ks)+bl(i,ks)-2.0_r8*qc(i,ks))
1234 ad_br(i,ks)=ad_br(i,ks)+0.5_r8*adfac2-adfac3
1235 ad_bl(i,ks)=ad_bl(i,ks)-0.5_r8*adfac2-adfac3
1236 ad_qc(i,ks)=ad_qc(i,ks)+2.0_r8*adfac3
1237!^ tl_cu=(0.5_r8+SIGN(0.5_r8, &
1238!^ & (1.0_r8-(WL(i,k)-z_w(i,j,ks-1))* &
1239!^ & Hz_inv(i,ks))))* &
1240!^ & ((tl_WL(i,k)-tl_z_w(i,j,ks-1))*Hz_inv(i,ks)+ &
1241!^ & (WL(i,k)-z_w(i,j,ks-1))*tl_Hz_inv(i,ks))
1242!^
1243 adfac=(0.5_r8+sign(0.5_r8, &
1244 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
1245 & hz_inv(i,ks))))*ad_cu
1246 adfac1=adfac*hz_inv(i,ks)
1247 ad_wl(i,k)=ad_wl(i,k)+adfac1
1248 ad_z_w(i,j,ks-1)=ad_z_w(i,j,ks-1)-adfac1
1249 ad_hz_inv(i,ks)=ad_hz_inv(i,ks)+ &
1250 & (wl(i,k)-z_w(i,j,ks-1))*adfac
1251 ad_cu=0.0_r8
1252 END DO
1253 END DO
1254!
1255! After this moment reconstruction is considered complete. The next
1256! stage is to compute vertical advective fluxes, FC. It is expected
1257! that sinking may occurs relatively fast, the algorithm is designed
1258! to be free of CFL criterion, which is achieved by allowing
1259! integration bounds for semi-Lagrangian advective flux to use as
1260! many grid boxes in upstream direction as necessary.
1261!
1262! In the two code segments below, WL is the z-coordinate of the
1263! departure point for grid box interface z_w with the same indices;
1264! FC is the finite volume flux; ksource(:,k) is index of vertical
1265! grid box which contains the departure point (restricted by N(ng)).
1266! During the search: also add in content of whole grid boxes
1267! participating in FC.
1268!
1269 cff=dtdays*abs(wbio(isink))
1270 DO k=1,n(ng)
1271 DO ks=k,n(ng)-1
1272 DO i=istr,iend
1273 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1274!^ tl_FC(i,k-1)=tl_FC(i,k-1)+tl_WR(i,ks)
1275!^
1276 ad_wr(i,ks)=ad_wr(i,ks)+ad_fc(i,k-1)
1277 END IF
1278 END DO
1279 END DO
1280 END DO
1281 DO k=1,n(ng)
1282 DO i=istr,iend
1283!^ tl_WR(i,k)=tl_Hz(i,j,k)*qc(i,k)+Hz(i,j,k)*tl_qc(i,k)
1284!^
1285 ad_hz(i,j,k)=ad_hz(i,j,k)+qc(i,k)*ad_wr(i,k)
1286 ad_qc(i,k)=ad_qc(i,k)+hz(i,j,k)*ad_wr(i,k)
1287 ad_wr(i,k)=0.0_r8
1288!^ tl_WL(i,k)=tl_z_w(i,j,k-1)+tl_cff
1289!^
1290 ad_z_w(i,j,k-1)=ad_z_w(i,j,k-1)+ad_wl(i,k)
1291 ad_cff=ad_cff+ad_wl(i,k)
1292 ad_wl(i,k)=0.0_r8
1293!^ tl_FC(i,k-1)=0.0_r8
1294!^
1295 ad_fc(i,k-1)=0.0_r8
1296 END DO
1297 END DO
1298!^ tl_cff=dtdays*SIGN(1.0_r8,Wbio(isink))*tl_Wbio(isink)
1299!^
1300 ad_wbio(isink)=ad_wbio(isink)+ &
1301 & dtdays*sign(1.0_r8,wbio(isink))*ad_cff
1302 ad_cff=0.0_r8
1303!
1304! Compute appropriate values of bR and bL.
1305!
1306! Copy concentration of biological particulates into scratch array
1307! "qc" (q-central, restrict it to be positive) which is hereafter
1308! interpreted as a set of grid-box averaged values for biogeochemical
1309! constituent concentration.
1310!
1311 DO k=1,n(ng)
1312 DO i=istr,iend
1313 qc(i,k)=bio(i,k,ibio)
1314 END DO
1315 END DO
1316!
1317 DO k=n(ng)-1,1,-1
1318 DO i=istr,iend
1319 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1320 END DO
1321 END DO
1322 DO k=2,n(ng)-1
1323 DO i=istr,iend
1324 dltr=hz(i,j,k)*fc(i,k)
1325 dltl=hz(i,j,k)*fc(i,k-1)
1326 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1327 cffr=cff*fc(i,k)
1328 cffl=cff*fc(i,k-1)
1329!
1330! Apply PPM monotonicity constraint to prevent oscillations within the
1331! grid box.
1332!
1333 IF ((dltr*dltl).le.0.0_r8) THEN
1334 dltr=0.0_r8
1335 dltl=0.0_r8
1336 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1337 dltr=cffl
1338 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1339 dltl=cffr
1340 END IF
1341!
1342! Compute right and left side values (bR,bL) of parabolic segments
1343! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1344!
1345! NOTE: Although each parabolic segment is monotonic within its grid
1346! box, monotonicity of the whole profile is not guaranteed,
1347! because bL(k+1)-bR(k) may still have different sign than
1348! qc(i,k+1)-qc(i,k). This possibility is excluded,
1349! after bL and bR are reconciled using WENO procedure.
1350!
1351 cff=(dltr-dltl)*hz_inv3(i,k)
1352 dltr=dltr-cff*hz(i,j,k+1)
1353 dltl=dltl+cff*hz(i,j,k-1)
1354 br(i,k)=qc(i,k)+dltr
1355 bl(i,k)=qc(i,k)-dltl
1356 wr(i,k)=(2.0_r8*dltr-dltl)**2
1357 wl(i,k)=(dltr-2.0_r8*dltl)**2
1358 END DO
1359 END DO
1360 cff=1.0e-14_r8
1361 DO k=2,n(ng)-2
1362 DO i=istr,iend
1363 dltl=max(cff,wl(i,k ))
1364 dltr=max(cff,wr(i,k+1))
1365 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1366 bl(i,k+1)=br(i,k)
1367 END DO
1368 END DO
1369 DO i=istr,iend
1370 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1371#if defined LINEAR_CONTINUATION
1372 bl(i,n(ng))=br(i,n(ng)-1)
1373 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1374#elif defined NEUMANN
1375 bl(i,n(ng))=br(i,n(ng)-1)
1376 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1377#else
1378 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1379 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1380 br(i,n(ng)-1)=qc(i,n(ng))
1381#endif
1382#if defined LINEAR_CONTINUATION
1383 br(i,1)=bl(i,2)
1384 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1385#elif defined NEUMANN
1386 br(i,1)=bl(i,2)
1387 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1388#else
1389 bl(i,2)=qc(i,1) ! bottom grid boxes are
1390 br(i,1)=qc(i,1) ! re-assumed to be
1391 bl(i,1)=qc(i,1) ! piecewise constant.
1392#endif
1393 END DO
1394!
1395! Apply monotonicity constraint again, since the reconciled interfacial
1396! values may cause a non-monotonic behavior of the parabolic segments
1397! inside the grid box.
1398!
1399 DO k=1,n(ng)
1400 DO i=istr,iend
1401 dltr=br(i,k)-qc(i,k)
1402 dltl=qc(i,k)-bl(i,k)
1403 cffr=2.0_r8*dltr
1404 cffl=2.0_r8*dltl
1405!^ tl_bL(i,k)=tl_qc(i,k)-tl_dltL
1406!^
1407 ad_qc(i,k)=ad_qc(i,k)+ad_bl(i,k)
1408 ad_dltl=ad_dltl-ad_bl(i,k)
1409 ad_bl(i,k)=0.0_r8
1410!^ tl_bR(i,k)=tl_qc(i,k)+tl_dltR
1411!^
1412 ad_qc(i,k)=ad_qc(i,k)+ad_br(i,k)
1413 ad_dltr=ad_dltr+ad_br(i,k)
1414 ad_br(i,k)=0.0_r8
1415 IF ((dltr*dltl).lt.0.0_r8) THEN
1416!^ tl_dltR=0.0_r8
1417!^
1418 ad_dltr=0.0_r8
1419!^ tl_dltL=0.0_r8
1420!^
1421 ad_dltl=0.0_r8
1422 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1423!^ tl_dltR=tl_cffL
1424!^
1425 ad_cffl=ad_cffl+ad_dltr
1426 ad_dltr=0.0_r8
1427 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1428!^ tl_dltL=tl_cffR
1429!^
1430 ad_cffr=ad_cffr+ad_dltl
1431 ad_dltl=0.0_r8
1432 END IF
1433!^ tl_cffL=2.0_r8*tl_dltL
1434!^
1435 ad_dltl=ad_dltl+2.0_r8*ad_cffl
1436 ad_cffl=0.0_r8
1437!^ tl_cffR=2.0_r8*tl_dltR
1438!^
1439 ad_dltr=ad_dltr+2.0_r8*ad_cffr
1440 ad_cffr=0.0_r8
1441!^ tl_dltL=tl_qc(i,k)-tl_bL(i,k)
1442!^
1443 ad_qc(i,k)=ad_qc(i,k)+ad_dltl
1444 ad_bl(i,k)=ad_bl(i,k)-ad_dltl
1445 ad_dltl=0.0_r8
1446!^ tl_dltR=tl_bR(i,k)-tl_qc(i,k)
1447!^
1448 ad_br(i,k)=ad_br(i,k)+ad_dltr
1449 ad_qc(i,k)=ad_qc(i,k)-ad_dltr
1450 ad_dltr=0.0_r8
1451 END DO
1452 END DO
1453 DO i=istr,iend
1454#if defined LINEAR_CONTINUATION
1455!^ tl_bR(i,1)=tl_bL(i,2)
1456!^
1457 ad_bl(i,2)=ad_bl(i,2)+ad_br(i,1)
1458 ad_br(i,1)=0.0_r8
1459!^ tl_bL(i,1)=2.0_r8*tl_qc(i,1)-tl_bR(i,1)
1460!^
1461 ad_qc(i,1)=ad_qc(i,1)+2.0_r8*ad_bl(i,1)
1462 ad_br(i,1)=ad_br(i,1)-ad_bl(i,1)
1463 ad_bl(i,1)=0.0_r8
1464#elif defined NEUMANN
1465!^ tl_bR(i,1)=tl_bL(i,2)
1466!^
1467 ad_bl(i,2)=ad_bl(i,2)+ad_br(i,1)
1468 ad_br(i,1)=0.0_r8
1469!^ tl_bL(i,1)=1.5_r8*tl_qc(i,1)-0.5_r8*tl_bR(i,1)
1470!^
1471 ad_qc(i,1)=ad_qc(i,1)+1.5_r8*ad_bl(i,1)
1472 ad_br(i,1)=ad_br(i,1)-0.5_r8*ad_bl(i,1)
1473 ad_bl(i,1)=0.0_r8
1474#else
1475!^ tl_bL(i,2)=tl_qc(i,1) ! bottom grid boxes are
1476!^ tl_bR(i,1)=tl_qc(i,1) ! re-assumed to be
1477!^ tl_bL(i,1)=tl_qc(i,1) ! piecewise constant.
1478!^
1479 ad_qc(i,1)=ad_qc(i,1)+ad_bl(i,1)+ &
1480 & ad_br(i,1)+ &
1481 & ad_bl(i,2)
1482 ad_bl(i,1)=0.0_r8
1483 ad_br(i,1)=0.0_r8
1484 ad_bl(i,2)=0.0_r8
1485#endif
1486#if defined LINEAR_CONTINUATION
1487!^ tl_bL(i,N(ng))=tl_bR(i,N(ng)-1)
1488!^
1489 ad_br(i,n(ng)-1)=ad_br(i,n(ng)-1)+ad_bl(i,n(ng))
1490 ad_bl(i,n(ng))=0.0_r8
1491!^ tl_bR(i,N(ng))=2.0_r8*tl_qc(i,N(ng))-tl_bL(i,N(ng))
1492!^
1493 ad_qc(i,n(ng))=ad_qc(i,n(ng))+2.0_r8*ad_br(i,n(ng))
1494 ad_bl(i,n(ng))=ad_bl(i,n(ng))-ad_br(i,n(ng))
1495 ad_br(i,n(ng))=0.0_r8
1496#elif defined NEUMANN
1497!^ tl_bL(i,N(ng))=tl_bR(i,N(ng)-1)
1498!^
1499 ad_br(i,n(ng)-1)=ad_br(i,n(ng)-1)+ad_bl(i,n(ng))
1500 ad_bl(i,n(ng))=0.0_r8
1501!^ tl_bR(i,N(ng))=1.5_r8*tl_qc(i,N(ng))-0.5_r8*tl_bL(i,N(ng))
1502!^
1503 ad_qc(i,n(ng))=ad_qc(i,n(ng))+1.5_r8*ad_br(i,n(ng))
1504 ad_bl(i,n(ng))=ad_bl(i,n(ng))-0.5_r8*ad_br(i,n(ng))
1505 ad_br(i,n(ng))=0.0_r8
1506#else
1507!^ tl_bR(i,N(ng))=tl_qc(i,N(ng)) ! default strictly monotonic
1508!^ tl_bL(i,N(ng))=tl_qc(i,N(ng)) ! conditions
1509!^ tl_bR(i,N(ng)-1)=tl_qc(i,N(ng))
1510!^
1511 ad_qc(i,n(ng))=ad_qc(i,n(ng))+ad_br(i,n(ng)-1)+ &
1512 & ad_bl(i,n(ng))+ &
1513 & ad_br(i,n(ng))
1514 ad_br(i,n(ng)-1)=0.0_r8
1515 ad_bl(i,n(ng))=0.0_r8
1516 ad_br(i,n(ng))=0.0_r8
1517#endif
1518 END DO
1519
1520!
1521! Compute WR and WL arrays appropriate for this part of the code.
1522!
1523! Copy concentration of biological particulates into scratch array
1524! "qc" (q-central, restrict it to be positive) which is hereafter
1525! interpreted as a set of grid-box averaged values for biogeochemical
1526! constituent concentration.
1527!
1528 DO k=1,n(ng)
1529 DO i=istr,iend
1530 qc(i,k)=bio(i,k,ibio)
1531 END DO
1532 END DO
1533!
1534 DO k=n(ng)-1,1,-1
1535 DO i=istr,iend
1536 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1537 END DO
1538 END DO
1539 DO k=2,n(ng)-1
1540 DO i=istr,iend
1541 dltr=hz(i,j,k)*fc(i,k)
1542 dltl=hz(i,j,k)*fc(i,k-1)
1543 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1544 cffr=cff*fc(i,k)
1545 cffl=cff*fc(i,k-1)
1546!
1547! Apply PPM monotonicity constraint to prevent oscillations within the
1548! grid box.
1549!
1550 IF ((dltr*dltl).le.0.0_r8) THEN
1551 dltr=0.0_r8
1552 dltl=0.0_r8
1553 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1554 dltr=cffl
1555 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1556 dltl=cffr
1557 END IF
1558!
1559! Compute right and left side values (bR,bL) of parabolic segments
1560! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1561!
1562! NOTE: Although each parabolic segment is monotonic within its grid
1563! box, monotonicity of the whole profile is not guaranteed,
1564! because bL(k+1)-bR(k) may still have different sign than
1565! qc(i,k+1)-qc(i,k). This possibility is excluded,
1566! after bL and bR are reconciled using WENO procedure.
1567!
1568 cff=(dltr-dltl)*hz_inv3(i,k)
1569 dltr=dltr-cff*hz(i,j,k+1)
1570 dltl=dltl+cff*hz(i,j,k-1)
1571 br(i,k)=qc(i,k)+dltr
1572 bl(i,k)=qc(i,k)-dltl
1573 wr(i,k)=(2.0_r8*dltr-dltl)**2
1574 wl(i,k)=(dltr-2.0_r8*dltl)**2
1575 END DO
1576 END DO
1577
1578 cff=1.0e-14_r8
1579 DO k=2,n(ng)-2
1580 DO i=istr,iend
1581 dltl=max(cff,wl(i,k ))
1582 dltr=max(cff,wr(i,k+1))
1583 br1(i,k)=br(i,k)
1584 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1585 bl1(i,k+1)=bl(i,k+1)
1586 bl(i,k+1)=br(i,k)
1587!^ tl_bL(i,k+1)=tl_bR(i,k)
1588!^
1589 ad_br(i,k)=ad_br(i,k)+ad_bl(i,k+1)
1590 ad_bl(i,k+1)=0.0_r8
1591!^ tl_bR(i,k)=(tl_dltR*bR1(i,k )+dltR*tl_bR(i,k )+ &
1592!^ & tl_dltL*bL1(i,k+1)+dltL*tl_bL(i,k+1))/ &
1593!^ & (dltR+dltL)- &
1594!^ & (tl_dltR+tl_dltL)*bR(i,k)/(dltR+dltL)
1595!^
1596 adfac=ad_br(i,k)/(dltr+dltl)
1597 adfac1=ad_br(i,k)*br(i,k)/(dltr+dltl)
1598 ad_dltr=ad_dltr+adfac*br1(i,k)
1599 ad_dltl=ad_dltl+adfac*bl1(i,k+1)
1600 ad_bl(i,k+1)=ad_bl(i,k+1)+dltl*adfac
1601 ad_dltr=ad_dltr-adfac1
1602 ad_dltl=ad_dltl-adfac1
1603 ad_br(i,k)=dltr*adfac
1604!^ tl_dltR=(0.5_r8-SIGN(0.5_r8,cff-WR(i,k+1)))* &
1605!^ & tl_WR(i,k+1)
1606!^
1607 ad_wr(i,k+1)=ad_wr(i,k+1)+ &
1608 & (0.5_r8-sign(0.5_r8,cff-wr(i,k+1)))* &
1609 & ad_dltr
1610 ad_dltr=0.0_r8
1611!^ tl_dltL=(0.5_r8-SIGN(0.5_r8,cff-WL(i,k )))* &
1612!^ & tl_WL(i,k )
1613 ad_wl(i,k )=ad_wl(i,k )+ &
1614 & (0.5_r8-sign(0.5_r8,cff-wl(i,k )))* &
1615 & ad_dltl
1616 ad_dltl=0.0_r8
1617 END DO
1618 END DO
1619
1620 DO k=2,n(ng)-1
1621 DO i=istr,iend
1622!
1623! Compute appropriate dltL and dltr.
1624!
1625 dltr=hz(i,j,k)*fc(i,k)
1626 dltl=hz(i,j,k)*fc(i,k-1)
1627 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1628 cffr=cff*fc(i,k)
1629 cffl=cff*fc(i,k-1)
1630!
1631! Apply PPM monotonicity constraint to prevent oscillations within the
1632! grid box.
1633!
1634 IF ((dltr*dltl).le.0.0_r8) THEN
1635 dltr=0.0_r8
1636 dltl=0.0_r8
1637 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1638 dltr=cffl
1639 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1640 dltl=cffr
1641 END IF
1642!
1643! Compute right and left side values (bR,bL) of parabolic segments
1644! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1645!
1646! NOTE: Although each parabolic segment is monotonic within its grid
1647! box, monotonicity of the whole profile is not guaranteed,
1648! because bL(k+1)-bR(k) may still have different sign than
1649! qc(i,k+1)-qc(i,k). This possibility is excluded,
1650! after bL and bR are reconciled using WENO procedure.
1651!
1652 cff=(dltr-dltl)*hz_inv3(i,k)
1653 dltr=dltr-cff*hz(i,j,k+1)
1654 dltl=dltl+cff*hz(i,j,k-1)
1655!^ tl_WL(i,k)=2.0_r8*(dltR-2.0_r8*dltL)* &
1656!^ & (tl_dltR-2.0_r8*tl_dltL)
1657!^
1658 adfac=ad_wl(i,k)*2.0_r8*(dltr-2.0_r8*dltl)
1659 ad_dltr=ad_dltr+adfac
1660 ad_dltl=ad_dltl-2.0_r8*adfac
1661 ad_wl(i,k)=0.0_r8
1662!^ tl_WR(i,k)=2.0_r8*(2.0_r8*dltR-dltL)* &
1663!^ & (2.0_r8*tl_dltR-tl_dltL)
1664!^
1665 adfac=ad_wr(i,k)*2.0_r8*(2.0_r8*dltr-dltl)
1666 ad_dltr=ad_dltr+2.0_r8*adfac
1667 ad_dltl=ad_dltl-adfac
1668 ad_wr(i,k)=0.0_r8
1669!^ tl_bL(i,k)=tl_qc(i,k)-tl_dltL
1670!^
1671 ad_qc(i,k)=ad_qc(i,k)+ad_bl(i,k)
1672 ad_dltl=ad_dltl-ad_bl(i,k)
1673 ad_bl(i,k)=0.0_r8
1674!^ tl_bR(i,k)=tl_qc(i,k)+tl_dltR
1675!^
1676 ad_qc(i,k)=ad_qc(i,k)+ad_br(i,k)
1677 ad_dltr=ad_dltr+ad_br(i,k)
1678 ad_br(i,k)=0.0_r8
1679
1680!
1681! Compute appropriate dltL and dltr.
1682!
1683 dltr=hz(i,j,k)*fc(i,k)
1684 dltl=hz(i,j,k)*fc(i,k-1)
1685 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1686 cffr=cff*fc(i,k)
1687 cffl=cff*fc(i,k-1)
1688!
1689! Apply PPM monotonicity constraint to prevent oscillations within the
1690! grid box.
1691!
1692 IF ((dltr*dltl).le.0.0_r8) THEN
1693 dltr=0.0_r8
1694 dltl=0.0_r8
1695 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1696 dltr=cffl
1697 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1698 dltl=cffr
1699 END IF
1700
1701 cff=(dltr-dltl)*hz_inv3(i,k)
1702!^ tl_dltL=tl_dltL+tl_cff*Hz(i,j,k-1)+cff*tl_Hz(i,j,k-1)
1703!^
1704 ad_cff=ad_cff+ad_dltl*hz(i,j,k-1)
1705 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)+cff*ad_dltl
1706!^ tl_dltR=tl_dltR-tl_cff*Hz(i,j,k+1)-cff*tl_Hz(i,j,k+1)
1707!^
1708 ad_cff=ad_cff-ad_dltr*hz(i,j,k+1)
1709 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)-cff*ad_dltr
1710!^ tl_cff=(tl_dltR-tl_dltL)*Hz_inv3(i,k)+ &
1711!^ & (dltR-dltL)*tl_Hz_inv3(i,k)
1712!^
1713 adfac=ad_cff*hz_inv3(i,k)
1714 ad_dltr=ad_dltr+adfac
1715 ad_dltl=ad_dltl-adfac
1716 ad_hz_inv3(i,k)=ad_hz_inv3(i,k)+(dltr-dltl)*ad_cff
1717 ad_cff=0.0_r8
1718!
1719! Compute appropriate dltL and dltr.
1720!
1721 dltr=hz(i,j,k)*fc(i,k)
1722 dltl=hz(i,j,k)*fc(i,k-1)
1723 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1724 cffr=cff*fc(i,k)
1725 cffl=cff*fc(i,k-1)
1726
1727 IF ((dltr*dltl).le.0.0_r8) THEN
1728!^ tl_dltR=0.0_r8
1729!^
1730 ad_dltr=0.0_r8
1731!^ tl_dltL=0.0_r8
1732!^
1733 ad_dltl=0.0_r8
1734 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1735!^ tl_dltR=tl_cffL
1736!^
1737 ad_cffl=ad_cffl+ad_dltr
1738 ad_dltr=0.0_r8
1739 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1740!^ tl_dltL=tl_cffR
1741!^
1742 ad_cffr=ad_cffr+ad_dltl
1743 ad_dltl=0.0_r8
1744 END IF
1745!^ tl_cffL=tl_cff*FC(i,k-1)+cff*tl_FC(i,k-1)
1746!^
1747 ad_cff=ad_cff+ad_cffl*fc(i,k-1)
1748 ad_fc(i,k-1)=ad_fc(i,k-1)+cff*ad_cffl
1749 ad_cffl=0.0_r8
1750!^ tl_cffR=tl_cff*FC(i,k)+cff*tl_FC(i,k)
1751!^
1752 ad_cff=ad_cff+ad_cffr*fc(i,k)
1753 ad_fc(i,k)=ad_fc(i,k)+cff*ad_cffr
1754 ad_cffr=0.0_r8
1755!^ tl_cff=tl_Hz(i,j,k-1)+2.0_r8*tl_Hz(i,j,k)+tl_Hz(i,j,k+1)
1756!^
1757 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)+ad_cff
1758 ad_hz(i,j,k)=ad_hz(i,j,k)+2.0_r8*ad_cff
1759 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)+ad_cff
1760 ad_cff=0.0_r8
1761!^ tl_dltL=tl_Hz(i,j,k)*FC(i,k-1)+Hz(i,j,k)*tl_FC(i,k-1)
1762!^
1763 ad_hz(i,j,k)=ad_hz(i,j,k)+ad_dltl*fc(i,k-1)
1764 ad_fc(i,k-1)=ad_fc(i,k-1)+ad_dltl*hz(i,j,k)
1765 ad_dltl=0.0_r8
1766!^ tl_dltR=tl_Hz(i,j,k)*FC(i,k)+Hz(i,j,k)*tl_FC(i,k)
1767!^
1768 ad_hz(i,j,k)=ad_hz(i,j,k)+ad_dltr*fc(i,k)
1769 ad_fc(i,k)=ad_fc(i,k)+ad_dltr*hz(i,j,k)
1770 ad_dltr=0.0_r8
1771 END DO
1772 END DO
1773 DO k=n(ng)-1,1,-1
1774 DO i=istr,iend
1775!^ tl_FC(i,k)=(tl_qc(i,k+1)-tl_qc(i,k))*Hz_inv2(i,k)+ &
1776!^ & (qc(i,k+1)-qc(i,k))*tl_Hz_inv2(i,k)
1777!^
1778 adfac=ad_fc(i,k)*hz_inv2(i,k)
1779 ad_qc(i,k+1)=ad_qc(i,k+1)+adfac
1780 ad_qc(i,k)=ad_qc(i,k)-adfac
1781 ad_hz_inv2(i,k)=ad_hz_inv2(i,k)+(qc(i,k+1)-qc(i,k))* &
1782 & ad_fc(i,k)
1783 ad_fc(i,k)=0.0_r8
1784 END DO
1785 END DO
1786 DO k=1,n(ng)
1787 DO i=istr,iend
1788!^ tl_qc(i,k)=tl_Bio(i,k,ibio)
1789!^
1790 ad_bio(i,k,ibio)=ad_bio(i,k,ibio)+ad_qc(i,k)
1791 ad_qc(i,k)=0.0_r8
1792 END DO
1793 END DO
1794
1795 END DO sink_loop
1796!
1797! Compute appropriate basic state arrays II.
1798!
1799! Determine Correction for negativity.
1800!
1801 DO k=1,n(ng)
1802 DO i=istr,iend
1803 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
1804 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
1805 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
1806 & max(0.0_r8,eps-bio_old(i,k,isdet))
1807!
1808! If correction needed, determine the largest pool to debit.
1809!
1810 IF (cff1.gt.0.0) THEN
1811 itrmx=idbio(1)
1812 cff=t(i,j,k,nstp,itrmx)
1813 DO ibio=idbio(2),idbio(nbt)
1814 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
1815 itrmx=ibio
1816 cff=t(i,j,k,nstp,ibio)
1817 END IF
1818 END DO
1819!
1820! Update new values.
1821!
1822 DO itrc=1,nbt
1823 ibio=idbio(itrc)
1824 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
1825 & cff1*(sign(0.5_r8, &
1826 & real(itrmx-ibio,r8)**2)+ &
1827 & sign(0.5_r8, &
1828 & -real(itrmx-ibio,r8)**2))
1829 END DO
1830 ELSE
1831 DO itrc=1,nbt
1832 ibio=idbio(itrc)
1833 bio(i,k,ibio)=bio_old(i,k,ibio)
1834 END DO
1835 END IF
1836 END DO
1837 END DO
1838!
1839!=======================================================================
1840! Start internal iterations to achieve convergence of the nonlinear
1841! backward-implicit solution.
1842!=======================================================================
1843!
1844 DO iteradj=1,iter
1845!
1846! Nutrient uptake by phytoplankton.
1847!
1848 cff1=dtdays*vm_no3(ng)
1849 DO k=1,n(ng)
1850 DO i=istr,iend
1851 cff=bio(i,k,iphyt)* &
1852 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
1853 & (k_no3(ng)+bio(i,k,ino3_))
1854 bio1(i,k,ino3_)=bio(i,k,ino3_)
1855 bio(i,k,ino3_)=bio(i,k,ino3_)/ &
1856 & (1.0_r8+cff)
1857 bio1(i,k,iphyt)=bio(i,k,iphyt)
1858 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
1859 & bio(i,k,ino3_)*cff
1860 END DO
1861 END DO
1862!
1863! Phytoplankton grazing by Zooplankton and mortality to Detritus
1864! (rate: PhyMR).
1865!
1866 cff1=dtdays*zoogr(ng)
1867 cff2=dtdays*phymr(ng)
1868 cff3=k_phy(ng)*k_phy(ng)
1869 DO k=1,n(ng)
1870 DO i=istr,iend
1871 cff=bio(i,k,izoop)*bio(i,k,iphyt)*cff1/ &
1872 & (cff3+bio(i,k,iphyt)*bio(i,k,iphyt))
1873 bio1(i,k,iphyt)=bio(i,k,iphyt)
1874 bio(i,k,iphyt)=bio(i,k,iphyt)/ &
1875 & (1.0_r8+cff+cff2)
1876 bio1(i,k,izoop)=bio(i,k,izoop)
1877 bio(i,k,izoop)=bio(i,k,izoop)+ &
1878 & bio(i,k,iphyt)*cff*(1.0_r8-zooga(ng))
1879 bio(i,k,isdet)=bio(i,k,isdet)+ &
1880 & bio(i,k,iphyt)* &
1881 & (cff2+cff*(zooga(ng)-zooec(ng)))
1882 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1883 & bio(i,k,iphyt)*cff*zooec(ng)
1884 END DO
1885 END DO
1886!
1887 IF (iteradj.ne.iter) THEN
1888!
1889! Zooplankton excretion to nutrients and mortality to Detritus.
1890!
1891 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
1892 cff2=dtdays*zoomr(ng)
1893 cff3=dtdays*zoomd(ng)
1894 DO k=1,n(ng)
1895 DO i=istr,iend
1896 bio(i,k,izoop)=bio(i,k,izoop)*cff1
1897 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1898 & bio(i,k,izoop)*cff2
1899 bio(i,k,isdet)=bio(i,k,isdet)+ &
1900 & bio(i,k,izoop)*cff3
1901 END DO
1902 END DO
1903!
1904! Detritus breakdown to nutrients.
1905!
1906 cff1=dtdays*detrr(ng)
1907 cff2=1.0_r8/(1.0_r8+cff1)
1908 DO k=1,n(ng)
1909 DO i=istr,iend
1910 bio(i,k,isdet)=bio(i,k,isdet)*cff2
1911 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
1912 & bio(i,k,isdet)*cff1
1913 END DO
1914 END DO
1915!
1916!-----------------------------------------------------------------------
1917! Vertical sinking terms.
1918!-----------------------------------------------------------------------
1919!
1920! Reconstruct vertical profile of selected biological constituents
1921! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
1922! grid box. Then, compute semi-Lagrangian flux due to sinking.
1923!
1924 DO isink=1,nsink
1925 ibio=idsink(isink)
1926!
1927! Copy concentration of biological particulates into scratch array
1928! "qc" (q-central, restrict it to be positive) which is hereafter
1929! interpreted as a set of grid-box averaged values for biogeochemical
1930! constituent concentration.
1931!
1932 DO k=1,n(ng)
1933 DO i=istr,iend
1934 qc(i,k)=bio(i,k,ibio)
1935 END DO
1936 END DO
1937!
1938 DO k=n(ng)-1,1,-1
1939 DO i=istr,iend
1940 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1941 END DO
1942 END DO
1943 DO k=2,n(ng)-1
1944 DO i=istr,iend
1945 dltr=hz(i,j,k)*fc(i,k)
1946 dltl=hz(i,j,k)*fc(i,k-1)
1947 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1948 cffr=cff*fc(i,k)
1949 cffl=cff*fc(i,k-1)
1950!
1951! Apply PPM monotonicity constraint to prevent oscillations within the
1952! grid box.
1953!
1954 IF ((dltr*dltl).le.0.0_r8) THEN
1955 dltr=0.0_r8
1956 dltl=0.0_r8
1957 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1958 dltr=cffl
1959 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1960 dltl=cffr
1961 END IF
1962!
1963! Compute right and left side values (bR,bL) of parabolic segments
1964! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1965!
1966! NOTE: Although each parabolic segment is monotonic within its grid
1967! box, monotonicity of the whole profile is not guaranteed,
1968! because bL(k+1)-bR(k) may still have different sign than
1969! qc(i,k+1)-qc(i,k). This possibility is excluded,
1970! after bL and bR are reconciled using WENO procedure.
1971!
1972 cff=(dltr-dltl)*hz_inv3(i,k)
1973 dltr=dltr-cff*hz(i,j,k+1)
1974 dltl=dltl+cff*hz(i,j,k-1)
1975 br(i,k)=qc(i,k)+dltr
1976 bl(i,k)=qc(i,k)-dltl
1977 wr(i,k)=(2.0_r8*dltr-dltl)**2
1978 wl(i,k)=(dltr-2.0_r8*dltl)**2
1979 END DO
1980 END DO
1981 cff=1.0e-14_r8
1982 DO k=2,n(ng)-2
1983 DO i=istr,iend
1984 dltl=max(cff,wl(i,k ))
1985 dltr=max(cff,wr(i,k+1))
1986 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1987 bl(i,k+1)=br(i,k)
1988 END DO
1989 END DO
1990 DO i=istr,iend
1991 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1992#if defined LINEAR_CONTINUATION
1993 bl(i,n(ng))=br(i,n(ng)-1)
1994 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1995#elif defined NEUMANN
1996 bl(i,n(ng))=br(i,n(ng)-1)
1997 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1998#else
1999 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
2000 bl(i,n(ng))=qc(i,n(ng)) ! conditions
2001 br(i,n(ng)-1)=qc(i,n(ng))
2002#endif
2003#if defined LINEAR_CONTINUATION
2004 br(i,1)=bl(i,2)
2005 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2006#elif defined NEUMANN
2007 br(i,1)=bl(i,2)
2008 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2009#else
2010 bl(i,2)=qc(i,1) ! bottom grid boxes are
2011 br(i,1)=qc(i,1) ! re-assumed to be
2012 bl(i,1)=qc(i,1) ! piecewise constant.
2013#endif
2014 END DO
2015!
2016! Apply monotonicity constraint again, since the reconciled interfacial
2017! values may cause a non-monotonic behavior of the parabolic segments
2018! inside the grid box.
2019!
2020 DO k=1,n(ng)
2021 DO i=istr,iend
2022 dltr=br(i,k)-qc(i,k)
2023 dltl=qc(i,k)-bl(i,k)
2024 cffr=2.0_r8*dltr
2025 cffl=2.0_r8*dltl
2026 IF ((dltr*dltl).lt.0.0_r8) THEN
2027 dltr=0.0_r8
2028 dltl=0.0_r8
2029 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2030 dltr=cffl
2031 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2032 dltl=cffr
2033 END IF
2034 br(i,k)=qc(i,k)+dltr
2035 bl(i,k)=qc(i,k)-dltl
2036 END DO
2037 END DO
2038!
2039! After this moment reconstruction is considered complete. The next
2040! stage is to compute vertical advective fluxes, FC. It is expected
2041! that sinking may occurs relatively fast, the algorithm is designed
2042! to be free of CFL criterion, which is achieved by allowing
2043! integration bounds for semi-Lagrangian advective flux to use as
2044! many grid boxes in upstream direction as necessary.
2045!
2046! In the two code segments below, WL is the z-coordinate of the
2047! departure point for grid box interface z_w with the same indices;
2048! FC is the finite volume flux; ksource(:,k) is index of vertical
2049! grid box which contains the departure point (restricted by N(ng)).
2050! During the search: also add in content of whole grid boxes
2051! participating in FC.
2052!
2053 cff=dtdays*abs(wbio(isink))
2054 DO k=1,n(ng)
2055 DO i=istr,iend
2056 fc(i,k-1)=0.0_r8
2057 wl(i,k)=z_w(i,j,k-1)+cff
2058 wr(i,k)=hz(i,j,k)*qc(i,k)
2059 ksource(i,k)=k
2060 END DO
2061 END DO
2062 DO k=1,n(ng)
2063 DO ks=k,n(ng)-1
2064 DO i=istr,iend
2065 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2066 ksource(i,k)=ks+1
2067 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2068 END IF
2069 END DO
2070 END DO
2071 END DO
2072!
2073! Finalize computation of flux: add fractional part.
2074!
2075 DO k=1,n(ng)
2076 DO i=istr,iend
2077 ks=ksource(i,k)
2078 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2079 fc(i,k-1)=fc(i,k-1)+ &
2080 & hz(i,j,ks)*cu* &
2081 & (bl(i,ks)+ &
2082 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2083 & (1.5_r8-cu)* &
2084 & (br(i,ks)+bl(i,ks)- &
2085 & 2.0_r8*qc(i,ks))))
2086 END DO
2087 END DO
2088 DO k=1,n(ng)
2089 DO i=istr,iend
2090 bio(i,k,ibio)=qc(i,k)+ &
2091 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2092 END DO
2093 END DO
2094 END DO
2095 END IF
2096 END DO
2097!
2098! End compute basic state arrays II.
2099!
2100! Adjoint Detritus breakdown to nutrients.
2101!
2102 cff1=dtdays*detrr(ng)
2103 cff2=1.0_r8/(1.0_r8+cff1)
2104 DO k=1,n(ng)
2105 DO i=istr,iend
2106!^ tl_Bio(i,k,iNO3_)=tl_Bio(i,k,iNO3_)+ &
2107!^ & tl_Bio(i,k,iSDet)*cff1
2108!^
2109 ad_bio(i,k,isdet)=ad_bio(i,k,isdet)+ &
2110 & cff1*ad_bio(i,k,ino3_)
2111!^ tl_Bio(i,k,iSDet)=tl_Bio(i,k,iSDet)*cff2
2112!^
2113 ad_bio(i,k,isdet)=cff2*ad_bio(i,k,isdet)
2114 END DO
2115 END DO
2116!
2117! Adjoint Zooplankton excretion to nutrients and mortality to Detritus.
2118!
2119 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
2120 cff2=dtdays*zoomr(ng)
2121 cff3=dtdays*zoomd(ng)
2122 DO k=1,n(ng)
2123 DO i=istr,iend
2124!^ tl_Bio(i,k,iSDet)=tl_Bio(i,k,iSDet)+ &
2125!^ & tl_Bio(i,k,iZoop)*cff3
2126!^
2127 ad_bio(i,k,izoop)=ad_bio(i,k,izoop)+ &
2128 & cff3*ad_bio(i,k,isdet)
2129!^ tl_Bio(i,k,iNO3_)=tl_Bio(i,k,iNO3_)+ &
2130!^ & tl_Bio(i,k,iZoop)*cff2
2131!^
2132 ad_bio(i,k,izoop)=ad_bio(i,k,izoop)+ &
2133 & cff2*ad_bio(i,k,ino3_)
2134!^ tl_Bio(i,k,iZoop)=tl_Bio(i,k,iZoop)*cff1
2135!^
2136 ad_bio(i,k,izoop)=cff1*ad_bio(i,k,izoop)
2137 END DO
2138 END DO
2139!
2140! Phytoplankton grazing by Zooplankton and mortality to Detritus
2141! (rate: PhyMR).
2142!
2143 cff1=dtdays*zoogr(ng)
2144 cff2=dtdays*phymr(ng)
2145 cff3=k_phy(ng)*k_phy(ng)
2146 DO k=1,n(ng)
2147 DO i=istr,iend
2148 cff=bio1(i,k,izoop)*bio1(i,k,iphyt)*cff1/ &
2149 & (cff3+bio1(i,k,iphyt)*bio1(i,k,iphyt))
2150!^ tl_Bio(i,k,iNO3_)=tl_Bio(i,k,iNO3_)+ &
2151!^ & tl_Bio(i,k,iPhyt)*cff*ZooEC(ng)+ &
2152!^ & Bio(i,k,iPhyt)*tl_cff*ZooEC(ng)
2153!^
2154 ad_cff=ad_cff+ &
2155 & bio(i,k,iphyt)*zooec(ng)*ad_bio(i,k,ino3_)
2156 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2157 & cff*zooec(ng)*ad_bio(i,k,ino3_)
2158!^ tl_Bio(i,k,iSDet)=tl_Bio(i,k,iSDet)+ &
2159!^ & tl_Bio(i,k,iPhyt)* &
2160!^ & (cff2+cff*(ZooGA(ng)-ZooEC(ng)))+ &
2161!^ & Bio(i,k,iPhyt)* &
2162!^ & tl_cff*(ZooGA(ng)-ZooEC(ng))
2163!^
2164 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2165 & (cff2+cff*(zooga(ng)-zooec(ng)))* &
2166 & ad_bio(i,k,isdet)
2167 ad_cff=ad_cff+ &
2168 & bio(i,k,iphyt)*(zooga(ng)-zooec(ng))* &
2169 & ad_bio(i,k,isdet)
2170!^ tl_Bio(i,k,iZoop)=tl_Bio(i,k,iZoop)+ &
2171!^ & tl_Bio(i,k,iPhyt)* &
2172!^ & cff*(1.0_r8-ZooGA(ng))+ &
2173!^ & Bio(i,k,iPhyt)* &
2174!^ & tl_cff*(1.0_r8-ZooGA(ng))
2175!^
2176 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2177 & cff*(1.0_r8-zooga(ng))* &
2178 & ad_bio(i,k,izoop)
2179 ad_cff=ad_cff+ &
2180 & bio(i,k,iphyt)*(1.0_r8-zooga(ng))* &
2181 & ad_bio(i,k,izoop)
2182!^ tl_Bio(i,k,iPhyt)=(tl_Bio(i,k,iPhyt)- &
2183!^ & tl_cff*Bio(i,k,iPhyt))/ &
2184!^ & (1.0_r8+cff+cff2)
2185!^
2186 adfac=ad_bio(i,k,iphyt)/(1.0_r8+cff+cff2)
2187 ad_cff=ad_cff-bio(i,k,iphyt)*adfac
2188 ad_bio(i,k,iphyt)=adfac
2189!^ tl_cff=((tl_Bio(i,k,iZoop)*Bio1(i,k,iPhyt)+ &
2190!^ & Bio1(i,k,iZoop)*tl_Bio(i,k,iPhyt))*cff1- &
2191!^ & 2.0_r8*Bio1(i,k,iPhyt)*tl_Bio(i,k,iPhyt)*cff)/ &
2192!^ & (cff3+Bio1(i,k,iPhyt)*Bio1(i,k,iPhyt))
2193!^
2194 adfac=ad_cff/(cff3+bio1(i,k,iphyt)*bio1(i,k,iphyt))
2195 adfac1=adfac*cff1
2196 ad_bio(i,k,izoop)=ad_bio(i,k,izoop)+ &
2197 & bio1(i,k,iphyt)*adfac1
2198 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2199 & bio1(i,k,izoop)*adfac1
2200 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)- &
2201 & 2.0_r8*bio1(i,k,iphyt)*cff*adfac
2202 ad_cff=0.0_r8
2203 END DO
2204 END DO
2205!
2206! Compute appropriate basic state arrays I.
2207!
2208! Determine Correction for negativity.
2209!
2210 DO k=1,n(ng)
2211 DO i=istr,iend
2212 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
2213 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
2214 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
2215 & max(0.0_r8,eps-bio_old(i,k,isdet))
2216!
2217! If correction needed, determine the largest pool to debit.
2218!
2219 IF (cff1.gt.0.0) THEN
2220 itrmx=idbio(1)
2221 cff=t(i,j,k,nstp,itrmx)
2222 DO ibio=idbio(2),idbio(nbt)
2223 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
2224 itrmx=ibio
2225 cff=t(i,j,k,nstp,ibio)
2226 END IF
2227 END DO
2228!
2229! Update new values.
2230!
2231 DO itrc=1,nbt
2232 ibio=idbio(itrc)
2233 bio(i,k,ibio)=max(eps,bio_old(i,k,ibio))- &
2234 & cff1*(sign(0.5_r8, &
2235 & real(itrmx-ibio,r8)**2)+ &
2236 & sign(0.5_r8, &
2237 & -real(itrmx-ibio,r8)**2))
2238 END DO
2239 ELSE
2240 DO itrc=1,nbt
2241 ibio=idbio(itrc)
2242 bio(i,k,ibio)=bio_old(i,k,ibio)
2243 END DO
2244 END IF
2245 END DO
2246 END DO
2247!
2248!=======================================================================
2249! Start internal iterations to achieve convergence of the nonlinear
2250! backward-implicit solution.
2251!=======================================================================
2252!
2253 DO iteradj=1,iter
2254!
2255! Nutrient uptake by phytoplankton.
2256!
2257 cff1=dtdays*vm_no3(ng)
2258 DO k=1,n(ng)
2259 DO i=istr,iend
2260 cff=bio(i,k,iphyt)* &
2261 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
2262 & (k_no3(ng)+bio(i,k,ino3_))
2263 bio1(i,k,ino3_)=bio(i,k,ino3_)
2264 bio(i,k,ino3_)=bio(i,k,ino3_)/ &
2265 & (1.0_r8+cff)
2266 bio1(i,k,iphyt)=bio(i,k,iphyt)
2267 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
2268 & bio(i,k,ino3_)*cff
2269 END DO
2270 END DO
2271!
2272 IF (iteradj.ne.iter) THEN
2273!
2274! Phytoplankton grazing by Zooplankton and mortality to Detritus
2275! (rate: PhyMR).
2276!
2277 cff1=dtdays*zoogr(ng)
2278 cff2=dtdays*phymr(ng)
2279 cff3=k_phy(ng)*k_phy(ng)
2280 DO k=1,n(ng)
2281 DO i=istr,iend
2282 cff=bio(i,k,izoop)*bio(i,k,iphyt)*cff1/ &
2283 & (cff3+bio(i,k,iphyt)*bio(i,k,iphyt))
2284 bio1(i,k,iphyt)=bio(i,k,iphyt)
2285 bio(i,k,iphyt)=bio(i,k,iphyt)/ &
2286 & (1.0_r8+cff+cff2)
2287 bio1(i,k,izoop)=bio(i,k,izoop)
2288 bio(i,k,izoop)=bio(i,k,izoop)+ &
2289 & bio(i,k,iphyt)*cff*(1.0_r8-zooga(ng))
2290 bio(i,k,isdet)=bio(i,k,isdet)+ &
2291 & bio(i,k,iphyt)* &
2292 & (cff2+cff*(zooga(ng)-zooec(ng)))
2293 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2294 & bio(i,k,iphyt)*cff*zooec(ng)
2295 END DO
2296 END DO
2297!
2298! Zooplankton excretion to nutrients and mortality to Detritus.
2299!
2300 cff1=1.0_r8/(1.0_r8+dtdays*(zoomr(ng)+zoomd(ng)))
2301 cff2=dtdays*zoomr(ng)
2302 cff3=dtdays*zoomd(ng)
2303 DO k=1,n(ng)
2304 DO i=istr,iend
2305 bio(i,k,izoop)=bio(i,k,izoop)*cff1
2306 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2307 & bio(i,k,izoop)*cff2
2308 bio(i,k,isdet)=bio(i,k,isdet)+ &
2309 & bio(i,k,izoop)*cff3
2310 END DO
2311 END DO
2312!
2313! Detritus breakdown to nutrients.
2314!
2315 cff1=dtdays*detrr(ng)
2316 cff2=1.0_r8/(1.0_r8+cff1)
2317 DO k=1,n(ng)
2318 DO i=istr,iend
2319 bio(i,k,isdet)=bio(i,k,isdet)*cff2
2320 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2321 & bio(i,k,isdet)*cff1
2322 END DO
2323 END DO
2324!
2325!-----------------------------------------------------------------------
2326! Vertical sinking terms.
2327!-----------------------------------------------------------------------
2328!
2329! Reconstruct vertical profile of selected biological constituents
2330! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
2331! grid box. Then, compute semi-Lagrangian flux due to sinking.
2332!
2333 DO isink=1,nsink
2334 ibio=idsink(isink)
2335!
2336! Copy concentration of biological particulates into scratch array
2337! "qc" (q-central, restrict it to be positive) which is hereafter
2338! interpreted as a set of grid-box averaged values for biogeochemical
2339! constituent concentration.
2340!
2341 DO k=1,n(ng)
2342 DO i=istr,iend
2343 qc(i,k)=bio(i,k,ibio)
2344 END DO
2345 END DO
2346!
2347 DO k=n(ng)-1,1,-1
2348 DO i=istr,iend
2349 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2350 END DO
2351 END DO
2352 DO k=2,n(ng)-1
2353 DO i=istr,iend
2354 dltr=hz(i,j,k)*fc(i,k)
2355 dltl=hz(i,j,k)*fc(i,k-1)
2356 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2357 cffr=cff*fc(i,k)
2358 cffl=cff*fc(i,k-1)
2359!
2360! Apply PPM monotonicity constraint to prevent oscillations within the
2361! grid box.
2362!
2363 IF ((dltr*dltl).le.0.0_r8) THEN
2364 dltr=0.0_r8
2365 dltl=0.0_r8
2366 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2367 dltr=cffl
2368 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2369 dltl=cffr
2370 END IF
2371!
2372! Compute right and left side values (bR,bL) of parabolic segments
2373! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
2374!
2375! NOTE: Although each parabolic segment is monotonic within its grid
2376! box, monotonicity of the whole profile is not guaranteed,
2377! because bL(k+1)-bR(k) may still have different sign than
2378! qc(i,k+1)-qc(i,k). This possibility is excluded,
2379! after bL and bR are reconciled using WENO procedure.
2380!
2381 cff=(dltr-dltl)*hz_inv3(i,k)
2382 dltr=dltr-cff*hz(i,j,k+1)
2383 dltl=dltl+cff*hz(i,j,k-1)
2384 br(i,k)=qc(i,k)+dltr
2385 bl(i,k)=qc(i,k)-dltl
2386 wr(i,k)=(2.0_r8*dltr-dltl)**2
2387 wl(i,k)=(dltr-2.0_r8*dltl)**2
2388 END DO
2389 END DO
2390 cff=1.0e-14_r8
2391 DO k=2,n(ng)-2
2392 DO i=istr,iend
2393 dltl=max(cff,wl(i,k ))
2394 dltr=max(cff,wr(i,k+1))
2395 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2396 bl(i,k+1)=br(i,k)
2397 END DO
2398 END DO
2399 DO i=istr,iend
2400 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
2401#if defined LINEAR_CONTINUATION
2402 bl(i,n(ng))=br(i,n(ng)-1)
2403 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
2404#elif defined NEUMANN
2405 bl(i,n(ng))=br(i,n(ng)-1)
2406 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
2407#else
2408 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
2409 bl(i,n(ng))=qc(i,n(ng)) ! conditions
2410 br(i,n(ng)-1)=qc(i,n(ng))
2411#endif
2412#if defined LINEAR_CONTINUATION
2413 br(i,1)=bl(i,2)
2414 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2415#elif defined NEUMANN
2416 br(i,1)=bl(i,2)
2417 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2418#else
2419 bl(i,2)=qc(i,1) ! bottom grid boxes are
2420 br(i,1)=qc(i,1) ! re-assumed to be
2421 bl(i,1)=qc(i,1) ! piecewise constant.
2422#endif
2423 END DO
2424!
2425! Apply monotonicity constraint again, since the reconciled interfacial
2426! values may cause a non-monotonic behavior of the parabolic segments
2427! inside the grid box.
2428!
2429 DO k=1,n(ng)
2430 DO i=istr,iend
2431 dltr=br(i,k)-qc(i,k)
2432 dltl=qc(i,k)-bl(i,k)
2433 cffr=2.0_r8*dltr
2434 cffl=2.0_r8*dltl
2435 IF ((dltr*dltl).lt.0.0_r8) THEN
2436 dltr=0.0_r8
2437 dltl=0.0_r8
2438 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2439 dltr=cffl
2440 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2441 dltl=cffr
2442 END IF
2443 br(i,k)=qc(i,k)+dltr
2444 bl(i,k)=qc(i,k)-dltl
2445 END DO
2446 END DO
2447!
2448! After this moment reconstruction is considered complete. The next
2449! stage is to compute vertical advective fluxes, FC. It is expected
2450! that sinking may occurs relatively fast, the algorithm is designed
2451! to be free of CFL criterion, which is achieved by allowing
2452! integration bounds for semi-Lagrangian advective flux to use as
2453! many grid boxes in upstream direction as necessary.
2454!
2455! In the two code segments below, WL is the z-coordinate of the
2456! departure point for grid box interface z_w with the same indices;
2457! FC is the finite volume flux; ksource(:,k) is index of vertical
2458! grid box which contains the departure point (restricted by N(ng)).
2459! During the search: also add in content of whole grid boxes
2460! participating in FC.
2461!
2462 cff=dtdays*abs(wbio(isink))
2463 DO k=1,n(ng)
2464 DO i=istr,iend
2465 fc(i,k-1)=0.0_r8
2466 wl(i,k)=z_w(i,j,k-1)+cff
2467 wr(i,k)=hz(i,j,k)*qc(i,k)
2468 ksource(i,k)=k
2469 END DO
2470 END DO
2471 DO k=1,n(ng)
2472 DO ks=k,n(ng)-1
2473 DO i=istr,iend
2474 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2475 ksource(i,k)=ks+1
2476 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2477 END IF
2478 END DO
2479 END DO
2480 END DO
2481!
2482! Finalize computation of flux: add fractional part.
2483!
2484 DO k=1,n(ng)
2485 DO i=istr,iend
2486 ks=ksource(i,k)
2487 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2488 fc(i,k-1)=fc(i,k-1)+ &
2489 & hz(i,j,ks)*cu* &
2490 & (bl(i,ks)+ &
2491 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2492 & (1.5_r8-cu)* &
2493 & (br(i,ks)+bl(i,ks)- &
2494 & 2.0_r8*qc(i,ks))))
2495 END DO
2496 END DO
2497 DO k=1,n(ng)
2498 DO i=istr,iend
2499 bio(i,k,ibio)=qc(i,k)+ &
2500 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2501 END DO
2502 END DO
2503 END DO
2504 END IF
2505 END DO
2506!
2507! End of compute basic state arrays I.
2508!
2509 cff1=dtdays*vm_no3(ng)
2510 DO k=1,n(ng)
2511 DO i=istr,iend
2512 cff=bio1(i,k,iphyt)* &
2513 & cff1*exp(k_ext(ng)*z_r(i,j,k))/ &
2514 & (k_no3(ng)+bio1(i,k,ino3_))
2515!^ tl_Bio(i,k,iPhyt)=tl_Bio(i,k,iPhyt)+ &
2516!^ & tl_Bio(i,k,iNO3_)*cff+ &
2517!^ & Bio(i,k,iNO3_)*tl_cff
2518!^
2519 ad_bio(i,k,ino3_)=ad_bio(i,k,ino3_)+ &
2520 & ad_bio(i,k,iphyt)*cff
2521 ad_cff=ad_cff+bio(i,k,ino3_)*ad_bio(i,k,iphyt)
2522!^ tl_Bio(i,k,iNO3_)=(tl_Bio(i,k,iNO3_)- &
2523!^ & tl_cff*Bio(i,k,iNO3_))/ &
2524!^ & (1.0_r8+cff)
2525!^
2526 adfac=ad_bio(i,k,ino3_)/(1.0_r8+cff)
2527 ad_cff=ad_cff-bio(i,k,ino3_)*adfac
2528 ad_bio(i,k,ino3_)=adfac
2529!^ tl_cff=(tl_Bio(i,k,iPhyt)* &
2530!^ & cff1*EXP(K_ext(ng)*z_r(i,j,k))- &
2531!^ & tl_Bio(i,k,iNO3_)*cff)/ &
2532!^ & (K_NO3(ng)+Bio1(i,k,iNO3_))+ &
2533!^ & K_ext(ng)*tl_z_r(i,j,k)*cff
2534!^
2535 adfac=ad_cff/(k_no3(ng)+bio1(i,k,ino3_))
2536 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2537 & adfac*cff1*exp(k_ext(ng)*z_r(i,j,k))
2538 ad_bio(i,k,ino3_)=ad_bio(i,k,ino3_)-adfac*cff
2539 ad_z_r(i,j,k)=ad_z_r(i,j,k)+k_ext(ng)*ad_cff*cff
2540 ad_cff=0.0_r8
2541 END DO
2542 END DO
2543
2544 END DO iter_loop
2545!
2546! Adjoint Correction for negativity.
2547!
2548! Extract biological variables from tracer arrays, place them into
2549! scratch arrays, and restrict their values to be positive definite.
2550! At input, all tracers (index nnew) from predictor step have
2551! transport units (m Tunits) since we do not have yet the new
2552! values for zeta and Hz. These are known after the 2D barotropic
2553! time-stepping. In this routine, this is not a problem because
2554! we only use index nstp in the right-hand-side equations.
2555!
2556 DO itrc=1,nbt
2557 ibio=idbio(itrc)
2558 DO k=1,n(ng)
2559 DO i=istr,iend
2560 bio_old(i,k,ibio)=t(i,j,k,nstp,ibio)
2561 END DO
2562 END DO
2563 END DO
2564!
2565 DO k=1,n(ng)
2566 DO i=istr,iend
2567 cff1=max(0.0_r8,eps-bio_old(i,k,ino3_))+ &
2568 & max(0.0_r8,eps-bio_old(i,k,iphyt))+ &
2569 & max(0.0_r8,eps-bio_old(i,k,izoop))+ &
2570 & max(0.0_r8,eps-bio_old(i,k,isdet))
2571!
2572! If correction needed, determine the largest pool to debit.
2573!
2574 IF (cff1.gt.0.0) THEN
2575 itrmx=idbio(1)
2576 cff=t(i,j,k,nstp,itrmx)
2577 DO ibio=idbio(2),idbio(nbt)
2578 IF (t(i,j,k,nstp,ibio).gt.cff) THEN
2579 itrmx=ibio
2580 cff=t(i,j,k,nstp,ibio)
2581 END IF
2582 END DO
2583!
2584 DO itrc=1,nbt
2585 ibio=idbio(itrc)
2586!^ tl_Bio(i,k,ibio)=(0.5_r8- &
2587!^ & SIGN(0.5_r8,eps-Bio_old(i,k,ibio)))* &
2588!^ & tl_Bio_old(i,k,ibio)- &
2589!^ & tl_cff1* &
2590!^ & (SIGN(0.5_r8, REAL(itrmx-ibio,r8)**2)+ &
2591!^ & SIGN(0.5_r8,-REAL(itrmx-ibio,r8)**2))
2592!^
2593 ad_bio_old(i,k,ibio)=ad_bio_old(i,k,ibio)+ &
2594 & (0.5_r8- &
2595 & sign(0.5_r8, &
2596 & eps-bio_old(i,k,ibio)))* &
2597 & ad_bio(i,k,ibio)
2598 ad_cff1=ad_cff1- &
2599 & ad_bio(i,k,ibio)* &
2600 & (sign(0.5_r8, real(itrmx-ibio,r8)**2)+ &
2601 & sign(0.5_r8,-real(itrmx-ibio,r8)**2))
2602 ad_bio(i,k,ibio)=0.0_r8
2603 END DO
2604 ELSE
2605 DO itrc=1,nbt
2606 ibio=idbio(itrc)
2607!^ tl_Bio(i,k,ibio)=tl_Bio_old(i,k,ibio)
2608!^
2609 ad_bio_old(i,k,ibio)=ad_bio_old(i,k,ibio)+ &
2610 & ad_bio(i,k,ibio)
2611 ad_bio(i,k,ibio)=0.0_r8
2612 END DO
2613 END IF
2614!^ tl_cff1=-(0.5_r8-SIGN(0.5_r8,Bio_old(i,k,iNO3_)-eps))* &
2615!^ & tl_Bio_old(i,k,iNO3_)- &
2616!^ & (0.5_r8-SIGN(0.5_r8,Bio_old(i,k,iPhyt)-eps))* &
2617!^ & tl_Bio_old(i,k,iPhyt)- &
2618!^ & (0.5_r8-SIGN(0.5_r8,Bio_old(i,k,iZoop)-eps))* &
2619!^ & tl_Bio_old(i,k,iZoop)- &
2620!^ & (0.5_r8-SIGN(0.5_r8,Bio_old(i,k,iSDet)-eps))* &
2621!^ & tl_Bio_old(i,k,iSDet)
2622!^
2623 ad_bio_old(i,k,ino3_)=ad_bio_old(i,k,ino3_)- &
2624 & (0.5_r8- &
2625 & sign(0.5_r8, &
2626 & bio_old(i,k,ino3_)-eps))*ad_cff1
2627 ad_bio_old(i,k,iphyt)=ad_bio_old(i,k,iphyt)- &
2628 & (0.5_r8- &
2629 & sign(0.5_r8, &
2630 & bio_old(i,k,iphyt)-eps))*ad_cff1
2631 ad_bio_old(i,k,izoop)=ad_bio_old(i,k,izoop)- &
2632 & (0.5_r8- &
2633 & sign(0.5_r8, &
2634 & bio_old(i,k,izoop)-eps))*ad_cff1
2635 ad_bio_old(i,k,isdet)=ad_bio_old(i,k,isdet)- &
2636 & (0.5_r8- &
2637 & sign(0.5_r8, &
2638 & bio_old(i,k,isdet)-eps))*ad_cff1
2639 ad_cff1=0.0_r8
2640 END DO
2641 END DO
2642!
2643 DO itrc=1,nbt
2644 ibio=idbio(itrc)
2645 DO k=1,n(ng)
2646 DO i=istr,iend
2647!^ tl_Bio_old(i,k,ibio)=tl_t(i,j,k,nstp,ibio)
2648!^
2649 ad_t(i,j,k,nstp,ibio)=ad_t(i,j,k,nstp,ibio)+ &
2650 & ad_bio_old(i,k,ibio)
2651 ad_bio_old(i,k,ibio)=0.0_r8
2652 END DO
2653 END DO
2654 END DO
2655!
2656! Adjoint inverse thickness to avoid repeated divisions.
2657!
2658 DO k=2,n(ng)-1
2659 DO i=istr,iend
2660!^ tl_Hz_inv3(i,k)=-Hz_inv3(i,k)*Hz_inv3(i,k)* &
2661!^ & (tl_Hz(i,j,k-1)+tl_Hz(i,j,k)+ &
2662!^ & tl_Hz(i,j,k+1))
2663!^
2664 adfac=-hz_inv3(i,k)*hz_inv3(i,k)*ad_hz_inv3(i,k)
2665 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)+adfac
2666 ad_hz(i,j,k )=ad_hz(i,j,k )+adfac
2667 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)+adfac
2668 ad_hz_inv3(i,k)=0.0_r8
2669 END DO
2670 END DO
2671 DO k=1,n(ng)-1
2672 DO i=istr,iend
2673!^ tl_Hz_inv2(i,k)=-Hz_inv2(i,k)*Hz_inv2(i,k)* &
2674!^ & (tl_Hz(i,j,k)+tl_Hz(i,j,k+1))
2675!^
2676 adfac=-hz_inv2(i,k)*hz_inv2(i,k)*ad_hz_inv2(i,k)
2677 ad_hz(i,j,k )=ad_hz(i,j,k )+adfac
2678 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)+adfac
2679 ad_hz_inv2(i,k)=0.0_r8
2680 END DO
2681 END DO
2682 DO k=1,n(ng)
2683 DO i=istr,iend
2684!^ tl_Hz_inv(i,k)=-Hz_inv(i,k)*Hz_inv(i,k)*tl_Hz(i,j,k)
2685!^
2686 ad_hz(i,j,k)=ad_hz(i,j,k)- &
2687 & hz_inv(i,k)*hz_inv(i,k)*ad_hz_inv(i,k)
2688 ad_hz_inv(i,k)=0.0_r8
2689 END DO
2690 END DO
2691
2692 END DO j_loop
2693
2694!
2695! Adjoint vertical sinking velocity vector.
2696!
2697!^ tl_Wbio(1)=tl_wDet(ng) ! Small detritus
2698!^
2699 ad_wdet(ng)=ad_wdet(ng)+ad_wbio(1)
2700 ad_wbio(1)=0.0_r8
2701
2702 RETURN
2703 END SUBROUTINE ad_npzd_franks_tile
2704
2705 END MODULE ad_biology_mod
ad_biology_mod
Definition ad_biology.F:26
ad_biology_mod::ad_npzd_franks_tile
subroutine ad_npzd_franks_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nnew, rmask, hz, ad_hz, z_r, ad_z_r, z_w, ad_z_w, t, ad_t)
Definition ad_npzd_Franks.h:97
ad_biology_mod::ad_biology
subroutine, public ad_biology(ng, tile)
Definition ad_npzd_Franks.h:32
mod_biology
Definition ecosim_mod.h:1
mod_biology::zooga
real(r8), dimension(:), allocatable zooga
Definition npzd_Franks_mod.h:67
mod_biology::zoomr
real(r8), dimension(:), allocatable zoomr
Definition fennel_mod.h:162
mod_biology::k_phy
real(r8), dimension(:), allocatable k_phy
Definition fennel_mod.h:135
mod_biology::detrr
real(r8), dimension(:), allocatable detrr
Definition npzd_Franks_mod.h:51
mod_biology::wdet
real(r8), dimension(:), allocatable wdet
Definition npzd_Franks_mod.h:57
mod_biology::isdet
integer isdet
Definition npzd_Franks_mod.h:42
mod_biology::zoomd
real(r8), dimension(:), allocatable zoomd
Definition npzd_Franks_mod.h:66
mod_biology::bioiter
integer, dimension(:), allocatable bioiter
Definition ecosim_mod.h:343
mod_biology::ino3_
integer ino3_
Definition ecosim_mod.h:277
mod_biology::zoogr
real(r8), dimension(:), allocatable zoogr
Definition fennel_mod.h:160
mod_biology::k_no3
real(r8), dimension(:), allocatable k_no3
Definition fennel_mod.h:133
mod_biology::phymr
real(r8), dimension(:), allocatable phymr
Definition fennel_mod.h:145
mod_biology::iphyt
integer iphyt
Definition fennel_mod.h:81
mod_biology::zooec
real(r8), dimension(:), allocatable zooec
Definition npzd_Franks_mod.h:68
mod_biology::idbio
integer, dimension(:), allocatable idbio
Definition ecosim_mod.h:256
mod_biology::ad_wdet
real(r8), dimension(:), allocatable ad_wdet
Definition npzd_Franks_mod.h:62
mod_biology::vm_no3
real(r8), dimension(:), allocatable vm_no3
Definition npzd_Franks_mod.h:56
mod_biology::izoop
integer izoop
Definition fennel_mod.h:82
mod_biology::k_ext
real(r8), dimension(:), allocatable k_ext
Definition npzd_Franks_mod.h:52
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ncparam
Definition mod_ncparam.F:2
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_param::nbt
integer nbt
Definition mod_param.F:509
mod_param::iadm
integer, parameter iadm
Definition mod_param.F:665
mod_param::nt
integer, dimension(:), allocatable nt
Definition mod_param.F:489
mod_scalars
Definition mod_scalars.F:2
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::sec2day
real(dp), parameter sec2day
Definition mod_scalars.F:817
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
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