1 | MODULE asmpco2bal |
---|
2 | !Contains subroutine to convert a global pco2 increment to independent DIC and alkalinity increments |
---|
3 | |
---|
4 | USE par_kind, ONLY: wp |
---|
5 | USE par_oce, ONLY: jpi, jpj, jpk |
---|
6 | USE dom_oce, ONLY: tmask |
---|
7 | #if defined key_hadocc |
---|
8 | USE had_bgc_const, ONLY: RHO_WATER_SI |
---|
9 | #else |
---|
10 | USE phycst, ONLY: RHO_WATER_SI => rau0 |
---|
11 | #endif |
---|
12 | |
---|
13 | IMPLICIT NONE |
---|
14 | PRIVATE |
---|
15 | |
---|
16 | PUBLIC asm_pco2_bal |
---|
17 | CONTAINS |
---|
18 | |
---|
19 | SUBROUTINE asm_pco2_bal(pco2inc, tco2_tracer_bkg, alk_tracer_bkg, & |
---|
20 | & tem_tracer_bkg, sal_tracer_bkg, & |
---|
21 | & tco2_tracer_inc, alk_tracer_inc) |
---|
22 | !Converts a global pco2 increment to independent DIC and alkalinity increments |
---|
23 | !Calculates the solution to the equation [DDIC, DAlk]=(DpCO2 / |Grad (pCO2)|)* Grad(pCO2); where Grad(pco2)=[dpCO2/dDIC, dpCO2/dalk] |
---|
24 | !Gradient code pco2_grad is a modification of assimcd1.mf77 by Ian Totterdell. |
---|
25 | |
---|
26 | !James While 06/09 |
---|
27 | |
---|
28 | IMPLICIT NONE |
---|
29 | |
---|
30 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pco2inc !pCO2 increment |
---|
31 | |
---|
32 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg, & !In: tCO2 (DIC) background |
---|
33 | alk_tracer_bkg !In alkalinity background |
---|
34 | |
---|
35 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg, & !In: temperature background |
---|
36 | sal_tracer_bkg !In salinity background |
---|
37 | |
---|
38 | REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: tco2_tracer_inc, & !Out: tCO2 (DIC) increment |
---|
39 | alk_tracer_inc !Out alkalinity |
---|
40 | |
---|
41 | !local variables |
---|
42 | REAL(wp), DIMENSION(jpi,jpj) :: Dpco2_tco2, Dpco2_alk |
---|
43 | REAL(wp) :: wgt_fct |
---|
44 | INTEGER :: ji,jj |
---|
45 | |
---|
46 | |
---|
47 | !calculate the gradients with respect to DIC and Alkalinity |
---|
48 | CALL pco2_grad(tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, Dpco2_tco2, Dpco2_alk) |
---|
49 | |
---|
50 | |
---|
51 | !loop over the ocean grid calculating the increments to tCO2 (DIC) and alkalinity |
---|
52 | DO jj=1, jpj |
---|
53 | DO ji=1, jpi |
---|
54 | |
---|
55 | wgt_fct=pco2inc(ji,jj)/( (Dpco2_tco2(ji,jj)*Dpco2_tco2(ji,jj)) + (Dpco2_alk(ji,jj)*Dpco2_alk(ji,jj)) ) |
---|
56 | tco2_tracer_inc(ji,jj)=wgt_fct * Dpco2_tco2(ji,jj) * tmask(ji,jj,1) |
---|
57 | alk_tracer_inc(ji,jj)=wgt_fct * Dpco2_alk(ji,jj) * tmask(ji,jj,1) |
---|
58 | |
---|
59 | ENDDO |
---|
60 | ENDDO |
---|
61 | |
---|
62 | |
---|
63 | |
---|
64 | END SUBROUTINE asm_pco2_bal |
---|
65 | |
---|
66 | |
---|
67 | SUBROUTINE pco2_grad( tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, DPDC, DPDA ) |
---|
68 | ! |
---|
69 | !-------M---------M---------M---------M---------M---------M---------M-- |
---|
70 | ! |
---|
71 | ! Ian J Totterdell |
---|
72 | ! 2007 August 08 |
---|
73 | ! |
---|
74 | ! Code for calculating derivatives of pCO2 wrt DIC and Alk |
---|
75 | ! |
---|
76 | ! J While |
---|
77 | ! June 2009 |
---|
78 | ! Modiefied for use in NEMO |
---|
79 | |
---|
80 | !-------M---------M---------M---------M---------M---------M---------M-- |
---|
81 | IMPLICIT NONE |
---|
82 | ! |
---|
83 | ! |
---|
84 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg, & !In: tCO2 (DIC) background |
---|
85 | alk_tracer_bkg !In alkalinity background |
---|
86 | |
---|
87 | REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg, & !In: temperature background |
---|
88 | sal_tracer_bkg !In salinity background |
---|
89 | ! |
---|
90 | REAL(wp), INTENT(OUT), DIMENSION(jpi,jpj) :: & |
---|
91 | DPDC, & |
---|
92 | DPDA |
---|
93 | ! |
---|
94 | ! |
---|
95 | ! Local variables. |
---|
96 | ! |
---|
97 | INTEGER :: & |
---|
98 | mxiter, & |
---|
99 | ji,jj, & |
---|
100 | iter |
---|
101 | ! |
---|
102 | REAL(wp) :: & |
---|
103 | a1_co2, a2_co2, a3_co2, & |
---|
104 | b1_co2, b2_co2, b3_co2, & |
---|
105 | a1_k1, a2_k1, a3_k1, & |
---|
106 | b1_k1, b2_k1, b3_k1, b4_k1, b5_k1, & |
---|
107 | a1_k2, a2_k2, a3_k2, & |
---|
108 | b1_k2, b2_k2, b3_k2, b4_k2, b5_k2, & |
---|
109 | a1_kb, a2_kb, a3_kb, a4_kb, a5_kb, & |
---|
110 | b1_kb, b2_kb, b3_kb, & |
---|
111 | c1_kb, c2_kb, c3_kb, d_kb, & |
---|
112 | a1_kw, a2_kw, a3_kw, & |
---|
113 | b1_kw, b2_kw, b3_kw, c_kw |
---|
114 | ! |
---|
115 | REAL(wp) :: & |
---|
116 | tk, logtk, ttk, & |
---|
117 | sal, sqsal, & |
---|
118 | kone, ktwo, kboric, kwater, & |
---|
119 | alkprime, qa,qb,qc, & |
---|
120 | fxa, xtrma, xtrmb, & |
---|
121 | dxda, dxdc, & |
---|
122 | dgda, dgdc |
---|
123 | ! |
---|
124 | REAL(wp), DIMENSION(jpi,jpj) :: & |
---|
125 | co2_sol, & |
---|
126 | rek, & |
---|
127 | rg1, & |
---|
128 | rg2, & |
---|
129 | rg3, & |
---|
130 | boron, & |
---|
131 | tco2, & |
---|
132 | alk_t, & |
---|
133 | xs, & |
---|
134 | x1, & |
---|
135 | x2, & |
---|
136 | x2p, & |
---|
137 | x3, & |
---|
138 | xx, & |
---|
139 | xf |
---|
140 | ! |
---|
141 | ! Set parameter and constant values. |
---|
142 | ! |
---|
143 | PARAMETER ( mxiter = 10 ) |
---|
144 | ! |
---|
145 | ! CO2 solubility, constants for the calculation of. |
---|
146 | ! From Weiss (1974) |
---|
147 | ! |
---|
148 | DATA & |
---|
149 | a1_co2 / -60.2409 /, a2_co2 / 93.4517 /, a3_co2 / 23.3585 /, & |
---|
150 | b1_co2 / 0.023517 /, b2_co2 / -0.023656 /, b3_co2 / 0.0047036 / |
---|
151 | ! |
---|
152 | ! kone, constants for the calculation of. |
---|
153 | ! From Roy et al (1993) |
---|
154 | ! |
---|
155 | DATA a1_k1 /-2307.1266/, a2_k1 / 2.83655 /, a3_k1 /-1.5529413/ & |
---|
156 | , b1_k1 / -4.0484 /, b2_k1 /-0.20760841/, b3_k1 /0.08468345/ & |
---|
157 | , b4_k1 /-0.00654208/, b5_k1 / -0.001005 / |
---|
158 | ! |
---|
159 | ! ktwo, constants for the calculation of. |
---|
160 | ! From Roy et al (1993) |
---|
161 | ! |
---|
162 | DATA a1_k2 /-3351.6106/, a2_k2 / -9.226508 /, a3_k2 /-0.2005743/ & |
---|
163 | , b1_k2 / -23.9722 /, b2_k2 /-0.106901773/,b3_k2 /0.1130822/ & |
---|
164 | , b4_k2 /-0.00846934/, b5_k2 / -0.001005 / |
---|
165 | ! |
---|
166 | ! kboric, constants for the calculation of. |
---|
167 | ! Using data from Dickson (1990); |
---|
168 | ! |
---|
169 | DATA a1_kb / -8966.90 /, a2_kb / -2890.53 /, a3_kb / -77.942 / & |
---|
170 | , a4_kb / 1.728 /, a5_kb /-0.0996/ & |
---|
171 | , b1_kb / 148.0248 /, b2_kb / 137.1942 /, b3_kb / 1.62142 / & |
---|
172 | , c1_kb / -24.4344 /, c2_kb / -25.085 /, c3_kb / -0.2474 / & |
---|
173 | , d_kb / 0.053105 / |
---|
174 | ! |
---|
175 | ! kwater, constants for the calculation of. |
---|
176 | ! From Millero (1995, p670), using composite data. |
---|
177 | ! |
---|
178 | DATA a1_kw / -13847.26 /, a2_kw / 148.9652 /, a3_kw / -23.6521 / & |
---|
179 | , b1_kw / 118.67 /, b2_kw / -5.977 /, b3_kw / 1.0495 / & |
---|
180 | , c_kw / -0.01615 / |
---|
181 | ! |
---|
182 | ! Executable code begins here. |
---|
183 | ! |
---|
184 | ! Calculate the carbon system constants. |
---|
185 | ! |
---|
186 | |
---|
187 | DO ji=1,jpi |
---|
188 | DO jj=1,jpj |
---|
189 | ! |
---|
190 | tk = 273.15 + MAX( -2.0, MIN( 40.0, tem_tracer_bkg(ji,jj) ) ) |
---|
191 | logtk = LOG(tk) |
---|
192 | ttk = MAX( 2.71, 0.01 * ( 273.15 + tem_tracer_bkg(ji,jj) ) ) |
---|
193 | sal = MAX( 0.0, MIN( 45.0, sal_tracer_bkg(ji,jj) ) ) |
---|
194 | sqsal = SQRT(sal) |
---|
195 | ! |
---|
196 | ! Calculate the CO2 solubility from Weiss (1974) |
---|
197 | ! Units are mol/(kg.soln)/atm . |
---|
198 | ! Check value (Zeebe and Wolf-G; tc=25,sal=35): ln(co2_sol) = -3.5617 |
---|
199 | ! |
---|
200 | co2_sol(ji,jj) = EXP( & |
---|
201 | a1_co2 + a2_co2 / ttk + a3_co2 * LOG(ttk) + & |
---|
202 | sal * ( b1_co2 + ttk * ( b2_co2 + b3_co2 * ttk ) ) ) |
---|
203 | ! |
---|
204 | ! Calculate the 1st constant of carbonic acid from Roy et al. (1993) |
---|
205 | ! Units are mol/(kg.soln) on the total hydrogen ion scale. |
---|
206 | ! K1 = (HCO3-).(H+)/(CO2(g)) |
---|
207 | ! Check value (Dickson and Goyet; tc=25,sal=35): ln(K1) = -13.4847 |
---|
208 | ! |
---|
209 | kone & |
---|
210 | = 1.0e-7 *( 1.0 + b5_k1*sal ) * EXP( a1_k1/tk + a2_k1 + a3_k1*logtk & |
---|
211 | + ( b1_k1/tk + b2_k1 )*sqsal + b3_k1*sal + b4_k1*sal*sqsal ) |
---|
212 | ! |
---|
213 | ! Calculate the 2nd constant of carbonic acid from Roy et al. (1993) |
---|
214 | ! Units are mol/(kg.soln) on the total hydrogen ion scale. |
---|
215 | ! K2 = (CO3--).(H+)/(HCO3-) |
---|
216 | ! Check value (Dickson and Goyet; tc=25,sal=35): ln(K2) = -20.5504 |
---|
217 | ! |
---|
218 | ktwo & |
---|
219 | = 1.0e-7 *( 1.0 + b5_k2*sal ) * EXP( a1_k2/tk + a2_k2 + a3_k2*logtk & |
---|
220 | + ( b1_k2/tk + b2_k2 )*sqsal + b3_k2*sal + b4_k2*sal*sqsal ) |
---|
221 | ! |
---|
222 | ! Calculate the equilibrium constant of boric acid from Dickson (1990). |
---|
223 | ! Units are mol/(kg.soln) on the total hydrogen ion scale. |
---|
224 | ! Kb = (B(OH)4-).(H+)/(B(OH)3) (with (Bt) = (B(OH)3) + (B(OH)4-) ). |
---|
225 | ! Check value (Dickson and Goyet; tc=25,sal=35): ln(Kb) = -19.7964 |
---|
226 | ! |
---|
227 | kboric & |
---|
228 | = 1.0e-7 * EXP( ( a1_kb + a2_kb*sqsal + a3_kb*sal + a4_kb*sal*sqsal & |
---|
229 | + a5_kb*sal*sal )/tk + ( b1_kb + b2_kb*sqsal + b3_kb*sal ) & |
---|
230 | + ( c1_kb + c2_kb*sqsal + c3_kb*sal )*logtk + d_kb*sqsal*tk ) |
---|
231 | ! |
---|
232 | ! Calculate the equilibrium constant of water from Millero (1994). |
---|
233 | ! Units are ( mol/(kg.soln) )^2 on the total hydrogen ion scale. |
---|
234 | ! Kw = (OH-).(H+) |
---|
235 | ! Check value (Dickson and Goyet; tc=25,sal=35): ln(Kw) = -30.434 |
---|
236 | ! |
---|
237 | kwater & |
---|
238 | = 1.0e-14 *EXP( ( a1_kw/tk + a2_kw + a3_kw*logtk ) & |
---|
239 | + ( b1_kw/tk + b2_kw + b3_kw*logtk )*sqsal + c_kw*sal ) |
---|
240 | ! |
---|
241 | ! Calculate the combinations (as in PPCO2). |
---|
242 | ! |
---|
243 | rek(ji,jj) = SQRT( kone / ktwo ) |
---|
244 | rg3(ji,jj) = SQRT( kone * ktwo ) |
---|
245 | rg1(ji,jj) = rg3(ji,jj) / kboric |
---|
246 | rg2(ji,jj) = kwater / rg3(ji,jj) |
---|
247 | ! |
---|
248 | ! Set Boron concentration according to the surface salinity |
---|
249 | ! (Boron accounts for about 10% of the total Alkalinity). |
---|
250 | ! Use relationship used by Peng et al, Tellus 39b 1987, |
---|
251 | ! originally from Culkin 1965: |
---|
252 | ! |
---|
253 | ! boron(ji,jj) = 4.16e-4 * sal / 35.0 ! Correct value but not in UM |
---|
254 | boron(ji,jj) = 4.106e-4 * sal / 35.0 |
---|
255 | ! |
---|
256 | ! tco2 and alk_t are got from the modelled tracers at the surface |
---|
257 | ! (converted to units of mol/kg). |
---|
258 | ! |
---|
259 | tco2(ji,jj) = MAX( 0.0, tco2_tracer_bkg(ji,jj) * 0.001 / RHO_WATER_SI ) |
---|
260 | alk_t(ji,jj) = MAX( 0.0, alk_tracer_bkg(ji,jj) * 0.001 / RHO_WATER_SI ) |
---|
261 | ! |
---|
262 | ! Set initial values for iterative method. |
---|
263 | ! |
---|
264 | xs(ji,jj) = 1.0 |
---|
265 | x2p(ji,jj) = 0.0 |
---|
266 | x3(ji,jj) = 0.0 |
---|
267 | xx(ji,jj) = 0.0 |
---|
268 | xf(ji,jj) = xs(ji,jj) |
---|
269 | ! |
---|
270 | ENDDO |
---|
271 | ENDDO |
---|
272 | |
---|
273 | ! |
---|
274 | ! Use the secant method of similar triangles to calculate the H+ ion |
---|
275 | ! concentration, [H] (though strictly xf refers to rg3/[H] ). |
---|
276 | ! To aid vectorisation this has been written to iterate a set |
---|
277 | ! number of times (mxiter), large enough to bring every example to |
---|
278 | ! within the 1e-5 fractional change criterion. |
---|
279 | ! |
---|
280 | ! The method follows that of Bacastow (1990) in separating the |
---|
281 | ! alkalinity into two parts: the major part is the carbonate alk |
---|
282 | ! which uses the [H] that we are trying to find, while the minor |
---|
283 | ! part (xw) is sum of the borate and water dissociation terms which |
---|
284 | ! use the previous approximation to [H]. |
---|
285 | ! This leaves an easy-to-solve quadratic. |
---|
286 | ! |
---|
287 | DO iter=1,mxiter |
---|
288 | ! |
---|
289 | |
---|
290 | DO ji=1,jpi |
---|
291 | DO jj=1,jpj |
---|
292 | ! |
---|
293 | IF((abs(xf(ji,jj)-xx(ji,jj))/xf(ji,jj))<1.e-5) CYCLE !prevent iteration collapsing to 0.0/0.0 |
---|
294 | |
---|
295 | |
---|
296 | x1(ji,jj) = x2p(ji,jj) |
---|
297 | x2(ji,jj) = x3(ji,jj) |
---|
298 | x2p(ji,jj) = xx(ji,jj) |
---|
299 | x3(ji,jj) = xf(ji,jj) |
---|
300 | |
---|
301 | IF (iter<=2) THEN |
---|
302 | xx(ji,jj) = xf(ji,jj) |
---|
303 | ELSE |
---|
304 | xx(ji,jj) = x2p(ji,jj) + ( x3(ji,jj) - x2p(ji,jj) )*( x1(ji,jj) - x2p(ji,jj) ) & |
---|
305 | /( x1(ji,jj) - x2(ji,jj) - x2p(ji,jj) + x3(ji,jj) ) |
---|
306 | ENDIF |
---|
307 | |
---|
308 | |
---|
309 | ! |
---|
310 | ! All the above has been setting up the similar triangles for the |
---|
311 | ! secant method and using it to find an improved interpolated |
---|
312 | ! estimate; it has been completely independent of the equations |
---|
313 | ! we are trying to solve. |
---|
314 | ! |
---|
315 | ! The next bit, however, IS dependent on the equations. |
---|
316 | ! alkprime is the carbonate alkalinity, which is calculated |
---|
317 | ! from the total alkalinity by subtracting the other, smaller, terms, |
---|
318 | ! which use the interpolated previous approximation to xf (xx). |
---|
319 | ! |
---|
320 | |
---|
321 | alkprime = alk_t(ji,jj) - rg2(ji,jj)*xx(ji,jj) + rg3(ji,jj)/xx(ji,jj) & |
---|
322 | - boron(ji,jj)/( 1.0 + rg1(ji,jj)/xx(ji,jj) ) |
---|
323 | |
---|
324 | ! |
---|
325 | ! Set up the quadratic terms. |
---|
326 | ! |
---|
327 | |
---|
328 | qa = 2.0*tco2(ji,jj) - alkprime |
---|
329 | qb = -rek(ji,jj)*( alkprime - tco2(ji,jj) ) |
---|
330 | qc = -alkprime |
---|
331 | |
---|
332 | xf(ji,jj) = ( -qb + SQRT( qb*qb - 4.0*qa*qc ) )*0.5/qa |
---|
333 | |
---|
334 | |
---|
335 | ENDDO |
---|
336 | ENDDO |
---|
337 | ! |
---|
338 | ENDDO |
---|
339 | |
---|
340 | ! |
---|
341 | DO ji=1,jpi |
---|
342 | DO jj=1,jpj |
---|
343 | ! |
---|
344 | xtrma = 1.0/( 1.0 + rek(ji,jj)*xf(ji,jj) + xf(ji,jj)*xf(ji,jj) ) |
---|
345 | xtrmb = rg2(ji,jj) + rg3(ji,jj)/(xf(ji,jj)*xf(ji,jj)) & |
---|
346 | + rg1(ji,jj)*boron(ji,jj)/(( rg1(ji,jj) + xf(ji,jj) )*( rg1(ji,jj) + xf(ji,jj) )) |
---|
347 | ! |
---|
348 | dxda = 1.0/( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) & |
---|
349 | + rek(ji,jj)*xf(ji,jj)*xf(ji,jj) )*xtrma*xtrma + xtrmb ) |
---|
350 | dxdc = -dxda*xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma |
---|
351 | ! |
---|
352 | dgdc = xtrma - dxdc*tco2(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma*xtrma |
---|
353 | dgda = ( 1.0 - dxda*( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) )*xtrma & |
---|
354 | + xtrmb ) )/(xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )) |
---|
355 | ! |
---|
356 | fxa = 1.0e-3*RHO_WATER_SI*co2_sol(ji,jj) |
---|
357 | DPDC(ji,jj) = dgdc/fxa |
---|
358 | DPDA(ji,jj) = dgda/fxa |
---|
359 | ! |
---|
360 | ENDDO |
---|
361 | ENDDO |
---|
362 | |
---|
363 | RETURN |
---|
364 | END SUBROUTINE pco2_grad |
---|
365 | |
---|
366 | END MODULE asmpco2bal |
---|