ROMS
Loading...
Searching...
No Matches
npzd_Powell.h
Go to the documentation of this file.
1 MODULE biology_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group Craig V. Lewis !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! Nutrient-Phytoplankton-Zooplankton-Detritus Model. !
11! !
12! This routine computes the biological sources and sinks and adds !
13! then the global biological fields. !
14! !
15! Reference: !
16! !
17! Powell, T.M., C.V.W. Lewis, E. Curchitser, D. Haidvogel, !
18! Q. Hermann and E. Dobbins, 2006: Results from a three- !
19! dimensional, nested biological-physical model of the !
20! California Current System: Comparisons with Statistics !
21! from Satellite Imagery, J. Geophys. Res. !
22! !
23! Adapted from code written originally by Craig V. Lewis. !
24! !
25!=======================================================================
26!
27 implicit none
28!
29 PRIVATE
30 PUBLIC :: biology
31!
32 CONTAINS
33!
34!***********************************************************************
35 SUBROUTINE biology (ng,tile)
36!***********************************************************************
37!
38 USE mod_param
39 USE mod_forces
40 USE mod_grid
41 USE mod_ncparam
42 USE mod_ocean
43 USE mod_stepping
44!
45! Imported variable declarations.
46!
47 integer, intent(in) :: ng, tile
48!
49! Local variable declarations.
50!
51 character (len=*), parameter :: MyFile = &
52 & __FILE__
53!
54#include "tile.h"
55!
56! Set header file name.
57!
58#ifdef DISTRIBUTE
59 IF (lbiofile(inlm)) THEN
60#else
61 IF (lbiofile(inlm).and.(tile.eq.0)) THEN
62#endif
63 lbiofile(inlm)=.false.
64 bioname(inlm)=myfile
65 END IF
66!
67#ifdef PROFILE
68 CALL wclock_on (ng, inlm, 15, __line__, myfile)
69#endif
70 CALL npzd_powell_tile (ng, tile, &
71 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
72 & imins, imaxs, jmins, jmaxs, &
73 & nstp(ng), nnew(ng), &
74#ifdef MASKING
75 & grid(ng) % rmask, &
76#endif
77 & grid(ng) % Hz, &
78 & grid(ng) % z_r, &
79 & grid(ng) % z_w, &
80 & forces(ng) % srflx, &
81 & ocean(ng) % t)
82
83#ifdef PROFILE
84 CALL wclock_off (ng, inlm, 15, __line__, myfile)
85#endif
86!
87 RETURN
88 END SUBROUTINE biology
89!
90!-----------------------------------------------------------------------
91 SUBROUTINE npzd_powell_tile (ng, tile, &
92 & LBi, UBi, LBj, UBj, UBk, UBt, &
93 & IminS, ImaxS, JminS, JmaxS, &
94 & nstp, nnew, &
95#ifdef MASKING
96 & rmask, &
97#endif
98 & Hz, z_r, z_w, &
99 & srflx, &
100 & t)
101!-----------------------------------------------------------------------
102!
103 USE mod_param
104 USE mod_biology
105 USE mod_ncparam
106 USE mod_scalars
107!
108! Imported variable declarations.
109!
110 integer, intent(in) :: ng, tile
111 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
112 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
113 integer, intent(in) :: nstp, nnew
114
115#ifdef ASSUMED_SHAPE
116# ifdef MASKING
117 real(r8), intent(in) :: rmask(LBi:,LBj:)
118# endif
119 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
120 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
121 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
122 real(r8), intent(in) :: srflx(LBi:,LBj:)
123 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
124#else
125# ifdef MASKING
126 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
127# endif
128 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
129 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
130 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
131 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
132 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
133#endif
134!
135! Local variable declarations.
136!
137 integer, parameter :: Nsink = 2
138
139 integer :: Iter, i, ibio, isink, itime, itrc, iTrcMax, j, k, ks
140
141 integer, dimension(Nsink) :: idsink
142
143 real(r8), parameter :: MinVal = 1.0e-6_r8
144
145 real(r8) :: Att, ExpAtt, Itop, PAR
146 real(r8) :: cff, cff1, cff2, cff3, cff4, dtdays
147 real(r8) :: cffL, cffR, cu, dltL, dltR
148
149 real(r8), dimension(Nsink) :: Wbio
150
151 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
152
153 real(r8), dimension(IminS:ImaxS) :: PARsur
154
155 real(r8), dimension(NT(ng),2) :: BioTrc
156
157 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
158 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
159
160 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
161
162 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
163 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
164 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
165 real(r8), dimension(IminS:ImaxS,N(ng)) :: Light
166 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
167 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
168 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
169 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
170 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
171
172#include "set_bounds.h"
173!
174!-----------------------------------------------------------------------
175! Add biological Source/Sink terms.
176!-----------------------------------------------------------------------
177!
178! Avoid computing source/sink terms if no biological iterations.
179!
180 IF (bioiter(ng).le.0) RETURN
181!
182! Set time-stepping size (days) according to the number of iterations.
183!
184 dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
185!
186! Set vertical sinking indentification vector.
187!
188 idsink(1)=iphyt ! Phytoplankton
189 idsink(2)=isdet ! Small detritus
190!
191! Set vertical sinking velocity vector in the same order as the
192! identification vector, IDSINK.
193!
194 wbio(1)=wphy(ng) ! Phytoplankton
195 wbio(2)=wdet(ng) ! Small detritus
196!
197! Compute inverse thickness to avoid repeated divisions.
198!
199 j_loop : DO j=jstr,jend
200 DO k=1,n(ng)
201 DO i=istr,iend
202 hz_inv(i,k)=1.0_r8/hz(i,j,k)
203 END DO
204 END DO
205 DO k=1,n(ng)-1
206 DO i=istr,iend
207 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
208 END DO
209 END DO
210 DO k=2,n(ng)-1
211 DO i=istr,iend
212 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
213 END DO
214 END DO
215!
216! Restrict biological tracer to be positive definite. If a negative
217! concentration is detected, nitrogen is drawn from the most abundant
218! pool to supplement the negative pools to a lower limit of MinVal
219! which is set to 1E-6 above.
220!
221 DO k=1,n(ng)
222 DO i=istr,iend
223!
224! At input, all tracers (index nnew) from predictor step have
225! transport units (m Tunits) since we do not have yet the new
226! values for zeta and Hz. These are known after the 2D barotropic
227! time-stepping.
228!
229 DO itrc=1,nbt
230 ibio=idbio(itrc)
231 biotrc(ibio,nstp)=t(i,j,k,nstp,ibio)
232 biotrc(ibio,nnew)=t(i,j,k,nnew,ibio)*hz_inv(i,k)
233 END DO
234!
235! Impose positive definite concentrations.
236!
237 cff2=0.0_r8
238 DO itime=1,2
239 cff1=0.0_r8
240 itrcmax=idbio(1)
241 DO itrc=1,nbt
242 ibio=idbio(itrc)
243 cff1=cff1+max(0.0_r8,minval-biotrc(ibio,itime))
244 IF (biotrc(ibio,itime).gt.biotrc(itrcmax,itime)) THEN
245 itrcmax=ibio
246 END IF
247 biotrc(ibio,itime)=max(minval,biotrc(ibio,itime))
248 END DO
249 IF (biotrc(itrcmax,itime).gt.cff1) THEN
250 biotrc(itrcmax,itime)=biotrc(itrcmax,itime)-cff1
251 END IF
252 END DO
253!
254! Load biological tracers into local arrays.
255!
256 DO itrc=1,nbt
257 ibio=idbio(itrc)
258 bio_old(i,k,ibio)=biotrc(ibio,nstp)
259 bio(i,k,ibio)=biotrc(ibio,nstp)
260 END DO
261 END DO
262 END DO
263!
264! Calculate surface Photosynthetically Available Radiation (PAR). The
265! net shortwave radiation is scaled back to Watts/m2 and multiplied by
266! the fraction that is photosynthetically available, PARfrac.
267!
268 DO i=istr,iend
269#ifdef CONST_PAR
270!
271! Specify constant surface irradiance a la Powell and Spitz.
272!
273 parsur(i)=158.075_r8
274#else
275 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
276#endif
277 END DO
278!
279!=======================================================================
280! Start internal iterations to achieve convergence of the nonlinear
281! backward-implicit solution.
282!=======================================================================
283!
284! During the iterative procedure a series of fractional time steps are
285! performed in a chained mode (splitting by different biological
286! conversion processes) in sequence of the main food chain. In all
287! stages the concentration of the component being consumed is treated
288! in a fully implicit manner, so the algorithm guarantees non-negative
289! values, no matter how strong the concentration of active consuming
290! component (Phytoplankton or Zooplankton). The overall algorithm,
291! as well as any stage of it, is formulated in conservative form
292! (except explicit sinking) in sense that the sum of concentration of
293! all components is conserved.
294!
295! In the implicit algorithm, we have for example (N: nutrient,
296! P: phytoplankton),
297!
298! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
299! {Michaelis-Menten}
300! below, we set
301! The N in the numerator of
302! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
303! as N(new)
304!
305! so the time-stepping of the equations becomes:
306!
307! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
308! consuming, divide by (1 + cff)
309! and
310!
311! P(new) = P(old) + cff * N(new) (2) when adding a source term,
312! growing, add (cff * source)
313!
314! Notice that if you substitute (1) in (2), you will get:
315!
316! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
317!
318! If you add (1) and (3), you get
319!
320! N(new) + P(new) = N(old) + P(old)
321!
322! implying conservation regardless how "cff" is computed. Therefore,
323! this scheme is unconditionally stable regardless of the conversion
324! rate. It does not generate negative values since the constituent
325! to be consumed is always treated implicitly. It is also biased
326! toward damping oscillations.
327!
328! The iterative loop below is to iterate toward an universal Backward-
329! Euler treatment of all terms. So if there are oscillations in the
330! system, they are only physical oscillations. These iterations,
331! however, do not improve the accuaracy of the solution.
332!
333 iter_loop: DO iter=1,bioiter(ng)
334!
335! Compute light attenuation as function of depth.
336!
337 DO i=istr,iend
338 par=parsur(i)
339 IF (parsur(i).gt.0.0_r8) THEN ! day time
340 DO k=n(ng),1,-1
341!
342! Compute average light attenuation for each grid cell. Here, AttSW is
343! the light attenuation due to seawater and AttPhy is the attenuation
344! due to phytoplankton (self-shading coefficient).
345!
346 att=(attsw(ng)+attphy(ng)*bio(i,k,iphyt))* &
347 & (z_w(i,j,k)-z_w(i,j,k-1))
348 expatt=exp(-att)
349 itop=par
350 par=itop*(1.0_r8-expatt)/att ! average at cell center
351 light(i,k)=par
352!
353! Light attenuation at the bottom of the grid cell. It is the starting
354! PAR value for the next (deeper) vertical grid cell.
355!
356 par=itop*expatt
357 END DO
358 ELSE ! night time
359 DO k=1,n(ng)
360 light(i,k)=0.0_r8
361 END DO
362 END IF
363 END DO
364!
365! Phytoplankton photosynthetic growth and nitrate uptake (Vm_NO3 rate).
366! The Michaelis-Menten curve is used to describe the change in uptake
367! rate as a function of nitrate concentration. Here, PhyIS is the
368! initial slope of the P-I curve and K_NO3 is the half saturation of
369! phytoplankton nitrate uptake.
370!
371 cff1=dtdays*vm_no3(ng)*phyis(ng)
372 cff2=vm_no3(ng)*vm_no3(ng)
373 cff3=phyis(ng)*phyis(ng)
374 DO k=1,n(ng)
375 DO i=istr,iend
376 cff4=1.0_r8/sqrt(cff2+cff3*light(i,k)*light(i,k))
377 cff=bio(i,k,iphyt)* &
378 & cff1*cff4*light(i,k)/ &
379 & (k_no3(ng)+bio(i,k,ino3_))
380 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff)
381 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
382 & bio(i,k,ino3_)*cff
383 END DO
384 END DO
385!
386! Grazing on phytoplankton by zooplankton (ZooGR rate) using the Ivlev
387! formulation (Ivlev, 1955) and lost of phytoplankton to the nitrate
388! pool as function of "sloppy feeding" and metabolic processes
389! (ZooEEN and ZooEED fractions).
390!
391 cff1=dtdays*zoogr(ng)
392 cff2=1.0_r8-zooeen(ng)-zooeed(ng)
393 DO k=1,n(ng)
394 DO i=istr,iend
395 cff=bio(i,k,izoop)* &
396 & cff1*(1.0_r8-exp(-ivlev(ng)*bio(i,k,iphyt)))/ &
397 & bio(i,k,iphyt)
398 bio(i,k,iphyt)=bio(i,k,iphyt)/(1.0_r8+cff)
399 bio(i,k,izoop)=bio(i,k,izoop)+ &
400 & bio(i,k,iphyt)*cff2*cff
401 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
402 & bio(i,k,iphyt)*zooeen(ng)*cff
403 bio(i,k,isdet)=bio(i,k,isdet)+ &
404 & bio(i,k,iphyt)*zooeed(ng)*cff
405 END DO
406 END DO
407!
408! Phytoplankton mortality to nutrients (PhyMRN rate) and detritus
409! (PhyMRD rate).
410!
411 cff3=dtdays*phymrd(ng)
412 cff2=dtdays*phymrn(ng)
413 cff1=1.0_r8/(1.0_r8+cff2+cff3)
414 DO k=1,n(ng)
415 DO i=istr,iend
416 bio(i,k,iphyt)=bio(i,k,iphyt)*cff1
417 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
418 & bio(i,k,iphyt)*cff2
419 bio(i,k,isdet)=bio(i,k,isdet)+ &
420 & bio(i,k,iphyt)*cff3
421 END DO
422 END DO
423!
424! Zooplankton mortality to nutrients (ZooMRN rate) and Detritus
425! (ZooMRD rate).
426!
427 cff3=dtdays*zoomrd(ng)
428 cff2=dtdays*zoomrn(ng)
429 cff1=1.0_r8/(1.0_r8+cff2+cff3)
430 DO k=1,n(ng)
431 DO i=istr,iend
432 bio(i,k,izoop)=bio(i,k,izoop)*cff1
433 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
434 & bio(i,k,izoop)*cff2
435 bio(i,k,isdet)=bio(i,k,isdet)+ &
436 & bio(i,k,izoop)*cff3
437 END DO
438 END DO
439!
440! Detritus breakdown to nutrients: remineralization (DetRR rate).
441!
442 cff2=dtdays*detrr(ng)
443 cff1=1.0_r8/(1.0_r8+cff2)
444 DO k=1,n(ng)
445 DO i=istr,iend
446 bio(i,k,isdet)=bio(i,k,isdet)*cff1
447 bio(i,k,ino3_)=bio(i,k,ino3_)+ &
448 & bio(i,k,isdet)*cff2
449 END DO
450 END DO
451!
452!-----------------------------------------------------------------------
453! Vertical sinking terms: Phytoplankton and Detritus
454!-----------------------------------------------------------------------
455!
456! Reconstruct vertical profile of selected biological constituents
457! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
458! grid box. Then, compute semi-Lagrangian flux due to sinking.
459!
460 sink_loop: DO isink=1,nsink
461 ibio=idsink(isink)
462!
463! Copy concentration of biological particulates into scratch array
464! "qc" (q-central, restrict it to be positive) which is hereafter
465! interpreted as a set of grid-box averaged values for biogeochemical
466! constituent concentration.
467!
468 DO k=1,n(ng)
469 DO i=istr,iend
470 qc(i,k)=bio(i,k,ibio)
471 END DO
472 END DO
473!
474 DO k=n(ng)-1,1,-1
475 DO i=istr,iend
476 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
477 END DO
478 END DO
479 DO k=2,n(ng)-1
480 DO i=istr,iend
481 dltr=hz(i,j,k)*fc(i,k)
482 dltl=hz(i,j,k)*fc(i,k-1)
483 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
484 cffr=cff*fc(i,k)
485 cffl=cff*fc(i,k-1)
486!
487! Apply PPM monotonicity constraint to prevent oscillations within the
488! grid box.
489!
490 IF ((dltr*dltl).le.0.0_r8) THEN
491 dltr=0.0_r8
492 dltl=0.0_r8
493 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
494 dltr=cffl
495 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
496 dltl=cffr
497 END IF
498!
499! Compute right and left side values (bR,bL) of parabolic segments
500! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
501!
502! NOTE: Although each parabolic segment is monotonic within its grid
503! box, monotonicity of the whole profile is not guaranteed,
504! because bL(k+1)-bR(k) may still have different sign than
505! qc(i,k+1)-qc(i,k). This possibility is excluded,
506! after bL and bR are reconciled using WENO procedure.
507!
508 cff=(dltr-dltl)*hz_inv3(i,k)
509 dltr=dltr-cff*hz(i,j,k+1)
510 dltl=dltl+cff*hz(i,j,k-1)
511 br(i,k)=qc(i,k)+dltr
512 bl(i,k)=qc(i,k)-dltl
513 wr(i,k)=(2.0_r8*dltr-dltl)**2
514 wl(i,k)=(dltr-2.0_r8*dltl)**2
515 END DO
516 END DO
517 cff=1.0e-14_r8
518 DO k=2,n(ng)-2
519 DO i=istr,iend
520 dltl=max(cff,wl(i,k ))
521 dltr=max(cff,wr(i,k+1))
522 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
523 bl(i,k+1)=br(i,k)
524 END DO
525 END DO
526 DO i=istr,iend
527 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
528#if defined LINEAR_CONTINUATION
529 bl(i,n(ng))=br(i,n(ng)-1)
530 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
531#elif defined NEUMANN
532 bl(i,n(ng))=br(i,n(ng)-1)
533 br(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
534#else
535 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
536 bl(i,n(ng))=qc(i,n(ng)) ! conditions
537 br(i,n(ng)-1)=qc(i,n(ng))
538#endif
539#if defined LINEAR_CONTINUATION
540 br(i,1)=bl(i,2)
541 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
542#elif defined NEUMANN
543 br(i,1)=bl(i,2)
544 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
545#else
546 bl(i,2)=qc(i,1) ! bottom grid boxes are
547 br(i,1)=qc(i,1) ! re-assumed to be
548 bl(i,1)=qc(i,1) ! piecewise constant.
549#endif
550 END DO
551!
552! Apply monotonicity constraint again, since the reconciled interfacial
553! values may cause a non-monotonic behavior of the parabolic segments
554! inside the grid box.
555!
556 DO k=1,n(ng)
557 DO i=istr,iend
558 dltr=br(i,k)-qc(i,k)
559 dltl=qc(i,k)-bl(i,k)
560 cffr=2.0_r8*dltr
561 cffl=2.0_r8*dltl
562 IF ((dltr*dltl).lt.0.0_r8) THEN
563 dltr=0.0_r8
564 dltl=0.0_r8
565 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
566 dltr=cffl
567 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
568 dltl=cffr
569 END IF
570 br(i,k)=qc(i,k)+dltr
571 bl(i,k)=qc(i,k)-dltl
572 END DO
573 END DO
574!
575! After this moment reconstruction is considered complete. The next
576! stage is to compute vertical advective fluxes, FC. It is expected
577! that sinking may occurs relatively fast, the algorithm is designed
578! to be free of CFL criterion, which is achieved by allowing
579! integration bounds for semi-Lagrangian advective flux to use as
580! many grid boxes in upstream direction as necessary.
581!
582! In the two code segments below, WL is the z-coordinate of the
583! departure point for grid box interface z_w with the same indices;
584! FC is the finite volume flux; ksource(:,k) is index of vertical
585! grid box which contains the departure point (restricted by N(ng)).
586! During the search: also add in content of whole grid boxes
587! participating in FC.
588!
589 cff=dtdays*abs(wbio(isink))
590 DO k=1,n(ng)
591 DO i=istr,iend
592 fc(i,k-1)=0.0_r8
593 wl(i,k)=z_w(i,j,k-1)+cff
594 wr(i,k)=hz(i,j,k)*qc(i,k)
595 ksource(i,k)=k
596 END DO
597 END DO
598 DO k=1,n(ng)
599 DO ks=k,n(ng)-1
600 DO i=istr,iend
601 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
602 ksource(i,k)=ks+1
603 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
604 END IF
605 END DO
606 END DO
607 END DO
608!
609! Finalize computation of flux: add fractional part.
610!
611 DO k=1,n(ng)
612 DO i=istr,iend
613 ks=ksource(i,k)
614 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
615 fc(i,k-1)=fc(i,k-1)+ &
616 & hz(i,j,ks)*cu* &
617 & (bl(i,ks)+ &
618 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
619 & (1.5_r8-cu)* &
620 & (br(i,ks)+bl(i,ks)- &
621 & 2.0_r8*qc(i,ks))))
622 END DO
623 END DO
624 DO k=1,n(ng)
625 DO i=istr,iend
626 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
627 END DO
628 END DO
629
630 END DO sink_loop
631 END DO iter_loop
632!
633!-----------------------------------------------------------------------
634! Update global tracer variables: Add increment due to BGC processes
635! to tracer array in time index "nnew". Index "nnew" is solution after
636! advection and mixing and has transport units (m Tunits) hence the
637! increment is multiplied by Hz. Notice that we need to subtract
638! original values "Bio_old" at the top of the routine to just account
639! for the concentractions affected by BGC processes. This also takes
640! into account any constraints (non-negative concentrations, carbon
641! concentration range) specified before entering BGC kernel. If "Bio"
642! were unchanged by BGC processes, the increment would be exactly
643! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
644! bounded >=0 so that we can preserve total inventory of nutrients
645! when advection causes tracer concentration to go negative.
646!-----------------------------------------------------------------------
647!
648 DO itrc=1,nbt
649 ibio=idbio(itrc)
650 DO k=1,n(ng)
651 DO i=istr,iend
652 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
653 t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*hz(i,j,k)
654 END DO
655 END DO
656 END DO
657
658 END DO j_loop
659!
660 RETURN
661 END SUBROUTINE npzd_powell_tile
662
663 END MODULE biology_mod
subroutine npzd_powell_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nnew, rmask, hz, z_r, z_w, srflx, t)
subroutine, public biology(ng, tile)
Definition ecosim.h:57
real(r8), dimension(:), allocatable parfrac
Definition fennel_mod.h:139
real(r8), dimension(:), allocatable phymrd
real(r8), dimension(:), allocatable zooeed
real(r8), dimension(:), allocatable detrr
real(r8), dimension(:), allocatable wdet
integer, dimension(:), allocatable bioiter
Definition ecosim_mod.h:343
real(r8), dimension(:), allocatable phyis
Definition fennel_mod.h:143
real(r8), dimension(:), allocatable attsw
Definition fennel_mod.h:125
integer ino3_
Definition ecosim_mod.h:277
real(r8), dimension(:), allocatable zoogr
Definition fennel_mod.h:160
real(r8), dimension(:), allocatable k_no3
Definition fennel_mod.h:133
integer iphyt
Definition fennel_mod.h:81
real(r8), dimension(:), allocatable ivlev
real(r8), dimension(:), allocatable attphy
integer, dimension(:), allocatable idbio
Definition ecosim_mod.h:256
real(r8), dimension(:), allocatable wphy
Definition fennel_mod.h:154
real(r8), dimension(:), allocatable vm_no3
real(r8), dimension(:), allocatable zoomrn
real(r8), dimension(:), allocatable phymrn
integer izoop
Definition fennel_mod.h:82
real(r8), dimension(:), allocatable zoomrd
real(r8), dimension(:), allocatable zooeen
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
integer, parameter inlm
Definition mod_param.F:662
integer, dimension(:), allocatable n
Definition mod_param.F:479
integer nbt
Definition mod_param.F:509
integer, dimension(:), allocatable nt
Definition mod_param.F:489
real(dp), dimension(:), allocatable dt
real(dp) cp
real(dp), parameter sec2day
real(dp) rho0
integer, dimension(:), allocatable nnew
integer, dimension(:), allocatable nstp
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3