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