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
|
---|