ROMS
Loading...
Searching...
No Matches
fennel.h
Go to the documentation of this file.
1 MODULE biology_mod
2!
3!git $Id$
4!=======================================================================
5! Copyright (c) 2002-2025 The ROMS Group !
6! Licensed under a MIT/X style license Hernan G. Arango !
7! See License_ROMS.md Katja Fennel !
8!========================================== Alexander F. Shchepetkin ===
9! !
10! This routine computes the biological sources and sinks for the !
11! Fennel et al. (2006) ecosystem model. Then, it adds those terms !
12! to the global biological fields. !
13! !
14! This model is loosely based on the model by Fasham et al. (1990) !
15! but it differs in many respects. The detailed equations of the !
16! nitrogen cycling component are given in Fennel et al. (2006). !
17! Nitrogen is the fundamental elemental currency in this model. !
18! This model was adapted from a code written originally by John !
19! Moisan and Emanule DiLorenzo. !
20! !
21! It is recommended to activate always the "BIO_SEDIMENT" option !
22! to ensure conservation of mass by converting the organic matter !
23! that is sinking out of the bottom most grid cell into inorganic !
24! nutrients (i.e., instantanaous remineralization at the water- !
25! sediment interface). Additionally, the "DENITRIFICATION" option !
26! can be activated. Hence, a fraction of the instantenous bottom !
27! remineralization is assumed to occur through the anearobic !
28! (denitrification) pathway and thus lost from the pool of !
29! biologically availalbe fixed nitrogen. See Fennel et al. (2006) !
30! for details. !
31! !
32! Additional options can be activated to enable simulation of !
33! inorganic carbon and dissolved oxygen. Accounting of inorganic !
34! carbon is activated by the "CARBON" option, and results in two !
35! additional biological tracer variables: DIC and alkalinity. !
36! See Fennel et al. (2008) for details. !
37! !
38! If the "pCO2_RZ" options is activated, in addition to "CARBON", !
39! the carbonate system routines by Zeebe and Wolf-Gladrow (2001) !
40! are used, while the OCMIP standard routines are the default. !
41! There are two different ways of treating alkalinity. It can be !
42! treated diagnostically (default), in this case alkalinity acts !
43! like a passive tracer that is not affected by changes in the !
44! concentration of nitrate or ammonium. However, if the option !
45! "TALK_NONCONSERV" is used, the alkalinity will be affected by !
46! sources and sinks in nitrate. See Fennel et al. (2008) for more !
47! details. !
48! !
49! If the "OXYGEN" option is activated, one additional biological !
50! tracer variable for dissolved oxygen. "OXYGEN" can be activated !
51! independently of the "CARBON" option. If "OCMIP_OXYGEN_SC" is !
52! used, in addition to "OXYGEN", the Schmidt number of oxygen in !
53! seawater will be computed using the formulation proposed by !
54! Keeling et al. (1998, Global Biogeochem. Cycles, 12, 141-163). !
55! Otherwise, the Wanninkhof (1992) formula will be used. See !
56! Fennel et al. (2013) for more details. !
57! !
58!***********************************************************************
59! UPDATE Sept 2020 !
60! !
61! In this version additional tracers, additional alkalinity !
62! fluxes, and an updated parameterization of air-sea O2 and CO2 !
63! fluxes are added as described in Laurent et al. (2017). !
64! !
65! If "PO4" is activated, one additional biological tracer is !
66! added representing phosphate. With this option phytoplankton !
67! growth can be limited by either nitrogen and phosphorus. This !
68! option was introduced in Laurent et al. (2012). !
69! !
70! If "RIVER_DON" is activated, an additional biological tracer !
71! (or 2 if "CARBON" is defined) is added representing non-sinking !
72! dissolved organic matter from rivers as described in Yu et al. !
73! (2015). !
74! !
75! If the "RW14_OXYGEN_SC" and/or "RW14_CO2_SC" options are used, !
76! the model will use Wanninkhof (2014) air-sea flux parameteri- !
77! zation. With the "TALK_NONCONSERV" option, alkalinity is !
78! affected by biological fluxes as described in Laurent et al. !
79! (2017). !
80! !
81! With the "PCO2AIR_DATA" option, atmospheric pCO2 uses the !
82! climatology of Laurent et al. (2017). The "PCO2AIR_SECULAR" !
83! option provides an alternative time-dependent atmospheric pCO2 !
84! evolution. If none of the 2 options are defined, atmospheric !
85! pCO2 is constant. !
86! !
87! !
88!***********************************************************************
89! References: !
90! !
91! Fennel, K., Wilkin, J., Levin, J., Moisan, J., O^Reilly, J., !
92! Haidvogel, D., 2006: Nitrogen cycling in the Mid Atlantic !
93! Bight and implications for the North Atlantic nitrogen !
94! budget: Results from a three-dimensional model. Global !
95! Biogeochemical Cycles 20, GB3007, doi:10.1029/2005GB002456. !
96! !
97! Fennel, K., Wilkin, J., Previdi, M., Najjar, R. 2008: !
98! Denitrification effects on air-sea CO2 flux in the coastal !
99! ocean: Simulations for the Northwest North Atlantic. !
100! Geophys. Res. Letters 35, L24608, doi:10.1029/2008GL036147. !
101! !
102! Fennel, K., Hu, J., Laurent, A., Marta-Almeida, M., Hetland, R. !
103! 2013: Sensitivity of Hypoxia Predictions for the Northern Gulf !
104! of Mexico to Sediment Oxygen Consumption and Model Nesting. J. !
105! Geophys. Res. Ocean 118 (2), 990-1002, doi:10.1002/jgrc.20077. !
106! !
107! Laurent, A., Fennel, K., Hu, J., Hetland, R. 2012: Simulating !
108! the Effects of Phosphorus Limitation in the Mississippi and !
109! Atchafalaya River Plumes. Biogeosciences, 9 (11), 4707-4723, !
110! doi:10.5194/bg-9-4707-2012. !
111! !
112! Yu, L., Fennel, K., Laurent, A., Murrell, M. C., Lehrter, J. C. !
113! 2015: Numerical Analysis of the Primary Processes Controlling !
114! Oxygen Dynamics on the Louisiana Shelf. Biogeosciences, 12 (7), !
115! 2063-2076, doi:10.5194/bg-12-2063-2015. !
116! !
117! Wanninkhof, R. 2014: Relationship between Wind Speed and Gas !
118! Exchange over the Ocean Revisited. Limnol. Oceanogr. Methods !
119! 12 (6), 351-362, doi:10.4319/lom.2014.12.351. !
120! !
121!=======================================================================
122!
123 implicit none
124!
125 PRIVATE
126 PUBLIC :: biology
127!
128 CONTAINS
129!
130!***********************************************************************
131 SUBROUTINE biology (ng,tile)
132!***********************************************************************
133!
134 USE mod_param
135#ifdef DIAGNOSTICS_BIO
136 USE mod_diags
137#endif
138 USE mod_forces
139 USE mod_grid
140 USE mod_ncparam
141 USE mod_ocean
142 USE mod_stepping
143!
144 implicit none
145!
146! Imported variable declarations.
147!
148 integer, intent(in) :: ng, tile
149!
150! Local variable declarations.
151!
152 character (len=*), parameter :: MyFile = &
153 & __FILE__
154!
155#include "tile.h"
156!
157! Set header file name.
158!
159#ifdef DISTRIBUTE
160 IF (lbiofile(inlm)) THEN
161#else
162 IF (lbiofile(inlm).and.(tile.eq.0)) THEN
163#endif
164 lbiofile(inlm)=.false.
165 bioname(inlm)=myfile
166 END IF
167!
168#ifdef PROFILE
169 CALL wclock_on (ng, inlm, 15, __line__, myfile)
170#endif
171 CALL fennel_tile (ng, tile, &
172 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
173 & imins, imaxs, jmins, jmaxs, &
174 & nstp(ng), nnew(ng), &
175#ifdef MASKING
176 & grid(ng) % rmask, &
177# ifdef WET_DRY
178 & grid(ng) % rmask_wet, &
179# ifdef DIAGNOSTICS_BIO
180 & grid(ng) % rmask_full, &
181# endif
182# endif
183#endif
184 & grid(ng) % Hz, &
185 & grid(ng) % z_r, &
186 & grid(ng) % z_w, &
187 & forces(ng) % srflx, &
188#if defined CARBON || defined OXYGEN
189# ifdef BULK_FLUXES
190 & forces(ng) % Uwind, &
191 & forces(ng) % Vwind, &
192# else
193 & forces(ng) % sustr, &
194 & forces(ng) % svstr, &
195# endif
196#endif
197#ifdef CARBON
198 & ocean(ng) % pH, &
199#endif
200#ifdef DIAGNOSTICS_BIO
201 & diags(ng) % DiaBio2d, &
202 & diags(ng) % DiaBio3d, &
203#endif
204 & ocean(ng) % t)
205
206#ifdef PROFILE
207 CALL wclock_off (ng, inlm, 15, __line__, myfile)
208#endif
209!
210 RETURN
211 END SUBROUTINE biology
212!
213!-----------------------------------------------------------------------
214 SUBROUTINE fennel_tile (ng, tile, &
215 & LBi, UBi, LBj, UBj, UBk, UBt, &
216 & IminS, ImaxS, JminS, JmaxS, &
217 & nstp, nnew, &
218#ifdef MASKING
219 & rmask, &
220# if defined WET_DRY
221 & rmask_wet, &
222# ifdef DIAGNOSTICS_BIO
223 & rmask_full, &
224# endif
225# endif
226#endif
227 & Hz, z_r, z_w, srflx, &
228#if defined CARBON || defined OXYGEN
229# ifdef BULK_FLUXES
230 & Uwind, Vwind, &
231# else
232 & sustr, svstr, &
233# endif
234#endif
235#ifdef CARBON
236 & pH, &
237#endif
238#ifdef DIAGNOSTICS_BIO
239 & DiaBio2d, DiaBio3d, &
240#endif
241 & t)
242!-----------------------------------------------------------------------
243!
244 USE mod_param
245 USE mod_biology
246 USE mod_ncparam
247 USE mod_scalars
248!
249 USE dateclock_mod, ONLY : caldate
250!
251! Imported variable declarations.
252!
253 integer, intent(in) :: ng, tile
254 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
255 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
256 integer, intent(in) :: nstp, nnew
257
258#ifdef ASSUMED_SHAPE
259# ifdef MASKING
260 real(r8), intent(in) :: rmask(LBi:,LBj:)
261# ifdef WET_DRY
262 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
263# ifdef DIAGNOSTICS_BIO
264 real(r8), intent(in) :: rmask_full(LBi:,LBj:)
265# endif
266# endif
267# endif
268 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
269 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
270 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
271 real(r8), intent(in) :: srflx(LBi:,LBj:)
272# if defined CARBON || defined OXYGEN
273# ifdef BULK_FLUXES
274 real(r8), intent(in) :: Uwind(LBi:,LBj:)
275 real(r8), intent(in) :: Vwind(LBi:,LBj:)
276# else
277 real(r8), intent(in) :: sustr(LBi:,LBj:)
278 real(r8), intent(in) :: svstr(LBi:,LBj:)
279# endif
280# endif
281# ifdef CARBON
282 real(r8), intent(inout) :: pH(LBi:,LBj:)
283# endif
284# ifdef DIAGNOSTICS_BIO
285 real(r8), intent(inout) :: DiaBio2d(LBi:,LBj:,:)
286 real(r8), intent(inout) :: DiaBio3d(LBi:,LBj:,:,:)
287# endif
288 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
289#else
290# ifdef MASKING
291 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
292# ifdef WET_DRY
293 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
294# ifdef DIAGNOSTICS_BIO
295 real(r8), intent(in) :: rmask_full(LBi:UBi,LBj:UBj)
296# endif
297# endif
298# endif
299 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
300 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
301 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
302 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
303# if defined CARBON || defined OXYGEN
304# ifdef BULK_FLUXES
305 real(r8), intent(in) :: Uwind(LBi:UBi,LBj:UBj)
306 real(r8), intent(in) :: Vwind(LBi:UBi,LBj:UBj)
307# else
308 real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
309 real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
310# endif
311# endif
312# ifdef CARBON
313 real(r8), intent(inout) :: pH(LBi:UBi,LBj:UBj)
314# endif
315# ifdef DIAGNOSTICS_BIO
316 real(r8), intent(inout) :: DiaBio2d(LBi:UBi,LBj:UBj,NDbio2d)
317 real(r8), intent(inout) :: DiaBio3d(LBi:UBi,LBj:UBj,UBk,NDbio3d)
318# endif
319 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
320#endif
321!
322! Local variable declarations.
323!
324#ifdef CARBON
325 integer, parameter :: Nsink = 6
326#else
327 integer, parameter :: Nsink = 4
328#endif
329
330 integer :: Iter, i, ibio, isink, itrc, ivar, j, k, ks
331
332 integer, dimension(Nsink) :: idsink
333
334 real(r8), parameter :: eps = 1.0e-20_r8
335
336#if defined CARBON || defined OXYGEN
337 real(r8) :: u10squ
338#endif
339#ifdef OXYGEN
340# if defined OCMIP_OXYGEN_SC
341!
342! Alternative formulation for Schmidt number coefficients (Sc will be
343! slightly smaller up to about 35C) using the formulation proposed by
344! Keeling et al. (1998, Global Biogeochem. Cycles, 12, 141-163).
345!
346 real(r8), parameter :: A_O2 = 1638.0_r8
347 real(r8), parameter :: B_O2 = 81.83_r8
348 real(r8), parameter :: C_O2 = 1.483_r8
349 real(r8), parameter :: D_O2 = 0.008004_r8
350 real(r8), parameter :: E_O2 = 0.0_r8
351
352# elif defined RW14_OXYGEN_SC
353!
354! Alternative formulation for Schmidt number coefficients using the
355! formulation of Wanninkhof (2014, L and O Methods, 12,351-362).
356!
357 real(r8), parameter :: A_O2 = 1920.4_r8
358 real(r8), parameter :: B_O2 = 135.6_r8
359 real(r8), parameter :: C_O2 = 5.2122_r8
360 real(r8), parameter :: D_O2 = 0.10939_r8
361 real(r8), parameter :: E_O2 = 0.00093777_r8
362
363# else
364!
365! Schmidt number coefficients using the formulation of
366! Wanninkhof (1992).
367!
368 real(r8), parameter :: A_O2 = 1953.4_r8
369 real(r8), parameter :: B_O2 = 128.0_r8
370 real(r8), parameter :: C_O2 = 3.9918_r8
371 real(r8), parameter :: D_O2 = 0.050091_r8
372 real(r8), parameter :: E_O2 = 0.0_r8
373#endif
374 real(r8), parameter :: OA0 = 2.00907_r8 ! Oxygen
375 real(r8), parameter :: OA1 = 3.22014_r8 ! saturation
376 real(r8), parameter :: OA2 = 4.05010_r8 ! coefficients
377 real(r8), parameter :: OA3 = 4.94457_r8
378 real(r8), parameter :: OA4 =-0.256847_r8
379 real(r8), parameter :: OA5 = 3.88767_r8
380 real(r8), parameter :: OB0 =-0.00624523_r8
381 real(r8), parameter :: OB1 =-0.00737614_r8
382 real(r8), parameter :: OB2 =-0.0103410_r8
383 real(r8), parameter :: OB3 =-0.00817083_r8
384 real(r8), parameter :: OC0 =-0.000000488682_r8
385 real(r8), parameter :: rOxNO3= 8.625_r8 ! 138/16
386 real(r8), parameter :: rOxNH4= 6.625_r8 ! 106/16
387 real(r8) :: l2mol = 1000.0_r8/22.3916_r8 ! liter to mol
388#endif
389#ifdef CARBON
390 integer :: year
391 integer, parameter :: DoNewton = 0 ! pCO2 solver
392
393# if defined RW14_CO2_SC
394 real(r8), parameter :: A_CO2 = 2116.8_r8 ! Schmidt number
395 real(r8), parameter :: B_CO2 = 136.25_r8 ! transfer coeff
396 real(r8), parameter :: C_CO2 = 4.7353_r8 ! according to
397 real(r8), parameter :: D_CO2 = 0.092307_r8 ! Wanninkhof (2014)
398 real(r8), parameter :: E_CO2 = 0.0007555_r8
399# else
400 real(r8), parameter :: A_CO2 = 2073.1_r8 ! Schmidt
401 real(r8), parameter :: B_CO2 = 125.62_r8 ! number
402 real(r8), parameter :: C_CO2 = 3.6276_r8 ! transfer
403 real(r8), parameter :: D_CO2 = 0.043219_r8 ! coefficients
404 real(r8), parameter :: E_CO2 = 0.0_r8
405#endif
406
407 real(r8), parameter :: A1 = -60.2409_r8 ! surface
408 real(r8), parameter :: A2 = 93.4517_r8 ! CO2
409 real(r8), parameter :: A3 = 23.3585_r8 ! solubility
410 real(r8), parameter :: B1 = 0.023517_r8 ! coefficients
411 real(r8), parameter :: B2 = -0.023656_r8
412 real(r8), parameter :: B3 = 0.0047036_r8
413
414 real(r8) :: pmonth ! months since Jan 1951
415 real(r8) :: pCO2air_secular
416 real(dp) :: yday
417
418 real(r8), parameter :: pi2 = 6.2831853071796_r8
419
420 real(r8), parameter :: D0 = 282.6_r8 ! coefficients
421 real(r8), parameter :: D1 = 0.125_r8 ! to calculate
422 real(r8), parameter :: D2 =-7.18_r8 ! secular trend in
423 real(r8), parameter :: D3 = 0.86_r8 ! atmospheric pCO2
424 real(r8), parameter :: D4 =-0.99_r8
425 real(r8), parameter :: D5 = 0.28_r8
426 real(r8), parameter :: D6 =-0.80_r8
427 real(r8), parameter :: D7 = 0.06_r8
428#endif
429
430 real(r8) :: Att, AttFac, ExpAtt, Itop, PAR
431 real(r8) :: Epp, L_NH4, L_NO3, LTOT, Vp
432#ifdef PO4
433 real(r8), parameter :: MinVal = 1.0e-6_r8
434
435 real(r8) :: L_PO4, LMIN, mu, cff6
436#endif
437 real(r8) :: Chl2C, dtdays, t_PPmax, inhNH4
438
439 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5
440#ifdef RIVER_DON
441 real(r8) :: cff7, cff8
442#endif
443 real(r8) :: fac1, fac2, fac3
444 real(r8) :: cffL, cffR, cu, dltL, dltR
445
446 real(r8) :: total_N
447
448#ifdef DIAGNOSTICS_BIO
449 real(r8) :: fiter
450#endif
451
452#ifdef OXYGEN
453 real(r8) :: SchmidtN_Ox, O2satu, O2_Flux
454 real(r8) :: TS, AA
455#endif
456
457#ifdef CARBON
458 real(r8) :: C_Flux_RemineL, C_Flux_RemineS, C_Flux_RemineR
459 real(r8) :: CO2_Flux, CO2_sol, SchmidtN, TempK
460#endif
461
462 real(r8) :: N_Flux_Assim
463 real(r8) :: N_Flux_CoagD, N_Flux_CoagP
464 real(r8) :: N_Flux_Egest
465 real(r8) :: N_Flux_NewProd, N_Flux_RegProd
466 real(r8) :: N_Flux_Nitrifi
467 real(r8) :: N_Flux_Pmortal, N_Flux_Zmortal
468 real(r8) :: N_Flux_RemineL, N_Flux_RemineS, N_Flux_RemineR
469 real(r8) :: N_Flux_Zexcret, N_Flux_Zmetabo
470
471 real(r8), dimension(Nsink) :: Wbio
472
473 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
474
475 real(r8), dimension(IminS:ImaxS) :: PARsur
476#ifdef CARBON
477 real(r8), dimension(IminS:ImaxS) :: pCO2
478#endif
479
480 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
481 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
482
483 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
484
485 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
486 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
487 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
488 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
489 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
490 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
491 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
492 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
493
494#include "set_bounds.h"
495#ifdef DIAGNOSTICS_BIO
496!
497!-----------------------------------------------------------------------
498! If appropriate, initialize time-averaged diagnostic arrays.
499!-----------------------------------------------------------------------
500!
501 IF (((iic(ng).gt.ntsdia(ng)).and. &
502 & (mod(iic(ng),ndia(ng)).eq.1)).or. &
503 & ((iic(ng).ge.ntsdia(ng)).and.(ndia(ng).eq.1)).or. &
504 & ((nrrec(ng).gt.0).and.(iic(ng).eq.ntstart(ng)))) THEN
505 DO ivar=1,ndbio2d
506 DO j=jstr,jend
507 DO i=istr,iend
508 diabio2d(i,j,ivar)=0.0_r8
509 END DO
510 END DO
511 END DO
512 DO ivar=1,ndbio3d
513 DO k=1,n(ng)
514 DO j=jstr,jend
515 DO i=istr,iend
516 diabio3d(i,j,k,ivar)=0.0_r8
517 END DO
518 END DO
519 END DO
520 END DO
521 END IF
522#endif
523!
524!-----------------------------------------------------------------------
525! Add biological Source/Sink terms.
526!-----------------------------------------------------------------------
527!
528! Avoid computing source/sink terms if no biological iterations.
529!
530 IF (bioiter(ng).le.0) RETURN
531!
532! Set time-stepping according to the number of iterations.
533!
534 dtdays=dt(ng)*sec2day/real(bioiter(ng),r8)
535#ifdef DIAGNOSTICS_BIO
536!
537! A factor to account for the number of iterations in accumulating
538! diagnostic rate variables.
539!
540 fiter=1.0_r8/real(bioiter(ng),r8)
541#endif
542!
543! Set vertical sinking indentification vector.
544!
545 idsink(1)=iphyt
546 idsink(2)=ichlo
547 idsink(3)=isden
548 idsink(4)=ilden
549#ifdef CARBON
550 idsink(5)=isdec
551 idsink(6)=ildec
552#endif
553!
554! Set vertical sinking velocity vector in the same order as the
555! identification vector, IDSINK.
556!
557 wbio(1)=wphy(ng) ! phytoplankton
558 wbio(2)=wphy(ng) ! chlorophyll
559 wbio(3)=wsdet(ng) ! small Nitrogen-detritus
560 wbio(4)=wldet(ng) ! large Nitrogen-detritus
561#ifdef CARBON
562 wbio(5)=wsdet(ng) ! small Carbon-detritus
563 wbio(6)=wldet(ng) ! large Carbon-detritus
564#endif
565!
566! Compute inverse thickness to avoid repeated divisions.
567!
568 j_loop : DO j=jstr,jend
569 DO k=1,n(ng)
570 DO i=istr,iend
571 hz_inv(i,k)=1.0_r8/hz(i,j,k)
572 END DO
573 END DO
574 DO k=1,n(ng)-1
575 DO i=istr,iend
576 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
577 END DO
578 END DO
579 DO k=2,n(ng)-1
580 DO i=istr,iend
581 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
582 END DO
583 END DO
584!
585! Extract biological variables from tracer arrays, place them into
586! scratch arrays, and restrict their values to be positive definite.
587! At input, all tracers (index nnew) from predictor step have
588! transport units (m Tunits) since we do not have yet the new
589! values for zeta and Hz. These are known after the 2D barotropic
590! time-stepping.
591!
592 DO itrc=1,nbt
593 ibio=idbio(itrc)
594 DO k=1,n(ng)
595 DO i=istr,iend
596 bio_old(i,k,ibio)=max(0.0_r8,t(i,j,k,nstp,ibio))
597 bio(i,k,ibio)=bio_old(i,k,ibio)
598 END DO
599 END DO
600 END DO
601#ifdef CARBON
602 DO k=1,n(ng)
603 DO i=istr,iend
604 bio_old(i,k,itic_)=min(bio_old(i,k,itic_),3000.0_r8)
605 bio_old(i,k,itic_)=max(bio_old(i,k,itic_),400.0_r8)
606 bio(i,k,itic_)=bio_old(i,k,itic_)
607 END DO
608 END DO
609#endif
610!
611! Extract potential temperature and salinity.
612!
613 DO k=1,n(ng)
614 DO i=istr,iend
615 bio(i,k,itemp)=min(t(i,j,k,nstp,itemp),35.0_r8)
616 bio(i,k,isalt)=max(t(i,j,k,nstp,isalt), 0.0_r8)
617 END DO
618 END DO
619!
620! Calculate surface Photosynthetically Available Radiation (PAR). The
621! net shortwave radiation is scaled back to Watts/m2 and multiplied by
622! the fraction that is photosynthetically available, PARfrac.
623!
624 DO i=istr,iend
625 parsur(i)=parfrac(ng)*srflx(i,j)*rho0*cp
626 END DO
627!
628!=======================================================================
629! Start internal iterations to achieve convergence of the nonlinear
630! backward-implicit solution.
631!=======================================================================
632!
633! During the iterative procedure a series of fractional time steps are
634! performed in a chained mode (splitting by different biological
635! conversion processes) in sequence of the main food chain. In all
636! stages the concentration of the component being consumed is treated
637! in fully implicit manner, so the algorithm guarantees non-negative
638! values, no matter how strong s the concentration of active consuming
639! component (Phytoplankton or Zooplankton). The overall algorithm,
640! as well as any stage of it, is formulated in conservative form
641! (except explicit sinking) in sense that the sum of concentration of
642! all components is conserved.
643!
644!
645! In the implicit algorithm, we have for example (N: nitrate,
646! P: phytoplankton),
647!
648! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
649! {Michaelis-Menten}
650! below, we set
651! The N in the numerator of
652! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
653! as N(new)
654!
655! so the time-stepping of the equations becomes:
656!
657! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
658! consuming, divide by (1 + cff)
659! and
660!
661! P(new) = P(old) + cff * N(new) (2) when adding a source term,
662! growing, add (cff * source)
663!
664! Notice that if you substitute (1) in (2), you will get:
665!
666! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
667!
668! If you add (1) and (3), you get
669!
670! N(new) + P(new) = N(old) + P(old)
671!
672! implying conservation regardless how "cff" is computed. Therefore,
673! this scheme is unconditionally stable regardless of the conversion
674! rate. It does not generate negative values since the constituent
675! to be consumed is always treated implicitly. It is also biased
676! toward damping oscillations.
677!
678! The iterative loop below is to iterate toward an universal Backward-
679! Euler treatment of all terms. So if there are oscillations in the
680! system, they are only physical oscillations. These iterations,
681! however, do not improve the accuaracy of the solution.
682!
683 iter_loop: DO iter=1,bioiter(ng)
684!
685!-----------------------------------------------------------------------
686! Light-limited computations.
687!-----------------------------------------------------------------------
688!
689! Compute attenuation coefficient based on the concentration of
690! chlorophyll-a within each grid box. Then, attenuate surface
691! photosynthetically available radiation (PARsur) down inot the
692! water column. Thus, PAR at certain depth depends on the whole
693! distribution of chlorophyll-a above.
694! To compute rate of maximum primary productivity (t_PPmax), one needs
695! PAR somewhat in the middle of the gridbox, so that attenuation "Att"
696! corresponds to half of the grid box height, while PAR is multiplied
697! by it twice: once to get it in the middle of grid-box and once the
698! compute on the lower grid-box interface.
699!
700 DO i=istr,iend
701 par=parsur(i)
702 attfac=0.0_r8
703 IF (parsur(i).gt.0.0_r8) THEN
704 DO k=n(ng),1,-1
705!
706! Compute average light attenuation for each grid cell. To include
707! other attenuation contributions like suspended sediment or CDOM
708! modify AttFac.
709!
710 att=(attsw(ng)+ &
711 & attchl(ng)*bio(i,k,ichlo)+ &
712 & attfac)* &
713 & (z_w(i,j,k)-z_w(i,j,k-1))
714 expatt=exp(-att)
715 itop=par
716 par=itop*(1.0_r8-expatt)/att ! average at cell center
717!
718! Compute Chlorophyll-a phytoplankton ratio, [mg Chla / (mg C)].
719!
720 cff=phycn(ng)*12.0_r8
721 chl2c=min(bio(i,k,ichlo)/(bio(i,k,iphyt)*cff+eps), &
722 & chl2c_m(ng))
723!
724! Temperature-limited and light-limited growth rate (Eppley, R.W.,
725! 1972, Fishery Bulletin, 70: 1063-1085; here 0.59=ln(2)*0.851).
726! Check value for Vp is 2.9124317 at 19.25 degC.
727!
728 vp=vp0(ng)*0.59_r8*(1.066_r8**bio(i,k,itemp))
729 fac1=par*phyis(ng)
730 epp=vp/sqrt(vp*vp+fac1*fac1)
731 t_ppmax=epp*fac1
732!
733! Nutrient-limitation terms (Parker 1993 Ecol Mod., 66, 113-120).
734!
735 cff1=bio(i,k,inh4_)*k_nh4(ng)
736 cff2=bio(i,k,ino3_)*k_no3(ng)
737 inhnh4=1.0_r8/(1.0_r8+cff1)
738 l_nh4=cff1/(1.0_r8+cff1)
739 l_no3=cff2*inhnh4/(1.0_r8+cff2)
740 ltot=l_no3+l_nh4
741#ifdef PO4
742 cff3=bio(i,k,ipo4_)*k_po4(ng)
743 l_po4=cff3/(1.0_r8+cff3)
744 lmin=min(ltot,l_po4)
745#endif
746!
747! Nitrate and ammonium uptake by Phytoplankton.
748!
749#ifdef PO4
750 mu=dtdays*t_ppmax*lmin
751 cff4=mu*bio(i,k,iphyt)*l_no3/ &
752 & max(minval,ltot)/max(minval,bio(i,k,ino3_))
753 cff5=mu*bio(i,k,iphyt)*l_nh4/ &
754 & max(minval,ltot)/max(minval,bio(i,k,inh4_))
755 cff6=r_p2n(ng)*mu*bio(i,k,iphyt)/ &
756 & max(minval,bio(i,k,ipo4_))
757#else
758 fac1=dtdays*t_ppmax
759 cff4=fac1*k_no3(ng)*inhnh4/(1.0_r8+cff2)*bio(i,k,iphyt)
760 cff5=fac1*k_nh4(ng)/(1.0_r8+cff1)*bio(i,k,iphyt)
761#endif
762 bio(i,k,ino3_)=bio(i,k,ino3_)/(1.0_r8+cff4)
763 bio(i,k,inh4_)=bio(i,k,inh4_)/(1.0_r8+cff5)
764#ifdef PO4
765 bio(i,k,ipo4_)=bio(i,k,ipo4_)/(1.0_r8+cff6)
766#endif
767 n_flux_newprod=bio(i,k,ino3_)*cff4
768 n_flux_regprod=bio(i,k,inh4_)*cff5
769 bio(i,k,iphyt)=bio(i,k,iphyt)+ &
770 & n_flux_newprod+n_flux_regprod
771!
772 bio(i,k,ichlo)=bio(i,k,ichlo)+ &
773#ifdef PO4
774 & (dtdays*t_ppmax*t_ppmax*lmin*lmin* &
775#else
776 & (dtdays*t_ppmax*t_ppmax*ltot*ltot* &
777#endif
778 & chl2c_m(ng)*bio(i,k,ichlo))/ &
779 & (phyis(ng)*max(chl2c,eps)*par+eps)
780#ifdef DIAGNOSTICS_BIO
781 diabio3d(i,j,k,ippro)=diabio3d(i,j,k,ippro)+ &
782# ifdef WET_DRY
783 & rmask_full(i,j)* &
784# endif
785 & (n_flux_newprod+n_flux_regprod)* &
786 & fiter
787 diabio3d(i,j,k,ino3u)=diabio3d(i,j,k,ino3u)+ &
788# ifdef WET_DRY
789 & rmask_full(i,j)* &
790# endif
791 & n_flux_newprod*fiter
792#endif
793#ifdef OXYGEN
794 bio(i,k,ioxyg)=bio(i,k,ioxyg)+ &
795 & n_flux_newprod*roxno3+ &
796 & n_flux_regprod*roxnh4
797#endif
798#ifdef CARBON
799!
800! Total inorganic carbon (CO2) uptake during phytoplankton growth.
801!
802 cff1=phycn(ng)*(n_flux_newprod+n_flux_regprod)
803 bio(i,k,itic_)=bio(i,k,itic_)-cff1
804# ifdef TALK_NONCONSERV
805!
806! Account for the uptake of NO3 on total alkalinity.
807!
808 bio(i,k,italk)=bio(i,k,italk)+n_flux_newprod- &
809 & n_flux_regprod
810# endif
811#endif
812!
813! The Nitrification of NH4 ==> NO3 is thought to occur only in dark and
814! only in aerobic water (see Olson, R. J., 1981, JMR: (39), 227-238.).
815!
816! NH4+ + 3/2 O2 ==> NO2- + H2O; via Nitrosomonas bacteria
817! NO2- + 1/2 O2 ==> NO3- ; via Nitrobacter bacteria
818!
819! Note that the entire process has a total loss of two moles of O2 per
820! mole of NH4. If we were to resolve NO2 profiles, this is where we
821! would change the code to split out the differential effects of the
822! two different bacteria types. If OXYGEN is defined, nitrification is
823! inhibited at low oxygen concentrations using a Michaelis-Menten term.
824!
825#ifdef OXYGEN
826 fac2=max(bio(i,k,ioxyg),0.0_r8) ! O2 max
827 fac3=max(fac2/(3.0_r8+fac2),0.0_r8) ! MM for O2 dependence
828 fac1=dtdays*nitrir(ng)*fac3
829#else
830 fac1=dtdays*nitrir(ng)
831#endif
832 cff1=(par-i_thnh4(ng))/ &
833 & (d_p5nh4(ng)+par-2.0_r8*i_thnh4(ng))
834 cff2=1.0_r8-max(0.0_r8,cff1)
835 cff3=fac1*cff2
836 bio(i,k,inh4_)=bio(i,k,inh4_)/(1.0_r8+cff3)
837 n_flux_nitrifi=bio(i,k,inh4_)*cff3
838 bio(i,k,ino3_)=bio(i,k,ino3_)+n_flux_nitrifi
839#ifdef DIAGNOSTICS_BIO
840 diabio3d(i,j,k,inifx)=diabio3d(i,j,k,inifx)+ &
841# ifdef WET_DRY
842 & rmask_full(i,j)* &
843# endif
844 & n_flux_nitrifi*fiter
845#endif
846#ifdef OXYGEN
847 bio(i,k,ioxyg)=bio(i,k,ioxyg)-2.0_r8*n_flux_nitrifi
848#endif
849#if defined CARBON && defined TALK_NONCONSERV
850 bio(i,k,italk)=bio(i,k,italk)-2.0_r8*n_flux_nitrifi
851#endif
852!
853! Light attenuation at the bottom of the grid cell. It is the starting
854! PAR value for the next (deeper) vertical grid cell.
855!
856 par=itop*expatt
857 END DO
858!
859! If PARsur=0, nitrification occurs at the maximum rate (NitriR).
860!
861 ELSE
862 cff3=dtdays*nitrir(ng)
863 DO k=n(ng),1,-1
864 bio(i,k,inh4_)=bio(i,k,inh4_)/(1.0_r8+cff3)
865 n_flux_nitrifi=bio(i,k,inh4_)*cff3
866 bio(i,k,ino3_)=bio(i,k,ino3_)+n_flux_nitrifi
867#ifdef DIAGNOSTICS_BIO
868 diabio3d(i,j,k,inifx)=diabio3d(i,j,k,inifx)+ &
869# ifdef WET_DRY
870 & rmask_full(i,j)* &
871# endif
872 & n_flux_nitrifi*fiter
873#endif
874#ifdef OXYGEN
875 bio(i,k,ioxyg)=bio(i,k,ioxyg)-2.0_r8*n_flux_nitrifi
876#endif
877#if defined CARBON && defined TALK_NONCONSERV
878 bio(i,k,italk)=bio(i,k,italk)-2.0_r8*n_flux_nitrifi
879#endif
880 END DO
881 END IF
882 END DO
883!
884!-----------------------------------------------------------------------
885! Phytoplankton grazing by zooplankton (rate: ZooGR), phytoplankton
886! assimilated to zooplankton (fraction: ZooAE_N) and egested to small
887! detritus, and phytoplankton mortality (rate: PhyMR) to small
888! detritus. [Landry 1993 L and O 38:468-472]
889!-----------------------------------------------------------------------
890!
891 fac1=dtdays*zoogr(ng)
892 cff2=dtdays*phymr(ng)
893 DO k=1,n(ng)
894 DO i=istr,iend
895!
896! Phytoplankton grazing by zooplankton.
897!
898 cff1=fac1*bio(i,k,izoop)*bio(i,k,iphyt)/ &
899 & (k_phy(ng)+bio(i,k,iphyt)*bio(i,k,iphyt))
900 cff3=1.0_r8/(1.0_r8+cff1)
901 bio(i,k,iphyt)=cff3*bio(i,k,iphyt)
902 bio(i,k,ichlo)=cff3*bio(i,k,ichlo)
903!
904! Phytoplankton assimilated to zooplankton and egested to small
905! detritus.
906!
907 n_flux_assim=cff1*bio(i,k,iphyt)*zooae_n(ng)
908 n_flux_egest=bio(i,k,iphyt)*cff1*(1.0_r8-zooae_n(ng))
909 bio(i,k,izoop)=bio(i,k,izoop)+ &
910 & n_flux_assim
911 bio(i,k,isden)=bio(i,k,isden)+ &
912 & n_flux_egest
913!
914! Phytoplankton mortality (limited by a phytoplankton minimum).
915!
916 n_flux_pmortal=cff2*max(bio(i,k,iphyt)-phymin(ng),0.0_r8)
917 bio(i,k,iphyt)=bio(i,k,iphyt)-n_flux_pmortal
918 bio(i,k,ichlo)=bio(i,k,ichlo)- &
919 & cff2*max(bio(i,k,ichlo)-chlmin(ng),0.0_r8)
920 bio(i,k,isden)=bio(i,k,isden)+ &
921 & n_flux_pmortal
922#ifdef CARBON
923 bio(i,k,isdec)=bio(i,k,isdec)+ &
924 & phycn(ng)*(n_flux_egest+n_flux_pmortal)+ &
925 & (phycn(ng)-zoocn(ng))*n_flux_assim
926#endif
927 END DO
928 END DO
929!
930!-----------------------------------------------------------------------
931! Zooplankton basal metabolism to NH4 (rate: ZooBM), zooplankton
932! mortality to small detritus (rate: ZooMR), zooplankton ingestion
933! related excretion (rate: ZooER).
934!-----------------------------------------------------------------------
935!
936 cff1=dtdays*zoobm(ng)
937 fac2=dtdays*zoomr(ng)
938 fac3=dtdays*zooer(ng)
939 DO k=1,n(ng)
940 DO i=istr,iend
941 fac1=fac3*bio(i,k,iphyt)*bio(i,k,iphyt)/ &
942 & (k_phy(ng)+bio(i,k,iphyt)*bio(i,k,iphyt))
943 cff2=fac2*bio(i,k,izoop)
944 cff3=fac1*zooae_n(ng)
945 bio(i,k,izoop)=bio(i,k,izoop)/ &
946 & (1.0_r8+cff2+cff3)
947!
948! Zooplankton mortality and excretion.
949!
950 n_flux_zmortal=cff2*bio(i,k,izoop)
951 n_flux_zexcret=cff3*bio(i,k,izoop)
952 bio(i,k,inh4_)=bio(i,k,inh4_)+n_flux_zexcret
953#ifdef PO4
954 bio(i,k,ipo4_)=bio(i,k,ipo4_)+r_p2n(ng)*n_flux_zexcret
955#endif
956 bio(i,k,isden)=bio(i,k,isden)+n_flux_zmortal
957!
958! Zooplankton basal metabolism (limited by a zooplankton minimum).
959!
960 n_flux_zmetabo=cff1*max(bio(i,k,izoop)-zoomin(ng),0.0_r8)
961 bio(i,k,izoop)=bio(i,k,izoop)-n_flux_zmetabo
962 bio(i,k,inh4_)=bio(i,k,inh4_)+n_flux_zmetabo
963#ifdef PO4
964 bio(i,k,ipo4_)=bio(i,k,ipo4_)+r_p2n(ng)*n_flux_zmetabo
965#endif
966#ifdef OXYGEN
967 bio(i,k,ioxyg)=bio(i,k,ioxyg)- &
968 & roxnh4*(n_flux_zmetabo+n_flux_zexcret)
969#endif
970#ifdef CARBON
971 bio(i,k,isdec)=bio(i,k,isdec)+ &
972 & zoocn(ng)*n_flux_zmortal
973 bio(i,k,itic_)=bio(i,k,itic_)+ &
974 & zoocn(ng)*(n_flux_zmetabo+n_flux_zexcret)
975#ifdef TALK_NONCONSERV
976 bio(i,k,italk)=bio(i,k,italk)+n_flux_zmetabo+ &
977 & n_flux_zexcret
978#endif
979#endif
980 END DO
981 END DO
982!
983!-----------------------------------------------------------------------
984! Coagulation of phytoplankton and small detritus to large detritus.
985!-----------------------------------------------------------------------
986!
987 fac1=dtdays*coagr(ng)
988 DO k=1,n(ng)
989 DO i=istr,iend
990 cff1=fac1*(bio(i,k,isden)+bio(i,k,iphyt))
991 cff2=1.0_r8/(1.0_r8+cff1)
992 bio(i,k,iphyt)=bio(i,k,iphyt)*cff2
993 bio(i,k,ichlo)=bio(i,k,ichlo)*cff2
994 bio(i,k,isden)=bio(i,k,isden)*cff2
995 n_flux_coagp=bio(i,k,iphyt)*cff1
996 n_flux_coagd=bio(i,k,isden)*cff1
997 bio(i,k,ilden)=bio(i,k,ilden)+ &
998 & n_flux_coagp+n_flux_coagd
999#ifdef CARBON
1000 bio(i,k,isdec)=bio(i,k,isdec)-phycn(ng)*n_flux_coagd
1001 bio(i,k,ildec)=bio(i,k,ildec)+ &
1002 & phycn(ng)*(n_flux_coagp+n_flux_coagd)
1003#endif
1004 END DO
1005 END DO
1006!
1007!-----------------------------------------------------------------------
1008! Detritus recycling to NH4, remineralization.
1009!-----------------------------------------------------------------------
1010!
1011#ifdef OXYGEN
1012 DO k=1,n(ng)
1013 DO i=istr,iend
1014 fac1=max(bio(i,k,ioxyg)-6.0_r8,0.0_r8) ! O2 off max
1015 fac2=max(fac1/(3.0_r8+fac1),0.0_r8) ! MM for O2 dependence
1016 cff1=dtdays*sderrn(ng)*fac2
1017 cff2=1.0_r8/(1.0_r8+cff1)
1018 cff3=dtdays*lderrn(ng)*fac2
1019 cff4=1.0_r8/(1.0_r8+cff3)
1020 bio(i,k,isden)=bio(i,k,isden)*cff2
1021 bio(i,k,ilden)=bio(i,k,ilden)*cff4
1022 n_flux_remines=bio(i,k,isden)*cff1
1023 n_flux_reminel=bio(i,k,ilden)*cff3
1024 bio(i,k,inh4_)=bio(i,k,inh4_)+ &
1025 & n_flux_remines+n_flux_reminel
1026# ifdef PO4
1027 bio(i,k,ipo4_)=bio(i,k,ipo4_)+r_p2n(ng) &
1028 & *(n_flux_remines+n_flux_reminel)
1029# endif
1030 bio(i,k,ioxyg)=bio(i,k,ioxyg)- &
1031 & (n_flux_remines+n_flux_reminel)*roxnh4
1032# if defined CARBON && defined TALK_NONCONSERV
1033 bio(i,k,italk)=bio(i,k,italk)+n_flux_remines+ &
1034 & n_flux_reminel
1035# endif
1036# ifdef RIVER_DON
1037 cff7=dtdays*rderrn(ng)*fac2
1038 cff8=1.0_r8/(1.0_r8+cff7)
1039 bio(i,k,irden)=bio(i,k,irden)*cff8
1040 n_flux_reminer=bio(i,k,irden)*cff7
1041 bio(i,k,inh4_)=bio(i,k,inh4_)+ &
1042 & n_flux_reminer
1043# ifdef PO4
1044 bio(i,k,ipo4_)=bio(i,k,ipo4_)+r_p2n(ng) &
1045 & *n_flux_reminer
1046# endif
1047 bio(i,k,ioxyg)=bio(i,k,ioxyg)-n_flux_reminer*roxnh4
1048# if defined CARBON && defined TALK_NONCONSERV
1049 bio(i,k,italk)=bio(i,k,italk)+n_flux_reminer
1050# endif
1051# endif
1052 END DO
1053 END DO
1054#else
1055 cff1=dtdays*sderrn(ng)
1056 cff2=1.0_r8/(1.0_r8+cff1)
1057 cff3=dtdays*lderrn(ng)
1058 cff4=1.0_r8/(1.0_r8+cff3)
1059# ifdef RIVER_DON
1060 cff7=dtdays*rderrn(ng)
1061 cff8=1.0_r8/(1.0_r8+cff7)
1062# endif
1063 DO k=1,n(ng)
1064 DO i=istr,iend
1065 bio(i,k,isden)=bio(i,k,isden)*cff2
1066 bio(i,k,ilden)=bio(i,k,ilden)*cff4
1067 n_flux_remines=bio(i,k,isden)*cff1
1068 n_flux_reminel=bio(i,k,ilden)*cff3
1069 bio(i,k,inh4_)=bio(i,k,inh4_)+ &
1070 & n_flux_remines+n_flux_reminel
1071# ifdef PO4
1072 bio(i,k,ipo4_)=bio(i,k,ipo4_)+r_p2n(ng) &
1073 & *(n_flux_remines+n_flux_reminel)
1074# endif
1075# if defined CARBON && defined TALK_NONCONSERV
1076 bio(i,k,italk)=bio(i,k,italk)+n_flux_remines+ &
1077 & n_flux_reminel
1078# endif
1079# ifdef RIVER_DON
1080 bio(i,k,irden)=bio(i,k,irden)*cff8
1081 n_flux_reminer=bio(i,k,irden)*cff7
1082 bio(i,k,inh4_)=bio(i,k,inh4_)+n_flux_reminer
1083# ifdef PO4
1084 bio(i,k,ipo4_)=bio(i,k,ipo4_)+r_p2n(ng) &
1085 & *n_flux_reminer
1086# endif
1087# if defined CARBON && defined TALK_NONCONSERV
1088 bio(i,k,italk)=bio(i,k,italk)+n_flux_reminer
1089# endif
1090# endif
1091 END DO
1092 END DO
1093#endif
1094#ifdef OXYGEN
1095!
1096!-----------------------------------------------------------------------
1097! Surface O2 gas exchange.
1098!-----------------------------------------------------------------------
1099!
1100! Compute surface O2 gas exchange.
1101!
1102 cff1=rho0*550.0_r8
1103# if defined RW14_OXYGEN_SC
1104 cff2=dtdays*0.251_r8*24.0_r8/100.0_r8
1105# else
1106 cff2=dtdays*0.31_r8*24.0_r8/100.0_r8
1107# endif
1108 k=n(ng)
1109 DO i=istr,iend
1110!
1111! Compute O2 transfer velocity : u10squared (u10 in m/s)
1112!
1113# ifdef BULK_FLUXES
1114 u10squ=uwind(i,j)*uwind(i,j)+vwind(i,j)*vwind(i,j)
1115# else
1116 u10squ=cff1*sqrt((0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2+ &
1117 & (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2)
1118# endif
1119 schmidtn_ox=a_o2-bio(i,k,itemp)*(b_o2-bio(i,k,itemp)*(c_o2- &
1120 & bio(i,k,itemp)*(d_o2- &
1121 & bio(i,k,itemp)*e_o2)))
1122 cff3=cff2*u10squ*sqrt(660.0_r8/schmidtn_ox)
1123!
1124! Calculate O2 saturation concentration using Garcia and Gordon
1125! L and O (1992) formula, (EXP(AA) is in ml/l).
1126!
1127 ts=log((298.15_r8-bio(i,k,itemp))/ &
1128 & (273.15_r8+bio(i,k,itemp)))
1129 aa=oa0+ts*(oa1+ts*(oa2+ts*(oa3+ts*(oa4+ts*oa5))))+ &
1130 & bio(i,k,isalt)*(ob0+ts*(ob1+ts*(ob2+ts*ob3)))+ &
1131 & oc0*bio(i,k,isalt)*bio(i,k,isalt)
1132!
1133! Convert from ml/l to mmol/m3.
1134!
1135 o2satu=l2mol*exp(aa)
1136!
1137! Add in O2 gas exchange.
1138!
1139 o2_flux=cff3*(o2satu-bio(i,k,ioxyg))
1140 bio(i,k,ioxyg)=bio(i,k,ioxyg)+ &
1141 & o2_flux*hz_inv(i,k)
1142# ifdef DIAGNOSTICS_BIO
1143 diabio2d(i,j,io2fx)=diabio2d(i,j,io2fx)+ &
1144# ifdef WET_DRY
1145 & rmask_full(i,j)* &
1146# endif
1147 & o2_flux*fiter
1148# endif
1149
1150 END DO
1151#endif
1152
1153#ifdef CARBON
1154!
1155!-----------------------------------------------------------------------
1156! Allow different remineralization rates for detrital C and detrital N.
1157!-----------------------------------------------------------------------
1158!
1159 cff1=dtdays*sderrc(ng)
1160 cff2=1.0_r8/(1.0_r8+cff1)
1161 cff3=dtdays*lderrc(ng)
1162 cff4=1.0_r8/(1.0_r8+cff3)
1163# ifdef RIVER_DON
1164 cff7=dtdays*rderrc(ng)
1165 cff8=1.0_r8/(1.0_r8+cff7)
1166# endif
1167 DO k=1,n(ng)
1168 DO i=istr,iend
1169 bio(i,k,isdec)=bio(i,k,isdec)*cff2
1170 bio(i,k,ildec)=bio(i,k,ildec)*cff4
1171 c_flux_remines=bio(i,k,isdec)*cff1
1172 c_flux_reminel=bio(i,k,ildec)*cff3
1173 bio(i,k,itic_)=bio(i,k,itic_)+ &
1174 & c_flux_remines+c_flux_reminel
1175# ifdef RIVER_DON
1176 bio(i,k,irdec)=bio(i,k,irdec)*cff8
1177 c_flux_reminer=bio(i,k,irdec)*cff7
1178 bio(i,k,itic_)=bio(i,k,itic_)+c_flux_reminer
1179# endif
1180 END DO
1181 END DO
1182# ifndef TALK_NONCONSERV
1183!
1184! Alkalinity is treated as a diagnostic variable. TAlk = f(S[PSU])
1185! following Brewer et al. (1986).
1186!
1187 DO k=1,n(ng)
1188 DO i=istr,iend
1189 bio(i,k,italk)=587.05_r8+50.56_r8*bio(i,k,isalt)
1190 END DO
1191 END DO
1192# endif
1193!
1194!-----------------------------------------------------------------------
1195! Surface CO2 gas exchange.
1196!-----------------------------------------------------------------------
1197!
1198! Compute equilibrium partial pressure inorganic carbon (ppmv) at the
1199! surface.
1200!
1201 k=n(ng)
1202# ifdef pCO2_RZ
1203 CALL pco2_water_rz (istr, iend, lbi, ubi, lbj, ubj, &
1204 & imins, imaxs, j, donewton, &
1205# ifdef MASKING
1206 & rmask, &
1207# endif
1208 & bio(imins:,k,itemp), bio(imins:,k,isalt), &
1209 & bio(imins:,k,itic_), bio(imins:,k,italk), &
1210 & ph, pco2)
1211# else
1212 CALL pco2_water (istr, iend, lbi, ubi, lbj, ubj, &
1213 & imins, imaxs, j, donewton, &
1214# ifdef MASKING
1215 & rmask, &
1216# endif
1217 & bio(imins:,k,itemp), bio(imins:,k,isalt), &
1218 & bio(imins:,k,itic_), bio(imins:,k,italk), &
1219 & 0.0_r8, 0.0_r8, ph, pco2)
1220# endif
1221!
1222! Compute surface CO2 gas exchange.
1223!
1224 cff1=rho0*550.0_r8
1225# if defined RW14_CO2_SC
1226 cff2=dtdays*0.251_r8*24.0_r8/100.0_r8
1227# else
1228 cff2=dtdays*0.31_r8*24.0_r8/100.0_r8
1229# endif
1230 DO i=istr,iend
1231!
1232! Compute CO2 transfer velocity : u10squared (u10 in m/s)
1233!
1234# ifdef BULK_FLUXES
1235 u10squ=uwind(i,j)**2+vwind(i,j)**2
1236# else
1237 u10squ=cff1*sqrt((0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2+ &
1238 & (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2)
1239# endif
1240 schmidtn=a_co2-bio(i,k,itemp)*(b_co2-bio(i,k,itemp)*(c_co2- &
1241 & bio(i,k,itemp)*(d_co2- &
1242 & bio(i,k,itemp)*e_co2)))
1243 cff3=cff2*u10squ*sqrt(660.0_r8/schmidtn)
1244!
1245! Calculate CO2 solubility [mol/(kg.atm)] using Weiss (1974) formula.
1246!
1247 tempk=0.01_r8*(bio(i,k,itemp)+273.15_r8)
1248 co2_sol=exp(a1+ &
1249 & a2/tempk+ &
1250 & a3*log(tempk)+ &
1251 & bio(i,k,isalt)*(b1+tempk*(b2+b3*tempk)))
1252!
1253! Add in CO2 gas exchange.
1254!
1255 CALL caldate (tdays(ng), yy_i=year, yd_dp=yday)
1256 pmonth=year-1951.0_r8+yday/365.0_r8
1257# if defined PCO2AIR_DATA
1258 pco2air_secular=380.464_r8+9.321_r8*sin(pi2*yday/365.25_r8+ &
1259 & 1.068_r8)
1260 co2_flux=cff3*co2_sol*(pco2air_secular-pco2(i))
1261# elif defined PCO2AIR_SECULAR
1262 pco2air_secular=d0+d1*pmonth*12.0_r8+ &
1263 & d2*sin(pi2*pmonth+d3)+ &
1264 & d4*sin(pi2*pmonth+d5)+ &
1265 & d6*sin(pi2*pmonth+d7)
1266 co2_flux=cff3*co2_sol*(pco2air_secular-pco2(i))
1267# else
1268 co2_flux=cff3*co2_sol*(pco2air(ng)-pco2(i))
1269# endif
1270 bio(i,k,itic_)=bio(i,k,itic_)+ &
1271 & co2_flux*hz_inv(i,k)
1272# ifdef DIAGNOSTICS_BIO
1273 diabio2d(i,j,icofx)=diabio2d(i,j,icofx)+ &
1274# ifdef WET_DRY
1275 & rmask_full(i,j)* &
1276# endif
1277 & co2_flux*fiter
1278 diabio2d(i,j,ipco2)=pco2(i)
1279# ifdef WET_DRY
1280 diabio2d(i,j,ipco2)=diabio2d(i,j,ipco2)*rmask_full(i,j)
1281# endif
1282# endif
1283 END DO
1284#endif
1285!
1286!-----------------------------------------------------------------------
1287! Vertical sinking terms.
1288!-----------------------------------------------------------------------
1289!
1290! Reconstruct vertical profile of selected biological constituents
1291! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
1292! grid box. Then, compute semi-Lagrangian flux due to sinking.
1293!
1294 sink_loop: DO isink=1,nsink
1295 ibio=idsink(isink)
1296!
1297! Copy concentration of biological particulates into scratch array
1298! "qc" (q-central, restrict it to be positive) which is hereafter
1299! interpreted as a set of grid-box averaged values for biogeochemical
1300! constituent concentration.
1301!
1302 DO k=1,n(ng)
1303 DO i=istr,iend
1304 qc(i,k)=bio(i,k,ibio)
1305 END DO
1306 END DO
1307!
1308 DO k=n(ng)-1,1,-1
1309 DO i=istr,iend
1310 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
1311 END DO
1312 END DO
1313 DO k=2,n(ng)-1
1314 DO i=istr,iend
1315 dltr=hz(i,j,k)*fc(i,k)
1316 dltl=hz(i,j,k)*fc(i,k-1)
1317 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
1318 cffr=cff*fc(i,k)
1319 cffl=cff*fc(i,k-1)
1320!
1321! Apply PPM monotonicity constraint to prevent oscillations within the
1322! grid box.
1323!
1324 IF ((dltr*dltl).le.0.0_r8) THEN
1325 dltr=0.0_r8
1326 dltl=0.0_r8
1327 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1328 dltr=cffl
1329 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1330 dltl=cffr
1331 END IF
1332!
1333! Compute right and left side values (bR,bL) of parabolic segments
1334! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1335!
1336! NOTE: Although each parabolic segment is monotonic within its grid
1337! box, monotonicity of the whole profile is not guaranteed,
1338! because bL(k+1)-bR(k) may still have different sign than
1339! qc(i,k+1)-qc(i,k). This possibility is excluded,
1340! after bL and bR are reconciled using WENO procedure.
1341!
1342 cff=(dltr-dltl)*hz_inv3(i,k)
1343 dltr=dltr-cff*hz(i,j,k+1)
1344 dltl=dltl+cff*hz(i,j,k-1)
1345 br(i,k)=qc(i,k)+dltr
1346 bl(i,k)=qc(i,k)-dltl
1347 wr(i,k)=(2.0_r8*dltr-dltl)**2
1348 wl(i,k)=(dltr-2.0_r8*dltl)**2
1349 END DO
1350 END DO
1351 cff=1.0e-14_r8
1352 DO k=2,n(ng)-2
1353 DO i=istr,iend
1354 dltl=max(cff,wl(i,k ))
1355 dltr=max(cff,wr(i,k+1))
1356 br(i,k)=(dltr*br(i,k)+dltl*bl(i,k+1))/(dltr+dltl)
1357 bl(i,k+1)=br(i,k)
1358 END DO
1359 END DO
1360 DO i=istr,iend
1361 fc(i,n(ng))=0.0_r8 ! NO-flux boundary condition
1362#if defined LINEAR_CONTINUATION
1363 bl(i,n(ng))=br(i,n(ng)-1)
1364 br(i,n(ng))=2.0_r8*qc(i,n(ng))-bl(i,n(ng))
1365#elif defined NEUMANN
1366 bl(i,n(ng))=br(i,n(ng)-1)
1367 br(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*bl(i,n(ng))
1368#else
1369 br(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
1370 bl(i,n(ng))=qc(i,n(ng)) ! conditions
1371 br(i,n(ng)-1)=qc(i,n(ng))
1372#endif
1373#if defined LINEAR_CONTINUATION
1374 br(i,1)=bl(i,2)
1375 bl(i,1)=2.0_r8*qc(i,1)-br(i,1)
1376#elif defined NEUMANN
1377 br(i,1)=bl(i,2)
1378 bl(i,1)=1.5_r8*qc(i,1)-0.5_r8*br(i,1)
1379#else
1380 bl(i,2)=qc(i,1) ! bottom grid boxes are
1381 br(i,1)=qc(i,1) ! re-assumed to be
1382 bl(i,1)=qc(i,1) ! piecewise constant.
1383#endif
1384 END DO
1385!
1386! Apply monotonicity constraint again, since the reconciled interfacial
1387! values may cause a non-monotonic behavior of the parabolic segments
1388! inside the grid box.
1389!
1390 DO k=1,n(ng)
1391 DO i=istr,iend
1392 dltr=br(i,k)-qc(i,k)
1393 dltl=qc(i,k)-bl(i,k)
1394 cffr=2.0_r8*dltr
1395 cffl=2.0_r8*dltl
1396 IF ((dltr*dltl).lt.0.0_r8) THEN
1397 dltr=0.0_r8
1398 dltl=0.0_r8
1399 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
1400 dltr=cffl
1401 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
1402 dltl=cffr
1403 END IF
1404 br(i,k)=qc(i,k)+dltr
1405 bl(i,k)=qc(i,k)-dltl
1406 END DO
1407 END DO
1408!
1409! After this moment reconstruction is considered complete. The next
1410! stage is to compute vertical advective fluxes, FC. It is expected
1411! that sinking may occurs relatively fast, the algorithm is designed
1412! to be free of CFL criterion, which is achieved by allowing
1413! integration bounds for semi-Lagrangian advective flux to use as
1414! many grid boxes in upstream direction as necessary.
1415!
1416! In the two code segments below, WL is the z-coordinate of the
1417! departure point for grid box interface z_w with the same indices;
1418! FC is the finite volume flux; ksource(:,k) is index of vertical
1419! grid box which contains the departure point (restricted by N(ng)).
1420! During the search: also add in content of whole grid boxes
1421! participating in FC.
1422!
1423 cff=dtdays*abs(wbio(isink))
1424 DO k=1,n(ng)
1425 DO i=istr,iend
1426 fc(i,k-1)=0.0_r8
1427 wl(i,k)=z_w(i,j,k-1)+cff
1428 wr(i,k)=hz(i,j,k)*qc(i,k)
1429 ksource(i,k)=k
1430 END DO
1431 END DO
1432 DO k=1,n(ng)
1433 DO ks=k,n(ng)-1
1434 DO i=istr,iend
1435 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
1436 ksource(i,k)=ks+1
1437 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
1438 END IF
1439 END DO
1440 END DO
1441 END DO
1442!
1443! Finalize computation of flux: add fractional part.
1444!
1445 DO k=1,n(ng)
1446 DO i=istr,iend
1447 ks=ksource(i,k)
1448 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
1449 fc(i,k-1)=fc(i,k-1)+ &
1450 & hz(i,j,ks)*cu* &
1451 & (bl(i,ks)+ &
1452 & cu*(0.5_r8*(br(i,ks)-bl(i,ks))- &
1453 & (1.5_r8-cu)* &
1454 & (br(i,ks)+bl(i,ks)- &
1455 & 2.0_r8*qc(i,ks))))
1456 END DO
1457 END DO
1458 DO k=1,n(ng)
1459 DO i=istr,iend
1460 bio(i,k,ibio)=qc(i,k)+(fc(i,k)-fc(i,k-1))*hz_inv(i,k)
1461 END DO
1462 END DO
1463
1464#ifdef BIO_SEDIMENT
1465!
1466! Particulate flux reaching the seafloor is remineralized and returned
1467! to the dissolved nitrate pool. Without this conversion, particulate
1468! material falls out of the system. This is a temporary fix to restore
1469! total nitrogen conservation. It will be replaced later by a
1470! parameterization that includes the time delay of remineralization
1471! and dissolved oxygen.
1472!
1473 cff2=4.0_r8/16.0_r8
1474# ifdef OXYGEN
1475 cff3=115.0_r8/16.0_r8
1476 cff4=106.0_r8/16.0_r8
1477# endif
1478 IF ((ibio.eq.iphyt).or. &
1479 & (ibio.eq.isden).or. &
1480 & (ibio.eq.ilden)) THEN
1481 DO i=istr,iend
1482 cff1=fc(i,0)*hz_inv(i,1)
1483# ifdef DENITRIFICATION
1484 bio(i,1,inh4_)=bio(i,1,inh4_)+cff1*cff2
1485# ifdef DIAGNOSTICS_BIO
1486 diabio2d(i,j,idnit)=diabio2d(i,j,idnit)+ &
1487# ifdef WET_DRY
1488 & rmask_full(i,j)* &
1489# endif
1490 & (1.0_r8-cff2)*cff1*hz(i,j,1)*fiter
1491# endif
1492# ifdef PO4
1493 bio(i,1,ipo4_)=bio(i,1,ipo4_)+cff1*r_p2n(ng)
1494# endif
1495# ifdef OXYGEN
1496 bio(i,1,ioxyg)=bio(i,1,ioxyg)-cff1*cff3
1497# endif
1498# else
1499 bio(i,1,inh4_)=bio(i,1,inh4_)+cff1
1500# ifdef PO4
1501 bio(i,1,ipo4_)=bio(i,1,ipo4_)+cff1*r_p2n(ng)
1502# endif
1503# ifdef OXYGEN
1504 bio(i,1,ioxyg)=bio(i,1,ioxyg)-cff1*cff4
1505# endif
1506# if defined CARBON && defined TALK_NONCONSERV
1507 bio(i,1,italk)=bio(i,1,italk)+cff1
1508# endif
1509# endif
1510 END DO
1511 END IF
1512# ifdef CARBON
1513# ifdef DENITRIFICATION
1514 cff3=12.0_r8
1515 cff4=0.74_r8
1516# endif
1517 IF ((ibio.eq.isdec).or. &
1518 & (ibio.eq.ildec))THEN
1519 DO i=istr,iend
1520 cff1=fc(i,0)*hz_inv(i,1)
1521 bio(i,1,itic_)=bio(i,1,itic_)+cff1
1522 END DO
1523 END IF
1524 IF (ibio.eq.iphyt)THEN
1525 DO i=istr,iend
1526 cff1=fc(i,0)*hz_inv(i,1)
1527 bio(i,1,itic_)=bio(i,1,itic_)+cff1*phycn(ng)
1528 END DO
1529 END IF
1530# endif
1531#endif
1532 END DO sink_loop
1533 END DO iter_loop
1534!
1535!-----------------------------------------------------------------------
1536! Update global tracer variables: Add increment due to BGC processes
1537! to tracer array in time index "nnew". Index "nnew" is solution after
1538! advection and mixing and has transport units (m Tunits) hence the
1539! increment is multiplied by Hz. Notice that we need to subtract
1540! original values "Bio_old" at the top of the routine to just account
1541! for the concentractions affected by BGC processes. This also takes
1542! into account any constraints (non-negative concentrations, carbon
1543! concentration range) specified before entering BGC kernel. If "Bio"
1544! were unchanged by BGC processes, the increment would be exactly
1545! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
1546! bounded >=0 so that we can preserve total inventory of N and
1547! C even when advection causes tracer concentration to go negative.
1548! (J. Wilkin and H. Arango, Apr 27, 2012)
1549!-----------------------------------------------------------------------
1550!
1551#ifdef CARBON
1552 DO k=1,n(ng)
1553 DO i=istr,iend
1554 bio(i,k,itic_)=min(bio(i,k,itic_),3000.0_r8)
1555 bio(i,k,itic_)=max(bio(i,k,itic_),400.0_r8)
1556 END DO
1557 END DO
1558#endif
1559 DO itrc=1,nbt
1560 ibio=idbio(itrc)
1561 DO k=1,n(ng)
1562 DO i=istr,iend
1563 cff=bio(i,k,ibio)-bio_old(i,k,ibio)
1564#ifdef MASKING
1565 cff=cff*rmask(i,j)
1566# ifdef WET_DRY
1567 cff=cff*rmask_wet(i,j)
1568# endif
1569#endif
1570 t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*hz(i,j,k)
1571 END DO
1572 END DO
1573 END DO
1574 END DO j_loop
1575!
1576 RETURN
1577 END SUBROUTINE fennel_tile
1578
1579#ifdef CARBON
1580# ifdef pCO2_RZ
1581 SUBROUTINE pco2_water_rz (Istr, Iend, &
1582 & LBi, UBi, LBj, UBj, IminS, ImaxS, &
1583 & j, DoNewton, &
1584# ifdef MASKING
1585 & rmask, &
1586# endif
1587 & T, S, TIC, TAlk, pH, pCO2)
1588!
1589!***********************************************************************
1590! !
1591! This routine computes equilibrium partial pressure of CO2 (pCO2) !
1592! in the surface seawater. !
1593! !
1594! On Input: !
1595! !
1596! Istr Starting tile index in the I-direction. !
1597! Iend Ending tile index in the I-direction. !
1598! LBi I-dimension lower bound. !
1599! UBi I-dimension upper bound. !
1600! LBj J-dimension lower bound. !
1601! UBj J-dimension upper bound. !
1602! IminS I-dimension lower bound for private arrays. !
1603! ImaxS I-dimension upper bound for private arrays. !
1604! j j-pipelined index. !
1605! DoNewton Iteration solver: !
1606! [0] Bracket and bisection. !
1607! [1] Newton-Raphson method. !
1608! rmask Land/Sea masking. !
1609! T Surface temperature (Celsius). !
1610! S Surface salinity (PSS). !
1611! TIC Total inorganic carbon (millimol/m3). !
1612! TAlk Total alkalinity (milli-equivalents/m3). !
1613! pH Best pH guess. !
1614! !
1615! On Output: !
1616! !
1617! pCO2 partial pressure of CO2 (ppmv). !
1618! !
1619! Check Value: (T=24, S=36.6, TIC=2040, TAlk=2390, PO4b=0, !
1620! SiO3=0, pH=8) !
1621! !
1622! pcO2= ppmv (DoNewton=0) !
1623! pCO2= ppmv (DoNewton=1) !
1624! !
1625! This subroutine was adapted by Katja Fennel (Nov 2005) from !
1626! Zeebe and Wolf-Gladrow (2001). !
1627! !
1628! Reference: !
1629! !
1630! Zeebe, R.E. and D. Wolf-Gladrow, 2005: CO2 in Seawater: !
1631! Equilibrium, kinetics, isotopes, Elsevier Oceanographic !
1632! Series, 65, pp 346. !
1633! !
1634!***********************************************************************
1635!
1636 USE mod_kinds
1637!
1638 implicit none
1639!
1640! Imported variable declarations.
1641!
1642 integer, intent(in) :: LBi, UBi, LBj, UBj, IminS, ImaxS
1643 integer, intent(in) :: Istr, Iend, j, DoNewton
1644!
1645# ifdef ASSUMED_SHAPE
1646# ifdef MASKING
1647 real(r8), intent(in) :: rmask(LBi:,LBj:)
1648# endif
1649 real(r8), intent(in) :: T(IminS:)
1650 real(r8), intent(in) :: S(IminS:)
1651 real(r8), intent(in) :: TIC(IminS:)
1652 real(r8), intent(in) :: TAlk(IminS:)
1653 real(r8), intent(inout) :: pH(LBi:,LBj:)
1654# else
1655# ifdef MASKING
1656 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1657# endif
1658 real(r8), intent(in) :: T(IminS:ImaxS)
1659 real(r8), intent(in) :: S(IminS:ImaxS)
1660 real(r8), intent(in) :: TIC(IminS:ImaxS)
1661 real(r8), intent(in) :: TAlk(IminS:ImaxS)
1662 real(r8), intent(inout) :: pH(LBi:UBi,LBj:UBj)
1663# endif
1664
1665 real(r8), intent(out) :: pCO2(IminS:ImaxS)
1666!
1667! Local variable declarations.
1668!
1669 integer, parameter :: InewtonMax = 10
1670 integer, parameter :: IbrackMax = 30
1671
1672 integer :: Hstep, Ibrack, Inewton, i
1673
1674 real(r8) :: Tk, centiTk, invTk, logTk
1675 real(r8) :: scl, sqrtS
1676 real(r8) :: borate, alk, dic
1677 real(r8) :: ff, K1, K2, K12, Kb, Kw
1678 real(r8) :: p5, p4, p3, p2, p1, p0
1679 real(r8) :: df, fn, fni(3), ftest
1680 real(r8) :: deltaX, invX, invX2, X, X2, X3
1681 real(r8) :: pH_guess, pH_hi, pH_lo
1682 real(r8) :: X_guess, X_hi, X_lo, X_mid
1683 real(r8) :: CO2star, Htotal, Htotal2
1684!
1685!=======================================================================
1686! Determine coefficients for surface carbon chemisty. If land/sea
1687! masking, compute only on water points.
1688!=======================================================================
1689!
1690 i_loop: DO i=istr,iend
1691# ifdef MASKING
1692 IF (rmask(i,j).gt.0.0_r8) THEN
1693# endif
1694 tk=t(i)+273.15_r8
1695 centitk=0.01_r8*tk
1696 invtk=1.0_r8/tk
1697 logtk=log(tk)
1698 sqrts=sqrt(s(i))
1699 scl=s(i)/1.80655_r8
1700
1701 alk= talk(i)*0.000001_r8
1702 dic = tic(i)*0.000001_r8
1703!
1704!-----------------------------------------------------------------------
1705! Correction term for non-ideality, ff=k0*(1-pH2O). Equation 13 with
1706! table 6 values from Weiss and Price (1980, Mar. Chem., 8, 347-359).
1707!-----------------------------------------------------------------------
1708!
1709 ff=exp(-162.8301_r8+ &
1710 & 218.2968_r8/centitk+ &
1711 & log(centitk)*90.9241_r8- &
1712 & centitk*centitk*1.47696_r8+ &
1713 & s(i)*(0.025695_r8- &
1714 & centitk*(0.025225_r8- &
1715 & centitk*0.0049867_r8)))
1716!
1717!-----------------------------------------------------------------------
1718! Compute first (K1) and second (K2) dissociation constant of carboinic
1719! acid:
1720!
1721! K1 = [H][HCO3]/[H2CO3]
1722! K2 = [H][CO3]/[HCO3]
1723!
1724! From Millero (1995; page 664) using Mehrbach et al. (1973) data on
1725! seawater scale.
1726!-----------------------------------------------------------------------
1727!
1728 k1=10.0_r8**(62.008_r8- &
1729 & invtk*3670.7_r8- &
1730 & logtk*9.7944_r8+ &
1731 & s(i)*(0.0118_r8- &
1732 & s(i)*0.000116_r8))
1733 k2=10.0_r8**(-4.777_r8- &
1734 & invtk*1394.7_r8+ &
1735 & s(i)*(0.0184_r8- &
1736 & s(i)*0.000118_r8))
1737!
1738!-----------------------------------------------------------------------
1739! Compute dissociation constant of boric acid, Kb=[H][BO2]/[HBO2].
1740! From Millero (1995; page 669) using data from Dickson (1990).
1741!-----------------------------------------------------------------------
1742!
1743 kb=exp(-invtk*(8966.90_r8+ &
1744 & sqrts*(2890.53_r8+ &
1745 & sqrts*(77.942_r8- &
1746 & sqrts*(1.728_r8- &
1747 & sqrts*0.0996_r8))))- &
1748 & logtk*(24.4344_r8+ &
1749 & sqrts*(25.085_r8+ &
1750 & sqrts*0.2474_r8))+ &
1751 & tk*(sqrts*0.053105_r8)+ &
1752 & 148.0248_r8+ &
1753 & sqrts*(137.1942_r8+ &
1754 & sqrts*1.62142_r8))
1755!
1756!-----------------------------------------------------------------------
1757! Compute ion product of whater, Kw = [H][OH].
1758! From Millero (1995; page 670) using composite data.
1759!-----------------------------------------------------------------------
1760!
1761 kw=exp(148.9652_r8- &
1762 & invtk*13847.26_r8- &
1763 & logtk*23.6521_r8- &
1764 & sqrts*(5.977_r8- &
1765 & invtk*118.67_r8- &
1766 & logtk*1.0495_r8)- &
1767 & s(i)*0.01615_r8)
1768!
1769!-----------------------------------------------------------------------
1770! Calculate concentrations for borate (Uppstrom, 1974).
1771!-----------------------------------------------------------------------
1772!
1773 borate=0.000232_r8*scl/10.811_r8
1774!
1775!=======================================================================
1776! Iteratively solver for computing hydrogen ions [H+] using either:
1777!
1778! (1) Newton-Raphson method with fixed number of iterations,
1779! use previous [H+] as first guess, or
1780! (2) bracket and bisection
1781!=======================================================================
1782!
1783! Solve for h in fifth-order polynomial. First calculate
1784! polynomial coefficients.
1785!
1786 k12 = k1*k2
1787
1788 p5 = -1.0_r8;
1789 p4 = -alk-kb-k1;
1790 p3 = dic*k1-alk*(kb+k1)+kb*borate+kw-kb*k1-k12
1791 p2 = dic*(kb*k1+2*k12)-alk*(kb*k1+k12)+kb*borate*k1 &
1792 & +(kw*kb+kw*k1-kb*k12)
1793 p1 = 2.0_r8*dic*kb*k12-alk*kb*k12+kb*borate*k12 &
1794 & +kw*kb*k1+kw*k12
1795 p0 = kw*kb*k12;
1796!
1797! Set first guess and brackets for [H+] solvers.
1798!
1799 ph_guess=ph(i,j) ! Newton-Raphson
1800 ph_hi=10.0_r8 ! high bracket/bisection
1801 ph_lo=5.0_r8 ! low bracket/bisection
1802!
1803! Convert to [H+].
1804!
1805 x_guess=10.0_r8**(-ph_guess)
1806 x_lo=10.0_r8**(-ph_hi)
1807 x_hi=10.0_r8**(-ph_lo)
1808 x_mid=0.5_r8*(x_lo+x_hi)
1809!
1810!-----------------------------------------------------------------------
1811! Newton-Raphson method.
1812!-----------------------------------------------------------------------
1813!
1814 IF (donewton.eq.1) THEN
1815 x=x_guess
1816!
1817 DO inewton=1,inewtonmax
1818!
1819! Evaluate f([H+]) = p5*x^5+...+p1*x+p0
1820!
1821 fn=((((p5*x+p4)*x+p3)*x+p2)*x+p1)*x+p0
1822!
1823! Evaluate derivative, df([H+])/dx:
1824!
1825! df= d(fn)/d(X)
1826!
1827 df=(((5*p5*x+4*p4)*x+3*p3)*x+2*p2)*x+p1
1828!
1829! Evaluate increment in [H+].
1830!
1831 deltax=-fn/df
1832!
1833! Update estimate of [H+].
1834!
1835 x=x+deltax
1836 END DO
1837!
1838!-----------------------------------------------------------------------
1839! Bracket and bisection method.
1840!-----------------------------------------------------------------------
1841!
1842 ELSE
1843!
1844! If first step, use Bracket and Bisection method with fixed, large
1845! number of iterations
1846!
1847 brack_it: DO ibrack=1,ibrackmax
1848 DO hstep=1,3
1849 IF (hstep.eq.1) x=x_hi
1850 IF (hstep.eq.2) x=x_lo
1851 IF (hstep.eq.3) x=x_mid
1852!
1853! Evaluate f([H+]) for bracketing and mid-value cases.
1854!
1855 fni(hstep)=((((p5*x+p4)*x+p3)*x+p2)*x+p1)*x+p0
1856 END DO
1857!
1858! Now, bracket solution within two of three.
1859!
1860 IF (fni(3).eq.0) THEN
1861 EXIT brack_it
1862 ELSE
1863 ftest=fni(1)/fni(3)
1864 IF (ftest.gt.0) THEN
1865 x_hi=x_mid
1866 ELSE
1867 x_lo=x_mid
1868 END IF
1869 x_mid=0.5_r8*(x_lo+x_hi)
1870 END IF
1871 END DO brack_it
1872!
1873! Last iteration gives value.
1874!
1875 x=x_mid
1876 END IF
1877!
1878!-----------------------------------------------------------------------
1879! Determine pCO2.
1880!-----------------------------------------------------------------------
1881!
1882! Total Hydrogen ion concentration, Htotal = [H+].
1883!
1884 htotal=x
1885 htotal2=htotal*htotal
1886!
1887! Calculate [CO2*] (mole/m3) as defined in DOE Methods Handbook 1994
1888! Version 2, ORNL/CDIAC-74, Dickson and Goyet, Eds. (Chapter 2,
1889! page 10, Eq A.49).
1890!
1891 co2star=dic*htotal2/(htotal2+k1*htotal+k1*k2)
1892!
1893! Save pH is used again outside this routine.
1894!
1895 ph(i,j)=-log10(htotal)
1896!
1897! Add two output arguments for storing pCO2surf.
1898!
1899 pco2(i)=co2star*1000000.0_r8/ff
1900
1901# ifdef MASKING
1902 ELSE
1903 ph(i,j)=0.0_r8
1904 pco2(i)=0.0_r8
1905 END IF
1906# endif
1907
1908 END DO i_loop
1909!
1910 RETURN
1911 END SUBROUTINE pco2_water_rz
1912# else
1913 SUBROUTINE pco2_water (Istr, Iend, &
1914 & LBi, UBi, LBj, UBj, &
1915 & IminS, ImaxS, j, DoNewton, &
1916# ifdef MASKING
1917 & rmask, &
1918# endif
1919 & T, S, TIC, TAlk, PO4b, SiO3, pH, pCO2)
1920!
1921!***********************************************************************
1922! !
1923! This routine computes equilibrium partial pressure of CO2 (pCO2) !
1924! in the surface seawater. !
1925! !
1926! On Input: !
1927! !
1928! Istr Starting tile index in the I-direction. !
1929! Iend Ending tile index in the I-direction. !
1930! LBi I-dimension lower bound. !
1931! UBi I-dimension upper bound. !
1932! LBj J-dimension lower bound. !
1933! UBj J-dimension upper bound. !
1934! IminS I-dimension lower bound for private arrays. !
1935! ImaxS I-dimension upper bound for private arrays. !
1936! j j-pipelined index. !
1937! DoNewton Iteration solver: !
1938! [0] Bracket and bisection. !
1939! [1] Newton-Raphson method. !
1940! rmask Land/Sea masking. !
1941! T Surface temperature (Celsius). !
1942! S Surface salinity (PSS). !
1943! TIC Total inorganic carbon (millimol/m3). !
1944! TAlk Total alkalinity (milli-equivalents/m3). !
1945! PO4b Inorganic phosphate (millimol/m3). !
1946! SiO3 Inorganic silicate (millimol/m3). !
1947! pH Best pH guess. !
1948! !
1949! On Output: !
1950! !
1951! pCO2 partial pressure of CO2 (ppmv). !
1952! !
1953! Check Value: (T=24, S=36.6, TIC=2040, TAlk=2390, PO4b=0, !
1954! SiO3=0, pH=8) !
1955! !
1956! pcO2=0.35074945E+03 ppmv (DoNewton=0) !
1957! pCO2=0.35073560E+03 ppmv (DoNewton=1) !
1958! !
1959! This subroutine was adapted by Mick Follows (Oct 1999) from OCMIP2 !
1960! code CO2CALC. Modified for ROMS by Hernan Arango (Nov 2003). !
1961! !
1962!***********************************************************************
1963!
1964 USE mod_kinds
1965!
1966 implicit none
1967!
1968! Imported variable declarations.
1969!
1970 integer, intent(in) :: LBi, UBi, LBj, UBj, IminS, ImaxS
1971 integer, intent(in) :: Istr, Iend, j, DoNewton
1972!
1973# ifdef ASSUMED_SHAPE
1974# ifdef MASKING
1975 real(r8), intent(in) :: rmask(LBi:,LBj:)
1976# endif
1977 real(r8), intent(in) :: T(IminS:)
1978 real(r8), intent(in) :: S(IminS:)
1979 real(r8), intent(in) :: TIC(IminS:)
1980 real(r8), intent(in) :: TAlk(IminS:)
1981 real(r8), intent(inout) :: pH(LBi:,LBj:)
1982# else
1983# ifdef MASKING
1984 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1985# endif
1986 real(r8), intent(in) :: T(IminS:ImaxS)
1987 real(r8), intent(in) :: S(IminS:ImaxS)
1988 real(r8), intent(in) :: TIC(IminS:ImaxS)
1989 real(r8), intent(in) :: TAlk(IminS:ImaxS)
1990 real(r8), intent(inout) :: pH(LBi:UBi,LBj:UBj)
1991# endif
1992 real(r8), intent(in) :: PO4b
1993 real(r8), intent(in) :: SiO3
1994
1995 real(r8), intent(out) :: pCO2(IminS:ImaxS)
1996!
1997! Local variable declarations.
1998!
1999 integer, parameter :: InewtonMax = 10
2000 integer, parameter :: IbrackMax = 30
2001
2002 integer :: Hstep, Ibrack, Inewton, i
2003
2004 real(r8) :: Tk, centiTk, invTk, logTk
2005 real(r8) :: SO4, scl, sqrtS, sqrtSO4
2006 real(r8) :: alk, dic, phos, sili
2007 real(r8) :: borate, sulfate, fluoride
2008 real(r8) :: ff, K1, K2, K1p, K2p, K3p, Kb, Kf, Ks, Ksi, Kw
2009 real(r8) :: K12, K12p, K123p, invKb, invKs, invKsi
2010 real(r8) :: A, A2, B, B2, C, dA, dB
2011 real(r8) :: df, fn, fni(3), ftest
2012 real(r8) :: deltaX, invX, invX2, X, X2, X3
2013 real(r8) :: pH_guess, pH_hi, pH_lo
2014 real(r8) :: X_guess, X_hi, X_lo, X_mid
2015 real(r8) :: CO2star, Htotal, Htotal2
2016!
2017!=======================================================================
2018! Determine coefficients for surface carbon chemisty. If land/sea
2019! masking, compute only on water points.
2020!=======================================================================
2021!
2022 i_loop: DO i=istr,iend
2023# ifdef MASKING
2024 IF (rmask(i,j).gt.0.0_r8) THEN
2025# endif
2026 tk=t(i)+273.15_r8
2027 centitk=0.01_r8*tk
2028 invtk=1.0_r8/tk
2029 logtk=log(tk)
2030 sqrts=sqrt(s(i))
2031 so4=19.924_r8*s(i)/(1000.0_r8-1.005_r8*s(i))
2032 sqrtso4=sqrt(so4)
2033 scl=s(i)/1.80655_r8
2034
2035 alk=talk(i)*0.000001_r8
2036 dic=tic(i)*0.000001_r8
2037 phos=po4b*0.000001_r8
2038 sili=sio3*0.000001_r8
2039!
2040!-----------------------------------------------------------------------
2041! Correction term for non-ideality, ff=k0*(1-pH2O). Equation 13 with
2042! table 6 values from Weiss and Price (1980, Mar. Chem., 8, 347-359).
2043!-----------------------------------------------------------------------
2044!
2045 ff=exp(-162.8301_r8+ &
2046 & 218.2968_r8/centitk+ &
2047 & log(centitk)*90.9241_r8- &
2048 & centitk*centitk*1.47696_r8+ &
2049 & s(i)*(0.025695_r8- &
2050 & centitk*(0.025225_r8- &
2051 & centitk*0.0049867_r8)))
2052!
2053!-----------------------------------------------------------------------
2054! Compute first (K1) and second (K2) dissociation constant of carboinic
2055! acid:
2056!
2057! K1 = [H][HCO3]/[H2CO3]
2058! K2 = [H][CO3]/[HCO3]
2059!
2060! From Millero (1995; page 664) using Mehrbach et al. (1973) data on
2061! seawater scale.
2062!-----------------------------------------------------------------------
2063!
2064 k1=10.0_r8**(62.008_r8- &
2065 & invtk*3670.7_r8- &
2066 & logtk*9.7944_r8+ &
2067 & s(i)*(0.0118_r8- &
2068 & s(i)*0.000116_r8))
2069 k2=10.0_r8**(-4.777_r8- &
2070 & invtk*1394.7_r8+ &
2071 & s(i)*(0.0184_r8- &
2072 & s(i)*0.000118_r8))
2073!
2074!-----------------------------------------------------------------------
2075! Compute dissociation constant of boric acid, Kb=[H][BO2]/[HBO2].
2076! From Millero (1995; page 669) using data from Dickson (1990).
2077!-----------------------------------------------------------------------
2078!
2079 kb=exp(-invtk*(8966.90_r8+ &
2080 & sqrts*(2890.53_r8+ &
2081 & sqrts*(77.942_r8- &
2082 & sqrts*(1.728_r8- &
2083 & sqrts*0.0996_r8))))- &
2084 & logtk*(24.4344_r8+ &
2085 & sqrts*(25.085_r8+ &
2086 & sqrts*0.2474_r8))+ &
2087 & tk*(sqrts*0.053105_r8)+ &
2088 & 148.0248_r8+ &
2089 & sqrts*(137.1942_r8+ &
2090 & sqrts*1.62142_r8))
2091!
2092!-----------------------------------------------------------------------
2093! Compute first (K1p), second (K2p), and third (K3p) dissociation
2094! constant of phosphoric acid:
2095!
2096! K1p = [H][H2PO4]/[H3PO4]
2097! K2p = [H][HPO4]/[H2PO4]
2098! K3p = [H][PO4]/[HPO4]
2099!
2100! From DOE (1994) equations 7.2.20, 7.2.23, and 7.2.26, respectively.
2101! With footnote using data from Millero (1974).
2102!-----------------------------------------------------------------------
2103!
2104 k1p=exp(115.525_r8- &
2105 & invtk*4576.752_r8- &
2106 & logtk*18.453_r8+ &
2107 & sqrts*(0.69171_r8-invtk*106.736_r8)- &
2108 & s(i)*(0.01844_r8+invtk*0.65643_r8))
2109 k2p=exp(172.0883_r8- &
2110 & invtk*8814.715_r8- &
2111 & logtk*27.927_r8+ &
2112 & sqrts*(1.3566_r8-invtk*160.340_r8)- &
2113 & s(i)*(0.05778_r8-invtk*0.37335_r8))
2114 k3p=exp(-18.141_r8- &
2115 & invtk*3070.75_r8+ &
2116 & sqrts*(2.81197_r8+invtk*17.27039_r8)- &
2117 & s(i)*(0.09984_r8+invtk*44.99486_r8))
2118!
2119!-----------------------------------------------------------------------
2120! Compute dissociation constant of silica, Ksi=[H][SiO(OH)3]/[Si(OH)4].
2121! From Millero (1995; page 671) using data from Yao and Millero (1995).
2122!-----------------------------------------------------------------------
2123!
2124 ksi=exp(117.385_r8- &
2125 & invtk*8904.2_r8- &
2126 & logtk*19.334_r8+ &
2127 & sqrtso4*(3.5913_r8-invtk*458.79_r8)- &
2128 & so4*(1.5998_r8-invtk*188.74_r8- &
2129 & so4*(0.07871_r8-invtk*12.1652_r8))+ &
2130 & log(1.0_r8-0.001005_r8*s(i)))
2131!
2132!-----------------------------------------------------------------------
2133! Compute ion product of whater, Kw = [H][OH].
2134! From Millero (1995; page 670) using composite data.
2135!-----------------------------------------------------------------------
2136!
2137 kw=exp(148.9652_r8- &
2138 & invtk*13847.26_r8- &
2139 & logtk*23.6521_r8- &
2140 & sqrts*(5.977_r8- &
2141 & invtk*118.67_r8- &
2142 & logtk*1.0495_r8)- &
2143 & s(i)*0.01615_r8)
2144!
2145!------------------------------------------------------------------------
2146! Compute salinity constant of hydrogen sulfate, Ks = [H][SO4]/[HSO4].
2147! From Dickson (1990, J. chem. Thermodynamics 22, 113)
2148!------------------------------------------------------------------------
2149!
2150 ks=exp(141.328_r8- &
2151 & invtk*4276.1_r8- &
2152 & logtk*23.093_r8+ &
2153 & sqrtso4*(324.57_r8-invtk*13856.0_r8-logtk*47.986_r8- &
2154 & so4*invtk*2698.0_r8)- &
2155 & so4*(771.54_r8-invtk*35474.0_r8-logtk*114.723_r8- &
2156 & so4*invtk*1776.0_r8)+ &
2157 & log(1.0_r8-0.001005_r8*s(i)))
2158!
2159!-----------------------------------------------------------------------
2160! Compute stability constant of hydrogen fluorid, Kf = [H][F]/[HF].
2161! From Dickson and Riley (1979) -- change pH scale to total.
2162!-----------------------------------------------------------------------
2163!
2164 kf=exp(-12.641_r8+ &
2165 & invtk*1590.2_r8+ &
2166 & sqrtso4*1.525_r8+ &
2167 & log(1.0_r8-0.001005_r8*s(i))+ &
2168 & log(1.0_r8+0.1400_r8*scl/(96.062_r8*ks)))
2169!
2170!-----------------------------------------------------------------------
2171! Calculate concentrations for borate (Uppstrom, 1974), sulfate (Morris
2172! and Riley, 1966), and fluoride (Riley, 1965).
2173!-----------------------------------------------------------------------
2174!
2175 borate=0.000232_r8*scl/10.811_r8
2176 sulfate=0.14_r8*scl/96.062_r8
2177 fluoride=0.000067_r8*scl/18.9984_r8
2178!
2179!=======================================================================
2180! Iteratively solver for computing hydrogen ions [H+] using either:
2181!
2182! (1) Newton-Raphson method with fixed number of iterations,
2183! use previous [H+] as first guess, or
2184! (2) bracket and bisection
2185!=======================================================================
2186!
2187! Set first guess and brackets for [H+] solvers.
2188!
2189 ph_guess=ph(i,j) ! Newton-Raphson
2190 ph_hi=10.0_r8 ! high bracket/bisection
2191 ph_lo=5.0_r8 ! low bracket/bisection
2192!
2193! Convert to [H+].
2194!
2195 x_guess=10.0_r8**(-ph_guess)
2196 x_lo=10.0_r8**(-ph_hi)
2197 x_hi=10.0_r8**(-ph_lo)
2198 x_mid=0.5_r8*(x_lo+x_hi)
2199!
2200!-----------------------------------------------------------------------
2201! Newton-Raphson method.
2202!-----------------------------------------------------------------------
2203!
2204 IF (donewton.eq.1) THEN
2205 x=x_guess
2206 k12=k1*k2
2207 k12p=k1p*k2p
2208 k123p=k12p*k3p
2209 invkb=1.0_r8/kb
2210 invks=1.0_r8/ks
2211 invksi=1.0_r8/ksi
2212!
2213 DO inewton=1,inewtonmax
2214!
2215! Set some common combinations of parameters used in the iterative [H+]
2216! solver.
2217!
2218 x2=x*x
2219 x3=x2*x
2220 invx=1.0_r8/x
2221 invx2=1.0_r8/x2
2222
2223 a=x*(k12p+x*(k1p+x))
2224 b=x*(k1+x)+k12
2225 c=1.0_r8/(1.0_r8+sulfate*invks)
2226
2227 a2=a*a
2228 b2=b*b
2229 da=x*(2.0_r8*k1p+3.0_r8*x)+k12p
2230 db=2.0_r8*x+k1
2231!
2232! Evaluate f([H+]):
2233!
2234! fn=HCO3+CO3+borate+OH+HPO4+2*PO4+H3PO4+silicate+Hfree+HSO4+HF-TALK
2235!
2236 fn=dic*k1*(x+2.0_r8*k2)/b+ &
2237 & borate/(1.0_r8+x*invkb)+ &
2238 & kw*invx+ &
2239 & phos*(k12p*x+2.0_r8*k123p-x3)/a+ &
2240 & sili/(1.0_r8+x*invksi)- &
2241 & x*c- &
2242 & sulfate/(1.0_r8+ks*invx*c)- &
2243 & fluoride/(1.0_r8+kf*invx)- &
2244 & alk
2245!
2246! Evaluate derivative, f(prime)([H+]):
2247!
2248! df= d(fn)/d(X)
2249!
2250 df=dic*k1*(b-db*(x+2.0_r8*k2))/b2- &
2251 & borate/(invkb*(1.0+x*invkb)**2)- &
2252 & kw*invx2+ &
2253 & phos*(a*(k12p-3.0_r8*x2)-da*(k12p*x+2.0_r8*k123p-x3))/a2-&
2254 & sili/(invksi*(1.0_r8+x*invksi)**2)+ &
2255 & c+ &
2256 & sulfate*ks*c*invx2/((1.0_r8+ks*invx*c)**2)+ &
2257 & fluoride*kf*invx2/((1.0_r8+kf*invx)**2)
2258!
2259! Evaluate increment in [H+].
2260!
2261 deltax=-fn/df
2262!
2263! Update estimate of [H+].
2264!
2265 x=x+deltax
2266 END DO
2267!
2268!-----------------------------------------------------------------------
2269! Bracket and bisection method.
2270!-----------------------------------------------------------------------
2271!
2272 ELSE
2273!
2274! If first step, use Bracket and Bisection method with fixed, large
2275! number of iterations
2276!
2277 k12=k1*k2
2278 k12p=k1p*k2p
2279 k123p=k12p*k3p
2280 invkb=1.0_r8/kb
2281 invks=1.0_r8/ks
2282 invksi=1.0_r8/ksi
2283!
2284 brack_it: DO ibrack=1,ibrackmax
2285 DO hstep=1,3
2286 IF (hstep.eq.1) x=x_hi
2287 IF (hstep.eq.2) x=x_lo
2288 IF (hstep.eq.3) x=x_mid
2289!
2290! Set some common combinations of parameters used in the iterative [H+]
2291! solver.
2292!
2293 x2=x*x
2294 x3=x2*x
2295 invx=1.0_r8/x
2296
2297 a=x*(k12p+x*(k1p+x))+k123p
2298 b=x*(k1+x)+k12
2299 c=1.0_r8/(1.0_r8+sulfate*invks)
2300
2301 a2=a*a
2302 b2=b*b
2303 da=x*(k1p*2.0_r8+3.0_r8*x2)+k12p
2304 db=2.0_r8*x+k1
2305!
2306! Evaluate f([H+]) for bracketing and mid-value cases.
2307!
2308 fni(hstep)=dic*(k1*x+2.0_r8*k12)/b+ &
2309 & borate/(1.0_r8+x*invkb)+ &
2310 & kw*invx+ &
2311 & phos*(k12p*x+2.0_r8*k123p-x3)/a+ &
2312 & sili/(1.0_r8+x*invksi)- &
2313 & x*c- &
2314 & sulfate/(1.0_r8+ks*invx*c)- &
2315 & fluoride/(1.0_r8+kf*invx)- &
2316 & alk
2317 END DO
2318!
2319! Now, bracket solution within two of three.
2320!
2321 IF (fni(3).eq.0.0_r8) THEN
2322 EXIT brack_it
2323 ELSE
2324 ftest=fni(1)/fni(3)
2325 IF (ftest.gt.0.0) THEN
2326 x_hi=x_mid
2327 ELSE
2328 x_lo=x_mid
2329 END IF
2330 x_mid=0.5_r8*(x_lo+x_hi)
2331 END IF
2332 END DO brack_it
2333!
2334! Last iteration gives value.
2335!
2336 x=x_mid
2337 END IF
2338!
2339!-----------------------------------------------------------------------
2340! Determine pCO2.
2341!-----------------------------------------------------------------------
2342!
2343! Total Hydrogen ion concentration, Htotal = [H+].
2344!
2345 htotal=x
2346 htotal2=htotal*htotal
2347!
2348! Calculate [CO2*] (mole/m3) as defined in DOE Methods Handbook 1994
2349! Version 2, ORNL/CDIAC-74, Dickson and Goyet, Eds. (Chapter 2,
2350! page 10, Eq A.49).
2351!
2352 co2star=dic*htotal2/(htotal2+k1*htotal+k1*k2)
2353!
2354! Save pH is used again outside this routine.
2355!
2356 ph(i,j)=-log10(htotal)
2357!
2358! Add two output arguments for storing pCO2surf.
2359!
2360 pco2(i)=co2star*1000000.0_r8/ff
2361
2362# ifdef MASKING
2363 ELSE
2364 ph(i,j)=0.0_r8
2365 pco2(i)=0.0_r8
2366 END IF
2367# endif
2368
2369 END DO i_loop
2370!
2371 RETURN
2372 END SUBROUTINE pco2_water
2373# endif
2374#endif
2375 END MODULE biology_mod
biology_mod
Definition biology.F:120
biology_mod::biology
subroutine, public biology(ng, tile)
Definition ecosim.h:57
biology_mod::pco2_water
subroutine pco2_water(istr, iend, lbi, ubi, lbj, ubj, imins, imaxs, j, donewton, rmask, t, s, tic, talk, po4b, sio3, ph, pco2)
Definition fennel.h:1920
biology_mod::pco2_water_rz
subroutine pco2_water_rz(istr, iend, lbi, ubi, lbj, ubj, imins, imaxs, j, donewton, rmask, t, s, tic, talk, ph, pco2)
Definition fennel.h:1588
biology_mod::fennel_tile
subroutine fennel_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nnew, rmask, rmask_wet, rmask_full, hz, z_r, z_w, srflx, uwind, vwind, sustr, svstr, ph, diabio2d, diabio3d, t)
Definition fennel.h:242
dateclock_mod
Definition dateclock.F:2
dateclock_mod::caldate
subroutine, public caldate(currenttime, yy_i, yd_i, mm_i, dd_i, h_i, m_i, s_i, yd_dp, dd_dp, h_dp, m_dp, s_dp)
Definition dateclock.F:76
mod_biology
Definition ecosim_mod.h:1
mod_biology::isden
integer isden
Definition fennel_mod.h:84
mod_biology::ipco2
integer ipco2
Definition fennel_mod.h:109
mod_biology::k_po4
real(r8), dimension(:), allocatable k_po4
Definition fennel_mod.h:134
mod_biology::idnit
integer idnit
Definition fennel_mod.h:108
mod_biology::sderrc
real(r8), dimension(:), allocatable sderrc
Definition fennel_mod.h:149
mod_biology::parfrac
real(r8), dimension(:), allocatable parfrac
Definition fennel_mod.h:139
mod_biology::zoomr
real(r8), dimension(:), allocatable zoomr
Definition fennel_mod.h:162
mod_biology::rderrc
real(r8), dimension(:), allocatable rderrc
Definition fennel_mod.h:151
mod_biology::wldet
real(r8), dimension(:), allocatable wldet
Definition fennel_mod.h:153
mod_biology::zoocn
real(r8), dimension(:), allocatable zoocn
Definition fennel_mod.h:158
mod_biology::vp0
real(r8), dimension(:), allocatable vp0
Definition fennel_mod.h:152
mod_biology::sderrn
real(r8), dimension(:), allocatable sderrn
Definition fennel_mod.h:148
mod_biology::k_phy
real(r8), dimension(:), allocatable k_phy
Definition fennel_mod.h:135
mod_biology::itic_
integer itic_
Definition fennel_mod.h:91
mod_biology::nitrir
real(r8), dimension(:), allocatable nitrir
Definition fennel_mod.h:138
mod_biology::r_p2n
real(r8), dimension(:), allocatable r_p2n
Definition fennel_mod.h:141
mod_biology::ilden
integer ilden
Definition fennel_mod.h:83
mod_biology::bioiter
integer, dimension(:), allocatable bioiter
Definition ecosim_mod.h:343
mod_biology::k_nh4
real(r8), dimension(:), allocatable k_nh4
Definition fennel_mod.h:132
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::chlmin
real(r8), dimension(:), allocatable chlmin
Definition fennel_mod.h:128
mod_biology::isdec
integer isdec
Definition fennel_mod.h:90
mod_biology::ino3_
integer ino3_
Definition ecosim_mod.h:277
mod_biology::irden
integer irden
Definition fennel_mod.h:86
mod_biology::chl2c_m
real(r8), dimension(:), allocatable chl2c_m
Definition fennel_mod.h:127
mod_biology::ipo4_
integer ipo4_
Definition ecosim_mod.h:279
mod_biology::zoogr
real(r8), dimension(:), allocatable zoogr
Definition fennel_mod.h:160
mod_biology::d_p5nh4
real(r8), dimension(:), allocatable d_p5nh4
Definition fennel_mod.h:130
mod_biology::rderrn
real(r8), dimension(:), allocatable rderrn
Definition fennel_mod.h:150
mod_biology::zoomin
real(r8), dimension(:), allocatable zoomin
Definition fennel_mod.h:161
mod_biology::italk
integer italk
Definition fennel_mod.h:92
mod_biology::k_no3
real(r8), dimension(:), allocatable k_no3
Definition fennel_mod.h:133
mod_biology::inifx
integer inifx
Definition fennel_mod.h:118
mod_biology::wsdet
real(r8), dimension(:), allocatable wsdet
Definition fennel_mod.h:155
mod_biology::phymr
real(r8), dimension(:), allocatable phymr
Definition fennel_mod.h:145
mod_biology::phymin
real(r8), dimension(:), allocatable phymin
Definition fennel_mod.h:144
mod_biology::iphyt
integer iphyt
Definition fennel_mod.h:81
mod_biology::i_thnh4
real(r8), dimension(:), allocatable i_thnh4
Definition fennel_mod.h:131
mod_biology::lderrn
real(r8), dimension(:), allocatable lderrn
Definition fennel_mod.h:136
mod_biology::phycn
real(r8), dimension(:), allocatable phycn
Definition fennel_mod.h:140
mod_biology::pco2air
real(r8), dimension(:), allocatable pco2air
Definition fennel_mod.h:163
mod_biology::coagr
real(r8), dimension(:), allocatable coagr
Definition fennel_mod.h:129
mod_biology::inh4_
integer inh4_
Definition ecosim_mod.h:278
mod_biology::idbio
integer, dimension(:), allocatable idbio
Definition ecosim_mod.h:256
mod_biology::ioxyg
integer ioxyg
Definition fennel_mod.h:98
mod_biology::wphy
real(r8), dimension(:), allocatable wphy
Definition fennel_mod.h:154
mod_biology::zooer
real(r8), dimension(:), allocatable zooer
Definition fennel_mod.h:159
mod_biology::ino3u
integer ino3u
Definition fennel_mod.h:117
mod_biology::ippro
integer ippro
Definition fennel_mod.h:116
mod_biology::lderrc
real(r8), dimension(:), allocatable lderrc
Definition fennel_mod.h:137
mod_biology::attchl
real(r8), dimension(:), allocatable attchl
Definition fennel_mod.h:126
mod_biology::izoop
integer izoop
Definition fennel_mod.h:82
mod_biology::zoobm
real(r8), dimension(:), allocatable zoobm
Definition fennel_mod.h:157
mod_biology::zooae_n
real(r8), dimension(:), allocatable zooae_n
Definition fennel_mod.h:156
mod_biology::ildec
integer ildec
Definition fennel_mod.h:89
mod_biology::icofx
integer icofx
Definition fennel_mod.h:107
mod_biology::io2fx
integer io2fx
Definition fennel_mod.h:110
mod_biology::irdec
integer irdec
Definition fennel_mod.h:94
mod_biology::ichlo
integer ichlo
Definition fennel_mod.h:80
mod_diags
Definition mod_diags.F:2
mod_diags::diags
type(t_diags), dimension(:), allocatable diags
Definition mod_diags.F:100
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_kinds
Definition mod_kinds.F:2
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::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_param::nbt
integer nbt
Definition mod_param.F:509
mod_param::nt
integer, dimension(:), allocatable nt
Definition mod_param.F:489
mod_scalars
Definition mod_scalars.F:2
mod_scalars::nrrec
integer, dimension(:), allocatable nrrec
Definition mod_scalars.F:1171
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
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::tdays
real(dp), dimension(:), allocatable tdays
Definition mod_scalars.F:261
mod_scalars::sec2day
real(dp), parameter sec2day
Definition mod_scalars.F:817
mod_scalars::isalt
integer isalt
Definition mod_scalars.F:127
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::ntstart
integer, dimension(:), allocatable ntstart
Definition mod_scalars.F:378
mod_scalars::ndia
integer, dimension(:), allocatable ndia
Definition mod_scalars.F:1074
mod_scalars::ntsdia
integer, dimension(:), allocatable ntsdia
Definition mod_scalars.F:1068
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