ROMS
Loading...
Searching...
No Matches
ad_npzd_Powell.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! 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 :: ad_biology
29!
30 CONTAINS
31!
32!***********************************************************************
33 SUBROUTINE ad_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(iadm)) THEN
58#else
59 IF (lbiofile(iadm).and.(tile.eq.0)) THEN
60#endif
61 lbiofile(iadm)=.false.
62 bioname(iadm)=myfile
63 END IF
64!
65#ifdef PROFILE
66 CALL wclock_on (ng, iadm, 15, __line__, myfile)
67#endif
68 CALL ad_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) % ad_Hz, &
77 & grid(ng) % z_r, &
78 & grid(ng) % ad_z_r, &
79 & grid(ng) % z_w, &
80 & grid(ng) % ad_z_w, &
81 & forces(ng) % srflx, &
82 & forces(ng) % ad_srflx, &
83 & ocean(ng) % t, &
84 & ocean(ng) % ad_t)
85
86#ifdef PROFILE
87 CALL wclock_off (ng, iadm, 15, __line__, myfile)
88#endif
89 RETURN
90 END SUBROUTINE ad_biology
91!
92!-----------------------------------------------------------------------
93 SUBROUTINE ad_npzd_powell_tile (ng, tile, &
94 & LBi, UBi, LBj, UBj, UBk, UBt, &
95 & IminS, ImaxS, JminS, JmaxS, &
96 & nstp, nnew, &
97#ifdef MASKING
98 & rmask, &
99#endif
100 & Hz, ad_Hz, z_r, ad_z_r, &
101 & z_w, ad_z_w, &
102 & srflx, ad_srflx, &
103 & t, ad_t)
104!-----------------------------------------------------------------------
105!
106 USE mod_param
107 USE mod_biology
108 USE mod_ncparam
109 USE mod_scalars
110!
111! Imported variable declarations.
112!
113 integer, intent(in) :: ng, tile
114 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
115 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
116 integer, intent(in) :: nstp, nnew
117
118#ifdef ASSUMED_SHAPE
119# ifdef MASKING
120 real(r8), intent(in) :: rmask(LBi:,LBj:)
121# endif
122 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
123 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
124 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
125 real(r8), intent(in) :: srflx(LBi:,LBj:)
126 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
127
128 real(r8), intent(inout) :: ad_Hz(LBi:,LBj:,:)
129 real(r8), intent(in) :: ad_z_r(LBi:,LBj:,:)
130 real(r8), intent(inout) :: ad_z_w(LBi:,LBj:,0:)
131 real(r8), intent(inout) :: ad_srflx(LBi:,LBj:)
132 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
133#else
134# ifdef MASKING
135 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
136# endif
137 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
138 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
139 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
140 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
141 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
142
143 real(r8), intent(inout) :: ad_Hz(LBi:UBi,LBj:UBj,UBk)
144 real(r8), intent(in) :: ad_z_r(LBi:UBi,LBj:UBj,UBk)
145 real(r8), intent(inout) :: ad_z_w(LBi:UBi,LBj:UBj,0:UBk)
146 real(r8), intent(inout) :: ad_srflx(LBi:UBi,LBj:UBj)
147 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,UBk,3,UBt)
148#endif
149!
150! Local variable declarations.
151!
152 integer, parameter :: Nsink = 2
153
154 integer :: Iter, i, ibio, isink, itime, itrc, iTrcMax, j, k, ks
155 integer :: Iteradj, kk
156
157 integer, dimension(Nsink) :: idsink
158
159 real(r8), parameter :: MinVal = 1.0e-6_r8
160
161 real(r8) :: Att, ExpAtt, Itop, PAR, PAR1
162 real(r8) :: ad_Att, ad_ExpAtt, ad_Itop, ad_PAR
163 real(r8) :: cff, cff1, cff2, cff3, cff4, dtdays
164 real(r8) :: ad_cff, ad_cff1, ad_cff4
165 real(r8) :: cffL, cffR, cu, dltL, dltR
166 real(r8) :: ad_cffL, ad_cffR, ad_cu, ad_dltL, ad_dltR
167 real(r8) :: fac, adfac, adfac1, adfac2, adfac3
168
169 real(r8), dimension(Nsink) :: Wbio
170 real(r8), dimension(Nsink) :: ad_Wbio
171
172 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
173
174 real(r8), dimension(IminS:ImaxS) :: PARsur
175 real(r8), dimension(IminS:ImaxS) :: ad_PARsur
176
177 real(r8), dimension(NT(ng),2) :: BioTrc
178 real(r8), dimension(NT(ng),2) :: BioTrc1
179 real(r8), dimension(NT(ng),2) :: ad_BioTrc
180 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
181 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio1
182 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
183
184 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: ad_Bio
185 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: ad_Bio_old
186
187 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
188 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_FC
189
190 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
191 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
192 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
193 real(r8), dimension(IminS:ImaxS,N(ng)) :: Light
194 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
195 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
196 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
197 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL1
198 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
199 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR1
200 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
201
202 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv
203 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv2
204 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hz_inv3
205 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Light
206 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_WL
207 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_WR
208 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_bL
209 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_bR
210 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_qc
211
212#include "set_bounds.h"
213!
214!-----------------------------------------------------------------------
215! Add biological Source/Sink terms.
216!-----------------------------------------------------------------------
217!
218! Avoid computing source/sink terms if no biological iterations.
219!
220 IF (bioiter(ng).le.0) RETURN
221!
222! Set time-stepping size (days) according to the number of iterations.
223!
224 dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
225!
226! Set vertical sinking indentification vector.
227!
228 idsink(1)=iphyt ! Phytoplankton
229 idsink(2)=isdet ! Small detritus
230!
231! Set vertical sinking velocity vector in the same order as the
232! identification vector, IDSINK.
233!
234 wbio(1)=wphy(ng) ! Phytoplankton
235 wbio(2)=wdet(ng) ! Small detritus
236!
237 ad_wbio(1)=0.0_r8
238 ad_wbio(2)=0.0_r8
239!
240 j_loop : DO j=jstr,jend
241!
242!-----------------------------------------------------------------------
243! Initialize adjoint private variables.
244!-----------------------------------------------------------------------
245!
246 ad_att=0.0_r8
247 ad_expatt=0.0_r8
248 ad_itop=0.0_r8
249 ad_par=0.0_r8
250 ad_dltl=0.0_r8
251 ad_dltr=0.0_r8
252 ad_cu=0.0_r8
253 ad_cff=0.0_r8
254 ad_cff1=0.0_r8
255 ad_cff4=0.0_r8
256 ad_cffl=0.0_r8
257 ad_cffr=0.0_r8
258 adfac=0.0_r8
259 adfac1=0.0_r8
260 adfac2=0.0_r8
261 adfac3=0.0_r8
262!
263 DO k=1,n(ng)
264 DO i=imins,imaxs
265 ad_hz_inv(i,k)=0.0_r8
266 ad_hz_inv2(i,k)=0.0_r8
267 ad_hz_inv3(i,k)=0.0_r8
268 ad_wl(i,k)=0.0_r8
269 ad_wr(i,k)=0.0_r8
270 ad_bl(i,k)=0.0_r8
271 ad_br(i,k)=0.0_r8
272 ad_qc(i,k)=0.0_r8
273 ad_light(i,k)=0.0_r8
274 END DO
275 END DO
276 DO itrc=1,nbt
277 ibio=idbio(itrc)
278 ad_biotrc(ibio,1)=0.0_r8
279 ad_biotrc(ibio,2)=0.0_r8
280 END DO
281 DO i=imins,imaxs
282 ad_parsur(i)=0.0_r8
283 END DO
284 DO k=0,n(ng)
285 DO i=imins,imaxs
286 ad_fc(i,k)=0.0_r8
287 END DO
288 END DO
289!
290! Clear ad_Bio and Bio arrays.
291!
292 DO itrc=1,nbt
293 ibio=idbio(itrc)
294 DO k=1,n(ng)
295 DO i=istr,iend
296 bio(i,k,ibio)=0.0_r8
297 bio1(i,k,ibio)=0.0_r8
298 ad_bio(i,k,ibio)=0.0_r8
299 ad_bio_old(i,k,ibio)=0.0_r8
300 END DO
301 END DO
302 END DO
303!
304! Compute inverse thickness to avoid repeated divisions.
305!
306 DO k=1,n(ng)
307 DO i=istr,iend
308 hz_inv(i,k)=1.0_r8/hz(i,j,k)
309 END DO
310 END DO
311 DO k=1,n(ng)-1
312 DO i=istr,iend
313 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
314 END DO
315 END DO
316 DO k=2,n(ng)-1
317 DO i=istr,iend
318 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
319 END DO
320 END DO
321!
322! Compute the required basic state arrays.
323!
324! Restrict biological tracer to be positive definite. If a negative
325! concentration is detected, nitrogen is drawn from the most abundant
326! pool to supplement the negative pools to a lower limit of MinVal
327! which is set to 1E-6 above.
328!
329 DO k=1,n(ng)
330 DO i=istr,iend
331!
332! At input, all tracers (index nnew) from predictor step have
333! transport units (m Tunits) since we do not have yet the new
334! values for zeta and Hz. These are known after the 2D barotropic
335! time-stepping.
336!
337! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
338! tracer times depth. However the basic state (nstp and nnew
339! indices) that is read from the forward file is in units of
340! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
341! use t instead of t*Hz_inv.
342!
343 DO itrc=1,nbt
344 ibio=idbio(itrc)
345!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
346!^
347 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
348!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
349!^
350 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
351 END DO
352!
353! Impose positive definite concentrations.
354!
355 cff2=0.0_r8
356 DO itime=1,2
357 cff1=0.0_r8
358 itrcmax=idbio(1)
359 DO itrc=1,nbt
360 ibio=idbio(itrc)
361 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
362 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
363 itrcmax=ibio
364 END IF
365 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
366 END DO
367 IF (biotrc(itrcmax,itime).gt.cff1) THEN
368 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
369 END IF
370 END DO
371!
372! Load biological tracers into local arrays.
373!
374 DO itrc=1,nbt
375 ibio=idbio(itrc)
376 bio_old(i,k,ibio)=biotrc(ibio,nstp)
377 bio(i,k,ibio)=biotrc(ibio,nstp)
378 END DO
379 END DO
380 END DO
381!
382! Calculate surface Photosynthetically Available Radiation (PAR). The
383! net shortwave radiation is scaled back to Watts/m2 and multiplied by
384! the fraction that is photosynthetically available, PARfrac.
385!
386 DO i=istr,iend
387#ifdef CONST_PAR
388!
389! Specify constant surface irradiance a la Powell and Spitz.
390!
391 parsur(i)=158.075_r8
392#else
393 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
394#endif
395 END DO
396!
397!=======================================================================
398! Start internal iterations to achieve convergence of the nonlinear
399! backward-implicit solution.
400!=======================================================================
401!
402! During the iterative procedure a series of fractional time steps are
403! performed in a chained mode (splitting by different biological
404! conversion processes) in sequence of the main food chain. In all
405! stages the concentration of the component being consumed is treated
406! in a fully implicit manner, so the algorithm guarantees non-negative
407! values, no matter how strong the concentration of active consuming
408! component (Phytoplankton or Zooplankton). The overall algorithm,
409! as well as any stage of it, is formulated in conservative form
410! (except explicit sinking) in sense that the sum of concentration of
411! all components is conserved.
412!
413! In the implicit algorithm, we have for example (N: nutrient,
414! P: phytoplankton),
415!
416! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
417! {Michaelis-Menten}
418! below, we set
419! The N in the numerator of
420! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
421! as N(new)
422!
423! so the time-stepping of the equations becomes:
424!
425! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
426! consuming, divide by (1 + cff)
427! and
428!
429! P(new) = P(old) + cff * N(new) (2) when adding a source term,
430! growing, add (cff * source)
431!
432! Notice that if you substitute (1) in (2), you will get:
433!
434! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
435!
436! If you add (1) and (3), you get
437!
438! N(new) + P(new) = N(old) + P(old)
439!
440! implying conservation regardless how "cff" is computed. Therefore,
441! this scheme is unconditionally stable regardless of the conversion
442! rate. It does not generate negative values since the constituent
443! to be consumed is always treated implicitly. It is also biased
444! toward damping oscillations.
445!
446! The iterative loop below is to iterate toward an universal Backward-
447! Euler treatment of all terms. So if there are oscillations in the
448! system, they are only physical oscillations. These iterations,
449! however, do not improve the accuaracy of the solution.
450!
451 iter_loop: DO iter=1,bioiter(ng)
452!
453! Compute light attenuation as function of depth.
454!
455 DO i=istr,iend
456 par=parsur(i)
457 IF (parsur(i).gt.0.0_r8) THEN ! day time
458 DO k=n(ng),1,-1
459!
460! Compute average light attenuation for each grid cell. Here, AttSW is
461! the light attenuation due to seawater and AttPhy is the attenuation
462! due to phytoplankton (self-shading coefficient).
463!
464 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
465 & (z_w(i,j,k)-z_w(i,j,k-1))
466 expatt=exp(-att)
467 itop=par
468 par=itop*(1.0_r8-expatt)/att ! average at cell center
469 light(i,k)=par
470!
471! Light attenuation at the bottom of the grid cell. It is the starting
472! PAR value for the next (deeper) vertical grid cell.
473!
474 par=itop*expatt
475 END DO
476 ELSE ! night time
477 DO k=1,n(ng)
478 light(i,k)=0.0_r8
479 END DO
480 END IF
481 END DO
482!
483! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
484! The Michaelis-Menten curve is used to describe the change in uptake
485! rate as a function of nitrate concentration. Here, PhyIS is the
486! initial slope of the P-I curve and K_NO3 is the half saturation of
487! phytoplankton nitrate uptake.
488!
489 cff1=dtdays*vm_no3(ng)*phyis(ng)
490 cff2=vm_no3(ng)*vm_no3(ng)
491 cff3=phyis(ng)*phyis(ng)
492 DO k=1,n(ng)
493 DO i=istr,iend
494 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
495 cff=bio(i,k,iphyt)* &
496 & cff1*cff4*light(i,k)/ &
497 & (k_no3(ng)+bio(i,k,ino3_))
498 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
499 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
500 & bio(i,k,ino3_)*cff
501 END DO
502 END DO
503!
504! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
505! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
506! pool as function of "sloppy feeding" and metabolic processes
507! (ZooEEN and ZooEED fractions).
508!
509 cff1=dtdays*zoogr(ng)
510 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
511 DO k=1,n(ng)
512 DO i=istr,iend
513 cff=bio(i,k,izoop)* &
514 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
515 & bio(i,k,iphyt)
516 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
517 bio(i,k,izoop)=bio(i,k,izoop)+ &
518 & bio(i,k,iphyt)*cff2*cff
519 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
520 & bio(i,k,iphyt)*zooeen(ng)*cff
521 bio(i,k,isdet)=bio(i,k,isdet)+ &
522 & bio(i,k,iphyt)*zooeed(ng)*cff
523 END DO
524 END DO
525!
526! Phytoplankton mortality to nutrients (PhyMRN rate) and detritus
527! (PhyMRD rate).
528!
529 cff3=dtdays*phymrd(ng)
530 cff2=dtdays*phymrn(ng)
531 cff1=1.0_r8/(1.0_r8+cff2+cff3)
532 DO k=1,n(ng)
533 DO i=istr,iend
534 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
535 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
536 & bio(i,k,iphyt)*cff2
537 bio(i,k,isdet)=bio(i,k,isdet)+ &
538 & bio(i,k,iphyt)*cff3
539 END DO
540 END DO
541!
542! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
543! (ZooMRD rate).
544!
545 cff3=dtdays*zoomrd(ng)
546 cff2=dtdays*zoomrn(ng)
547 cff1=1.0_r8/(1.0_r8+cff2+cff3)
548 DO k=1,n(ng)
549 DO i=istr,iend
550 bio(i,k,izoop)=bio(i,k,izoop)*cff1
551 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
552 & bio(i,k,izoop)*cff2
553 bio(i,k,isdet)=bio(i,k,isdet)+ &
554 & bio(i,k,izoop)*cff3
555 END DO
556 END DO
557!
558! Detritus breakdown to nutrients: remineralization (DetRR rate).
559!
560 cff2=dtdays*detrr(ng)
561 cff1=1.0_r8/(1.0_r8+cff2)
562 DO k=1,n(ng)
563 DO i=istr,iend
564 bio(i,k,isdet)=bio(i,k,isdet)*cff1
565 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
566 & bio(i,k,isdet)*cff2
567 END DO
568 END DO
569!
570!-----------------------------------------------------------------------
571! Vertical sinking terms: Phytoplankton and Detritus
572!-----------------------------------------------------------------------
573!
574! Reconstruct vertical profile of selected biological constituents
575! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
576! grid box. Then, compute semi-Lagrangian flux due to sinking.
577!
578 sink_loop: DO isink=1,nsink
579 ibio=idsink(isink)
580!
581! Copy concentration of biological particulates into scratch array
582! "qc" (q-central, restrict it to be positive) which is hereafter
583! interpreted as a set of grid-box averaged values for biogeochemical
584! constituent concentration.
585!
586 DO k=1,n(ng)
587 DO i=istr,iend
588 qc(i,k)=bio(i,k,ibio)
589 END DO
590 END DO
591!
592 DO k=n(ng)-1,1,-1
593 DO i=istr,iend
594 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
595 END DO
596 END DO
597 DO k=2,n(ng)-1
598 DO i=istr,iend
599 dltr=hz(i,j,k)*fc(i,k)
600 dltl=hz(i,j,k)*fc(i,k-1)
601 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
602 cffr=cff*fc(i,k)
603 cffl=cff*fc(i,k-1)
604!
605! Apply PPM monotonicity constraint to prevent oscillations within the
606! grid box.
607!
608 IF ((dltr*dltl).le.0.0_r8) THEN
609 dltr=0.0_r8
610 dltl=0.0_r8
611 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
612 dltr=cffl
613 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
614 dltl=cffr
615 END IF
616!
617! Compute right and left side values (bR,bL) of parabolic segments
618! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
619!
620! NOTE: Although each parabolic segment is monotonic within its grid
621! box, monotonicity of the whole profile is not guaranteed,
622! because bL(k+1)-bR(k) may still have different sign than
623! qc(i,k+1)-qc(i,k). This possibility is excluded,
624! after bL and bR are reconciled using WENO procedure.
625!
626 cff=(dltr-dltl)*hz_inv3(i,k)
627 dltr=dltr-cff*hz(i,j,k+1)
628 dltl=dltl+cff*hz(i,j,k-1)
629 br(i,k)=qc(i,k)+dltr
630 bl(i,k)=qc(i,k)-dltl
631 wr(i,k)=(2.0_r8*dltr-dltl)**2
632 wl(i,k)=(dltr-2.0_r8*dltl)**2
633 END DO
634 END DO
635 cff=1.0e-14_r8
636 DO k=2,n(ng)-2
637 DO i=istr,iend
638 dltl=max(cff,wl(i,k ))
639 dltr=max(cff,wr(i,k+1))
640 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
641 bl(i,k+1)=br(i,k)
642 END DO
643 END DO
644 DO i=istr,iend
645 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
646#if defined LINEAR_CONTINUATION
647 bl(i,n(ng))=br(i,n(ng)-1)
648 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
649#elif defined NEUMANN
650 bl(i,n(ng))=br(i,n(ng)-1)
651 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
652#else
653 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
654 bl(i,n(ng))=qc(i,n(ng)) ! conditions
655 br(i,n(ng)-1)=qc(i,n(ng))
656#endif
657#if defined LINEAR_CONTINUATION
658 br(i,1)=bl(i,2)
659 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
660#elif defined NEUMANN
661 br(i,1)=bl(i,2)
662 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
663#else
664 bl(i,2)=qc(i,1) ! bottom grid boxes are
665 br(i,1)=qc(i,1) ! re-assumed to be
666 bl(i,1)=qc(i,1) ! piecewise constant.
667#endif
668 END DO
669!
670! Apply monotonicity constraint again, since the reconciled interfacial
671! values may cause a non-monotonic behavior of the parabolic segments
672! inside the grid box.
673!
674 DO k=1,n(ng)
675 DO i=istr,iend
676 dltr=br(i,k)-qc(i,k)
677 dltl=qc(i,k)-bl(i,k)
678 cffr=2.0_r8*dltr
679 cffl=2.0_r8*dltl
680 IF ((dltr*dltl).lt.0.0_r8) THEN
681 dltr=0.0_r8
682 dltl=0.0_r8
683 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
684 dltr=cffl
685 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
686 dltl=cffr
687 END IF
688 br(i,k)=qc(i,k)+dltr
689 bl(i,k)=qc(i,k)-dltl
690 END DO
691 END DO
692!
693! After this moment reconstruction is considered complete. The next
694! stage is to compute vertical advective fluxes, FC. It is expected
695! that sinking may occurs relatively fast, the algorithm is designed
696! to be free of CFL criterion, which is achieved by allowing
697! integration bounds for semi-Lagrangian advective flux to use as
698! many grid boxes in upstream direction as necessary.
699!
700! In the two code segments below, WL is the z-coordinate of the
701! departure point for grid box interface z_w with the same indices;
702! FC is the finite volume flux; ksource(:,k) is index of vertical
703! grid box which contains the departure point (restricted by N(ng)).
704! During the search: also add in content of whole grid boxes
705! participating in FC.
706!
707 cff=dtdays*abs(wbio(isink))
708 DO k=1,n(ng)
709 DO i=istr,iend
710 fc(i,k-1)=0.0_r8
711 wl(i,k)=z_w(i,j,k-1)+cff
712 wr(i,k)=hz(i,j,k)*qc(i,k)
713 ksource(i,k)=k
714 END DO
715 END DO
716 DO k=1,n(ng)
717 DO ks=k,n(ng)-1
718 DO i=istr,iend
719 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
720 ksource(i,k)=ks+1
721 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
722 END IF
723 END DO
724 END DO
725 END DO
726!
727! Finalize computation of flux: add fractional part.
728!
729 DO k=1,n(ng)
730 DO i=istr,iend
731 ks=ksource(i,k)
732 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
733 fc(i,k-1)=fc(i,k-1)+ &
734 & hz(i,j,ks)*cu* &
735 & (bl(i,ks)+ &
736 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
737 & (1.5_r8-cu)* &
738 & (br(i,ks)+bl(i,ks)- &
739 & 2.0_r8*qc(i,ks))))
740 END DO
741 END DO
742 DO k=1,n(ng)
743 DO i=istr,iend
744 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
745 END DO
746 END DO
747
748 END DO sink_loop
749 END DO iter_loop
750!
751!-----------------------------------------------------------------------
752! Update global tracer variables: Add increment due to BGC processes
753! to tracer array in time index "nnew". Index "nnew" is solution after
754! advection and mixing and has transport units (m Tunits) hence the
755! increment is multiplied by Hz. Notice that we need to subtract
756! original values "Bio_old" at the top of the routine to just account
757! for the concentractions affected by BGC processes. This also takes
758! into account any constraints (non-negative concentrations, carbon
759! concentration range) specified before entering BGC kernel. If "Bio"
760! were unchanged by BGC processes, the increment would be exactly
761! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
762! bounded >=0 so that we can preserve total inventory of nutrients
763! when advection causes tracer concentration to go negative.
764!-----------------------------------------------------------------------
765!
766 DO itrc=1,nbt
767 ibio=idbio(itrc)
768 DO k=1,n(ng)
769 DO i=istr,iend
770 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
771!^ tl_t(i,j,k,nnew,ibio)=tl_t(i,j,k,nnew,ibio)+ &
772!^ & tl_cff*Hz(i,j,k)+cff*tl_Hz(i,j,k)
773!^
774 ad_hz(i,j,k)=ad_hz(i,j,k)+cff*ad_t(i,j,k,nnew,ibio)
775 ad_cff=ad_cff+hz(i,j,k)*ad_t(i,j,k,nnew,ibio)
776!^ tl_cff=tl_Bio(i,k,ibio)-tl_Bio_old(i,k,ibio)
777!^
778 ad_bio_old(i,k,ibio)=ad_bio_old(i,k,ibio)-ad_cff
779 ad_bio(i,k,ibio)=ad_bio(i,k,ibio)+ad_cff
780 END DO
781 END DO
782 END DO
783!
784!=======================================================================
785! Start internal iterations to achieve convergence of the nonlinear
786! backward-implicit solution.
787!=======================================================================
788!
789 iter_loop1: DO iter=bioiter(ng),1,-1
790!
791!-----------------------------------------------------------------------
792! Adjoint vertical sinking terms.
793!-----------------------------------------------------------------------
794!
795! Reconstruct vertical profile of selected biological constituents
796! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
797! grid box. Then, compute semi-Lagrangian flux due to sinking.
798!
799! Compute appropriate basic state arrays III.
800!
801 DO k=1,n(ng)
802 DO i=istr,iend
803!
804! At input, all tracers (index nnew) from predictor step have
805! transport units (m Tunits) since we do not have yet the new
806! values for zeta and Hz. These are known after the 2D barotropic
807! time-stepping.
808!
809! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
810! tracer times depth. However the basic state (nstp and nnew
811! indices) that is read from the forward file is in units of
812! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
813! use t instead of t*Hz_inv.
814!
815 DO itrc=1,nbt
816 ibio=idbio(itrc)
817!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
818!^
819 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
820!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
821!^
822 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
823 END DO
824!
825! Impose positive definite concentrations.
826!
827 cff2=0.0_r8
828 DO itime=1,2
829 cff1=0.0_r8
830 itrcmax=idbio(1)
831 DO itrc=1,nbt
832 ibio=idbio(itrc)
833 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
834 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
835 itrcmax=ibio
836 END IF
837 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
838 END DO
839 IF (biotrc(itrcmax,itime).gt.cff1) THEN
840 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
841 END IF
842 END DO
843!
844! Load biological tracers into local arrays.
845!
846 DO itrc=1,nbt
847 ibio=idbio(itrc)
848 bio_old(i,k,ibio)=biotrc(ibio,nstp)
849 bio(i,k,ibio)=biotrc(ibio,nstp)
850 END DO
851 END DO
852 END DO
853!
854! Calculate surface Photosynthetically Available Radiation (PAR). The
855! net shortwave radiation is scaled back to Watts/m2 and multiplied by
856! the fraction that is photosynthetically available, PARfrac.
857!
858 DO i=istr,iend
859#ifdef CONST_PAR
860!
861! Specify constant surface irradiance a la Powell and Spitz.
862!
863 parsur(i)=158.075_r8
864#else
865 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
866#endif
867 END DO
868!
869!=======================================================================
870! Start internal iterations to achieve convergence of the nonlinear
871! backward-implicit solution.
872!=======================================================================
873!
874 DO iteradj=1,iter
875!
876! Compute light attenuation as function of depth.
877!
878 DO i=istr,iend
879 par=parsur(i)
880 IF (parsur(i).gt.0.0_r8) THEN ! day time
881 DO k=n(ng),1,-1
882!
883! Compute average light attenuation for each grid cell. Here, AttSW is
884! the light attenuation due to seawater and AttPhy is the attenuation
885! due to phytoplankton (self-shading coefficient).
886!
887 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
888 & (z_w(i,j,k)-z_w(i,j,k-1))
889 expatt=exp(-att)
890 itop=par
891 par=itop*(1.0_r8-expatt)/att ! average at cell center
892 light(i,k)=par
893!
894! Light attenuation at the bottom of the grid cell. It is the starting
895! PAR value for the next (deeper) vertical grid cell.
896!
897 par=itop*expatt
898 END DO
899 ELSE ! night time
900 DO k=1,n(ng)
901 light(i,k)=0.0_r8
902 END DO
903 END IF
904 END DO
905!
906! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
907! The Michaelis-Menten curve is used to describe the change in uptake
908! rate as a function of nitrate concentration. Here, PhyIS is the
909! initial slope of the P-I curve and K_NO3 is the half saturation of
910! phytoplankton nitrate uptake.
911!
912 cff1=dtdays*vm_no3(ng)*phyis(ng)
913 cff2=vm_no3(ng)*vm_no3(ng)
914 cff3=phyis(ng)*phyis(ng)
915 DO k=1,n(ng)
916 DO i=istr,iend
917 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
918 cff=bio(i,k,iphyt)* &
919 & cff1*cff4*light(i,k)/ &
920 & (k_no3(ng)+bio(i,k,ino3_))
921 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
922 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
923 & bio(i,k,ino3_)*cff
924 END DO
925 END DO
926!
927! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
928! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
929! pool as function of "sloppy feeding" and metabolic processes
930! (ZooEEN and ZooEED fractions).
931!
932 cff1=dtdays*zoogr(ng)
933 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
934 DO k=1,n(ng)
935 DO i=istr,iend
936 cff=bio(i,k,izoop)* &
937 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
938 & bio(i,k,iphyt)
939 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
940 bio(i,k,izoop)=bio(i,k,izoop)+ &
941 & bio(i,k,iphyt)*cff2*cff
942 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
943 & bio(i,k,iphyt)*zooeen(ng)*cff
944 bio(i,k,isdet)=bio(i,k,isdet)+ &
945 & bio(i,k,iphyt)*zooeed(ng)*cff
946 END DO
947 END DO
948!
949! Phytoplankton mortality to nutrients (PhyMRN rate) and detritus
950! (PhyMRD rate).
951!
952 cff3=dtdays*phymrd(ng)
953 cff2=dtdays*phymrn(ng)
954 cff1=1.0_r8/(1.0_r8+cff2+cff3)
955 DO k=1,n(ng)
956 DO i=istr,iend
957 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
958 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
959 & bio(i,k,iphyt)*cff2
960 bio(i,k,isdet)=bio(i,k,isdet)+ &
961 & bio(i,k,iphyt)*cff3
962 END DO
963 END DO
964!
965! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
966! (ZooMRD rate).
967!
968 cff3=dtdays*zoomrd(ng)
969 cff2=dtdays*zoomrn(ng)
970 cff1=1.0_r8/(1.0_r8+cff2+cff3)
971 DO k=1,n(ng)
972 DO i=istr,iend
973 bio(i,k,izoop)=bio(i,k,izoop)*cff1
974 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
975 & bio(i,k,izoop)*cff2
976 bio(i,k,isdet)=bio(i,k,isdet)+ &
977 & bio(i,k,izoop)*cff3
978 END DO
979 END DO
980!
981! Detritus breakdown to nutrients: remineralization (DetRR rate).
982!
983 cff2=dtdays*detrr(ng)
984 cff1=1.0_r8/(1.0_r8+cff2)
985 DO k=1,n(ng)
986 DO i=istr,iend
987 bio(i,k,isdet)=bio(i,k,isdet)*cff1
988 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
989 & bio(i,k,isdet)*cff2
990 END DO
991 END DO
992!
993 IF (iteradj.ne.iter) THEN
994!
995!-----------------------------------------------------------------------
996! Vertical sinking terms: Phytoplankton and Detritus
997!-----------------------------------------------------------------------
998!
999! Reconstruct vertical profile of selected biological constituents
1000! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
1001! grid box. Then, compute semi-Lagrangian flux due to sinking.
1002!
1003 DO isink=1,nsink
1004 ibio=idsink(isink)
1005!
1006! Copy concentration of biological particulates into scratch array
1007! "qc" (q-central, restrict it to be positive) which is hereafter
1008! interpreted as a set of grid-box averaged values for biogeochemical
1009! constituent concentration.
1010!
1011 DO k=1,n(ng)
1012 DO i=istr,iend
1013 qc(i,k)=bio(i,k,ibio)
1014 END DO
1015 END DO
1016!
1017 DO k=n(ng)-1,1,-1
1018 DO i=istr,iend
1019 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1020 END DO
1021 END DO
1022 DO k=2,n(ng)-1
1023 DO i=istr,iend
1024 dltr=hz(i,j,k)*fc(i,k)
1025 dltl=hz(i,j,k)*fc(i,k-1)
1026 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1027 cffr=cff*fc(i,k)
1028 cffl=cff*fc(i,k-1)
1029!
1030! Apply PPM monotonicity constraint to prevent oscillations within the
1031! grid box.
1032!
1033 IF ((dltr*dltl).le.0.0_r8) THEN
1034 dltr=0.0_r8
1035 dltl=0.0_r8
1036 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1037 dltr=cffl
1038 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1039 dltl=cffr
1040 END IF
1041!
1042! Compute right and left side values (bR,bL) of parabolic segments
1043! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1044!
1045! NOTE: Although each parabolic segment is monotonic within its grid
1046! box, monotonicity of the whole profile is not guaranteed,
1047! because bL(k+1)-bR(k) may still have different sign than
1048! qc(i,k+1)-qc(i,k). This possibility is excluded,
1049! after bL and bR are reconciled using WENO procedure.
1050!
1051 cff=(dltr-dltl)*hz_inv3(i,k)
1052 dltr=dltr-cff*hz(i,j,k+1)
1053 dltl=dltl+cff*hz(i,j,k-1)
1054 br(i,k)=qc(i,k)+dltr
1055 bl(i,k)=qc(i,k)-dltl
1056 wr(i,k)=(2.0_r8*dltr-dltl)**2
1057 wl(i,k)=(dltr-2.0_r8*dltl)**2
1058 END DO
1059 END DO
1060 cff=1.0e-14_r8
1061 DO k=2,n(ng)-2
1062 DO i=istr,iend
1063 dltl=max(cff,wl(i,k ))
1064 dltr=max(cff,wr(i,k+1))
1065 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1066 bl(i,k+1)=br(i,k)
1067 END DO
1068 END DO
1069 DO i=istr,iend
1070 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1071#if defined LINEAR_CONTINUATION
1072 bl(i,n(ng))=br(i,n(ng)-1)
1073 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1074#elif defined NEUMANN
1075 bl(i,n(ng))=br(i,n(ng)-1)
1076 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1077#else
1078 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1079 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1080 br(i,n(ng)-1)=qc(i,n(ng))
1081#endif
1082#if defined LINEAR_CONTINUATION
1083 br(i,1)=bl(i,2)
1084 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1085#elif defined NEUMANN
1086 br(i,1)=bl(i,2)
1087 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1088#else
1089 bl(i,2)=qc(i,1) ! bottom grid boxes are
1090 br(i,1)=qc(i,1) ! re-assumed to be
1091 bl(i,1)=qc(i,1) ! piecewise constant.
1092#endif
1093 END DO
1094!
1095! Apply monotonicity constraint again, since the reconciled interfacial
1096! values may cause a non-monotonic behavior of the parabolic segments
1097! inside the grid box.
1098!
1099 DO k=1,n(ng)
1100 DO i=istr,iend
1101 dltr=br(i,k)-qc(i,k)
1102 dltl=qc(i,k)-bl(i,k)
1103 cffr=2.0_r8*dltr
1104 cffl=2.0_r8*dltl
1105 IF ((dltr*dltl).lt.0.0_r8) THEN
1106 dltr=0.0_r8
1107 dltl=0.0_r8
1108 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1109 dltr=cffl
1110 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1111 dltl=cffr
1112 END IF
1113 br(i,k)=qc(i,k)+dltr
1114 bl(i,k)=qc(i,k)-dltl
1115 END DO
1116 END DO
1117!
1118! After this moment reconstruction is considered complete. The next
1119! stage is to compute vertical advective fluxes, FC. It is expected
1120! that sinking may occurs relatively fast, the algorithm is designed
1121! to be free of CFL criterion, which is achieved by allowing
1122! integration bounds for semi-Lagrangian advective flux to use as
1123! many grid boxes in upstream direction as necessary.
1124!
1125! In the two code segments below, WL is the z-coordinate of the
1126! departure point for grid box interface z_w with the same indices;
1127! FC is the finite volume flux; ksource(:,k) is index of vertical
1128! grid box which contains the departure point (restricted by N(ng)).
1129! During the search: also add in content of whole grid boxes
1130! participating in FC.
1131!
1132 cff=dtdays*abs(wbio(isink))
1133 DO k=1,n(ng)
1134 DO i=istr,iend
1135 fc(i,k-1)=0.0_r8
1136 wl(i,k)=z_w(i,j,k-1)+cff
1137 wr(i,k)=hz(i,j,k)*qc(i,k)
1138 ksource(i,k)=k
1139 END DO
1140 END DO
1141 DO k=1,n(ng)
1142 DO ks=k,n(ng)-1
1143 DO i=istr,iend
1144 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1145 ksource(i,k)=ks+1
1146 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1147 END IF
1148 END DO
1149 END DO
1150 END DO
1151!
1152! Finalize computation of flux: add fractional part.
1153!
1154 DO k=1,n(ng)
1155 DO i=istr,iend
1156 ks=ksource(i,k)
1157 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1158 fc(i,k-1)=fc(i,k-1)+ &
1159 & hz(i,j,ks)*cu* &
1160 & (bl(i,ks)+ &
1161 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1162 & (1.5_r8-cu)* &
1163 & (br(i,ks)+bl(i,ks)- &
1164 & 2.0_r8*qc(i,ks))))
1165 END DO
1166 END DO
1167 DO k=1,n(ng)
1168 DO i=istr,iend
1169 bio(i,k,ibio)=qc(i,k)+ &
1170 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1171 END DO
1172 END DO
1173 END DO
1174 END IF
1175 END DO
1176!
1177! End of compute basic state arrays III.
1178!
1179 sink_loop1: DO isink=1,nsink
1180 ibio=idsink(isink)
1181
1182!
1183! Compute required flux arrays.
1184!
1185! Copy concentration of biological particulates into scratch array
1186! "qc" (q-central, restrict it to be positive) which is hereafter
1187! interpreted as a set of grid-box averaged values for biogeochemical
1188! constituent concentration.
1189!
1190 DO k=1,n(ng)
1191 DO i=istr,iend
1192 qc(i,k)=bio(i,k,ibio)
1193 END DO
1194 END DO
1195!
1196 DO k=n(ng)-1,1,-1
1197 DO i=istr,iend
1198 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1199 END DO
1200 END DO
1201 DO k=2,n(ng)-1
1202 DO i=istr,iend
1203 dltr=hz(i,j,k)*fc(i,k)
1204 dltl=hz(i,j,k)*fc(i,k-1)
1205 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1206 cffr=cff*fc(i,k)
1207 cffl=cff*fc(i,k-1)
1208!
1209! Apply PPM monotonicity constraint to prevent oscillations within the
1210! grid box.
1211!
1212 IF ((dltr*dltl).le.0.0_r8) THEN
1213 dltr=0.0_r8
1214 dltl=0.0_r8
1215 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1216 dltr=cffl
1217 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1218 dltl=cffr
1219 END IF
1220!
1221! Compute right and left side values (bR,bL) of parabolic segments
1222! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1223!
1224! NOTE: Although each parabolic segment is monotonic within its grid
1225! box, monotonicity of the whole profile is not guaranteed,
1226! because bL(k+1)-bR(k) may still have different sign than
1227! qc(i,k+1)-qc(i,k). This possibility is excluded,
1228! after bL and bR are reconciled using WENO procedure.
1229!
1230 cff=(dltr-dltl)*hz_inv3(i,k)
1231 dltr=dltr-cff*hz(i,j,k+1)
1232 dltl=dltl+cff*hz(i,j,k-1)
1233 br(i,k)=qc(i,k)+dltr
1234 bl(i,k)=qc(i,k)-dltl
1235 wr(i,k)=(2.0_r8*dltr-dltl)**2
1236 wl(i,k)=(dltr-2.0_r8*dltl)**2
1237 END DO
1238 END DO
1239 cff=1.0e-14_r8
1240 DO k=2,n(ng)-2
1241 DO i=istr,iend
1242 dltl=max(cff,wl(i,k ))
1243 dltr=max(cff,wr(i,k+1))
1244 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1245 bl(i,k+1)=br(i,k)
1246 END DO
1247 END DO
1248 DO i=istr,iend
1249 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1250#if defined LINEAR_CONTINUATION
1251 bl(i,n(ng))=br(i,n(ng)-1)
1252 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1253#elif defined NEUMANN
1254 bl(i,n(ng))=br(i,n(ng)-1)
1255 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1256#else
1257 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1258 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1259 br(i,n(ng)-1)=qc(i,n(ng))
1260#endif
1261#if defined LINEAR_CONTINUATION
1262 br(i,1)=bl(i,2)
1263 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1264#elif defined NEUMANN
1265 br(i,1)=bl(i,2)
1266 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1267#else
1268 bl(i,2)=qc(i,1) ! bottom grid boxes are
1269 br(i,1)=qc(i,1) ! re-assumed to be
1270 bl(i,1)=qc(i,1) ! piecewise constant.
1271#endif
1272 END DO
1273!
1274! Apply monotonicity constraint again, since the reconciled interfacial
1275! values may cause a non-monotonic behavior of the parabolic segments
1276! inside the grid box.
1277!
1278 DO k=1,n(ng)
1279 DO i=istr,iend
1280 dltr=br(i,k)-qc(i,k)
1281 dltl=qc(i,k)-bl(i,k)
1282 cffr=2.0_r8*dltr
1283 cffl=2.0_r8*dltl
1284 IF ((dltr*dltl).lt.0.0_r8) THEN
1285 dltr=0.0_r8
1286 dltl=0.0_r8
1287 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1288 dltr=cffl
1289 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1290 dltl=cffr
1291 END IF
1292 br(i,k)=qc(i,k)+dltr
1293 bl(i,k)=qc(i,k)-dltl
1294 END DO
1295 END DO
1296!
1297! After this moment reconstruction is considered complete. The next
1298! stage is to compute vertical advective fluxes, FC. It is expected
1299! that sinking may occurs relatively fast, the algorithm is designed
1300! to be free of CFL criterion, which is achieved by allowing
1301! integration bounds for semi-Lagrangian advective flux to use as
1302! many grid boxes in upstream direction as necessary.
1303!
1304! In the two code segments below, WL is the z-coordinate of the
1305! departure point for grid box interface z_w with the same indices;
1306! FC is the finite volume flux; ksource(:,k) is index of vertical
1307! grid box which contains the departure point (restricted by N(ng)).
1308! During the search: also add in content of whole grid boxes
1309! participating in FC.
1310!
1311 cff=dtdays*abs(wbio(isink))
1312 DO k=1,n(ng)
1313 DO i=istr,iend
1314 fc(i,k-1)=0.0_r8
1315 wl(i,k)=z_w(i,j,k-1)+cff
1316 wr(i,k)=hz(i,j,k)*qc(i,k)
1317 ksource(i,k)=k
1318 END DO
1319 END DO
1320 DO k=1,n(ng)
1321 DO ks=k,n(ng)-1
1322 DO i=istr,iend
1323 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1324 ksource(i,k)=ks+1
1325 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1326 END IF
1327 END DO
1328 END DO
1329 END DO
1330!
1331! Finalize computation of flux: add fractional part.
1332!
1333 DO k=1,n(ng)
1334 DO i=istr,iend
1335 ks=ksource(i,k)
1336 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1337 fc(i,k-1)=fc(i,k-1)+ &
1338 & hz(i,j,ks)*cu* &
1339 & (bl(i,ks)+ &
1340 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1341 & (1.5_r8-cu)* &
1342 & (br(i,ks)+bl(i,ks)- &
1343 & 2.0_r8*qc(i,ks))))
1344 END DO
1345 END DO
1346 DO k=1,n(ng)
1347 DO i=istr,iend
1348!^ tl_Bio(i,k,ibio)=tl_qc(i,k)+ &
1349!^ & (tl_FC(i,k)-tl_FC(i,k-1))*Hz_inv(i,k)+ &
1350!^ & (FC(i,k)-FC(i,k-1))*tl_Hz_inv(i,k)
1351!^
1352 ad_qc(i,k)=ad_qc(i,k)+ad_bio(i,k,ibio)
1353 ad_fc(i,k)=ad_fc(i,k)+hz_inv(i,k)*ad_bio(i,k,ibio)
1354 ad_fc(i,k-1)=ad_fc(i,k-1)-hz_inv(i,k)*ad_bio(i,k,ibio)
1355 ad_hz_inv(i,k)=ad_hz_inv(i,k)+ &
1356 & (fc(i,k)-fc(i,k-1))*ad_bio(i,k,ibio)
1357 ad_bio(i,k,ibio)=0.0_r8
1358 END DO
1359 END DO
1360!
1361! Adjoint of final computation of flux: add fractional part.
1362!
1363 DO k=1,n(ng)
1364 DO i=istr,iend
1365 ks=ksource(i,k)
1366 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1367!^ tl_FC(i,k-1)=tl_FC(i,k-1)+ &
1368!^ & (tl_Hz(i,j,ks)*cu+Hz(i,j,ks)*tl_cu)* &
1369!^ & (bL(i,ks)+ &
1370!^ & cu*(0.5_r8*(bR(i,ks)-bL(i,ks))- &
1371!^ & (1.5_r8-cu)* &
1372!^ & (bR(i,ks)+bL(i,ks)- &
1373!^ & 2.0_r8*qc(i,ks))))+ &
1374!^ & Hz(i,j,ks)*cu* &
1375!^ & (tl_bL(i,ks)+ &
1376!^ & tl_cu*(0.5_r8*(bR(i,ks)-bL(i,ks))- &
1377!^ & (1.5_r8-cu)* &
1378!^ & (bR(i,ks)+bL(i,ks)- &
1379!^ & 2.0_r8*qc(i,ks)))+ &
1380!^ & cu*(0.5_r8*(tl_bR(i,ks)-tl_bL(i,ks))+ &
1381!^ & tl_cu* &
1382!^ & (bR(i,ks)+bL(i,ks)-2.0_r8*qc(i,ks))- &
1383!^ & (1.5_r8-cu)* &
1384!^ & (tl_bR(i,ks)+tl_bL(i,ks)- &
1385!^ & 2.0_r8*tl_qc(i,ks))))
1386!^
1387 adfac=(bl(i,ks)+ &
1388 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1389 & (1.5_r8-cu)* &
1390 & (br(i,ks)+bl(i,ks)- &
1391 & 2.0_r8*qc(i,ks))))*ad_fc(i,k-1)
1392 adfac1=hz(i,j,ks)*cu*ad_fc(i,k-1)
1393 adfac2=adfac1*cu
1394 adfac3=adfac2*(1.5_r8-cu)
1395 ad_hz(i,j,ks)=ad_hz(i,j,ks)+cu*adfac
1396 ad_cu=ad_cu+hz(i,j,ks)*adfac
1397 ad_bl(i,ks)=ad_bl(i,ks)+adfac1
1398 ad_cu=ad_cu+ &
1399 & adfac1*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1400 & (1.5_r8-cu)* &
1401 & (br(i,ks)+bl(i,ks)- &
1402 & 2.0_r8*qc(i,ks)))+ &
1403 & adfac2*(br(i,ks)+bl(i,ks)-2.0_r8*qc(i,ks))
1404 ad_br(i,ks)=ad_br(i,ks)+0.5_r8*adfac2-adfac3
1405 ad_bl(i,ks)=ad_bl(i,ks)-0.5_r8*adfac2-adfac3
1406 ad_qc(i,ks)=ad_qc(i,ks)+2.0_r8*adfac3
1407!^ tl_cu=(0.5_r8+SIGN(0.5_r8, &
1408!^ & (1.0_r8-(WL(i,k)-z_w(i,j,ks-1))* &
1409!^ & Hz_inv(i,ks))))* &
1410!^ & ((tl_WL(i,k)-tl_z_w(i,j,ks-1))*Hz_inv(i,ks)+ &
1411!^ & (WL(i,k)-z_w(i,j,ks-1))*tl_Hz_inv(i,ks))
1412!^
1413 adfac=(0.5_r8+sign(0.5_r8, &
1414 & (1.0_r8-(wl(i,k)-z_w(i,j,ks-1))* &
1415 & hz_inv(i,ks))))*ad_cu
1416 adfac1=adfac*hz_inv(i,ks)
1417 ad_wl(i,k)=ad_wl(i,k)+adfac1
1418 ad_z_w(i,j,ks-1)=ad_z_w(i,j,ks-1)-adfac1
1419 ad_hz_inv(i,ks)=ad_hz_inv(i,ks)+ &
1420 & (wl(i,k)-z_w(i,j,ks-1))*adfac
1421 ad_cu=0.0_r8
1422 END DO
1423 END DO
1424!
1425! After this moment reconstruction is considered complete. The next
1426! stage is to compute vertical advective fluxes, FC. It is expected
1427! that sinking may occurs relatively fast, the algorithm is designed
1428! to be free of CFL criterion, which is achieved by allowing
1429! integration bounds for semi-Lagrangian advective flux to use as
1430! many grid boxes in upstream direction as necessary.
1431!
1432! In the two code segments below, WL is the z-coordinate of the
1433! departure point for grid box interface z_w with the same indices;
1434! FC is the finite volume flux; ksource(:,k) is index of vertical
1435! grid box which contains the departure point (restricted by N(ng)).
1436! During the search: also add in content of whole grid boxes
1437! participating in FC.
1438!
1439 cff=dtdays*abs(wbio(isink))
1440 DO k=1,n(ng)
1441 DO ks=k,n(ng)-1
1442 DO i=istr,iend
1443 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1444!^ tl_FC(i,k-1)=tl_FC(i,k-1)+tl_WR(i,ks)
1445!^
1446 ad_wr(i,ks)=ad_wr(i,ks)+ad_fc(i,k-1)
1447 END IF
1448 END DO
1449 END DO
1450 END DO
1451 DO k=1,n(ng)
1452 DO i=istr,iend
1453!^ tl_WR(i,k)=tl_Hz(i,j,k)*qc(i,k)+Hz(i,j,k)*tl_qc(i,k)
1454!^
1455 ad_hz(i,j,k)=ad_hz(i,j,k)+qc(i,k)*ad_wr(i,k)
1456 ad_qc(i,k)=ad_qc(i,k)+hz(i,j,k)*ad_wr(i,k)
1457 ad_wr(i,k)=0.0_r8
1458!^ tl_WL(i,k)=tl_z_w(i,j,k-1)+tl_cff
1459!^
1460 ad_z_w(i,j,k-1)=ad_z_w(i,j,k-1)+ad_wl(i,k)
1461 ad_cff=ad_cff+ad_wl(i,k)
1462 ad_wl(i,k)=0.0_r8
1463!^ tl_FC(i,k-1)=0.0_r8
1464!^
1465 ad_fc(i,k-1)=0.0_r8
1466 END DO
1467 END DO
1468!^ tl_cff=dtdays*SIGN(1.0_r8,Wbio(isink))*tl_Wbio(isink)
1469!^
1470 ad_wbio(isink)=ad_wbio(isink)+ &
1471 & dtdays*sign(1.0_r8,wbio(isink))*ad_cff
1472 ad_cff=0.0_r8
1473!
1474! Compute appropriate values of bR and bL.
1475!
1476! Copy concentration of biological particulates into scratch array
1477! "qc" (q-central, restrict it to be positive) which is hereafter
1478! interpreted as a set of grid-box averaged values for biogeochemical
1479! constituent concentration.
1480!
1481 DO k=1,n(ng)
1482 DO i=istr,iend
1483 qc(i,k)=bio(i,k,ibio)
1484 END DO
1485 END DO
1486!
1487 DO k=n(ng)-1,1,-1
1488 DO i=istr,iend
1489 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1490 END DO
1491 END DO
1492 DO k=2,n(ng)-1
1493 DO i=istr,iend
1494 dltr=hz(i,j,k)*fc(i,k)
1495 dltl=hz(i,j,k)*fc(i,k-1)
1496 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1497 cffr=cff*fc(i,k)
1498 cffl=cff*fc(i,k-1)
1499!
1500! Apply PPM monotonicity constraint to prevent oscillations within the
1501! grid box.
1502!
1503 IF ((dltr*dltl).le.0.0_r8) THEN
1504 dltr=0.0_r8
1505 dltl=0.0_r8
1506 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1507 dltr=cffl
1508 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1509 dltl=cffr
1510 END IF
1511!
1512! Compute right and left side values (bR,bL) of parabolic segments
1513! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1514!
1515! NOTE: Although each parabolic segment is monotonic within its grid
1516! box, monotonicity of the whole profile is not guaranteed,
1517! because bL(k+1)-bR(k) may still have different sign than
1518! qc(i,k+1)-qc(i,k). This possibility is excluded,
1519! after bL and bR are reconciled using WENO procedure.
1520!
1521 cff=(dltr-dltl)*hz_inv3(i,k)
1522 dltr=dltr-cff*hz(i,j,k+1)
1523 dltl=dltl+cff*hz(i,j,k-1)
1524 br(i,k)=qc(i,k)+dltr
1525 bl(i,k)=qc(i,k)-dltl
1526 wr(i,k)=(2.0_r8*dltr-dltl)**2
1527 wl(i,k)=(dltr-2.0_r8*dltl)**2
1528 END DO
1529 END DO
1530 cff=1.0e-14_r8
1531 DO k=2,n(ng)-2
1532 DO i=istr,iend
1533 dltl=max(cff,wl(i,k ))
1534 dltr=max(cff,wr(i,k+1))
1535 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1536 bl(i,k+1)=br(i,k)
1537 END DO
1538 END DO
1539 DO i=istr,iend
1540 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1541#if defined LINEAR_CONTINUATION
1542 bl(i,n(ng))=br(i,n(ng)-1)
1543 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1544#elif defined NEUMANN
1545 bl(i,n(ng))=br(i,n(ng)-1)
1546 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1547#else
1548 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1549 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1550 br(i,n(ng)-1)=qc(i,n(ng))
1551#endif
1552#if defined LINEAR_CONTINUATION
1553 br(i,1)=bl(i,2)
1554 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1555#elif defined NEUMANN
1556 br(i,1)=bl(i,2)
1557 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1558#else
1559 bl(i,2)=qc(i,1) ! bottom grid boxes are
1560 br(i,1)=qc(i,1) ! re-assumed to be
1561 bl(i,1)=qc(i,1) ! piecewise constant.
1562#endif
1563 END DO
1564!
1565! Apply monotonicity constraint again, since the reconciled interfacial
1566! values may cause a non-monotonic behavior of the parabolic segments
1567! inside the grid box.
1568!
1569 DO k=1,n(ng)
1570 DO i=istr,iend
1571 dltr=br(i,k)-qc(i,k)
1572 dltl=qc(i,k)-bl(i,k)
1573 cffr=2.0_r8*dltr
1574 cffl=2.0_r8*dltl
1575!^ tl_bL(i,k)=tl_qc(i,k)-tl_dltL
1576!^
1577 ad_qc(i,k)=ad_qc(i,k)+ad_bl(i,k)
1578 ad_dltl=ad_dltl-ad_bl(i,k)
1579 ad_bl(i,k)=0.0_r8
1580!^ tl_bR(i,k)=tl_qc(i,k)+tl_dltR
1581!^
1582 ad_qc(i,k)=ad_qc(i,k)+ad_br(i,k)
1583 ad_dltr=ad_dltr+ad_br(i,k)
1584 ad_br(i,k)=0.0_r8
1585 IF ((dltr*dltl).lt.0.0_r8) THEN
1586!^ tl_dltR=0.0_r8
1587!^
1588 ad_dltr=0.0_r8
1589!^ tl_dltL=0.0_r8
1590!^
1591 ad_dltl=0.0_r8
1592 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1593!^ tl_dltR=tl_cffL
1594!^
1595 ad_cffl=ad_cffl+ad_dltr
1596 ad_dltr=0.0_r8
1597 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1598!^ tl_dltL=tl_cffR
1599!^
1600 ad_cffr=ad_cffr+ad_dltl
1601 ad_dltl=0.0_r8
1602 END IF
1603!^ tl_cffL=2.0_r8*tl_dltL
1604!^
1605 ad_dltl=ad_dltl+2.0_r8*ad_cffl
1606 ad_cffl=0.0_r8
1607!^ tl_cffR=2.0_r8*tl_dltR
1608!^
1609 ad_dltr=ad_dltr+2.0_r8*ad_cffr
1610 ad_cffr=0.0_r8
1611!^ tl_dltL=tl_qc(i,k)-tl_bL(i,k)
1612!^
1613 ad_qc(i,k)=ad_qc(i,k)+ad_dltl
1614 ad_bl(i,k)=ad_bl(i,k)-ad_dltl
1615 ad_dltl=0.0_r8
1616!^ tl_dltR=tl_bR(i,k)-tl_qc(i,k)
1617!^
1618 ad_br(i,k)=ad_br(i,k)+ad_dltr
1619 ad_qc(i,k)=ad_qc(i,k)-ad_dltr
1620 ad_dltr=0.0_r8
1621 END DO
1622 END DO
1623 DO i=istr,iend
1624#if defined LINEAR_CONTINUATION
1625!^ tl_bR(i,1)=tl_bL(i,2)
1626!^
1627 ad_bl(i,2)=ad_bl(i,2)+ad_br(i,1)
1628 ad_br(i,1)=0.0_r8
1629!^ tl_bL(i,1)=2.0_r8*tl_qc(i,1)-tl_bR(i,1)
1630!^
1631 ad_qc(i,1)=ad_qc(i,1)+2.0_r8*ad_bl(i,1)
1632 ad_br(i,1)=ad_br(i,1)-ad_bl(i,1)
1633 ad_bl(i,1)=0.0_r8
1634#elif defined NEUMANN
1635!^ tl_bR(i,1)=tl_bL(i,2)
1636!^
1637 ad_bl(i,2)=ad_bl(i,2)+ad_br(i,1)
1638 ad_br(i,1)=0.0_r8
1639!^ tl_bL(i,1)=1.5_r8*tl_qc(i,1)-0.5_r8*tl_bR(i,1)
1640!^
1641 ad_qc(i,1)=ad_qc(i,1)+1.5_r8*ad_bl(i,1)
1642 ad_br(i,1)=ad_br(i,1)-0.5_r8*ad_bl(i,1)
1643 ad_bl(i,1)=0.0_r8
1644#else
1645!^ tl_bL(i,2)=tl_qc(i,1) ! bottom grid boxes are
1646!^ tl_bR(i,1)=tl_qc(i,1) ! re-assumed to be
1647!^ tl_bL(i,1)=tl_qc(i,1) ! piecewise constant.
1648!^
1649 ad_qc(i,1)=ad_qc(i,1)+ad_bl(i,1)+ &
1650 & ad_br(i,1)+ &
1651 & ad_bl(i,2)
1652 ad_bl(i,1)=0.0_r8
1653 ad_br(i,1)=0.0_r8
1654 ad_bl(i,2)=0.0_r8
1655#endif
1656#if defined LINEAR_CONTINUATION
1657!^ tl_bL(i,N(ng))=tl_bR(i,N(ng)-1)
1658!^
1659 ad_br(i,n(ng)-1)=ad_br(i,n(ng)-1)+ad_bl(i,n(ng))
1660 ad_bl(i,n(ng))=0.0_r8
1661!^ tl_bR(i,N(ng))=2.0_r8*tl_qc(i,N(ng))-tl_bL(i,N(ng))
1662!^
1663 ad_qc(i,n(ng))=ad_qc(i,n(ng))+2.0_r8*ad_br(i,n(ng))
1664 ad_bl(i,n(ng))=ad_bl(i,n(ng))-ad_br(i,n(ng))
1665 ad_br(i,n(ng))=0.0_r8
1666#elif defined NEUMANN
1667!^ tl_bL(i,N(ng))=tl_bR(i,N(ng)-1)
1668!^
1669 ad_br(i,n(ng)-1)=ad_br(i,n(ng)-1)+ad_bl(i,n(ng))
1670 ad_bl(i,n(ng))=0.0_r8
1671!^ tl_bR(i,N(ng))=1.5_r8*tl_qc(i,N(ng))-0.5_r8*tl_bL(i,N(ng))
1672!^
1673 ad_qc(i,n(ng))=ad_qc(i,n(ng))+1.5_r8*ad_br(i,n(ng))
1674 ad_bl(i,n(ng))=ad_bl(i,n(ng))-0.5_r8*ad_br(i,n(ng))
1675 ad_br(i,n(ng))=0.0_r8
1676#else
1677!^ tl_bR(i,N(ng))=tl_qc(i,N(ng)) ! default strictly monotonic
1678!^ tl_bL(i,N(ng))=tl_qc(i,N(ng)) ! conditions
1679!^ tl_bR(i,N(ng)-1)=tl_qc(i,N(ng))
1680!^
1681 ad_qc(i,n(ng))=ad_qc(i,n(ng))+ad_br(i,n(ng)-1)+ &
1682 & ad_bl(i,n(ng))+ &
1683 & ad_br(i,n(ng))
1684 ad_br(i,n(ng)-1)=0.0_r8
1685 ad_bl(i,n(ng))=0.0_r8
1686 ad_br(i,n(ng))=0.0_r8
1687#endif
1688 END DO
1689!
1690! Compute WR and WL arrays appropriate for this part of the code.
1691!
1692! Copy concentration of biological particulates into scratch array
1693! "qc" (q-central, restrict it to be positive) which is hereafter
1694! interpreted as a set of grid-box averaged values for biogeochemical
1695! constituent concentration.
1696!
1697 DO k=1,n(ng)
1698 DO i=istr,iend
1699 qc(i,k)=bio(i,k,ibio)
1700 END DO
1701 END DO
1702!
1703 DO k=n(ng)-1,1,-1
1704 DO i=istr,iend
1705 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1706 END DO
1707 END DO
1708 DO k=2,n(ng)-1
1709 DO i=istr,iend
1710 dltr=hz(i,j,k)*fc(i,k)
1711 dltl=hz(i,j,k)*fc(i,k-1)
1712 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1713 cffr=cff*fc(i,k)
1714 cffl=cff*fc(i,k-1)
1715!
1716! Apply PPM monotonicity constraint to prevent oscillations within the
1717! grid box.
1718!
1719 IF ((dltr*dltl).le.0.0_r8) THEN
1720 dltr=0.0_r8
1721 dltl=0.0_r8
1722 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1723 dltr=cffl
1724 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1725 dltl=cffr
1726 END IF
1727!
1728! Compute right and left side values (bR,bL) of parabolic segments
1729! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1730!
1731! NOTE: Although each parabolic segment is monotonic within its grid
1732! box, monotonicity of the whole profile is not guaranteed,
1733! because bL(k+1)-bR(k) may still have different sign than
1734! qc(i,k+1)-qc(i,k). This possibility is excluded,
1735! after bL and bR are reconciled using WENO procedure.
1736!
1737 cff=(dltr-dltl)*hz_inv3(i,k)
1738 dltr=dltr-cff*hz(i,j,k+1)
1739 dltl=dltl+cff*hz(i,j,k-1)
1740 br(i,k)=qc(i,k)+dltr
1741 bl(i,k)=qc(i,k)-dltl
1742 wr(i,k)=(2.0_r8*dltr-dltl)**2
1743 wl(i,k)=(dltr-2.0_r8*dltl)**2
1744 END DO
1745 END DO
1746
1747 cff=1.0e-14_r8
1748 DO k=2,n(ng)-2
1749 DO i=istr,iend
1750 dltl=max(cff,wl(i,k ))
1751 dltr=max(cff,wr(i,k+1))
1752 br1(i,k)=br(i,k)
1753 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1754 bl1(i,k+1)=bl(i,k+1)
1755 bl(i,k+1)=br(i,k)
1756!^ tl_bL(i,k+1)=tl_bR(i,k)
1757!^
1758 ad_br(i,k)=ad_br(i,k)+ad_bl(i,k+1)
1759 ad_bl(i,k+1)=0.0_r8
1760!^ tl_bR(i,k)=(tl_dltR*bR1(i,k)+dltR*tl_bR(i,k)+ &
1761!^ & tl_dltL*bL1(i,k+1)+dltL*tl_bL(i,k+1))/ &
1762!^ & (dltR+dltL)- &
1763!^ & (tl_dltR+tl_dltL)*bR(i,k)/(dltR+dltL)
1764!^
1765 adfac=ad_br(i,k)/(dltr+dltl)
1766 adfac1=ad_br(i,k)*br(i,k)/(dltr+dltl)
1767 ad_dltr=ad_dltr+adfac*br1(i,k)
1768 ad_dltl=ad_dltl+adfac*bl1(i,k+1)
1769 ad_bl(i,k+1)=ad_bl(i,k+1)+dltl*adfac
1770 ad_dltr=ad_dltr-adfac1
1771 ad_dltl=ad_dltl-adfac1
1772 ad_br(i,k)=dltr*adfac
1773!^ tl_dltR=(0.5_r8-SIGN(0.5_r8,cff-WR(i,k+1)))* &
1774!^ & tl_WR(i,k+1)
1775!^
1776 ad_wr(i,k+1)=ad_wr(i,k+1)+ &
1777 & (0.5_r8-sign(0.5_r8,cff-wr(i,k+1)))* &
1778 & ad_dltr
1779 ad_dltr=0.0_r8
1780!^ tl_dltL=(0.5_r8-SIGN(0.5_r8,cff-WL(i,k )))* &
1781!^ & tl_WL(i,k )
1782!^
1783 ad_wl(i,k )=ad_wl(i,k )+ &
1784 & (0.5_r8-sign(0.5_r8,cff-wl(i,k )))* &
1785 & ad_dltl
1786 ad_dltl=0.0_r8
1787 END DO
1788 END DO
1789
1790 DO k=2,n(ng)-1
1791 DO i=istr,iend
1792!
1793! Compute appropriate dltL and dltr.
1794!
1795 dltr=hz(i,j,k)*fc(i,k)
1796 dltl=hz(i,j,k)*fc(i,k-1)
1797 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1798 cffr=cff*fc(i,k)
1799 cffl=cff*fc(i,k-1)
1800!
1801! Apply PPM monotonicity constraint to prevent oscillations within the
1802! grid box.
1803!
1804 IF ((dltr*dltl).le.0.0_r8) THEN
1805 dltr=0.0_r8
1806 dltl=0.0_r8
1807 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1808 dltr=cffl
1809 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1810 dltl=cffr
1811 END IF
1812!
1813! Compute right and left side values (bR,bL) of parabolic segments
1814! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1815!
1816! NOTE: Although each parabolic segment is monotonic within its grid
1817! box, monotonicity of the whole profile is not guaranteed,
1818! because bL(k+1)-bR(k) may still have different sign than
1819! qc(i,k+1)-qc(i,k). This possibility is excluded,
1820! after bL and bR are reconciled using WENO procedure.
1821!
1822 cff=(dltr-dltl)*hz_inv3(i,k)
1823 dltr=dltr-cff*hz(i,j,k+1)
1824 dltl=dltl+cff*hz(i,j,k-1)
1825!^ tl_WL(i,k)=2.0_r8*(dltR-2.0_r8*dltL)* &
1826!^ & (tl_dltR-2.0_r8*tl_dltL)
1827!^
1828 adfac=ad_wl(i,k)*2.0_r8*(dltr-2.0_r8*dltl)
1829 ad_dltr=ad_dltr+adfac
1830 ad_dltl=ad_dltl-2.0_r8*adfac
1831 ad_wl(i,k)=0.0_r8
1832
1833!^ tl_WR(i,k)=2.0_r8*(2.0_r8*dltR-dltL)* &
1834!^ & (2.0_r8*tl_dltR-tl_dltL)
1835!^
1836 adfac=ad_wr(i,k)*2.0_r8*(2.0_r8*dltr-dltl)
1837 ad_dltr=ad_dltr+2.0_r8*adfac
1838 ad_dltl=ad_dltl-adfac
1839 ad_wr(i,k)=0.0_r8
1840!^ tl_bL(i,k)=tl_qc(i,k)-tl_dltL
1841!^
1842 ad_qc(i,k)=ad_qc(i,k)+ad_bl(i,k)
1843 ad_dltl=ad_dltl-ad_bl(i,k)
1844 ad_bl(i,k)=0.0_r8
1845!^ tl_bR(i,k)=tl_qc(i,k)+tl_dltR
1846!^
1847 ad_qc(i,k)=ad_qc(i,k)+ad_br(i,k)
1848 ad_dltr=ad_dltr+ad_br(i,k)
1849 ad_br(i,k)=0.0_r8
1850
1851!
1852! Compute appropriate dltL and dltr.
1853!
1854 dltr=hz(i,j,k)*fc(i,k)
1855 dltl=hz(i,j,k)*fc(i,k-1)
1856 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1857 cffr=cff*fc(i,k)
1858 cffl=cff*fc(i,k-1)
1859!
1860! Apply PPM monotonicity constraint to prevent oscillations within the
1861! grid box.
1862!
1863 IF ((dltr*dltl).le.0.0_r8) THEN
1864 dltr=0.0_r8
1865 dltl=0.0_r8
1866 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1867 dltr=cffl
1868 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1869 dltl=cffr
1870 END IF
1871
1872 cff=(dltr-dltl)*hz_inv3(i,k)
1873!^ tl_dltL=tl_dltL+tl_cff*Hz(i,j,k-1)+cff*tl_Hz(i,j,k-1)
1874!^
1875 ad_cff=ad_cff+ad_dltl*hz(i,j,k-1)
1876 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)+cff*ad_dltl
1877!^ tl_dltR=tl_dltR-tl_cff*Hz(i,j,k+1)-cff*tl_Hz(i,j,k+1)
1878!^
1879 ad_cff=ad_cff-ad_dltr*hz(i,j,k+1)
1880 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)-cff*ad_dltr
1881!^ tl_cff=(tl_dltR-tl_dltL)*Hz_inv3(i,k)+ &
1882!^ & (dltR-dltL)*tl_Hz_inv3(i,k)
1883!^
1884 adfac=ad_cff*hz_inv3(i,k)
1885 ad_dltr=ad_dltr+adfac
1886 ad_dltl=ad_dltl-adfac
1887 ad_hz_inv3(i,k)=ad_hz_inv3(i,k)+(dltr-dltl)*ad_cff
1888 ad_cff=0.0_r8
1889!
1890! Compute appropriate dltL and dltr.
1891!
1892 dltr=hz(i,j,k)*fc(i,k)
1893 dltl=hz(i,j,k)*fc(i,k-1)
1894 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1895 cffr=cff*fc(i,k)
1896 cffl=cff*fc(i,k-1)
1897
1898 IF ((dltr*dltl).le.0.0_r8) THEN
1899!^ tl_dltR=0.0_r8
1900!^
1901 ad_dltr=0.0_r8
1902!^ tl_dltL=0.0_r8
1903!^
1904 ad_dltl=0.0_r8
1905 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1906!^ tl_dltR=tl_cffL
1907!^
1908 ad_cffl=ad_cffl+ad_dltr
1909 ad_dltr=0.0_r8
1910 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1911!^ tl_dltL=tl_cffR
1912!^
1913 ad_cffr=ad_cffr+ad_dltl
1914 ad_dltl=0.0_r8
1915 END IF
1916!^ tl_cffL=tl_cff*FC(i,k-1)+cff*tl_FC(i,k-1)
1917!^
1918 ad_cff=ad_cff+ad_cffl*fc(i,k-1)
1919 ad_fc(i,k-1)=ad_fc(i,k-1)+cff*ad_cffl
1920 ad_cffl=0.0_r8
1921!^ tl_cffR=tl_cff*FC(i,k)+cff*tl_FC(i,k)
1922!^
1923 ad_cff=ad_cff+ad_cffr*fc(i,k)
1924 ad_fc(i,k)=ad_fc(i,k)+cff*ad_cffr
1925 ad_cffr=0.0_r8
1926!^ tl_cff=tl_Hz(i,j,k-1)+2.0_r8*tl_Hz(i,j,k)+tl_Hz(i,j,k+1)
1927!^
1928 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)+ad_cff
1929 ad_hz(i,j,k)=ad_hz(i,j,k)+2.0_r8*ad_cff
1930 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)+ad_cff
1931 ad_cff=0.0_r8
1932!^ tl_dltL=tl_Hz(i,j,k)*FC(i,k-1)+Hz(i,j,k)*tl_FC(i,k-1)
1933!^
1934 ad_hz(i,j,k)=ad_hz(i,j,k)+ad_dltl*fc(i,k-1)
1935 ad_fc(i,k-1)=ad_fc(i,k-1)+ad_dltl*hz(i,j,k)
1936 ad_dltl=0.0_r8
1937!^ tl_dltR=tl_Hz(i,j,k)*FC(i,k)+Hz(i,j,k)*tl_FC(i,k)
1938!^
1939 ad_hz(i,j,k)=ad_hz(i,j,k)+ad_dltr*fc(i,k)
1940 ad_fc(i,k)=ad_fc(i,k)+ad_dltr*hz(i,j,k)
1941 ad_dltr=0.0_r8
1942 END DO
1943 END DO
1944 DO k=n(ng)-1,1,-1
1945 DO i=istr,iend
1946!^ tl_FC(i,k)=(tl_qc(i,k+1)-tl_qc(i,k))*Hz_inv2(i,k)+ &
1947!^ & (qc(i,k+1)-qc(i,k))*tl_Hz_inv2(i,k)
1948!^
1949 adfac=ad_fc(i,k)*hz_inv2(i,k)
1950 ad_qc(i,k+1)=ad_qc(i,k+1)+adfac
1951 ad_qc(i,k)=ad_qc(i,k)-adfac
1952 ad_hz_inv2(i,k)=ad_hz_inv2(i,k)+(qc(i,k+1)-qc(i,k))* &
1953 & ad_fc(i,k)
1954 ad_fc(i,k)=0.0_r8
1955 END DO
1956 END DO
1957 DO k=1,n(ng)
1958 DO i=istr,iend
1959!^ tl_qc(i,k)=tl_Bio(i,k,ibio)
1960!^
1961 ad_bio(i,k,ibio)=ad_bio(i,k,ibio)+ad_qc(i,k)
1962 ad_qc(i,k)=0.0_r8
1963 END DO
1964 END DO
1965
1966 END DO sink_loop1
1967
1968!
1969! Compute appropriate basic state arrays II.
1970!
1971 DO k=1,n(ng)
1972 DO i=istr,iend
1973!
1974! At input, all tracers (index nnew) from predictor step have
1975! transport units (m Tunits) since we do not have yet the new
1976! values for zeta and Hz. These are known after the 2D barotropic
1977! time-stepping.
1978!
1979! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
1980! tracer times depth. However the basic state (nstp and nnew
1981! indices) that is read from the forward file is in units of
1982! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
1983! use t instead of t*Hz_inv.
1984!
1985 DO itrc=1,nbt
1986 ibio=idbio(itrc)
1987!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1988!^
1989 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
1990!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
1991!^
1992 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
1993 END DO
1994!
1995! Impose positive definite concentrations.
1996!
1997 cff2=0.0_r8
1998 DO itime=1,2
1999 cff1=0.0_r8
2000 itrcmax=idbio(1)
2001 DO itrc=1,nbt
2002 ibio=idbio(itrc)
2003 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
2004 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
2005 itrcmax=ibio
2006 END IF
2007 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2008 END DO
2009 IF (biotrc(itrcmax,itime).gt.cff1) THEN
2010 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
2011 END IF
2012 END DO
2013!
2014! Load biological tracers into local arrays.
2015!
2016 DO itrc=1,nbt
2017 ibio=idbio(itrc)
2018 bio_old(i,k,ibio)=biotrc(ibio,nstp)
2019 bio(i,k,ibio)=biotrc(ibio,nstp)
2020 END DO
2021 END DO
2022 END DO
2023!
2024! Calculate surface Photosynthetically Available Radiation (PAR). The
2025! net shortwave radiation is scaled back to Watts/m2 and multiplied by
2026! the fraction that is photosynthetically available, PARfrac.
2027!
2028 DO i=istr,iend
2029#ifdef CONST_PAR
2030!
2031! Specify constant surface irradiance a la Powell and Spitz.
2032!
2033 parsur(i)=158.075_r8
2034#else
2035 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
2036#endif
2037 END DO
2038!
2039!=======================================================================
2040! Start internal iterations to achieve convergence of the nonlinear
2041! backward-implicit solution.
2042!=======================================================================
2043!
2044 DO iteradj=1,iter
2045!
2046! Compute light attenuation as function of depth.
2047!
2048 DO i=istr,iend
2049 par=parsur(i)
2050 IF (parsur(i).gt.0.0_r8) THEN ! day time
2051 DO k=n(ng),1,-1
2052!
2053! Compute average light attenuation for each grid cell. Here, AttSW is
2054! the light attenuation due to seawater and AttPhy is the attenuation
2055! due to phytoplankton (self-shading coefficient).
2056!
2057 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
2058 & (z_w(i,j,k)-z_w(i,j,k-1))
2059 expatt=exp(-att)
2060 itop=par
2061 par=itop*(1.0_r8-expatt)/att ! average at cell center
2062 light(i,k)=par
2063!
2064! Light attenuation at the bottom of the grid cell. It is the starting
2065! PAR value for the next (deeper) vertical grid cell.
2066!
2067 par=itop*expatt
2068 END DO
2069 ELSE ! night time
2070 DO k=1,n(ng)
2071 light(i,k)=0.0_r8
2072 END DO
2073 END IF
2074 END DO
2075!
2076! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
2077! The Michaelis-Menten curve is used to describe the change in uptake
2078! rate as a function of nitrate concentration. Here, PhyIS is the
2079! initial slope of the P-I curve and K_NO3 is the half saturation of
2080! phytoplankton nitrate uptake.
2081!
2082 cff1=dtdays*vm_no3(ng)*phyis(ng)
2083 cff2=vm_no3(ng)*vm_no3(ng)
2084 cff3=phyis(ng)*phyis(ng)
2085 DO k=1,n(ng)
2086 DO i=istr,iend
2087 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
2088 cff=bio(i,k,iphyt)* &
2089 & cff1*cff4*light(i,k)/ &
2090 & (k_no3(ng)+bio(i,k,ino3_))
2091 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
2092 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
2093 & bio(i,k,ino3_)*cff
2094 END DO
2095 END DO
2096!
2097! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
2098! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
2099! pool as function of "sloppy feeding" and metabolic processes
2100! (ZooEEN and ZooEED fractions).
2101!
2102 cff1=dtdays*zoogr(ng)
2103 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
2104 DO k=1,n(ng)
2105 DO i=istr,iend
2106 cff=bio(i,k,izoop)* &
2107 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
2108 & bio(i,k,iphyt)
2109 bio1(i,k,iphyt)=bio(i,k,iphyt)
2110 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
2111 bio1(i,k,izoop)=bio(i,k,izoop)
2112 bio(i,k,izoop)=bio(i,k,izoop)+ &
2113 & bio(i,k,iphyt)*cff2*cff
2114 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2115 & bio(i,k,iphyt)*zooeen(ng)*cff
2116 bio(i,k,isdet)=bio(i,k,isdet)+ &
2117 & bio(i,k,iphyt)*zooeed(ng)*cff
2118 END DO
2119 END DO
2120!
2121 IF (iteradj.ne.iter) THEN
2122!
2123! Phytoplankton mortality to nutrients (PhyMRN rate) and detritus
2124! (PhyMRD rate).
2125!
2126 cff3=dtdays*phymrd(ng)
2127 cff2=dtdays*phymrn(ng)
2128 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2129 DO k=1,n(ng)
2130 DO i=istr,iend
2131 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
2132 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2133 & bio(i,k,iphyt)*cff2
2134 bio(i,k,isdet)=bio(i,k,isdet)+ &
2135 & bio(i,k,iphyt)*cff3
2136 END DO
2137 END DO
2138!
2139! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
2140! (ZooMRD rate).
2141!
2142 cff3=dtdays*zoomrd(ng)
2143 cff2=dtdays*zoomrn(ng)
2144 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2145 DO k=1,n(ng)
2146 DO i=istr,iend
2147 bio(i,k,izoop)=bio(i,k,izoop)*cff1
2148 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2149 & bio(i,k,izoop)*cff2
2150 bio(i,k,isdet)=bio(i,k,isdet)+ &
2151 & bio(i,k,izoop)*cff3
2152 END DO
2153 END DO
2154!
2155! Detritus breakdown to nutrients: remineralization (DetRR rate).
2156!
2157 cff2=dtdays*detrr(ng)
2158 cff1=1.0_r8/(1.0_r8+cff2)
2159 DO k=1,n(ng)
2160 DO i=istr,iend
2161 bio(i,k,isdet)=bio(i,k,isdet)*cff1
2162 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2163 & bio(i,k,isdet)*cff2
2164 END DO
2165 END DO
2166!
2167!-----------------------------------------------------------------------
2168! Vertical sinking terms: Phytoplankton and Detritus
2169!-----------------------------------------------------------------------
2170!
2171! Reconstruct vertical profile of selected biological constituents
2172! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
2173! grid box. Then, compute semi-Lagrangian flux due to sinking.
2174!
2175 DO isink=1,nsink
2176 ibio=idsink(isink)
2177!
2178! Copy concentration of biological particulates into scratch array
2179! "qc" (q-central, restrict it to be positive) which is hereafter
2180! interpreted as a set of grid-box averaged values for biogeochemical
2181! constituent concentration.
2182!
2183 DO k=1,n(ng)
2184 DO i=istr,iend
2185 qc(i,k)=bio(i,k,ibio)
2186 END DO
2187 END DO
2188!
2189 DO k=n(ng)-1,1,-1
2190 DO i=istr,iend
2191 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2192 END DO
2193 END DO
2194 DO k=2,n(ng)-1
2195 DO i=istr,iend
2196 dltr=hz(i,j,k)*fc(i,k)
2197 dltl=hz(i,j,k)*fc(i,k-1)
2198 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2199 cffr=cff*fc(i,k)
2200 cffl=cff*fc(i,k-1)
2201!
2202! Apply PPM monotonicity constraint to prevent oscillations within the
2203! grid box.
2204!
2205 IF ((dltr*dltl).le.0.0_r8) THEN
2206 dltr=0.0_r8
2207 dltl=0.0_r8
2208 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2209 dltr=cffl
2210 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2211 dltl=cffr
2212 END IF
2213!
2214! Compute right and left side values (bR,bL) of parabolic segments
2215! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
2216!
2217! NOTE: Although each parabolic segment is monotonic within its grid
2218! box, monotonicity of the whole profile is not guaranteed,
2219! because bL(k+1)-bR(k) may still have different sign than
2220! qc(i,k+1)-qc(i,k). This possibility is excluded,
2221! after bL and bR are reconciled using WENO procedure.
2222!
2223 cff=(dltr-dltl)*hz_inv3(i,k)
2224 dltr=dltr-cff*hz(i,j,k+1)
2225 dltl=dltl+cff*hz(i,j,k-1)
2226 br(i,k)=qc(i,k)+dltr
2227 bl(i,k)=qc(i,k)-dltl
2228 wr(i,k)=(2.0_r8*dltr-dltl)**2
2229 wl(i,k)=(dltr-2.0_r8*dltl)**2
2230 END DO
2231 END DO
2232 cff=1.0e-14_r8
2233 DO k=2,n(ng)-2
2234 DO i=istr,iend
2235 dltl=max(cff,wl(i,k ))
2236 dltr=max(cff,wr(i,k+1))
2237 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2238 bl(i,k+1)=br(i,k)
2239 END DO
2240 END DO
2241 DO i=istr,iend
2242 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
2243#if defined LINEAR_CONTINUATION
2244 bl(i,n(ng))=br(i,n(ng)-1)
2245 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
2246#elif defined NEUMANN
2247 bl(i,n(ng))=br(i,n(ng)-1)
2248 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
2249#else
2250 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
2251 bl(i,n(ng))=qc(i,n(ng)) ! conditions
2252 br(i,n(ng)-1)=qc(i,n(ng))
2253#endif
2254#if defined LINEAR_CONTINUATION
2255 br(i,1)=bl(i,2)
2256 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2257#elif defined NEUMANN
2258 br(i,1)=bl(i,2)
2259 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2260#else
2261 bl(i,2)=qc(i,1) ! bottom grid boxes are
2262 br(i,1)=qc(i,1) ! re-assumed to be
2263 bl(i,1)=qc(i,1) ! piecewise constant.
2264#endif
2265 END DO
2266!
2267! Apply monotonicity constraint again, since the reconciled interfacial
2268! values may cause a non-monotonic behavior of the parabolic segments
2269! inside the grid box.
2270!
2271 DO k=1,n(ng)
2272 DO i=istr,iend
2273 dltr=br(i,k)-qc(i,k)
2274 dltl=qc(i,k)-bl(i,k)
2275 cffr=2.0_r8*dltr
2276 cffl=2.0_r8*dltl
2277 IF ((dltr*dltl).lt.0.0_r8) THEN
2278 dltr=0.0_r8
2279 dltl=0.0_r8
2280 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2281 dltr=cffl
2282 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2283 dltl=cffr
2284 END IF
2285 br(i,k)=qc(i,k)+dltr
2286 bl(i,k)=qc(i,k)-dltl
2287 END DO
2288 END DO
2289!
2290! After this moment reconstruction is considered complete. The next
2291! stage is to compute vertical advective fluxes, FC. It is expected
2292! that sinking may occurs relatively fast, the algorithm is designed
2293! to be free of CFL criterion, which is achieved by allowing
2294! integration bounds for semi-Lagrangian advective flux to use as
2295! many grid boxes in upstream direction as necessary.
2296!
2297! In the two code segments below, WL is the z-coordinate of the
2298! departure point for grid box interface z_w with the same indices;
2299! FC is the finite volume flux; ksource(:,k) is index of vertical
2300! grid box which contains the departure point (restricted by N(ng)).
2301! During the search: also add in content of whole grid boxes
2302! participating in FC.
2303!
2304 cff=dtdays*abs(wbio(isink))
2305 DO k=1,n(ng)
2306 DO i=istr,iend
2307 fc(i,k-1)=0.0_r8
2308 wl(i,k)=z_w(i,j,k-1)+cff
2309 wr(i,k)=hz(i,j,k)*qc(i,k)
2310 ksource(i,k)=k
2311 END DO
2312 END DO
2313 DO k=1,n(ng)
2314 DO ks=k,n(ng)-1
2315 DO i=istr,iend
2316 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2317 ksource(i,k)=ks+1
2318 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2319 END IF
2320 END DO
2321 END DO
2322 END DO
2323!
2324! Finalize computation of flux: add fractional part.
2325!
2326 DO k=1,n(ng)
2327 DO i=istr,iend
2328 ks=ksource(i,k)
2329 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2330 fc(i,k-1)=fc(i,k-1)+ &
2331 & hz(i,j,ks)*cu* &
2332 & (bl(i,ks)+ &
2333 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2334 & (1.5_r8-cu)* &
2335 & (br(i,ks)+bl(i,ks)- &
2336 & 2.0_r8*qc(i,ks))))
2337 END DO
2338 END DO
2339 DO k=1,n(ng)
2340 DO i=istr,iend
2341 bio(i,k,ibio)=qc(i,k)+ &
2342 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2343 END DO
2344 END DO
2345 END DO
2346 END IF
2347 END DO
2348!
2349! End of compute basic state arrays II.
2350!
2351! Adjoint detritus breakdown to nutrients: remineralization
2352! (DetRR rate).
2353!
2354 cff2=dtdays*detrr(ng)
2355 cff1=1.0_r8/(1.0_r8+cff2)
2356 DO k=1,n(ng)
2357 DO i=istr,iend
2358!^ tl_Bio(i,k,iNO3_)=tl_Bio(i,k,iNO3_)+ &
2359!^ & tl_Bio(i,k,iSDet)*cff2
2360!^
2361 ad_bio(i,k,isdet)=ad_bio(i,k,isdet)+ &
2362 & cff2*ad_bio(i,k,ino3_)
2363!^ tl_Bio(i,k,iSDet)=tl_Bio(i,k,iSDet)*cff1
2364!^
2365 ad_bio(i,k,isdet)=ad_bio(i,k,isdet)*cff1
2366 END DO
2367 END DO
2368!
2369! Adjoint Zooplankton mortality to nutrients (ZooMRN rate) and
2370! Detritus (ZooMRD rate).
2371!
2372 cff3=dtdays*zoomrd(ng)
2373 cff2=dtdays*zoomrn(ng)
2374 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2375 DO k=1,n(ng)
2376 DO i=istr,iend
2377!^ tl_Bio(i,k,iSDet)=tl_Bio(i,k,iSDet)+ &
2378!^ & tl_Bio(i,k,iZoop)*cff3
2379!^
2380 ad_bio(i,k,izoop)=ad_bio(i,k,izoop)+ &
2381 & cff3*ad_bio(i,k,isdet)
2382!^ tl_Bio(i,k,iNO3_)=tl_Bio(i,k,iNO3_)+ &
2383!^ & tl_Bio(i,k,iZoop)*cff2
2384!^
2385 ad_bio(i,k,izoop)=ad_bio(i,k,izoop)+ &
2386 & cff2*ad_bio(i,k,ino3_)
2387!^ tl_Bio(i,k,iZoop)=tl_Bio(i,k,iZoop)*cff1
2388!^
2389 ad_bio(i,k,izoop)=ad_bio(i,k,izoop)*cff1
2390 END DO
2391 END DO
2392!
2393! Adjoint Phytoplankton mortality to nutrients (PhyMRN rate) and
2394! detritus (PhyMRD rate).
2395!
2396 cff3=dtdays*phymrd(ng)
2397 cff2=dtdays*phymrn(ng)
2398 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2399 DO k=1,n(ng)
2400 DO i=istr,iend
2401!^ tl_Bio(i,k,iSDet)=tl_Bio(i,k,iSDet)+ &
2402!^ & tl_Bio(i,k,iPhyt)*cff3
2403!^
2404 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2405 & cff3*ad_bio(i,k,isdet)
2406!^ tl_Bio(i,k,iNO3_)=tl_Bio(i,k,iNO3_)+ &
2407!^ & tl_Bio(i,k,iPhyt)*cff2
2408!^
2409 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2410 & cff2*ad_bio(i,k,ino3_)
2411!^ tl_Bio(i,k,iPhyt)=tl_Bio(i,k,iPhyt)*cff1
2412!^
2413 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)*cff1
2414 END DO
2415 END DO
2416!
2417! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
2418! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
2419! pool as function of "sloppy feeding" and metabolic processes
2420! (ZooEEN and ZooEED fractions).
2421!
2422 cff1=dtdays*zoogr(ng)
2423 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
2424 DO k=1,n(ng)
2425 DO i=istr,iend
2426 cff=bio1(i,k,izoop)* &
2427 & cff1*(1.0_r8-exp(-ivlev(ng)*bio1(i,k,iphyt)))/ &
2428 & bio1(i,k,iphyt)
2429!^ tl_Bio(i,k,iSDet)=tl_Bio(i,k,iSDet)+ &
2430!^ & ZooEED(ng)*(tl_Bio(i,k,iPhyt)*cff+ &
2431!^ & Bio(i,k,iPhyt)*tl_cff)
2432!^
2433 ad_cff=ad_cff+zooeed(ng)*bio(i,k,iphyt)*ad_bio(i,k,isdet)
2434 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2435 & zooeed(ng)*cff*ad_bio(i,k,isdet)
2436!^ tl_Bio(i,k,iNO3_)=tl_Bio(i,k,iNO3_)+ &
2437!^ & ZooEEN(ng)*(tl_Bio(i,k,iPhyt)*cff+ &
2438!^ & Bio(i,k,iPhyt)*tl_cff)
2439!^
2440 ad_cff=ad_cff+ &
2441 & zooeen(ng)*bio(i,k,iphyt)*ad_bio(i,k,ino3_)
2442 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2443 & zooeen(ng)*cff*ad_bio(i,k,ino3_)
2444!^ tl_Bio(i,k,iZoop)=tl_Bio(i,k,iZoop)+ &
2445!^ & cff2*(tl_Bio(i,k,iPhyt)*cff+ &
2446!^ & Bio(i,k,iPhyt)*tl_cff)
2447!^
2448 ad_cff=ad_cff+ &
2449 & cff2*bio(i,k,iphyt)*ad_bio(i,k,izoop)
2450 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2451 & cff2*cff*ad_bio(i,k,izoop)
2452!^ tl_Bio(i,k,iPhyt)=(tl_Bio(i,k,iPhyt)- &
2453!^ & tl_cff*Bio(i,k,iPhyt))/ &
2454!^ & (1.0_r8+cff)
2455!^
2456 adfac=ad_bio(i,k,iphyt)/(1.0_r8+cff)
2457 ad_cff=ad_cff-bio(i,k,iphyt)*adfac
2458 ad_bio(i,k,iphyt)=adfac
2459!^ tl_cff=(tl_Bio(i,k,iZoop)* &
2460!^ & cff1*(1.0_r8-EXP(-Ivlev(ng)*Bio1(i,k,iPhyt)))+ &
2461!^ & Bio1(i,k,iZoop)*Ivlev(ng)*tl_Bio(i,k,iPhyt)*cff1* &
2462!^ & EXP(-Ivlev(ng)*Bio1(i,k,iPhyt))- &
2463!^ & tl_Bio(i,k,iPhyt)*cff)/ &
2464!^ & Bio1(i,k,iPhyt)
2465!^
2466 fac=exp(-ivlev(ng)*bio1(i,k,iphyt))
2467 adfac=ad_cff/bio1(i,k,iphyt)
2468 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)-adfac*cff
2469 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2470 & bio1(i,k,izoop)*ivlev(ng)* &
2471 & adfac*cff1*fac
2472 ad_bio(i,k,izoop)=ad_bio(i,k,izoop)+ &
2473 & adfac*cff1*(1.0_r8-fac)
2474 ad_cff=0.0_r8
2475 END DO
2476 END DO
2477!
2478! Compute appropriate basic state arrays I.
2479!
2480 DO k=1,n(ng)
2481 DO i=istr,iend
2482!
2483! At input, all tracers (index nnew) from predictor step have
2484! transport units (m Tunits) since we do not have yet the new
2485! values for zeta and Hz. These are known after the 2D barotropic
2486! time-stepping.
2487!
2488! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
2489! tracer times depth. However the basic state (nstp and nnew
2490! indices) that is read from the forward file is in units of
2491! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
2492! use t instead of t*Hz_inv.
2493!
2494 DO itrc=1,nbt
2495 ibio=idbio(itrc)
2496!^ BioTrc(ibio,nstp)=t(i,j,k,nstp,ibio)
2497!^
2498 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
2499!^ BioTrc(ibio,nnew)=t(i,j,k,nnew,ibio)*Hz_inv(i,k)
2500!^
2501 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)
2502 END DO
2503!
2504! Impose positive definite concentrations.
2505!
2506 cff2=0.0_r8
2507 DO itime=1,2
2508 cff1=0.0_r8
2509 itrcmax=idbio(1)
2510 DO itrc=1,nbt
2511 ibio=idbio(itrc)
2512 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
2513 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
2514 itrcmax=ibio
2515 END IF
2516 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
2517 END DO
2518 IF (biotrc(itrcmax,itime).gt.cff1) THEN
2519 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
2520 END IF
2521 END DO
2522!
2523! Load biological tracers into local arrays.
2524!
2525 DO itrc=1,nbt
2526 ibio=idbio(itrc)
2527 bio_old(i,k,ibio)=biotrc(ibio,nstp)
2528 bio(i,k,ibio)=biotrc(ibio,nstp)
2529 END DO
2530 END DO
2531 END DO
2532!
2533! Calculate surface Photosynthetically Available Radiation (PAR). The
2534! net shortwave radiation is scaled back to Watts/m2 and multiplied by
2535! the fraction that is photosynthetically available, PARfrac.
2536!
2537 DO i=istr,iend
2538#ifdef CONST_PAR
2539!
2540! Specify constant surface irradiance a la Powell and Spitz.
2541!
2542 parsur(i)=158.075_r8
2543#else
2544 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
2545#endif
2546 END DO
2547!
2548!=======================================================================
2549! Start internal iterations to achieve convergence of the nonlinear
2550! backward-implicit solution.
2551!=======================================================================
2552!
2553 DO iteradj=1,iter
2554!
2555! Compute light attenuation as function of depth.
2556!
2557 DO i=istr,iend
2558 par=parsur(i)
2559 IF (parsur(i).gt.0.0_r8) THEN ! day time
2560 DO k=n(ng),1,-1
2561!
2562! Compute average light attenuation for each grid cell. Here, AttSW is
2563! the light attenuation due to seawater and AttPhy is the attenuation
2564! due to phytoplankton (self-shading coefficient).
2565!
2566 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
2567 & (z_w(i,j,k)-z_w(i,j,k-1))
2568 expatt=exp(-att)
2569 itop=par
2570 par=itop*(1.0_r8-expatt)/att ! average at cell center
2571 light(i,k)=par
2572!
2573! Light attenuation at the bottom of the grid cell. It is the starting
2574! PAR value for the next (deeper) vertical grid cell.
2575!
2576 par=itop*expatt
2577 END DO
2578 ELSE ! night time
2579 DO k=1,n(ng)
2580 light(i,k)=0.0_r8
2581 END DO
2582 END IF
2583 END DO
2584!
2585! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
2586! The Michaelis-Menten curve is used to describe the change in uptake
2587! rate as a function of nitrate concentration. Here, PhyIS is the
2588! initial slope of the P-I curve and K_NO3 is the half saturation of
2589! phytoplankton nitrate uptake.
2590!
2591 cff1=dtdays*vm_no3(ng)*phyis(ng)
2592 cff2=vm_no3(ng)*vm_no3(ng)
2593 cff3=phyis(ng)*phyis(ng)
2594 DO k=1,n(ng)
2595 DO i=istr,iend
2596 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
2597 cff=bio(i,k,iphyt)* &
2598 & cff1*cff4*light(i,k)/ &
2599 & (k_no3(ng)+bio(i,k,ino3_))
2600 bio1(i,k,ino3_)=bio(i,k,ino3_)
2601 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
2602 bio1(i,k,iphyt)=bio(i,k,iphyt)
2603 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
2604 & bio(i,k,ino3_)*cff
2605 END DO
2606 END DO
2607!
2608 IF (iteradj.ne.iter) THEN
2609!
2610! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
2611! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
2612! pool as function of "sloppy feeding" and metabolic processes
2613! (ZooEEN and ZooEED fractions).
2614!
2615 cff1=dtdays*zoogr(ng)
2616 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
2617 DO k=1,n(ng)
2618 DO i=istr,iend
2619 cff=bio(i,k,izoop)* &
2620 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
2621 & bio(i,k,iphyt)
2622 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
2623 bio(i,k,izoop)=bio(i,k,izoop)+ &
2624 & bio(i,k,iphyt)*cff2*cff
2625 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2626 & bio(i,k,iphyt)*zooeen(ng)*cff
2627 bio(i,k,isdet)=bio(i,k,isdet)+ &
2628 & bio(i,k,iphyt)*zooeed(ng)*cff
2629 END DO
2630 END DO
2631!
2632! Phytoplankton mortality to nutrients (PhyMRN rate) and detritus
2633! (PhyMRD rate).
2634!
2635 cff3=dtdays*phymrd(ng)
2636 cff2=dtdays*phymrn(ng)
2637 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2638 DO k=1,n(ng)
2639 DO i=istr,iend
2640 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
2641 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2642 & bio(i,k,iphyt)*cff2
2643 bio(i,k,isdet)=bio(i,k,isdet)+ &
2644 & bio(i,k,iphyt)*cff3
2645 END DO
2646 END DO
2647!
2648! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
2649! (ZooMRD rate).
2650!
2651 cff3=dtdays*zoomrd(ng)
2652 cff2=dtdays*zoomrn(ng)
2653 cff1=1.0_r8/(1.0_r8+cff2+cff3)
2654 DO k=1,n(ng)
2655 DO i=istr,iend
2656 bio(i,k,izoop)=bio(i,k,izoop)*cff1
2657 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2658 & bio(i,k,izoop)*cff2
2659 bio(i,k,isdet)=bio(i,k,isdet)+ &
2660 & bio(i,k,izoop)*cff3
2661 END DO
2662 END DO
2663!
2664! Detritus breakdown to nutrients: remineralization (DetRR rate).
2665!
2666 cff2=dtdays*detrr(ng)
2667 cff1=1.0_r8/(1.0_r8+cff2)
2668 DO k=1,n(ng)
2669 DO i=istr,iend
2670 bio(i,k,isdet)=bio(i,k,isdet)*cff1
2671 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
2672 & bio(i,k,isdet)*cff2
2673 END DO
2674 END DO
2675!
2676!-----------------------------------------------------------------------
2677! Vertical sinking terms: Phytoplankton and Detritus
2678!-----------------------------------------------------------------------
2679!
2680! Reconstruct vertical profile of selected biological constituents
2681! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
2682! grid box. Then, compute semi-Lagrangian flux due to sinking.
2683!
2684 DO isink=1,nsink
2685 ibio=idsink(isink)
2686!
2687! Copy concentration of biological particulates into scratch array
2688! "qc" (q-central, restrict it to be positive) which is hereafter
2689! interpreted as a set of grid-box averaged values for biogeochemical
2690! constituent concentration.
2691!
2692 DO k=1,n(ng)
2693 DO i=istr,iend
2694 qc(i,k)=bio(i,k,ibio)
2695 END DO
2696 END DO
2697!
2698 DO k=n(ng)-1,1,-1
2699 DO i=istr,iend
2700 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
2701 END DO
2702 END DO
2703 DO k=2,n(ng)-1
2704 DO i=istr,iend
2705 dltr=hz(i,j,k)*fc(i,k)
2706 dltl=hz(i,j,k)*fc(i,k-1)
2707 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
2708 cffr=cff*fc(i,k)
2709 cffl=cff*fc(i,k-1)
2710!
2711! Apply PPM monotonicity constraint to prevent oscillations within the
2712! grid box.
2713!
2714 IF ((dltr*dltl).le.0.0_r8) THEN
2715 dltr=0.0_r8
2716 dltl=0.0_r8
2717 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2718 dltr=cffl
2719 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2720 dltl=cffr
2721 END IF
2722!
2723! Compute right and left side values (bR,bL) of parabolic segments
2724! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
2725!
2726! NOTE: Although each parabolic segment is monotonic within its grid
2727! box, monotonicity of the whole profile is not guaranteed,
2728! because bL(k+1)-bR(k) may still have different sign than
2729! qc(i,k+1)-qc(i,k). This possibility is excluded,
2730! after bL and bR are reconciled using WENO procedure.
2731!
2732 cff=(dltr-dltl)*hz_inv3(i,k)
2733 dltr=dltr-cff*hz(i,j,k+1)
2734 dltl=dltl+cff*hz(i,j,k-1)
2735 br(i,k)=qc(i,k)+dltr
2736 bl(i,k)=qc(i,k)-dltl
2737 wr(i,k)=(2.0_r8*dltr-dltl)**2
2738 wl(i,k)=(dltr-2.0_r8*dltl)**2
2739 END DO
2740 END DO
2741 cff=1.0e-14_r8
2742 DO k=2,n(ng)-2
2743 DO i=istr,iend
2744 dltl=max(cff,wl(i,k ))
2745 dltr=max(cff,wr(i,k+1))
2746 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
2747 bl(i,k+1)=br(i,k)
2748 END DO
2749 END DO
2750 DO i=istr,iend
2751 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
2752#if defined LINEAR_CONTINUATION
2753 bl(i,n(ng))=br(i,n(ng)-1)
2754 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
2755#elif defined NEUMANN
2756 bl(i,n(ng))=br(i,n(ng)-1)
2757 br(i,n(ng))=1.5*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
2758#else
2759 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
2760 bl(i,n(ng))=qc(i,n(ng)) ! conditions
2761 br(i,n(ng)-1)=qc(i,n(ng))
2762#endif
2763#if defined LINEAR_CONTINUATION
2764 br(i,1)=bl(i,2)
2765 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
2766#elif defined NEUMANN
2767 br(i,1)=bl(i,2)
2768 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
2769#else
2770 bl(i,2)=qc(i,1) ! bottom grid boxes are
2771 br(i,1)=qc(i,1) ! re-assumed to be
2772 bl(i,1)=qc(i,1) ! piecewise constant.
2773#endif
2774 END DO
2775!
2776! Apply monotonicity constraint again, since the reconciled interfacial
2777! values may cause a non-monotonic behavior of the parabolic segments
2778! inside the grid box.
2779!
2780 DO k=1,n(ng)
2781 DO i=istr,iend
2782 dltr=br(i,k)-qc(i,k)
2783 dltl=qc(i,k)-bl(i,k)
2784 cffr=2.0_r8*dltr
2785 cffl=2.0_r8*dltl
2786 IF ((dltr*dltl).lt.0.0_r8) THEN
2787 dltr=0.0_r8
2788 dltl=0.0_r8
2789 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
2790 dltr=cffl
2791 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
2792 dltl=cffr
2793 END IF
2794 br(i,k)=qc(i,k)+dltr
2795 bl(i,k)=qc(i,k)-dltl
2796 END DO
2797 END DO
2798!
2799! After this moment reconstruction is considered complete. The next
2800! stage is to compute vertical advective fluxes, FC. It is expected
2801! that sinking may occurs relatively fast, the algorithm is designed
2802! to be free of CFL criterion, which is achieved by allowing
2803! integration bounds for semi-Lagrangian advective flux to use as
2804! many grid boxes in upstream direction as necessary.
2805!
2806! In the two code segments below, WL is the z-coordinate of the
2807! departure point for grid box interface z_w with the same indices;
2808! FC is the finite volume flux; ksource(:,k) is index of vertical
2809! grid box which contains the departure point (restricted by N(ng)).
2810! During the search: also add in content of whole grid boxes
2811! participating in FC.
2812!
2813 cff=dtdays*abs(wbio(isink))
2814 DO k=1,n(ng)
2815 DO i=istr,iend
2816 fc(i,k-1)=0.0_r8
2817 wl(i,k)=z_w(i,j,k-1)+cff
2818 wr(i,k)=hz(i,j,k)*qc(i,k)
2819 ksource(i,k)=k
2820 END DO
2821 END DO
2822 DO k=1,n(ng)
2823 DO ks=k,n(ng)-1
2824 DO i=istr,iend
2825 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
2826 ksource(i,k)=ks+1
2827 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
2828 END IF
2829 END DO
2830 END DO
2831 END DO
2832!
2833! Finalize computation of flux: add fractional part.
2834!
2835 DO k=1,n(ng)
2836 DO i=istr,iend
2837 ks=ksource(i,k)
2838 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
2839 fc(i,k-1)=fc(i,k-1)+ &
2840 & hz(i,j,ks)*cu* &
2841 & (bl(i,ks)+ &
2842 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
2843 & (1.5_r8-cu)* &
2844 & (br(i,ks)+bl(i,ks)- &
2845 & 2.0_r8*qc(i,ks))))
2846 END DO
2847 END DO
2848 DO k=1,n(ng)
2849 DO i=istr,iend
2850 bio(i,k,ibio)=qc(i,k)+ &
2851 & (fc(i,k)-fc(i,k-1))*hz_inv(i,k)
2852 END DO
2853 END DO
2854 END DO
2855 END IF
2856 END DO
2857!
2858! End of compute basic state arrays I.
2859!
2860! Adjoint Phytoplankton photosynthetic growth and nitrate uptake
2861! (Vm_NO3 rate). The Michaelis-Menten curve is used to describe the
2862! change in uptake rate as a function of nitrate concentration.
2863! Here, PhyIS is the initial slope of the P-I curve and K_NO3 is the
2864! half saturation of phytoplankton nitrate uptake.
2865!
2866 cff1=dtdays*vm_no3(ng)*phyis(ng)
2867 cff2=vm_no3(ng)*vm_no3(ng)
2868 cff3=phyis(ng)*phyis(ng)
2869 DO k=1,n(ng)
2870 DO i=istr,iend
2871 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
2872 cff=bio1(i,k,iphyt)* &
2873 & cff1*cff4*light(i,k)/ &
2874 & (k_no3(ng)+bio1(i,k,ino3_))
2875!^ tl_Bio(i,k,iPhyt)=tl_Bio(i,k,iPhyt)+ &
2876!^ & tl_Bio(i,k,iNO3_)*cff+ &
2877!^ & Bio(i,k,iNO3_)*tl_cff
2878!^
2879 ad_cff=ad_cff+bio(i,k,ino3_)*ad_bio(i,k,iphyt)
2880 ad_bio(i,k,ino3_)=ad_bio(i,k,ino3_)+ &
2881 & cff*ad_bio(i,k,iphyt)
2882!^ tl_Bio(i,k,iNO3_)=(tl_Bio(i,k,iNO3_)- &
2883!^ & tl_cff*Bio(i,k,iNO3_))/ &
2884!^ & (1.0_r8+cff)
2885!^
2886 adfac=ad_bio(i,k,ino3_)/(1.0_r8+cff)
2887 ad_cff=ad_cff-bio(i,k,ino3_)*adfac
2888 ad_bio(i,k,ino3_)=adfac
2889!^ tl_cff=(tl_Bio(i,k,iPhyt)*cff1*cff4*Light(i,k)+ &
2890!^ & Bio1(i,k,iPhyt)*cff1* &
2891!^ & (tl_cff4*Light(i,k)+cff4*tl_Light(i,k))- &
2892!^ & tl_Bio(i,k,iNO3_)*cff)/ &
2893!^ & (K_NO3(ng)+Bio1(i,k,iNO3_))
2894!^
2895 adfac=ad_cff/(k_no3(ng)+bio1(i,k,ino3_))
2896 adfac1=adfac*bio1(i,k,iphyt)*cff1
2897 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2898 & cff1*cff4*light(i,k)*adfac
2899 ad_cff4=adfac1*light(i,k)
2900 ad_light(i,k)=ad_light(i,k)+ &
2901 & adfac1*cff4
2902 ad_bio(i,k,ino3_)=ad_bio(i,k,ino3_)- &
2903 & adfac*cff
2904 ad_cff=0.0_r8
2905!^ tl_cff4=-cff3*tl_Light(i,k)*Light(i,k)*cff4*cff4*cff4
2906!^
2907 ad_light(i,k)=ad_light(i,k)- &
2908 & cff3*light(i,k)* &
2909 & cff4*cff4*cff4*ad_cff4
2910 ad_cff4=0.0_r8
2911 END DO
2912 END DO
2913!
2914! Compute adjoint light attenuation as function of depth.
2915!
2916 DO i=istr,iend
2917 par=parsur(i)
2918 IF (parsur(i).gt.0.0_r8) THEN ! day time
2919 DO k=1,n(ng)
2920!
2921! Compute the basic state PAR appropriate for each level.
2922!
2923 par=parsur(i)
2924 DO kk=n(ng),k,-1
2925!
2926! Compute average light attenuation for each grid cell. Here, AttSW is
2927! the light attenuation due to seawater and AttPhy is the attenuation
2928! due to phytoplankton (self-shading coefficient).
2929!
2930 att=(attsw(ng)+attphy(ng)*bio1(i,kk,iphyt))* &
2931 & (z_w(i,j,kk)-z_w(i,j,kk-1))
2932 expatt=exp(-att)
2933 itop=par
2934 par=itop*(1.0_r8-expatt)/att ! average at cell center
2935 par1=par
2936!
2937! Light attenuation at the bottom of the grid cell. It is the starting
2938! PAR value for the next (deeper) vertical grid cell.
2939!
2940 par=itop*expatt
2941 END DO
2942!
2943! Adjoint of light attenuation at the bottom of the grid cell. It is
2944! the starting PAR value for the next (deeper) vertical grid cell.
2945!
2946!^ tl_PAR=tl_Itop*ExpAtt+Itop*tl_ExpAtt
2947!^
2948 ad_expatt=ad_expatt+itop*ad_par
2949 ad_itop=ad_itop+expatt*ad_par
2950 ad_par=0.0_r8
2951!
2952! Adjoint of compute average light attenuation for each grid cell.
2953! Here, AttSW is the light attenuation due to seawater and AttPhy is
2954! the attenuation due to phytoplankton (self-shading coefficient).
2955!
2956!^ tl_Light(i,k)=tl_PAR
2957!^
2958 ad_par=ad_par+ad_light(i,k)
2959 ad_light(i,k)=0.0_r8
2960!^ tl_PAR=(-tl_Att*PAR1+tl_Itop*(1.0_r8-ExpAtt)- &
2961!^ & Itop*tl_ExpAtt)/Att
2962!^
2963 adfac=ad_par/att
2964 ad_att=ad_att-par1*adfac
2965 ad_expatt=ad_expatt-itop*adfac
2966 ad_itop=ad_itop+(1.0_r8-expatt)*adfac
2967 ad_par=0.0_r8
2968!^ tl_Itop=tl_PAR
2969!^
2970 ad_par=ad_par+ad_itop
2971 ad_itop=0.0_r8
2972!^ tl_ExpAtt=-ExpAtt*tl_Att
2973!^
2974 ad_att=ad_att-expatt*ad_expatt
2975 ad_expatt=0.0_r8
2976!^ tl_Att=AttPhy(ng)*tl_Bio(i,k,iPhyt)* &
2977!^ & (z_w(i,j,k)-z_w(i,j,k-1))+ &
2978!^ & (AttSW(ng)+AttPhy(ng)*Bio1(i,k,iPhyt))* &
2979!^ & (tl_z_w(i,j,k)-tl_z_w(i,j,k-1))
2980!^
2981 adfac=(attsw(ng)+attphy(ng)*bio1(i,k,iphyt))*ad_att
2982 ad_bio(i,k,iphyt)=ad_bio(i,k,iphyt)+ &
2983 & attphy(ng)*(z_w(i,j,k)-z_w(i,j,k-1))* &
2984 & ad_att
2985 ad_z_w(i,j,k-1)=ad_z_w(i,j,k-1)-adfac
2986 ad_z_w(i,j,k )=ad_z_w(i,j,k )+adfac
2987 ad_att=0.0_r8
2988 END DO
2989 ELSE ! night time
2990 DO k=1,n(ng)
2991!^ tl_Light(i,k)=0.0_r8
2992!^
2993 ad_light(i,k)=0.0_r8
2994 END DO
2995 END IF
2996!^ tl_PAR=tl_PARsur(i)
2997!^
2998 ad_parsur(i)=ad_parsur(i)+ad_par
2999 ad_par=0.0_r8
3000 END DO
3001
3002 END DO iter_loop1
3003!
3004! Calculate adjoint surface Photosynthetically Available Radiation
3005! (PAR). The net shortwave radiation is scaled back to Watts/m2
3006! and multiplied by the fraction that is photosynthetically
3007! available, PARfrac.
3008!
3009 DO i=istr,iend
3010#ifdef CONST_PAR
3011!
3012! Specify constant surface irradiance a la Powell and Spitz.
3013!
3014!^ tl_PARsur(i)=0.0_r8
3015!^
3016!! ad_PARsur(i)=0.0_r8
3017#else
3018!^ tl_PARsur(i)=(tl_PARfrac(ng)*srflx(i,j)+ &
3019!^ & PARfrac(ng)*tl_srflx(i,j))*rho0*Cp
3020!^
3021 adfac=rho0*cp*ad_parsur(i)
3022 ad_srflx(i,j)=ad_srflx(i,j)+parfrac(ng)*adfac
3023 ad_parfrac(ng)=ad_parfrac(ng)+srflx(i,j)*adfac
3024!! ad_PARsur(i)=0.0_r8
3025#endif
3026 END DO
3027!
3028! Restrict biological tracer to be positive definite. If a negative
3029! concentration is detected, nitrogen is drawn from the most abundant
3030! pool to supplement the negative pools to a lower limit of MinVal
3031! which is set to 1E-6 above.
3032!
3033 DO k=1,n(ng)
3034 DO i=istr,iend
3035!
3036! Adjoint load biological tracers into local arrays.
3037!
3038 DO itrc=1,nbt
3039 ibio=idbio(itrc)
3040!^ tl_Bio(i,k,ibio)=tl_BioTrc(ibio,nstp)
3041!^
3042 ad_biotrc(ibio,nstp)=ad_biotrc(ibio,nstp)+ &
3043 & ad_bio(i,k,ibio)
3044 ad_bio(i,k,ibio)=0.0_r8
3045!^ tl_Bio_old(i,k,ibio)=tl_BioTrc(ibio,nstp)
3046!^
3047 ad_biotrc(ibio,nstp)=ad_biotrc(ibio,nstp)+ &
3048 & ad_bio_old(i,k,ibio)
3049 ad_bio_old(i,k,ibio)=0.0_r8
3050 END DO
3051!
3052! Adjoint positive definite concentrations.
3053!
3054 cff2=0.0_r8
3055 DO itime=1,2
3056 cff1=0.0_r8
3057 itrcmax=idbio(1)
3058 DO itrc=1,nbt
3059 ibio=idbio(itrc)
3060!
3061! The basic state (nstp and nnew indices) that is read from the
3062! forward file is in units of tracer. Since BioTrc(ibio,:) is in
3063! tracer units, we simply use t instead of t*Hz_inv.
3064!
3065 biotrc(ibio,itime)=t(i,j,k,itime,ibio)
3066 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
3067 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
3068 itrcmax=ibio
3069 END IF
3070 biotrc1(ibio,itime)=biotrc(ibio,itime)
3071 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
3072 END DO
3073 IF (biotrc(itrcmax,itime).gt.cff1) THEN
3074!^ tl_BioTrc(iTrcMax,itime)=tl_BioTrc(iTrcMax,itime)- &
3075!^ & tl_cff1
3076!^
3077 ad_cff1=-ad_biotrc(itrcmax,itime)
3078 END IF
3079
3080 cff1=0.0_r8
3081 DO itrc=1,nbt
3082 ibio=idbio(itrc)
3083!
3084! The basic state (nstp and nnew indices) that is read from the
3085! forward file is in units of tracer. Since BioTrc(ibio,:) is in
3086! tracer units, we simply use t instead of t*Hz_inv.
3087!
3088 biotrc(ibio,itime)=t(i,j,k,itime,ibio)
3089 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
3090!^ tl_BioTrc(ibio,itime)=(0.5_r8- &
3091!^ & SIGN(0.5_r8, &
3092!^ & MinVal- &
3093!^ & BioTrc1(ibio,itime)))* &
3094!^ & tl_BioTrc(ibio,itime)
3095!^
3096 ad_biotrc(ibio,itime)=(0.5_r8- &
3097 & sign(0.5_r8, &
3098 & minval- &
3099 & biotrc1(ibio,itime)))* &
3100 & ad_biotrc(ibio,itime)
3101!^ tl_cff1=tl_cff1- &
3102!^ & (0.5_r8-SIGN(0.5_r8, &
3103!^ & BioTrc(ibio,itime)-MinVal))* &
3104!^ & tl_BioTrc(ibio,itime)
3105!^
3106 ad_biotrc(ibio,itime)=ad_biotrc(ibio,itime)- &
3107 & (0.5_r8-sign(0.5_r8, &
3108 & biotrc(ibio,itime)- &
3109 & minval))*ad_cff1
3110 END DO
3111 ad_cff1=0.0_r8
3112 END DO
3113!
3114! At input, all tracers (index nnew) from predictor step have
3115! transport units (m Tunits) since we do not have yet the new
3116! values for zeta and Hz. These are known after the 2D barotropic
3117! time-stepping.
3118!
3119! NOTE: In the following code, t(:,:,:,nnew,:) should be in units of
3120! tracer times depth. However the basic state (nstp and nnew
3121! indices) that is read from the forward file is in units of
3122! tracer. Since BioTrc(ibio,nnew) is in tracer units, we simply
3123! use t instead of t*Hz_inv.
3124!
3125 DO itrc=1,nbt
3126 ibio=idbio(itrc)
3127!^ tl_BioTrc(ibio,nnew)=tl_t(i,j,k,nnew,ibio)* &
3128!^ & Hz_inv(i,k)+ &
3129!^ & t(i,j,k,nnew,ibio)*Hz(i,j,k)* &
3130!^ & tl_Hz_inv(i,k)
3131!^
3132 ad_hz_inv(i,k)=ad_hz_inv(i,k)+ &
3133 & t(i,j,k,nnew,ibio)*hz(i,j,k)* &
3134 & ad_biotrc(ibio,nnew)
3135 ad_t(i,j,k,nnew,ibio)=ad_t(i,j,k,nnew,ibio)+ &
3136 & hz_inv(i,k)*ad_biotrc(ibio,nnew)
3137 ad_biotrc(ibio,nnew)=0.0_r8
3138!^ tl_BioTrc(ibio,nstp)=tl_t(i,j,k,nstp,ibio)
3139!^
3140 ad_t(i,j,k,nstp,ibio)=ad_t(i,j,k,nstp,ibio)+ &
3141 & ad_biotrc(ibio,nstp)
3142 ad_biotrc(ibio,nstp)=0.0_r8
3143 END DO
3144 END DO
3145 END DO
3146!
3147! Adjoint inverse thickness to avoid repeated divisions.
3148!
3149 DO k=2,n(ng)-1
3150 DO i=istr,iend
3151!^ tl_Hz_inv3(i,k)=-Hz_inv3(i,k)*Hz_inv3(i,k)* &
3152!^ & (tl_Hz(i,j,k-1)+tl_Hz(i,j,k)+ &
3153!^ & tl_Hz(i,j,k+1))
3154!^
3155 adfac=hz_inv3(i,k)*hz_inv3(i,k)*ad_hz_inv3(i,k)
3156 ad_hz(i,j,k-1)=ad_hz(i,j,k-1)-adfac
3157 ad_hz(i,j,k )=ad_hz(i,j,k )-adfac
3158 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)-adfac
3159 ad_hz_inv3(i,k)=0.0_r8
3160 END DO
3161 END DO
3162 DO k=1,n(ng)-1
3163 DO i=istr,iend
3164!^ tl_Hz_inv2(i,k)=-Hz_inv2(i,k)*Hz_inv2(i,k)* &
3165!^ & (tl_Hz(i,j,k)+tl_Hz(i,j,k+1))
3166!^
3167 adfac=hz_inv2(i,k)*hz_inv2(i,k)*ad_hz_inv2(i,k)
3168 ad_hz(i,j,k )=ad_hz(i,j,k )-adfac
3169 ad_hz(i,j,k+1)=ad_hz(i,j,k+1)-adfac
3170 ad_hz_inv2(i,k)=0.0_r8
3171 END DO
3172 END DO
3173 DO k=1,n(ng)
3174 DO i=istr,iend
3175!^ tl_Hz_inv(i,k)=-Hz_inv(i,k)*Hz_inv(i,k)*tl_Hz(i,j,k)
3176!^
3177 ad_hz(i,j,k)=ad_hz(i,j,k)- &
3178 & hz_inv(i,k)*hz_inv(i,k)*ad_hz_inv(i,k)
3179 ad_hz_inv(i,k)=0.0_r8
3180 END DO
3181 END DO
3182
3183 END DO j_loop
3184!
3185! Set adjoint vertical sinking velocity vector in the same order as the
3186! identification vector, IDSINK.
3187!
3188!^ tl_Wbio(2)=tl_wDet(ng) ! Small detritus
3189!^
3190 ad_wdet(ng)=ad_wdet(ng)+ad_wbio(2)
3191 ad_wbio(2)=0.0_r8
3192!^ tl_Wbio(1)=tl_wPhy(ng) ! Phytoplankton
3193!^
3194 ad_wphy(ng)=ad_wphy(ng)+ad_wbio(1)
3195 ad_wbio(1)=0.0_r8
3196
3197 RETURN
3198 END SUBROUTINE ad_npzd_powell_tile
3199
3200 END MODULE ad_biology_mod
ad_biology_mod
Definition ad_biology.F:26
ad_biology_mod::ad_npzd_powell_tile
subroutine ad_npzd_powell_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, srflx, ad_srflx, t, ad_t)
Definition ad_npzd_Powell.h:104
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::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::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::ad_parfrac
real(r8), dimension(:), allocatable ad_parfrac
Definition npzd_iron_mod.h:99
mod_biology::wphy
real(r8), dimension(:), allocatable wphy
Definition fennel_mod.h:154
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::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_biology::ad_wphy
real(r8), dimension(:), allocatable ad_wphy
Definition npzd_iron_mod.h:117
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::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::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
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