Ticket #665: fennel.h

File fennel.h, 82.7 KB (added by m.hadfield, 11 years ago)

Corrected version of fennel.h

Line 
1 SUBROUTINE biology (ng,tile)
2!
3!svn $Id$
4!***********************************************************************
5! Copyright (c) 2002-2015 The ROMS/TOMS Group !
6! Licensed under a MIT/X style license Hernan G. Arango !
7! See License_ROMS.txt Katja Fennel !
8!****************************************** Alexander F. Shchepetkin ***
9! !
10! This routine computes the biological sources and sinks for the !
11! Fennel et at. (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's (1992) formula will be used. !
56! !
57! References: !
58! !
59! Fennel, K., Wilkin, J., Levin, J., Moisan, J., O'Reilly, J., !
60! Haidvogel, D., 2006: Nitrogen cycling in the Mid Atlantic !
61! Bight and implications for the North Atlantic nitrogen !
62! budget: Results from a three-dimensional model. Global !
63! Biogeochemical Cycles 20, GB3007, doi:10.1029/2005GB002456. !
64! !
65! Fennel, K., Wilkin, J., Previdi, M., Najjar, R. 2008: !
66! Denitrification effects on air-sea CO2 flux in the coastal !
67! ocean: Simulations for the Northwest North Atlantic. !
68! Geophys. Res. Letters 35, L24608, doi:10.1029/2008GL036147. !
69! !
70!***********************************************************************
71!
72 USE mod_param
73#ifdef DIAGNOSTICS_BIO
74 USE mod_diags
75#endif
76 USE mod_forces
77 USE mod_grid
78 USE mod_ncparam
79 USE mod_ocean
80 USE mod_stepping
81!
82 implicit none
83!
84! Imported variable declarations.
85!
86 integer, intent(in) :: ng, tile
87!
88! Local variable declarations.
89!
90#include "tile.h"
91!
92! Set header file name.
93!
94#ifdef DISTRIBUTE
95 IF (Lbiofile(iNLM)) THEN
96#else
97 IF (Lbiofile(iNLM).and.(tile.eq.0)) THEN
98#endif
99 Lbiofile(iNLM)=.FALSE.
100 BIONAME(iNLM)=__FILE__
101 END IF
102!
103#ifdef PROFILE
104 CALL wclock_on (ng, iNLM, 15)
105#endif
106 CALL biology_tile (ng, tile, &
107 & LBi, UBi, LBj, UBj, N(ng), NT(ng), &
108 & IminS, ImaxS, JminS, JmaxS, &
109 & nstp(ng), nnew(ng), &
110#ifdef MASKING
111 & GRID(ng) % rmask, &
112# if defined WET_DRY && defined DIAGNOSTICS_BIO
113 & GRID(ng) % rmask_full, &
114# endif
115#endif
116 & GRID(ng) % Hz, &
117 & GRID(ng) % z_r, &
118 & GRID(ng) % z_w, &
119 & FORCES(ng) % srflx, &
120#if defined CARBON || defined OXYGEN
121# ifdef BULK_FLUXES
122 & FORCES(ng) % Uwind, &
123 & FORCES(ng) % Vwind, &
124# else
125 & FORCES(ng) % sustr, &
126 & FORCES(ng) % svstr, &
127# endif
128#endif
129#ifdef CARBON
130 & OCEAN(ng) % pH, &
131#endif
132#ifdef DIAGNOSTICS_BIO
133 & DIAGS(ng) % DiaBio2d, &
134 & DIAGS(ng) % DiaBio3d, &
135#endif
136 & OCEAN(ng) % t)
137
138#ifdef PROFILE
139 CALL wclock_off (ng, iNLM, 15)
140#endif
141
142 RETURN
143 END SUBROUTINE biology
144!
145!-----------------------------------------------------------------------
146 SUBROUTINE biology_tile (ng, tile, &
147 & LBi, UBi, LBj, UBj, UBk, UBt, &
148 & IminS, ImaxS, JminS, JmaxS, &
149 & nstp, nnew, &
150#ifdef MASKING
151 & rmask, &
152# if defined WET_DRY && defined DIAGNOSTICS_BIO
153 & rmask_full, &
154# endif
155#endif
156 & Hz, z_r, z_w, srflx, &
157#if defined CARBON || defined OXYGEN
158# ifdef BULK_FLUXES
159 & Uwind, Vwind, &
160# else
161 & sustr, svstr, &
162# endif
163#endif
164#ifdef CARBON
165 & pH, &
166#endif
167#ifdef DIAGNOSTICS_BIO
168 & DiaBio2d, DiaBio3d, &
169#endif
170 & t)
171!-----------------------------------------------------------------------
172!
173 USE mod_param
174 USE mod_biology
175 USE mod_ncparam
176 USE mod_scalars
177!
178! Imported variable declarations.
179!
180 integer, intent(in) :: ng, tile
181 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk, UBt
182 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
183 integer, intent(in) :: nstp, nnew
184
185#ifdef ASSUMED_SHAPE
186# ifdef MASKING
187 real(r8), intent(in) :: rmask(LBi:,LBj:)
188# if defined WET_DRY && defined DIAGNOSTICS_BIO
189 real(r8), intent(in) :: rmask_full(LBi:,LBj:)
190# endif
191# endif
192 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
193 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
194 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
195 real(r8), intent(in) :: srflx(LBi:,LBj:)
196# if defined CARBON || defined OXYGEN
197# ifdef BULK_FLUXES
198 real(r8), intent(in) :: Uwind(LBi:,LBj:)
199 real(r8), intent(in) :: Vwind(LBi:,LBj:)
200# else
201 real(r8), intent(in) :: sustr(LBi:,LBj:)
202 real(r8), intent(in) :: svstr(LBi:,LBj:)
203# endif
204# endif
205# ifdef CARBON
206 real(r8), intent(inout) :: pH(LBi:,LBj:)
207# endif
208# ifdef DIAGNOSTICS_BIO
209 real(r8), intent(inout) :: DiaBio2d(LBi:,LBj:,:)
210 real(r8), intent(inout) :: DiaBio3d(LBi:,LBj:,:,:)
211# endif
212 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
213#else
214# ifdef MASKING
215 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
216# if defined WET_DRY && defined DIAGNOSTICS_BIO
217 real(r8), intent(in) :: rmask_full(LBi:UBi,LBj:UBj)
218# endif
219# endif
220 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,UBk)
221 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,UBk)
222 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:UBk)
223 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
224# if defined CARBON || defined OXYGEN
225# ifdef BULK_FLUXES
226 real(r8), intent(in) :: Uwind(LBi:UBi,LBj:UBj)
227 real(r8), intent(in) :: Vwind(LBi:UBi,LBj:UBj)
228# else
229 real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
230 real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
231# endif
232# endif
233# ifdef CARBON
234 real(r8), intent(inout) :: pH(LBi:UBi,LBj:UBj)
235# endif
236# ifdef DIAGNOSTICS_BIO
237 real(r8), intent(inout) :: DiaBio2d(LBi:UBi,LBj:UBj,NDbio2d)
238 real(r8), intent(inout) :: DiaBio3d(LBi:UBi,LBj:UBj,UBk,NDbio3d)
239# endif
240 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,UBk,3,UBt)
241#endif
242!
243! Local variable declarations.
244!
245#if defined CARBON || defined OXYGEN
246 real(r8) :: u10squ
247#endif
248#ifdef CARBON
249 integer, parameter :: Nsink = 6
250#else
251 integer, parameter :: Nsink = 4
252#endif
253
254 integer :: Iter, i, ibio, isink, itrc, ivar, j, k, ks
255
256 integer, dimension(Nsink) :: idsink
257
258 real(r8), parameter :: eps = 1.0e-20_r8
259
260#ifdef OXYGEN
261 real(r8), parameter :: OA0 = 2.00907_r8 ! Oxygen
262 real(r8), parameter :: OA1 = 3.22014_r8 ! saturation
263 real(r8), parameter :: OA2 = 4.05010_r8 ! coefficients
264 real(r8), parameter :: OA3 = 4.94457_r8
265 real(r8), parameter :: OA4 =-0.256847_r8
266 real(r8), parameter :: OA5 = 3.88767_r8
267 real(r8), parameter :: OB0 =-0.00624523_r8
268 real(r8), parameter :: OB1 =-0.00737614_r8
269 real(r8), parameter :: OB2 =-0.0103410_r8
270 real(r8), parameter :: OB3 =-0.00817083_r8
271 real(r8), parameter :: OC0 =-0.000000488682_r8
272 real(r8), parameter :: rOxNO3= 8.625_r8 ! 138/16
273 real(r8), parameter :: rOxNH4= 6.625_r8 ! 106/16
274 real(r8) :: l2mol = 1000.0_r8/22.3916_r8 ! liter to mol
275#endif
276#ifdef CARBON
277 integer :: iday, month, year
278
279 integer, parameter :: DoNewton = 0 ! pCO2 solver
280
281 real(r8), parameter :: Acoef = 2073.1_r8 ! Schmidt
282 real(r8), parameter :: Bcoef = 125.62_r8 ! number
283 real(r8), parameter :: Ccoef = 3.6276_r8 ! transfer
284 real(r8), parameter :: Dcoef = 0.043219_r8 ! coefficients
285
286 real(r8), parameter :: A1 = -60.2409_r8 ! surface
287 real(r8), parameter :: A2 = 93.4517_r8 ! CO2
288 real(r8), parameter :: A3 = 23.3585_r8 ! solubility
289 real(r8), parameter :: B1 = 0.023517_r8 ! coefficients
290 real(r8), parameter :: B2 = -0.023656_r8
291 real(r8), parameter :: B3 = 0.0047036_r8
292
293 real(r8) :: pmonth ! months since Jan 1951
294 real(r8) :: pCO2air_secular
295 real(r8) :: yday, hour
296
297 real(r8), parameter :: pi2 = 6.2831853071796_r8
298
299 real(r8), parameter :: D0 = 282.6_r8 ! coefficients
300 real(r8), parameter :: D1 = 0.125_r8 ! to calculate
301 real(r8), parameter :: D2 =-7.18_r8 ! secular trend in
302 real(r8), parameter :: D3 = 0.86_r8 ! atmospheric pCO2
303 real(r8), parameter :: D4 =-0.99_r8
304 real(r8), parameter :: D5 = 0.28_r8
305 real(r8), parameter :: D6 =-0.80_r8
306 real(r8), parameter :: D7 = 0.06_r8
307#endif
308
309 real(r8) :: Att, AttFac, ExpAtt, Itop, PAR
310 real(r8) :: Epp, L_NH4, L_NO3, LTOT, Vp
311 real(r8) :: Chl2C, dtdays, t_PPmax, inhNH4
312
313 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5
314 real(r8) :: fac1, fac2, fac3
315 real(r8) :: cffL, cffR, cu, dltL, dltR
316
317 real(r8) :: total_N
318
319#ifdef DIAGNOSTICS_BIO
320 real(r8) :: fiter
321#endif
322
323#ifdef OXYGEN
324 real(r8) :: SchmidtN_Ox, O2satu, O2_Flux
325 real(r8) :: TS, AA
326#endif
327
328#ifdef CARBON
329 real(r8) :: C_Flux_Remine
330 real(r8) :: CO2_Flux, CO2_sol, SchmidtN, TempK
331#endif
332
333 real(r8) :: N_Flux_Assim
334 real(r8) :: N_Flux_CoagD, N_Flux_CoagP
335 real(r8) :: N_Flux_Egest
336 real(r8) :: N_Flux_NewProd, N_Flux_RegProd
337 real(r8) :: N_Flux_Nitrifi
338 real(r8) :: N_Flux_Pmortal, N_Flux_Zmortal
339 real(r8) :: N_Flux_Remine
340 real(r8) :: N_Flux_Zexcret, N_Flux_Zmetabo
341
342 real(r8), dimension(Nsink) :: Wbio
343
344 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
345
346 real(r8), dimension(IminS:ImaxS) :: PARsur
347#ifdef CARBON
348 real(r8), dimension(IminS:ImaxS) :: pCO2
349#endif
350
351 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio
352 real(r8), dimension(IminS:ImaxS,N(ng),NT(ng)) :: Bio_old
353
354 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
355
356 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
357 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
358 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
359 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
360 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
361 real(r8), dimension(IminS:ImaxS,N(ng)) :: bL
362 real(r8), dimension(IminS:ImaxS,N(ng)) :: bR
363 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
364
365#include "set_bounds.h"
366#ifdef DIAGNOSTICS_BIO
367!
368!-----------------------------------------------------------------------
369! If appropriate, initialize time-averaged diagnostic arrays.
370!-----------------------------------------------------------------------
371!
372 IF (((iic(ng).gt.ntsDIA(ng)).and. &
373 & (MOD(iic(ng),nDIA(ng)).eq.1)).or. &
374 & ((iic(ng).ge.ntsDIA(ng)).and.(nDIA(ng).eq.1)).or. &
375 & ((nrrec(ng).gt.0).and.(iic(ng).eq.ntstart(ng)))) THEN
376 DO ivar=1,NDbio2d
377 DO j=Jstr,Jend
378 DO i=Istr,Iend
379 DiaBio2d(i,j,ivar)=0.0_r8
380 END DO
381 END DO
382 END DO
383 DO ivar=1,NDbio3d
384 DO k=1,N(ng)
385 DO j=Jstr,Jend
386 DO i=Istr,Iend
387 DiaBio3d(i,j,k,ivar)=0.0_r8
388 END DO
389 END DO
390 END DO
391 END DO
392 END IF
393#endif
394!
395!-----------------------------------------------------------------------
396! Add biological Source/Sink terms.
397!-----------------------------------------------------------------------
398!
399! Avoid computing source/sink terms if no biological iterations.
400!
401 IF (BioIter(ng).le.0) RETURN
402!
403! Set time-stepping according to the number of iterations.
404!
405 dtdays=dt(ng)*sec2day/REAL(BioIter(ng),r8)
406#ifdef DIAGNOSTICS_BIO
407!
408! A factor to account for the number of iterations in accumulating
409! diagnostic rate variables.
410!
411 fiter=1.0_r8/REAL(BioIter(ng),r8)
412#endif
413!
414! Set vertical sinking indentification vector.
415!
416 idsink(1)=iPhyt
417 idsink(2)=iChlo
418 idsink(3)=iSDeN
419 idsink(4)=iLDeN
420#ifdef CARBON
421 idsink(5)=iSDeC
422 idsink(6)=iLDeC
423#endif
424!
425! Set vertical sinking velocity vector in the same order as the
426! identification vector, IDSINK.
427!
428 Wbio(1)=wPhy(ng) ! phytoplankton
429 Wbio(2)=wPhy(ng) ! chlorophyll
430 Wbio(3)=wSDet(ng) ! small Nitrogen-detritus
431 Wbio(4)=wLDet(ng) ! large Nitrogen-detritus
432#ifdef CARBON
433 Wbio(5)=wSDet(ng) ! small Carbon-detritus
434 Wbio(6)=wLDet(ng) ! large Carbon-detritus
435#endif
436!
437! Compute inverse thickness to avoid repeated divisions.
438!
439 J_LOOP : DO j=Jstr,Jend
440 DO k=1,N(ng)
441 DO i=Istr,Iend
442 Hz_inv(i,k)=1.0_r8/Hz(i,j,k)
443 END DO
444 END DO
445 DO k=1,N(ng)-1
446 DO i=Istr,Iend
447 Hz_inv2(i,k)=1.0_r8/(Hz(i,j,k)+Hz(i,j,k+1))
448 END DO
449 END DO
450 DO k=2,N(ng)-1
451 DO i=Istr,Iend
452 Hz_inv3(i,k)=1.0_r8/(Hz(i,j,k-1)+Hz(i,j,k)+Hz(i,j,k+1))
453 END DO
454 END DO
455!
456! Extract biological variables from tracer arrays, place them into
457! scratch arrays, and restrict their values to be positive definite.
458! At input, all tracers (index nnew) from predictor step have
459! transport units (m Tunits) since we do not have yet the new
460! values for zeta and Hz. These are known after the 2D barotropic
461! time-stepping.
462!
463 DO itrc=1,NBT
464 ibio=idbio(itrc)
465 DO k=1,N(ng)
466 DO i=Istr,Iend
467 Bio_old(i,k,ibio)=MAX(0.0_r8,t(i,j,k,nstp,ibio))
468 Bio(i,k,ibio)=Bio_old(i,k,ibio)
469 END DO
470 END DO
471 END DO
472#ifdef CARBON
473 DO k=1,N(ng)
474 DO i=Istr,Iend
475 Bio_old(i,k,iTIC_)=MIN(Bio_old(i,k,iTIC_),3000.0_r8)
476 Bio_old(i,k,iTIC_)=MAX(Bio_old(i,k,iTIC_),400.0_r8)
477 Bio(i,k,iTIC_)=Bio_old(i,k,iTIC_)
478 END DO
479 END DO
480#endif
481!
482! Extract potential temperature and salinity.
483!
484 DO k=1,N(ng)
485 DO i=Istr,Iend
486 Bio(i,k,itemp)=MIN(t(i,j,k,nstp,itemp),35.0_r8)
487 Bio(i,k,isalt)=MAX(t(i,j,k,nstp,isalt), 0.0_r8)
488 END DO
489 END DO
490!
491! Calculate surface Photosynthetically Available Radiation (PAR). The
492! net shortwave radiation is scaled back to Watts/m2 and multiplied by
493! the fraction that is photosynthetically available, PARfrac.
494!
495 DO i=Istr,Iend
496 PARsur(i)=PARfrac(ng)*srflx(i,j)*rho0*Cp
497 END DO
498!
499!=======================================================================
500! Start internal iterations to achieve convergence of the nonlinear
501! backward-implicit solution.
502!=======================================================================
503!
504! During the iterative procedure a series of fractional time steps are
505! performed in a chained mode (splitting by different biological
506! conversion processes) in sequence of the main food chain. In all
507! stages the concentration of the component being consumed is treated
508! in fully implicit manner, so the algorithm guarantees non-negative
509! values, no matter how strong s the concentration of active consuming
510! component (Phytoplankton or Zooplankton). The overall algorithm,
511! as well as any stage of it, is formulated in conservative form
512! (except explicit sinking) in sense that the sum of concentration of
513! all components is conserved.
514!
515!
516! In the implicit algorithm, we have for example (N: nitrate,
517! P: phytoplankton),
518!
519! N(new) = N(old) - uptake * P(old) uptake = mu * N / (Kn + N)
520! {Michaelis-Menten}
521! below, we set
522! The N in the numerator of
523! cff = mu * P(old) / (Kn + N(old)) uptake is treated implicitly
524! as N(new)
525!
526! so the time-stepping of the equations becomes:
527!
528! N(new) = N(old) / (1 + cff) (1) when substracting a sink term,
529! consuming, divide by (1 + cff)
530! and
531!
532! P(new) = P(old) + cff * N(new) (2) when adding a source term,
533! growing, add (cff * source)
534!
535! Notice that if you substitute (1) in (2), you will get:
536!
537! P(new) = P(old) + cff * N(old) / (1 + cff) (3)
538!
539! If you add (1) and (3), you get
540!
541! N(new) + P(new) = N(old) + P(old)
542!
543! implying conservation regardless how "cff" is computed. Therefore,
544! this scheme is unconditionally stable regardless of the conversion
545! rate. It does not generate negative values since the constituent
546! to be consumed is always treated implicitly. It is also biased
547! toward damping oscillations.
548!
549! The iterative loop below is to iterate toward an universal Backward-
550! Euler treatment of all terms. So if there are oscillations in the
551! system, they are only physical oscillations. These iterations,
552! however, do not improve the accuaracy of the solution.
553!
554 ITER_LOOP: DO Iter=1,BioIter(ng)
555!
556!-----------------------------------------------------------------------
557! Light-limited computations.
558!-----------------------------------------------------------------------
559!
560! Compute attenuation coefficient based on the concentration of
561! chlorophyll-a within each grid box. Then, attenuate surface
562! photosynthetically available radiation (PARsur) down inot the
563! water column. Thus, PAR at certain depth depends on the whole
564! distribution of chlorophyll-a above.
565! To compute rate of maximum primary productivity (t_PPmax), one needs
566! PAR somewhat in the middle of the gridbox, so that attenuation "Att"
567! corresponds to half of the grid box height, while PAR is multiplied
568! by it twice: once to get it in the middle of grid-box and once the
569! compute on the lower grid-box interface.
570!
571 DO i=Istr,Iend
572 PAR=PARsur(i)
573 AttFac=0.0_r8
574 IF (PARsur(i).gt.0.0_r8) THEN
575 DO k=N(ng),1,-1
576!
577! Compute average light attenuation for each grid cell. To include
578! other attenuation contributions like suspended sediment or CDOM
579! modify AttFac.
580!
581 Att=(AttSW(ng)+ &
582 & AttChl(ng)*Bio(i,k,iChlo)+ &
583 & AttFac)* &
584 & (z_w(i,j,k)-z_w(i,j,k-1))
585 ExpAtt=EXP(-Att)
586 Itop=PAR
587 PAR=Itop*(1.0_r8-ExpAtt)/Att ! average at cell center
588!
589! Compute Chlorophyll-a phytoplankton ratio, [mg Chla / (mg C)].
590!
591 cff=PhyCN(ng)*12.0_r8
592 Chl2C=MIN(Bio(i,k,iChlo)/(Bio(i,k,iPhyt)*cff+eps), &
593 & Chl2C_m(ng))
594!
595! Temperature-limited and light-limited growth rate (Eppley, R.W.,
596! 1972, Fishery Bulletin, 70: 1063-1085; here 0.59=ln(2)*0.851).
597! Check value for Vp is 2.9124317 at 19.25 degC.
598!
599 Vp=Vp0(ng)*0.59_r8*(1.066_r8**Bio(i,k,itemp))
600 fac1=PAR*PhyIS(ng)
601 Epp=Vp/SQRT(Vp*Vp+fac1*fac1)
602 t_PPmax=Epp*fac1
603!
604! Nutrient-limitation terms (Parker 1993 Ecol Mod., 66, 113-120).
605!
606 cff1=Bio(i,k,iNH4_)*K_NH4(ng)
607 cff2=Bio(i,k,iNO3_)*K_NO3(ng)
608 inhNH4=1.0_r8/(1.0_r8+cff1)
609 L_NH4=cff1/(1.0_r8+cff1)
610 L_NO3=cff2*inhNH4/(1.0_r8+cff2)
611 LTOT=L_NO3+L_NH4
612!
613! Nitrate and ammonium uptake by Phytoplankton.
614!
615 fac1=dtdays*t_PPmax
616 cff4=fac1*K_NO3(ng)*inhNH4/(1.0_r8+cff2)*Bio(i,k,iPhyt)
617 cff5=fac1*K_NH4(ng)/(1.0_r8+cff1)*Bio(i,k,iPhyt)
618 Bio(i,k,iNO3_)=Bio(i,k,iNO3_)/(1.0_r8+cff4)
619 Bio(i,k,iNH4_)=Bio(i,k,iNH4_)/(1.0_r8+cff5)
620 N_Flux_NewProd=Bio(i,k,iNO3_)*cff4
621 N_Flux_RegProd=Bio(i,k,iNH4_)*cff5
622 Bio(i,k,iPhyt)=Bio(i,k,iPhyt)+ &
623 & N_Flux_NewProd+N_Flux_RegProd
624!
625 Bio(i,k,iChlo)=Bio(i,k,iChlo)+ &
626 & (dtdays*t_PPmax*t_PPmax*LTOT*LTOT* &
627 & Chl2C_m(ng)*Bio(i,k,iChlo))/ &
628 & (PhyIS(ng)*MAX(Chl2C,eps)*PAR+eps)
629#ifdef DIAGNOSTICS_BIO
630 DiaBio3d(i,j,k,iPPro)=DiaBio3d(i,j,k,iPPro)+ &
631# ifdef WET_DRY
632 & rmask_full(i,j)* &
633# endif
634 & (N_Flux_NewProd+N_Flux_RegProd)* &
635 & fiter
636 DiaBio3d(i,j,k,iNO3u)=DiaBio3d(i,j,k,iNO3u)+ &
637# ifdef WET_DRY
638 & rmask_full(i,j)* &
639# endif
640 & N_Flux_NewProd*fiter
641#endif
642#ifdef OXYGEN
643 Bio(i,k,iOxyg)=Bio(i,k,iOxyg)+ &
644 & N_Flux_NewProd*rOxNO3+ &
645 & N_Flux_RegProd*rOxNH4
646#endif
647#ifdef CARBON
648!
649! Total inorganic carbon (CO2) uptake during phytoplankton growth.
650!
651 cff1=PhyCN(ng)*(N_Flux_NewProd+N_Flux_RegProd)
652 Bio(i,k,iTIC_)=Bio(i,k,iTIC_)-cff1
653# ifdef TALK_NONCONSERV
654!
655! Account for the uptake of NO3 on total alkalinity.
656!
657 Bio(i,k,iTAlk)=Bio(i,k,iTAlk)+N_Flux_NewProd
658# endif
659#endif
660!
661! The Nitrification of NH4 ==> NO3 is thought to occur only in dark and
662! only in aerobic water (see Olson, R. J., 1981, JMR: (39), 227-238.).
663!
664! NH4+ + 3/2 O2 ==> NO2- + H2O; via Nitrosomonas bacteria
665! NO2- + 1/2 O2 ==> NO3- ; via Nitrobacter bacteria
666!
667! Note that the entire process has a total loss of two moles of O2 per
668! mole of NH4. If we were to resolve NO2 profiles, this is where we
669! would change the code to split out the differential effects of the
670! two different bacteria types. If OXYGEN is defined, nitrification is
671! inhibited at low oxygen concentrations using a Michaelis-Menten term.
672!
673#ifdef OXYGEN
674 fac2=MAX(Bio(i,k,iOxyg),0.0_r8) ! O2 max
675 fac3=MAX(fac2/(3.0_r8+fac2),0.0_r8) ! MM for O2 dependence
676 fac1=dtdays*NitriR(ng)*fac3
677#else
678 fac1=dtdays*NitriR(ng)
679#endif
680 cff1=(PAR-I_thNH4(ng))/ &
681 & (D_p5NH4(ng)+PAR-2.0_r8*I_thNH4(ng))
682 cff2=1.0_r8-MAX(0.0_r8,cff1)
683 cff3=fac1*cff2
684 Bio(i,k,iNH4_)=Bio(i,k,iNH4_)/(1.0_r8+cff3)
685 N_Flux_Nitrifi=Bio(i,k,iNH4_)*cff3
686 Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+N_Flux_Nitrifi
687#ifdef OXYGEN
688 Bio(i,k,iOxyg)=Bio(i,k,iOxyg)-2.0_r8*N_Flux_Nitrifi
689#endif
690#if defined CARBON && defined TALK_NONCONSERV
691 Bio(i,k,iTAlk)=Bio(i,k,iTAlk)-N_Flux_Nitrifi
692#endif
693!
694! Light attenuation at the bottom of the grid cell. It is the starting
695! PAR value for the next (deeper) vertical grid cell.
696!
697 PAR=Itop*ExpAtt
698 END DO
699!
700! If PARsur=0, nitrification occurs at the maximum rate (NitriR).
701!
702 ELSE
703 cff3=dtdays*NitriR(ng)
704 DO k=N(ng),1,-1
705 Bio(i,k,iNH4_)=Bio(i,k,iNH4_)/(1.0_r8+cff3)
706 N_Flux_Nitrifi=Bio(i,k,iNH4_)*cff3
707 Bio(i,k,iNO3_)=Bio(i,k,iNO3_)+N_Flux_Nitrifi
708#ifdef OXYGEN
709 Bio(i,k,iOxyg)=Bio(i,k,iOxyg)-2.0_r8*N_Flux_Nitrifi
710#endif
711#if defined CARBON && defined TALK_NONCONSERV
712 Bio(i,k,iTAlk)=Bio(i,k,iTAlk)-N_Flux_Nitrifi
713#endif
714 END DO
715 END IF
716 END DO
717!
718!-----------------------------------------------------------------------
719! Phytoplankton grazing by zooplankton (rate: ZooGR), phytoplankton
720! assimilated to zooplankton (fraction: ZooAE_N) and egested to small
721! detritus, and phytoplankton mortality (rate: PhyMR) to small
722! detritus. [Landry 1993 L&O 38:468-472]
723!-----------------------------------------------------------------------
724!
725 fac1=dtdays*ZooGR(ng)
726 cff2=dtdays*PhyMR(ng)
727 DO k=1,N(ng)
728 DO i=Istr,Iend
729!
730! Phytoplankton grazing by zooplankton.
731!
732 cff1=fac1*Bio(i,k,iZoop)*Bio(i,k,iPhyt)/ &
733 & (K_Phy(ng)+Bio(i,k,iPhyt)*Bio(i,k,iPhyt))
734 cff3=1.0_r8/(1.0_r8+cff1)
735 Bio(i,k,iPhyt)=cff3*Bio(i,k,iPhyt)
736 Bio(i,k,iChlo)=cff3*Bio(i,k,iChlo)
737!
738! Phytoplankton assimilated to zooplankton and egested to small
739! detritus.
740!
741 N_Flux_Assim=cff1*Bio(i,k,iPhyt)*ZooAE_N(ng)
742 N_Flux_Egest=Bio(i,k,iPhyt)*cff1*(1.0_r8-ZooAE_N(ng))
743 Bio(i,k,iZoop)=Bio(i,k,iZoop)+ &
744 & N_Flux_Assim
745 Bio(i,k,iSDeN)=Bio(i,k,iSDeN)+ &
746 & N_Flux_Egest
747!
748! Phytoplankton mortality (limited by a phytoplankton minimum).
749!
750 N_Flux_Pmortal=cff2*MAX(Bio(i,k,iPhyt)-PhyMin(ng),0.0_r8)
751 Bio(i,k,iPhyt)=Bio(i,k,iPhyt)-N_Flux_Pmortal
752 Bio(i,k,iChlo)=Bio(i,k,iChlo)- &
753 & cff2*MAX(Bio(i,k,iChlo)-ChlMin(ng),0.0_r8)
754 Bio(i,k,iSDeN)=Bio(i,k,iSDeN)+ &
755 & N_Flux_Pmortal
756#ifdef CARBON
757 Bio(i,k,iSDeC)=Bio(i,k,iSDeC)+ &
758 & PhyCN(ng)*(N_Flux_Egest+N_Flux_Pmortal)+ &
759 & (PhyCN(ng)-ZooCN(ng))*N_Flux_Assim
760#endif
761 END DO
762 END DO
763!
764!-----------------------------------------------------------------------
765! Zooplankton basal metabolism to NH4 (rate: ZooBM), zooplankton
766! mortality to small detritus (rate: ZooMR), zooplankton ingestion
767! related excretion (rate: ZooER).
768!-----------------------------------------------------------------------
769!
770 cff1=dtdays*ZooBM(ng)
771 fac2=dtdays*ZooMR(ng)
772 fac3=dtdays*ZooER(ng)
773 DO k=1,N(ng)
774 DO i=Istr,Iend
775 fac1=fac3*Bio(i,k,iPhyt)*Bio(i,k,iPhyt)/ &
776 & (K_Phy(ng)+Bio(i,k,iPhyt)*Bio(i,k,iPhyt))
777 cff2=fac2*Bio(i,k,iZoop)
778 cff3=fac1*ZooAE_N(ng)
779 Bio(i,k,iZoop)=Bio(i,k,iZoop)/ &
780 & (1.0_r8+cff2+cff3)
781!
782! Zooplankton mortality and excretion.
783!
784 N_Flux_Zmortal=cff2*Bio(i,k,iZoop)
785 N_Flux_Zexcret=cff3*Bio(i,k,iZoop)
786 Bio(i,k,iNH4_)=Bio(i,k,iNH4_)+N_Flux_Zexcret
787 Bio(i,k,iSDeN)=Bio(i,k,iSDeN)+N_Flux_Zmortal
788!
789! Zooplankton basal metabolism (limited by a zooplankton minimum).
790!
791 N_Flux_Zmetabo=cff1*MAX(Bio(i,k,iZoop)-ZooMin(ng),0.0_r8)
792 Bio(i,k,iZoop)=Bio(i,k,iZoop)-N_Flux_Zmetabo
793 Bio(i,k,iNH4_)=Bio(i,k,iNH4_)+N_Flux_Zmetabo
794#ifdef OXYGEN
795 Bio(i,k,iOxyg)=Bio(i,k,iOxyg)- &
796 & rOxNH4*(N_Flux_Zmetabo+N_Flux_Zexcret)
797#endif
798#ifdef CARBON
799 Bio(i,k,iSDeC)=Bio(i,k,iSDeC)+ &
800 & ZooCN(ng)*N_Flux_Zmortal
801 Bio(i,k,iTIC_)=Bio(i,k,iTIC_)+ &
802 & ZooCN(ng)*(N_Flux_Zmetabo+N_Flux_Zexcret)
803#endif
804 END DO
805 END DO
806!
807!-----------------------------------------------------------------------
808! Coagulation of phytoplankton and small detritus to large detritus.
809!-----------------------------------------------------------------------
810!
811 fac1=dtdays*CoagR(ng)
812 DO k=1,N(ng)
813 DO i=Istr,Iend
814 cff1=fac1*(Bio(i,k,iSDeN)+Bio(i,k,iPhyt))
815 cff2=1.0_r8/(1.0_r8+cff1)
816 Bio(i,k,iPhyt)=Bio(i,k,iPhyt)*cff2
817 Bio(i,k,iChlo)=Bio(i,k,iChlo)*cff2
818 Bio(i,k,iSDeN)=Bio(i,k,iSDeN)*cff2
819 N_Flux_CoagP=Bio(i,k,iPhyt)*cff1
820 N_Flux_CoagD=Bio(i,k,iSDeN)*cff1
821 Bio(i,k,iLDeN)=Bio(i,k,iLDeN)+ &
822 & N_Flux_CoagP+N_Flux_CoagD
823#ifdef CARBON
824 Bio(i,k,iSDeC)=Bio(i,k,iSDeC)-PhyCN(ng)*N_Flux_CoagD
825 Bio(i,k,iLDeC)=Bio(i,k,iLDeC)+ &
826 & PhyCN(ng)*(N_Flux_CoagP+N_Flux_CoagD)
827#endif
828 END DO
829 END DO
830!
831!-----------------------------------------------------------------------
832! Detritus recycling to NH4, remineralization.
833!-----------------------------------------------------------------------
834!
835#ifdef OXYGEN
836 DO k=1,N(ng)
837 DO i=Istr,Iend
838 fac1=MAX(Bio(i,k,iOxyg)-6.0_r8,0.0_r8) ! O2 off max
839 fac2=MAX(fac1/(3.0_r8+fac1),0.0_r8) ! MM for O2 dependence
840 cff1=dtdays*SDeRRN(ng)*fac2
841 cff2=1.0_r8/(1.0_r8+cff1)
842 cff3=dtdays*LDeRRN(ng)*fac2
843 cff4=1.0_r8/(1.0_r8+cff3)
844 Bio(i,k,iSDeN)=Bio(i,k,iSDeN)*cff2
845 Bio(i,k,iLDeN)=Bio(i,k,iLDeN)*cff4
846 N_Flux_Remine=Bio(i,k,iSDeN)*cff1+Bio(i,k,iLDeN)*cff3
847 Bio(i,k,iNH4_)=Bio(i,k,iNH4_)+N_Flux_Remine
848 Bio(i,k,iOxyg)=Bio(i,k,iOxyg)-N_Flux_Remine*rOxNH4
849 END DO
850 END DO
851#else
852 cff1=dtdays*SDeRRN(ng)
853 cff2=1.0_r8/(1.0_r8+cff1)
854 cff3=dtdays*LDeRRN(ng)
855 cff4=1.0_r8/(1.0_r8+cff3)
856 DO k=1,N(ng)
857 DO i=Istr,Iend
858 Bio(i,k,iSDeN)=Bio(i,k,iSDeN)*cff2
859 Bio(i,k,iLDeN)=Bio(i,k,iLDeN)*cff4
860 N_Flux_Remine=Bio(i,k,iSDeN)*cff1+Bio(i,k,iLDeN)*cff3
861 Bio(i,k,iNH4_)=Bio(i,k,iNH4_)+N_Flux_Remine
862 END DO
863 END DO
864#endif
865#ifdef OXYGEN
866!
867!-----------------------------------------------------------------------
868! Surface O2 gas exchange.
869!-----------------------------------------------------------------------
870!
871! Compute surface O2 gas exchange.
872!
873 cff1=rho0*550.0_r8
874 cff2=dtdays*0.31_r8*24.0_r8/100.0_r8
875 k=N(ng)
876 DO i=Istr,Iend
877!
878! Compute O2 transfer velocity : u10squared (u10 in m/s)
879!
880# ifdef BULK_FLUXES
881 u10squ=Uwind(i,j)*Uwind(i,j)+Vwind(i,j)*Vwind(i,j)
882# else
883 u10squ=cff1*SQRT((0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2+ &
884 & (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2)
885# endif
886# ifdef OCMIP_OXYGEN_SC
887!
888! Alternative formulation for Schmidt number (Sc will be slightly
889! smaller up to about 35 C): Compute the Schmidt number of oxygen
890! in seawater using the formulation proposed by Keeling et al.
891! (1998, Global Biogeochem. Cycles, 12, 141-163). Input temperature
892! in Celsius.
893!
894 SchmidtN_Ox=1638.0_r8- &
895 & Bio(i,k,itemp)*(81.83_r8- &
896 & Bio(i,k,itemp)* &
897 & (1.483_r8- &
898 & Bio(i,k,itemp)*0.008004_r8))
899# else
900!
901! Calculate the Schmidt number for O2 in sea water (Wanninkhof, 1992).
902!
903 SchmidtN_Ox=1953.4_r8- &
904 & Bio(i,k,itemp)*(128.0_r8- &
905 & Bio(i,k,itemp)* &
906 & (3.9918_r8- &
907 & Bio(i,k,itemp)*0.050091_r8))
908# endif
909
910 cff3=cff2*u10squ*SQRT(660.0_r8/SchmidtN_Ox)
911!
912! Calculate O2 saturation concentration using Garcia and Gordon
913! L&O (1992) formula, (EXP(AA) is in ml/l).
914!
915 TS=LOG((298.15_r8-Bio(i,k,itemp))/ &
916 & (273.15_r8+Bio(i,k,itemp)))
917 AA=OA0+TS*(OA1+TS*(OA2+TS*(OA3+TS*(OA4+TS*OA5))))+ &
918 & Bio(i,k,isalt)*(OB0+TS*(OB1+TS*(OB2+TS*OB3)))+ &
919 & OC0*Bio(i,k,isalt)*Bio(i,k,isalt)
920!
921! Convert from ml/l to mmol/m3.
922!
923 O2satu=l2mol*EXP(AA)
924!
925! Add in O2 gas exchange.
926!
927 O2_Flux=cff3*(O2satu-Bio(i,k,iOxyg))
928 Bio(i,k,iOxyg)=Bio(i,k,iOxyg)+ &
929 & O2_Flux*Hz_inv(i,k)
930# ifdef DIAGNOSTICS_BIO
931 DiaBio2d(i,j,iO2fx)=DiaBio2d(i,j,iO2fx)+ &
932# ifdef WET_DRY
933 & rmask_full(i,j)* &
934# endif
935 & O2_Flux*fiter
936# endif
937
938 END DO
939#endif
940
941#ifdef CARBON
942!
943!-----------------------------------------------------------------------
944! Allow different remineralization rates for detrital C and detrital N.
945!-----------------------------------------------------------------------
946!
947 cff1=dtdays*SDeRRC(ng)
948 cff2=1.0_r8/(1.0_r8+cff1)
949 cff3=dtdays*LDeRRC(ng)
950 cff4=1.0_r8/(1.0_r8+cff3)
951 DO k=1,N(ng)
952 DO i=Istr,Iend
953 Bio(i,k,iSDeC)=Bio(i,k,iSDeC)*cff2
954 Bio(i,k,iLDeC)=Bio(i,k,iLDeC)*cff4
955 C_Flux_Remine=Bio(i,k,iSDeC)*cff1+Bio(i,k,iLDeC)*cff3
956 Bio(i,k,iTIC_)=Bio(i,k,iTIC_)+C_Flux_Remine
957 END DO
958 END DO
959!
960!-----------------------------------------------------------------------
961! Surface CO2 gas exchange.
962!-----------------------------------------------------------------------
963!
964! Compute equilibrium partial pressure inorganic carbon (ppmv) at the
965! surface.
966!
967 k=N(ng)
968# ifdef pCO2_RZ
969 CALL pCO2_water_RZ (Istr, Iend, LBi, UBi, LBj, UBj, &
970 & IminS, ImaxS, j, DoNewton, &
971# ifdef MASKING
972 & rmask, &
973# endif
974 & Bio(IminS:,k,itemp), Bio(IminS:,k,isalt), &
975 & Bio(IminS:,k,iTIC_), Bio(IminS:,k,iTAlk), &
976 & pH, pCO2)
977# else
978 CALL pCO2_water (Istr, Iend, LBi, UBi, LBj, UBj, &
979 & IminS, ImaxS, j, DoNewton, &
980# ifdef MASKING
981 & rmask, &
982# endif
983 & Bio(IminS:,k,itemp), Bio(IminS:,k,isalt), &
984 & Bio(IminS:,k,iTIC_), Bio(IminS:,k,iTAlk), &
985 & 0.0_r8, 0.0_r8, pH, pCO2)
986# endif
987!
988! Compute surface CO2 gas exchange.
989!
990 cff1=rho0*550.0_r8
991 cff2=dtdays*0.31_r8*24.0_r8/100.0_r8
992 DO i=Istr,Iend
993!
994! Compute CO2 transfer velocity : u10squared (u10 in m/s)
995!
996# ifdef BULK_FLUXES
997 u10squ=Uwind(i,j)**2+Vwind(i,j)**2
998# else
999 u10squ=cff1*SQRT((0.5_r8*(sustr(i,j)+sustr(i+1,j)))**2+ &
1000 & (0.5_r8*(svstr(i,j)+svstr(i,j+1)))**2)
1001# endif
1002 SchmidtN=Acoef- &
1003 & Bio(i,k,itemp)*(Bcoef- &
1004 & Bio(i,k,itemp)*(Ccoef- &
1005 & Bio(i,k,itemp)*Dcoef))
1006 cff3=cff2*u10squ*SQRT(660.0_r8/SchmidtN)
1007!
1008! Calculate CO2 solubility [mol/(kg.atm)] using Weiss (1974) formula.
1009!
1010 TempK=0.01_r8*(Bio(i,k,itemp)+273.15_r8)
1011 CO2_sol=EXP(A1+ &
1012 & A2/TempK+ &
1013 & A3*LOG(TempK)+ &
1014 & Bio(i,k,isalt)*(B1+TempK*(B2+B3*TempK)))
1015!
1016! Add in CO2 gas exchange.
1017!
1018 CALL caldate (r_date, tdays(ng), year, yday, month, iday, &
1019 & hour)
1020 pmonth=2003.0_r8-1951.0_r8+yday/365.0_r8
1021!! pCO2air_secular=D0+D1*pmonth*12.0_r8+ &
1022!! & D2*SIN(pi2*pmonth+D3)+ &
1023!! & D4*SIN(pi2*pmonth+D5)+ &
1024!! & D6*SIN(pi2*pmonth+D7)
1025!! CO2_Flux=cff3*CO2_sol*(pCO2air_secular-pCO2(i))
1026 CO2_Flux=cff3*CO2_sol*(pCO2air(ng)-pCO2(i))
1027 Bio(i,k,iTIC_)=Bio(i,k,iTIC_)+ &
1028 & CO2_Flux*Hz_inv(i,k)
1029# ifdef DIAGNOSTICS_BIO
1030 DiaBio2d(i,j,iCOfx)=DiaBio2d(i,j,iCOfx)+ &
1031# ifdef WET_DRY
1032 & rmask_full(i,j)* &
1033# endif
1034 & CO2_Flux*fiter
1035 DiaBio2d(i,j,ipCO2)=pCO2(i)
1036# ifdef WET_DRY
1037 DiaBio2d(i,j,ipCO2)=DiaBio2d(i,j,ipCO2)*rmask_full(i,j)
1038# endif
1039# endif
1040 END DO
1041#endif
1042!
1043!-----------------------------------------------------------------------
1044! Vertical sinking terms.
1045!-----------------------------------------------------------------------
1046!
1047! Reconstruct vertical profile of selected biological constituents
1048! "Bio(:,:,isink)" in terms of a set of parabolic segments within each
1049! grid box. Then, compute semi-Lagrangian flux due to sinking.
1050!
1051 SINK_LOOP: DO isink=1,Nsink
1052 ibio=idsink(isink)
1053!
1054! Copy concentration of biological particulates into scratch array
1055! "qc" (q-central, restrict it to be positive) which is hereafter
1056! interpreted as a set of grid-box averaged values for biogeochemical
1057! constituent concentration.
1058!
1059 DO k=1,N(ng)
1060 DO i=Istr,Iend
1061 qc(i,k)=Bio(i,k,ibio)
1062 END DO
1063 END DO
1064!
1065 DO k=N(ng)-1,1,-1
1066 DO i=Istr,Iend
1067 FC(i,k)=(qc(i,k+1)-qc(i,k))*Hz_inv2(i,k)
1068 END DO
1069 END DO
1070 DO k=2,N(ng)-1
1071 DO i=Istr,Iend
1072 dltR=Hz(i,j,k)*FC(i,k)
1073 dltL=Hz(i,j,k)*FC(i,k-1)
1074 cff=Hz(i,j,k-1)+2.0_r8*Hz(i,j,k)+Hz(i,j,k+1)
1075 cffR=cff*FC(i,k)
1076 cffL=cff*FC(i,k-1)
1077!
1078! Apply PPM monotonicity constraint to prevent oscillations within the
1079! grid box.
1080!
1081 IF ((dltR*dltL).le.0.0_r8) THEN
1082 dltR=0.0_r8
1083 dltL=0.0_r8
1084 ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN
1085 dltR=cffL
1086 ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN
1087 dltL=cffR
1088 END IF
1089!
1090! Compute right and left side values (bR,bL) of parabolic segments
1091! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
1092!
1093! NOTE: Although each parabolic segment is monotonic within its grid
1094! box, monotonicity of the whole profile is not guaranteed,
1095! because bL(k+1)-bR(k) may still have different sign than
1096! qc(i,k+1)-qc(i,k). This possibility is excluded,
1097! after bL and bR are reconciled using WENO procedure.
1098!
1099 cff=(dltR-dltL)*Hz_inv3(i,k)
1100 dltR=dltR-cff*Hz(i,j,k+1)
1101 dltL=dltL+cff*Hz(i,j,k-1)
1102 bR(i,k)=qc(i,k)+dltR
1103 bL(i,k)=qc(i,k)-dltL
1104 WR(i,k)=(2.0_r8*dltR-dltL)**2
1105 WL(i,k)=(dltR-2.0_r8*dltL)**2
1106 END DO
1107 END DO
1108 cff=1.0E-14_r8
1109 DO k=2,N(ng)-2
1110 DO i=Istr,Iend
1111 dltL=MAX(cff,WL(i,k ))
1112 dltR=MAX(cff,WR(i,k+1))
1113 bR(i,k)=(dltR*bR(i,k)+dltL*bL(i,k+1))/(dltR+dltL)
1114 bL(i,k+1)=bR(i,k)
1115 END DO
1116 END DO
1117 DO i=Istr,Iend
1118 FC(i,N(ng))=0.0_r8 ! NO-flux boundary condition
1119#if defined LINEAR_CONTINUATION
1120 bL(i,N(ng))=bR(i,N(ng)-1)
1121 bR(i,N(ng))=2.0_r8*qc(i,N(ng))-bL(i,N(ng))
1122#elif defined NEUMANN
1123 bL(i,N(ng))=bR(i,N(ng)-1)
1124 bR(i,N(ng))=1.5_r8*qc(i,N(ng))-0.5_r8*bL(i,N(ng))
1125#else
1126 bR(i,N(ng))=qc(i,N(ng)) ! default strictly monotonic
1127 bL(i,N(ng))=qc(i,N(ng)) ! conditions
1128 bR(i,N(ng)-1)=qc(i,N(ng))
1129#endif
1130#if defined LINEAR_CONTINUATION
1131 bR(i,1)=bL(i,2)
1132 bL(i,1)=2.0_r8*qc(i,1)-bR(i,1)
1133#elif defined NEUMANN
1134 bR(i,1)=bL(i,2)
1135 bL(i,1)=1.5_r8*qc(i,1)-0.5_r8*bR(i,1)
1136#else
1137 bL(i,2)=qc(i,1) ! bottom grid boxes are
1138 bR(i,1)=qc(i,1) ! re-assumed to be
1139 bL(i,1)=qc(i,1) ! piecewise constant.
1140#endif
1141 END DO
1142!
1143! Apply monotonicity constraint again, since the reconciled interfacial
1144! values may cause a non-monotonic behavior of the parabolic segments
1145! inside the grid box.
1146!
1147 DO k=1,N(ng)
1148 DO i=Istr,Iend
1149 dltR=bR(i,k)-qc(i,k)
1150 dltL=qc(i,k)-bL(i,k)
1151 cffR=2.0_r8*dltR
1152 cffL=2.0_r8*dltL
1153 IF ((dltR*dltL).lt.0.0_r8) THEN
1154 dltR=0.0_r8
1155 dltL=0.0_r8
1156 ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN
1157 dltR=cffL
1158 ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN
1159 dltL=cffR
1160 END IF
1161 bR(i,k)=qc(i,k)+dltR
1162 bL(i,k)=qc(i,k)-dltL
1163 END DO
1164 END DO
1165!
1166! After this moment reconstruction is considered complete. The next
1167! stage is to compute vertical advective fluxes, FC. It is expected
1168! that sinking may occurs relatively fast, the algorithm is designed
1169! to be free of CFL criterion, which is achieved by allowing
1170! integration bounds for semi-Lagrangian advective flux to use as
1171! many grid boxes in upstream direction as necessary.
1172!
1173! In the two code segments below, WL is the z-coordinate of the
1174! departure point for grid box interface z_w with the same indices;
1175! FC is the finite volume flux; ksource(:,k) is index of vertical
1176! grid box which contains the departure point (restricted by N(ng)).
1177! During the search: also add in content of whole grid boxes
1178! participating in FC.
1179!
1180 cff=dtdays*ABS(Wbio(isink))
1181 DO k=1,N(ng)
1182 DO i=Istr,Iend
1183 FC(i,k-1)=0.0_r8
1184 WL(i,k)=z_w(i,j,k-1)+cff
1185 WR(i,k)=Hz(i,j,k)*qc(i,k)
1186 ksource(i,k)=k
1187 END DO
1188 END DO
1189 DO k=1,N(ng)
1190 DO ks=k,N(ng)-1
1191 DO i=Istr,Iend
1192 IF (WL(i,k).gt.z_w(i,j,ks)) THEN
1193 ksource(i,k)=ks+1
1194 FC(i,k-1)=FC(i,k-1)+WR(i,ks)
1195 END IF
1196 END DO
1197 END DO
1198 END DO
1199!
1200! Finalize computation of flux: add fractional part.
1201!
1202 DO k=1,N(ng)
1203 DO i=Istr,Iend
1204 ks=ksource(i,k)
1205 cu=MIN(1.0_r8,(WL(i,k)-z_w(i,j,ks-1))*Hz_inv(i,ks))
1206 FC(i,k-1)=FC(i,k-1)+ &
1207 & Hz(i,j,ks)*cu* &
1208 & (bL(i,ks)+ &
1209 & cu*(0.5_r8*(bR(i,ks)-bL(i,ks))- &
1210 & (1.5_r8-cu)* &
1211 & (bR(i,ks)+bL(i,ks)- &
1212 & 2.0_r8*qc(i,ks))))
1213 END DO
1214 END DO
1215 DO k=1,N(ng)
1216 DO i=Istr,Iend
1217 Bio(i,k,ibio)=qc(i,k)+(FC(i,k)-FC(i,k-1))*Hz_inv(i,k)
1218 END DO
1219 END DO
1220
1221#ifdef BIO_SEDIMENT
1222!
1223! Particulate flux reaching the seafloor is remineralized and returned
1224! to the dissolved nitrate pool. Without this conversion, particulate
1225! material falls out of the system. This is a temporary fix to restore
1226! total nitrogen conservation. It will be replaced later by a
1227! parameterization that includes the time delay of remineralization
1228! and dissolved oxygen.
1229!
1230 cff2=4.0_r8/16.0_r8
1231# ifdef OXYGEN
1232 cff3=115.0_r8/16.0_r8
1233 cff4=106.0_r8/16.0_r8
1234# endif
1235 IF ((ibio.eq.iPhyt).or. &
1236 & (ibio.eq.iSDeN).or. &
1237 & (ibio.eq.iLDeN)) THEN
1238 DO i=Istr,Iend
1239 cff1=FC(i,0)*Hz_inv(i,1)
1240# ifdef DENITRIFICATION
1241 Bio(i,1,iNH4_)=Bio(i,1,iNH4_)+cff1*cff2
1242# ifdef DIAGNOSTICS_BIO
1243 DiaBio2d(i,j,iDNIT)=DiaBio2d(i,j,iDNIT)+ &
1244# ifdef WET_DRY
1245 & rmask_full(i,j)* &
1246# endif
1247 & (1.0_r8-cff2)*cff1*Hz(i,j,1)*fiter
1248# endif
1249# ifdef OXYGEN
1250 Bio(i,1,iOxyg)=Bio(i,1,iOxyg)-cff1*cff3
1251# endif
1252# else
1253 Bio(i,1,iNH4_)=Bio(i,1,iNH4_)+cff1
1254# ifdef OXYGEN
1255 Bio(i,1,iOxyg)=Bio(i,1,iOxyg)-cff1*cff4
1256# endif
1257# endif
1258 END DO
1259 END IF
1260# ifdef CARBON
1261# ifdef DENITRIFICATION
1262 cff3=12.0_r8
1263 cff4=0.74_r8
1264# endif
1265 IF ((ibio.eq.iSDeC).or. &
1266 & (ibio.eq.iLDeC))THEN
1267 DO i=Istr,Iend
1268 cff1=FC(i,0)*Hz_inv(i,1)
1269 Bio(i,1,iTIC_)=Bio(i,1,iTIC_)+cff1
1270 END DO
1271 END IF
1272 IF (ibio.eq.iPhyt)THEN
1273 DO i=Istr,Iend
1274 cff1=FC(i,0)*Hz_inv(i,1)
1275 Bio(i,1,iTIC_)=Bio(i,1,iTIC_)+cff1*PhyCN(ng)
1276 END DO
1277 END IF
1278# endif
1279#endif
1280 END DO SINK_LOOP
1281 END DO ITER_LOOP
1282!
1283!-----------------------------------------------------------------------
1284! Update global tracer variables: Add increment due to BGC processes
1285! to tracer array in time index "nnew". Index "nnew" is solution after
1286! advection and mixing and has transport units (m Tunits) hence the
1287! increment is multiplied by Hz. Notice that we need to subtract
1288! original values "Bio_old" at the top of the routine to just account
1289! for the concentractions affected by BGC processes. This also takes
1290! into account any constraints (non-negative concentrations, carbon
1291! concentration range) specified before entering BGC kernel. If "Bio"
1292! were unchanged by BGC processes, the increment would be exactly
1293! zero. Notice that final tracer values, t(:,:,:,nnew,:) are not
1294! bounded >=0 so that we can preserve total inventory of N and
1295! C even when advection causes tracer concentration to go negative.
1296! (J. Wilkin and H. Arango, Apr 27, 2012)
1297!-----------------------------------------------------------------------
1298!
1299 DO itrc=1,NBT
1300 ibio=idbio(itrc)
1301 DO k=1,N(ng)
1302 DO i=Istr,Iend
1303 cff=Bio(i,k,ibio)-Bio_old(i,k,ibio)
1304 t(i,j,k,nnew,ibio)=t(i,j,k,nnew,ibio)+cff*Hz(i,j,k)
1305 END DO
1306 END DO
1307 END DO
1308 END DO J_LOOP
1309
1310 RETURN
1311 END SUBROUTINE biology_tile
1312
1313#ifdef CARBON
1314# ifdef pCO2_RZ
1315 SUBROUTINE pCO2_water_RZ (Istr, Iend, &
1316 & LBi, UBi, LBj, UBj, IminS, ImaxS, &
1317 & j, DoNewton, &
1318# ifdef MASKING
1319 & rmask, &
1320# endif
1321 & T, S, TIC, TAlk, pH, pCO2)
1322!
1323!***********************************************************************
1324! !
1325! This routine computes equilibrium partial pressure of CO2 (pCO2) !
1326! in the surface seawater. !
1327! !
1328! On Input: !
1329! !
1330! Istr Starting tile index in the I-direction. !
1331! Iend Ending tile index in the I-direction. !
1332! LBi I-dimension lower bound. !
1333! UBi I-dimension upper bound. !
1334! LBj J-dimension lower bound. !
1335! UBj J-dimension upper bound. !
1336! IminS I-dimension lower bound for private arrays. !
1337! ImaxS I-dimension upper bound for private arrays. !
1338! j j-pipelined index. !
1339! DoNewton Iteration solver: !
1340! [0] Bracket and bisection. !
1341! [1] Newton-Raphson method. !
1342! rmask Land/Sea masking. !
1343! T Surface temperature (Celsius). !
1344! S Surface salinity (PSS). !
1345! TIC Total inorganic carbon (millimol/m3). !
1346! TAlk Total alkalinity (milli-equivalents/m3). !
1347! pH Best pH guess. !
1348! !
1349! On Output: !
1350! !
1351! pCO2 partial pressure of CO2 (ppmv). !
1352! !
1353! Check Value: (T=24, S=36.6, TIC=2040, TAlk=2390, PO4=0, !
1354! SiO3=0, pH=8) !
1355! !
1356! pcO2= ppmv (DoNewton=0) !
1357! pCO2= ppmv (DoNewton=1) !
1358! !
1359! This subroutine was adapted by Katja Fennel (Nov 2005) from !
1360! Zeebe and Wolf-Gladrow (2001). !
1361! !
1362! Reference: !
1363! !
1364! Zeebe, R.E. and D. Wolf-Gladrow, 2005: CO2 in Seawater: !
1365! Equilibrium, kinetics, isotopes, Elsevier Oceanographic !
1366! Series, 65, pp 346. !
1367! !
1368!***********************************************************************
1369!
1370 USE mod_kinds
1371!
1372 implicit none
1373!
1374! Imported variable declarations.
1375!
1376 integer, intent(in) :: LBi, UBi, LBj, UBj, IminS, ImaxS
1377 integer, intent(in) :: Istr, Iend, j, DoNewton
1378!
1379# ifdef ASSUMED_SHAPE
1380# ifdef MASKING
1381 real(r8), intent(in) :: rmask(LBi:,LBj:)
1382# endif
1383 real(r8), intent(in) :: T(IminS:)
1384 real(r8), intent(in) :: S(IminS:)
1385 real(r8), intent(in) :: TIC(IminS:)
1386 real(r8), intent(in) :: TAlk(IminS:)
1387 real(r8), intent(inout) :: pH(LBi:,LBj:)
1388# else
1389# ifdef MASKING
1390 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1391# endif
1392 real(r8), intent(in) :: T(IminS:ImaxS)
1393 real(r8), intent(in) :: S(IminS:ImaxS)
1394 real(r8), intent(in) :: TIC(IminS:ImaxS)
1395 real(r8), intent(in) :: TAlk(IminS:ImaxS)
1396 real(r8), intent(inout) :: pH(LBi:UBi,LBj:UBj)
1397# endif
1398
1399 real(r8), intent(out) :: pCO2(IminS:ImaxS)
1400!
1401! Local variable declarations.
1402!
1403 integer, parameter :: InewtonMax = 10
1404 integer, parameter :: IbrackMax = 30
1405
1406 integer :: Hstep, Ibrack, Inewton, i
1407
1408 real(r8) :: Tk, centiTk, invTk, logTk
1409 real(r8) :: scl, sqrtS
1410 real(r8) :: borate, alk, dic
1411 real(r8) :: ff, K1, K2, K12, Kb, Kw
1412 real(r8) :: p5, p4, p3, p2, p1, p0
1413 real(r8) :: df, fn, fni(3), ftest
1414 real(r8) :: deltaX, invX, invX2, X, X2, X3
1415 real(r8) :: pH_guess, pH_hi, pH_lo
1416 real(r8) :: X_guess, X_hi, X_lo, X_mid
1417 real(r8) :: CO2star, Htotal, Htotal2
1418!
1419!=======================================================================
1420! Determine coefficients for surface carbon chemisty. If land/sea
1421! masking, compute only on water points.
1422!=======================================================================
1423!
1424 I_LOOP: DO i=Istr,Iend
1425# ifdef MASKING
1426 IF (rmask(i,j).gt.0.0_r8) THEN
1427# endif
1428 Tk=T(i)+273.15_r8
1429 centiTk=0.01_r8*Tk
1430 invTk=1.0_r8/Tk
1431 logTk=LOG(Tk)
1432 sqrtS=SQRT(S(i))
1433 scl=S(i)/1.80655_r8
1434
1435 alk= TAlk(i)*0.000001_r8
1436 dic = TIC(i)*0.000001_r8
1437!
1438!-----------------------------------------------------------------------
1439! Correction term for non-ideality, ff=k0*(1-pH2O). Equation 13 with
1440! table 6 values from Weiss and Price (1980, Mar. Chem., 8, 347-359).
1441!-----------------------------------------------------------------------
1442!
1443 ff=EXP(-162.8301_r8+ &
1444 & 218.2968_r8/centiTk+ &
1445 & LOG(centiTk)*90.9241_r8- &
1446 & centiTk*centiTk*1.47696_r8+ &
1447 & S(i)*(0.025695_r8- &
1448 & centiTk*(0.025225_r8- &
1449 & centiTk*0.0049867_r8)))
1450!
1451!-----------------------------------------------------------------------
1452! Compute first (K1) and second (K2) dissociation constant of carboinic
1453! acid:
1454!
1455! K1 = [H][HCO3]/[H2CO3]
1456! K2 = [H][CO3]/[HCO3]
1457!
1458! From Millero (1995; page 664) using Mehrbach et al. (1973) data on
1459! seawater scale.
1460!-----------------------------------------------------------------------
1461!
1462 K1=10.0_r8**(62.008_r8- &
1463 & invTk*3670.7_r8- &
1464 & logTk*9.7944_r8+ &
1465 & S(i)*(0.0118_r8- &
1466 & S(i)*0.000116_r8))
1467 K2=10.0_r8**(-4.777_r8- &
1468 & invTk*1394.7_r8+ &
1469 & S(i)*(0.0184_r8- &
1470 & S(i)*0.000118_r8))
1471!
1472!-----------------------------------------------------------------------
1473! Compute dissociation constant of boric acid, Kb=[H][BO2]/[HBO2].
1474! From Millero (1995; page 669) using data from Dickson (1990).
1475!-----------------------------------------------------------------------
1476!
1477 Kb=EXP(-invTk*(8966.90_r8+ &
1478 & sqrtS*(2890.53_r8+ &
1479 & sqrtS*(77.942_r8- &
1480 & sqrtS*(1.728_r8- &
1481 & sqrtS*0.0996_r8))))- &
1482 & logTk*(24.4344_r8+ &
1483 & sqrtS*(25.085_r8+ &
1484 & sqrtS*0.2474_r8))+ &
1485 & Tk*(sqrtS*0.053105_r8)+ &
1486 & 148.0248_r8+ &
1487 & sqrtS*(137.1942_r8+ &
1488 & sqrtS*1.62142_r8))
1489!
1490!-----------------------------------------------------------------------
1491! Compute ion product of whater, Kw = [H][OH].
1492! From Millero (1995; page 670) using composite data.
1493!-----------------------------------------------------------------------
1494!
1495 Kw=EXP(148.9652_r8- &
1496 & invTk*13847.26_r8- &
1497 & logTk*23.6521_r8- &
1498 & sqrtS*(5.977_r8- &
1499 & invTk*118.67_r8- &
1500 & logTk*1.0495_r8)- &
1501 & S(i)*0.01615_r8)
1502!
1503!-----------------------------------------------------------------------
1504! Calculate concentrations for borate (Uppstrom, 1974).
1505!-----------------------------------------------------------------------
1506!
1507 borate=0.000232_r8*scl/10.811_r8
1508!
1509!=======================================================================
1510! Iteratively solver for computing hydrogen ions [H+] using either:
1511!
1512! (1) Newton-Raphson method with fixed number of iterations,
1513! use previous [H+] as first guess, or
1514! (2) bracket and bisection
1515!=======================================================================
1516!
1517! Solve for h in fifth-order polynomial. First calculate
1518! polynomial coefficients.
1519!
1520 K12 = K1*K2
1521
1522 p5 = -1.0_r8;
1523 p4 = -alk-Kb-K1;
1524 p3 = dic*K1-alk*(Kb+K1)+Kb*borate+Kw-Kb*K1-K12
1525 p2 = dic*(Kb*K1+2*K12)-alk*(Kb*K1+K12)+Kb*borate*K1 &
1526 & +(Kw*Kb+Kw*K1-Kb*K12)
1527 p1 = 2.0_r8*dic*Kb*K12-alk*Kb*K12+Kb*borate*K12 &
1528 & +Kw*Kb*K1+Kw*K12
1529 p0 = Kw*Kb*K12;
1530!
1531! Set first guess and brackets for [H+] solvers.
1532!
1533 pH_guess=pH(i,j) ! Newton-Raphson
1534 pH_hi=10.0_r8 ! high bracket/bisection
1535 pH_lo=5.0_r8 ! low bracket/bisection
1536!
1537! Convert to [H+].
1538!
1539 X_guess=10.0_r8**(-pH_guess)
1540 X_lo=10.0_r8**(-pH_hi)
1541 X_hi=10.0_r8**(-pH_lo)
1542 X_mid=0.5_r8*(X_lo+X_hi)
1543!
1544!-----------------------------------------------------------------------
1545! Newton-Raphson method.
1546!-----------------------------------------------------------------------
1547!
1548 IF (DoNewton.eq.1) THEN
1549 X=X_guess
1550!
1551 DO Inewton=1,InewtonMax
1552!
1553! Evaluate f([H+]) = p5*x^5+...+p1*x+p0
1554!
1555 fn=((((p5*X+p4)*X+p3)*X+p2)*X+p1)*X+p0
1556!
1557! Evaluate derivative, df([H+])/dx:
1558!
1559! df= d(fn)/d(X)
1560!
1561 df=(((5*p5*X+4*p4)*X+3*p3)*X+2*p2)*X+p1
1562!
1563! Evaluate increment in [H+].
1564!
1565 deltaX=-fn/df
1566!
1567! Update estimate of [H+].
1568!
1569 X=X+deltaX
1570 END DO
1571!
1572!-----------------------------------------------------------------------
1573! Bracket and bisection method.
1574!-----------------------------------------------------------------------
1575!
1576 ELSE
1577!
1578! If first step, use Bracket and Bisection method with fixed, large
1579! number of iterations
1580!
1581 BRACK_IT: DO Ibrack=1,IbrackMax
1582 DO Hstep=1,3
1583 IF (Hstep.eq.1) X=X_hi
1584 IF (Hstep.eq.2) X=X_lo
1585 IF (Hstep.eq.3) X=X_mid
1586!
1587! Evaluate f([H+]) for bracketing and mid-value cases.
1588!
1589 fni(Hstep)=((((p5*X+p4)*X+p3)*X+p2)*X+p1)*X+p0
1590 END DO
1591!
1592! Now, bracket solution within two of three.
1593!
1594 IF (fni(3).eq.0) THEN
1595 EXIT BRACK_IT
1596 ELSE
1597 ftest=fni(1)/fni(3)
1598 IF (ftest.gt.0) THEN
1599 X_hi=X_mid
1600 ELSE
1601 X_lo=X_mid
1602 END IF
1603 X_mid=0.5_r8*(X_lo+X_hi)
1604 END IF
1605 END DO BRACK_IT
1606!
1607! Last iteration gives value.
1608!
1609 X=X_mid
1610 END IF
1611!
1612!-----------------------------------------------------------------------
1613! Determine pCO2.
1614!-----------------------------------------------------------------------
1615!
1616! Total Hydrogen ion concentration, Htotal = [H+].
1617!
1618 Htotal=X
1619 Htotal2=Htotal*Htotal
1620!
1621! Calculate [CO2*] (mole/m3) as defined in DOE Methods Handbook 1994
1622! Version 2, ORNL/CDIAC-74, Dickson and Goyet, Eds. (Chapter 2,
1623! page 10, Eq A.49).
1624!
1625 CO2star=dic*Htotal2/(Htotal2+K1*Htotal+K1*K2)
1626!
1627! Save pH is used again outside this routine.
1628!
1629 pH(i,j)=-LOG10(Htotal)
1630!
1631! Add two output arguments for storing pCO2surf.
1632!
1633 pCO2(i)=CO2star*1000000.0_r8/ff
1634
1635# ifdef MASKING
1636 ELSE
1637 pH(i,j)=0.0_r8
1638 pCO2(i)=0.0_r8
1639 END IF
1640# endif
1641
1642 END DO I_LOOP
1643
1644 RETURN
1645 END SUBROUTINE pCO2_water_RZ
1646# else
1647 SUBROUTINE pCO2_water (Istr, Iend, &
1648 & LBi, UBi, LBj, UBj, &
1649 & IminS, ImaxS, j, DoNewton, &
1650# ifdef MASKING
1651 & rmask, &
1652# endif
1653 & T, S, TIC, TAlk, PO4, SiO3, pH, pCO2)
1654!
1655!***********************************************************************
1656! !
1657! This routine computes equilibrium partial pressure of CO2 (pCO2) !
1658! in the surface seawater. !
1659! !
1660! On Input: !
1661! !
1662! Istr Starting tile index in the I-direction. !
1663! Iend Ending tile index in the I-direction. !
1664! LBi I-dimension lower bound. !
1665! UBi I-dimension upper bound. !
1666! LBj J-dimension lower bound. !
1667! UBj J-dimension upper bound. !
1668! IminS I-dimension lower bound for private arrays. !
1669! ImaxS I-dimension upper bound for private arrays. !
1670! j j-pipelined index. !
1671! DoNewton Iteration solver: !
1672! [0] Bracket and bisection. !
1673! [1] Newton-Raphson method. !
1674! rmask Land/Sea masking. !
1675! T Surface temperature (Celsius). !
1676! S Surface salinity (PSS). !
1677! TIC Total inorganic carbon (millimol/m3). !
1678! TAlk Total alkalinity (milli-equivalents/m3). !
1679! PO4 Inorganic phosphate (millimol/m3). !
1680! SiO3 Inorganic silicate (millimol/m3). !
1681! pH Best pH guess. !
1682! !
1683! On Output: !
1684! !
1685! pCO2 partial pressure of CO2 (ppmv). !
1686! !
1687! Check Value: (T=24, S=36.6, TIC=2040, TAlk=2390, PO4=0, !
1688! SiO3=0, pH=8) !
1689! !
1690! pcO2=0.35074945E+03 ppmv (DoNewton=0) !
1691! pCO2=0.35073560E+03 ppmv (DoNewton=1) !
1692! !
1693! This subroutine was adapted by Mick Follows (Oct 1999) from OCMIP2 !
1694! code CO2CALC. Modified for ROMS by Hernan Arango (Nov 2003). !
1695! !
1696!***********************************************************************
1697!
1698 USE mod_kinds
1699!
1700 implicit none
1701!
1702! Imported variable declarations.
1703!
1704 integer, intent(in) :: LBi, UBi, LBj, UBj, IminS, ImaxS
1705 integer, intent(in) :: Istr, Iend, j, DoNewton
1706!
1707# ifdef ASSUMED_SHAPE
1708# ifdef MASKING
1709 real(r8), intent(in) :: rmask(LBi:,LBj:)
1710# endif
1711 real(r8), intent(in) :: T(IminS:)
1712 real(r8), intent(in) :: S(IminS:)
1713 real(r8), intent(in) :: TIC(IminS:)
1714 real(r8), intent(in) :: TAlk(IminS:)
1715 real(r8), intent(inout) :: pH(LBi:,LBj:)
1716# else
1717# ifdef MASKING
1718 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1719# endif
1720 real(r8), intent(in) :: T(IminS:ImaxS)
1721 real(r8), intent(in) :: S(IminS:ImaxS)
1722 real(r8), intent(in) :: TIC(IminS:ImaxS)
1723 real(r8), intent(in) :: TAlk(IminS:ImaxS)
1724 real(r8), intent(inout) :: pH(LBi:UBi,LBj:UBj)
1725# endif
1726 real(r8), intent(in) :: PO4
1727 real(r8), intent(in) :: SiO3
1728
1729 real(r8), intent(out) :: pCO2(IminS:ImaxS)
1730!
1731! Local variable declarations.
1732!
1733 integer, parameter :: InewtonMax = 10
1734 integer, parameter :: IbrackMax = 30
1735
1736 integer :: Hstep, Ibrack, Inewton, i
1737
1738 real(r8) :: Tk, centiTk, invTk, logTk
1739 real(r8) :: SO4, scl, sqrtS, sqrtSO4
1740 real(r8) :: alk, dic, phos, sili
1741 real(r8) :: borate, sulfate, fluoride
1742 real(r8) :: ff, K1, K2, K1p, K2p, K3p, Kb, Kf, Ks, Ksi, Kw
1743 real(r8) :: K12, K12p, K123p, invKb, invKs, invKsi
1744 real(r8) :: A, A2, B, B2, C, dA, dB
1745 real(r8) :: df, fn, fni(3), ftest
1746 real(r8) :: deltaX, invX, invX2, X, X2, X3
1747 real(r8) :: pH_guess, pH_hi, pH_lo
1748 real(r8) :: X_guess, X_hi, X_lo, X_mid
1749 real(r8) :: CO2star, Htotal, Htotal2
1750!
1751!=======================================================================
1752! Determine coefficients for surface carbon chemisty. If land/sea
1753! masking, compute only on water points.
1754!=======================================================================
1755!
1756 I_LOOP: DO i=Istr,Iend
1757# ifdef MASKING
1758 IF (rmask(i,j).gt.0.0_r8) THEN
1759# endif
1760 Tk=T(i)+273.15_r8
1761 centiTk=0.01_r8*Tk
1762 invTk=1.0_r8/Tk
1763 logTk=LOG(Tk)
1764 sqrtS=SQRT(S(i))
1765 SO4=19.924_r8*S(i)/(1000.0_r8-1.005_r8*S(i))
1766 sqrtSO4=SQRT(SO4)
1767 scl=S(i)/1.80655_r8
1768
1769 alk=TAlk(i)*0.000001_r8
1770 dic=TIC(i)*0.000001_r8
1771 phos=PO4*0.000001_r8
1772 sili=SiO3*0.000001_r8
1773!
1774!-----------------------------------------------------------------------
1775! Correction term for non-ideality, ff=k0*(1-pH2O). Equation 13 with
1776! table 6 values from Weiss and Price (1980, Mar. Chem., 8, 347-359).
1777!-----------------------------------------------------------------------
1778!
1779 ff=EXP(-162.8301_r8+ &
1780 & 218.2968_r8/centiTk+ &
1781 & LOG(centiTk)*90.9241_r8- &
1782 & centiTk*centiTk*1.47696_r8+ &
1783 & S(i)*(0.025695_r8- &
1784 & centiTk*(0.025225_r8- &
1785 & centiTk*0.0049867_r8)))
1786!
1787!-----------------------------------------------------------------------
1788! Compute first (K1) and second (K2) dissociation constant of carboinic
1789! acid:
1790!
1791! K1 = [H][HCO3]/[H2CO3]
1792! K2 = [H][CO3]/[HCO3]
1793!
1794! From Millero (1995; page 664) using Mehrbach et al. (1973) data on
1795! seawater scale.
1796!-----------------------------------------------------------------------
1797!
1798 K1=10.0_r8**(62.008_r8- &
1799 & invTk*3670.7_r8- &
1800 & logTk*9.7944_r8+ &
1801 & S(i)*(0.0118_r8- &
1802 & S(i)*0.000116_r8))
1803 K2=10.0_r8**(-4.777_r8- &
1804 & invTk*1394.7_r8+ &
1805 & S(i)*(0.0184_r8- &
1806 & S(i)*0.000118_r8))
1807!
1808!-----------------------------------------------------------------------
1809! Compute dissociation constant of boric acid, Kb=[H][BO2]/[HBO2].
1810! From Millero (1995; page 669) using data from Dickson (1990).
1811!-----------------------------------------------------------------------
1812!
1813 Kb=EXP(-invTk*(8966.90_r8+ &
1814 & sqrtS*(2890.53_r8+ &
1815 & sqrtS*(77.942_r8- &
1816 & sqrtS*(1.728_r8- &
1817 & sqrtS*0.0996_r8))))- &
1818 & logTk*(24.4344_r8+ &
1819 & sqrtS*(25.085_r8+ &
1820 & sqrtS*0.2474_r8))+ &
1821 & Tk*(sqrtS*0.053105_r8)+ &
1822 & 148.0248_r8+ &
1823 & sqrtS*(137.1942_r8+ &
1824 & sqrtS*1.62142_r8))
1825!
1826!-----------------------------------------------------------------------
1827! Compute first (K1p), second (K2p), and third (K3p) dissociation
1828! constant of phosphoric acid:
1829!
1830! K1p = [H][H2PO4]/[H3PO4]
1831! K2p = [H][HPO4]/[H2PO4]
1832! K3p = [H][PO4]/[HPO4]
1833!
1834! From DOE (1994) equations 7.2.20, 7.2.23, and 7.2.26, respectively.
1835! With footnote using data from Millero (1974).
1836!-----------------------------------------------------------------------
1837!
1838 K1p=EXP(115.525_r8- &
1839 & invTk*4576.752_r8- &
1840 & logTk*18.453_r8+ &
1841 & sqrtS*(0.69171_r8-invTk*106.736_r8)- &
1842 & S(i)*(0.01844_r8+invTk*0.65643_r8))
1843 K2p=EXP(172.0883_r8- &
1844 & invTk*8814.715_r8- &
1845 & logTk*27.927_r8+ &
1846 & sqrtS*(1.3566_r8-invTk*160.340_r8)- &
1847 & S(i)*(0.05778_r8-invTk*0.37335_r8))
1848 K3p=EXP(-18.141_r8- &
1849 & invTk*3070.75_r8+ &
1850 & sqrtS*(2.81197_r8+invTk*17.27039_r8)- &
1851 & S(i)*(0.09984_r8+invTk*44.99486_r8))
1852!
1853!-----------------------------------------------------------------------
1854! Compute dissociation constant of silica, Ksi=[H][SiO(OH)3]/[Si(OH)4].
1855! From Millero (1995; page 671) using data from Yao and Millero (1995).
1856!-----------------------------------------------------------------------
1857!
1858 Ksi=EXP(117.385_r8- &
1859 & invTk*8904.2_r8- &
1860 & logTk*19.334_r8+ &
1861 & sqrtSO4*(3.5913_r8-invTk*458.79_r8)- &
1862 & SO4*(1.5998_r8-invTk*188.74_r8- &
1863 & SO4*(0.07871_r8-invTk*12.1652_r8))+ &
1864 & LOG(1.0_r8-0.001005_r8*S(i)))
1865!
1866!-----------------------------------------------------------------------
1867! Compute ion product of whater, Kw = [H][OH].
1868! From Millero (1995; page 670) using composite data.
1869!-----------------------------------------------------------------------
1870!
1871 Kw=EXP(148.9652_r8- &
1872 & invTk*13847.26_r8- &
1873 & logTk*23.6521_r8- &
1874 & sqrtS*(5.977_r8- &
1875 & invTk*118.67_r8- &
1876 & logTk*1.0495_r8)- &
1877 & S(i)*0.01615_r8)
1878!
1879!------------------------------------------------------------------------
1880! Compute salinity constant of hydrogen sulfate, Ks = [H][SO4]/[HSO4].
1881! From Dickson (1990, J. chem. Thermodynamics 22, 113)
1882!------------------------------------------------------------------------
1883!
1884 Ks=EXP(141.328_r8- &
1885 & invTk*4276.1_r8- &
1886 & logTk*23.093_r8+ &
1887 & sqrtSO4*(324.57_r8-invTk*13856.0_r8-logTk*47.986_r8- &
1888 & SO4*invTk*2698.0_r8)- &
1889 & SO4*(771.54_r8-invTk*35474.0_r8-logTk*114.723_r8- &
1890 & SO4*invTk*1776.0_r8)+ &
1891 & LOG(1.0_r8-0.001005_r8*S(i)))
1892!
1893!-----------------------------------------------------------------------
1894! Compute stability constant of hydrogen fluorid, Kf = [H][F]/[HF].
1895! From Dickson and Riley (1979) -- change pH scale to total.
1896!-----------------------------------------------------------------------
1897!
1898 Kf=EXP(-12.641_r8+ &
1899 & invTk*1590.2_r8+ &
1900 & sqrtSO4*1.525_r8+ &
1901 & LOG(1.0_r8-0.001005_r8*S(i))+ &
1902 & LOG(1.0_r8+0.1400_r8*scl/(96.062_r8*Ks)))
1903!
1904!-----------------------------------------------------------------------
1905! Calculate concentrations for borate (Uppstrom, 1974), sulfate (Morris
1906! and Riley, 1966), and fluoride (Riley, 1965).
1907!-----------------------------------------------------------------------
1908!
1909 borate=0.000232_r8*scl/10.811_r8
1910 sulfate=0.14_r8*scl/96.062_r8
1911 fluoride=0.000067_r8*scl/18.9984_r8
1912!
1913!=======================================================================
1914! Iteratively solver for computing hydrogen ions [H+] using either:
1915!
1916! (1) Newton-Raphson method with fixed number of iterations,
1917! use previous [H+] as first guess, or
1918! (2) bracket and bisection
1919!=======================================================================
1920!
1921! Set first guess and brackets for [H+] solvers.
1922!
1923 pH_guess=pH(i,j) ! Newton-Raphson
1924 pH_hi=10.0_r8 ! high bracket/bisection
1925 pH_lo=5.0_r8 ! low bracket/bisection
1926!
1927! Convert to [H+].
1928!
1929 X_guess=10.0_r8**(-pH_guess)
1930 X_lo=10.0_r8**(-pH_hi)
1931 X_hi=10.0_r8**(-pH_lo)
1932 X_mid=0.5_r8*(X_lo+X_hi)
1933!
1934!-----------------------------------------------------------------------
1935! Newton-Raphson method.
1936!-----------------------------------------------------------------------
1937!
1938 IF (DoNewton.eq.1) THEN
1939 X=X_guess
1940 K12=K1*K2
1941 K12p=K1p*K2p
1942 K123p=K12p*K3p
1943 invKb=1.0_r8/Kb
1944 invKs=1.0_r8/Ks
1945 invKsi=1.0_r8/Ksi
1946!
1947 DO Inewton=1,InewtonMax
1948!
1949! Set some common combinations of parameters used in the iterative [H+]
1950! solver.
1951!
1952 X2=X*X
1953 X3=X2*X
1954 invX=1.0_r8/X
1955 invX2=1.0_r8/X2
1956
1957 A=X*(K12p+X*(K1p+X))
1958 B=X*(K1+X)+K12
1959 C=1.0_r8/(1.0_r8+sulfate*invKs)
1960
1961 A2=A*A
1962 B2=B*B
1963 dA=X*(2.0_r8*K1p+3.0_r8*X)+K12p
1964 dB=2.0_r8*X+K1
1965!
1966! Evaluate f([H+]):
1967!
1968! fn=HCO3+CO3+borate+OH+HPO4+2*PO4+H3PO4+silicate+Hfree+HSO4+HF-TALK
1969!
1970 fn=dic*K1*(X+2.0_r8*K2)/B+ &
1971 & borate/(1.0_r8+X*invKb)+ &
1972 & Kw*invX+ &
1973 & phos*(K12p*X+2.0_r8*K123p-X3)/A+ &
1974 & sili/(1.0_r8+X*invKsi)- &
1975 & X*C- &
1976 & sulfate/(1.0_r8+Ks*invX*C)- &
1977 & fluoride/(1.0_r8+Kf*invX)- &
1978 & alk
1979!
1980! Evaluate derivative, f(prime)([H+]):
1981!
1982! df= d(fn)/d(X)
1983!
1984 df=dic*K1*(B-dB*(X+2.0_r8*K2))/B2- &
1985 & borate/(invKb*(1.0+X*invKb)**2)- &
1986 & Kw*invX2+ &
1987 & phos*(A*(K12p-3.0_r8*X2)-dA*(K12p*X+2.0_r8*K123p-X3))/A2-&
1988 & sili/(invKsi*(1.0_r8+X*invKsi)**2)+ &
1989 & C+ &
1990 & sulfate*Ks*C*invX2/((1.0_r8+Ks*invX*C)**2)+ &
1991 & fluoride*Kf*invX2/((1.0_r8+Kf*invX)**2)
1992!
1993! Evaluate increment in [H+].
1994!
1995 deltaX=-fn/df
1996!
1997! Update estimate of [H+].
1998!
1999 X=X+deltaX
2000 END DO
2001!
2002!-----------------------------------------------------------------------
2003! Bracket and bisection method.
2004!-----------------------------------------------------------------------
2005!
2006 ELSE
2007!
2008! If first step, use Bracket and Bisection method with fixed, large
2009! number of iterations
2010!
2011 K12=K1*K2
2012 K12p=K1p*K2p
2013 K123p=K12p*K3p
2014 invKb=1.0_r8/Kb
2015 invKs=1.0_r8/Ks
2016 invKsi=1.0_r8/Ksi
2017!
2018 BRACK_IT: DO Ibrack=1,IbrackMax
2019 DO Hstep=1,3
2020 IF (Hstep.eq.1) X=X_hi
2021 IF (Hstep.eq.2) X=X_lo
2022 IF (Hstep.eq.3) X=X_mid
2023!
2024! Set some common combinations of parameters used in the iterative [H+]
2025! solver.
2026!
2027 X2=X*X
2028 X3=X2*X
2029 invX=1.0_r8/X
2030
2031 A=X*(K12p+X*(K1p+X))+K123p
2032 B=X*(K1+X)+K12
2033 C=1.0_r8/(1.0_r8+sulfate*invKs)
2034
2035 A2=A*A
2036 B2=B*B
2037 dA=X*(K1p*2.0_r8+3.0_r8*X2)+K12p
2038 dB=2.0_r8*X+K1
2039!
2040! Evaluate f([H+]) for bracketing and mid-value cases.
2041!
2042 fni(Hstep)=dic*(K1*X+2.0_r8*K12)/B+ &
2043 & borate/(1.0_r8+X*invKb)+ &
2044 & Kw*invX+ &
2045 & phos*(K12p*X+2.0_r8*K123p-X3)/A+ &
2046 & sili/(1.0_r8+X*invKsi)- &
2047 & X*C- &
2048 & sulfate/(1.0_r8+Ks*invX*C)- &
2049 & fluoride/(1.0_r8+Kf*invX)- &
2050 & alk
2051 END DO
2052!
2053! Now, bracket solution within two of three.
2054!
2055 IF (fni(3).eq.0.0_r8) THEN
2056 EXIT BRACK_IT
2057 ELSE
2058 ftest=fni(1)/fni(3)
2059 IF (ftest.gt.0.0) THEN
2060 X_hi=X_mid
2061 ELSE
2062 X_lo=X_mid
2063 END IF
2064 X_mid=0.5_r8*(X_lo+X_hi)
2065 END IF
2066 END DO BRACK_IT
2067!
2068! Last iteration gives value.
2069!
2070 X=X_mid
2071 END IF
2072!
2073!-----------------------------------------------------------------------
2074! Determine pCO2.
2075!-----------------------------------------------------------------------
2076!
2077! Total Hydrogen ion concentration, Htotal = [H+].
2078!
2079 Htotal=X
2080 Htotal2=Htotal*Htotal
2081!
2082! Calculate [CO2*] (mole/m3) as defined in DOE Methods Handbook 1994
2083! Version 2, ORNL/CDIAC-74, Dickson and Goyet, Eds. (Chapter 2,
2084! page 10, Eq A.49).
2085!
2086 CO2star=dic*Htotal2/(Htotal2+K1*Htotal+K1*K2)
2087!
2088! Save pH is used again outside this routine.
2089!
2090 pH(i,j)=-LOG10(Htotal)
2091!
2092! Add two output arguments for storing pCO2surf.
2093!
2094 pCO2(i)=CO2star*1000000.0_r8/ff
2095
2096# ifdef MASKING
2097 ELSE
2098 pH(i,j)=0.0_r8
2099 pCO2(i)=0.0_r8
2100 END IF
2101# endif
2102
2103 END DO I_LOOP
2104
2105 RETURN
2106 END SUBROUTINE pCO2_water
2107# endif
2108#endif