1 | MODULE mocsy_mainmod |
---|
2 | |
---|
3 | USE in_out_manager ! I/O manager |
---|
4 | |
---|
5 | CONTAINS |
---|
6 | |
---|
7 | ! ---------------------------------------------------------------------- |
---|
8 | ! SW_ADTG |
---|
9 | ! ---------------------------------------------------------------------- |
---|
10 | ! |
---|
11 | !> \file sw_adtg.f90 |
---|
12 | !! \BRIEF |
---|
13 | !> Module with sw_adtg function - compute adiabatic temp. gradient from S,T,P |
---|
14 | !> Function to calculate adiabatic temperature gradient as per UNESCO 1983 routines. |
---|
15 | FUNCTION sw_adtg (s,t,p) |
---|
16 | |
---|
17 | ! ================================================================== |
---|
18 | ! Calculates adiabatic temperature gradient as per UNESCO 1983 routines. |
---|
19 | ! Armin Koehl akoehl@ucsd.edu |
---|
20 | ! ================================================================== |
---|
21 | USE mocsy_singledouble |
---|
22 | IMPLICIT NONE |
---|
23 | !> salinity [psu (PSU-78)] |
---|
24 | REAL(kind=wp) :: s |
---|
25 | !> temperature [degree C (IPTS-68)] |
---|
26 | REAL(kind=wp) :: t |
---|
27 | !> pressure [db] |
---|
28 | REAL(kind=wp) :: p |
---|
29 | |
---|
30 | REAL(kind=wp) :: a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2 |
---|
31 | REAL(kind=wp) :: sref |
---|
32 | |
---|
33 | REAL(kind=wp) :: sw_adtg |
---|
34 | |
---|
35 | sref = 35.d0 |
---|
36 | a0 = 3.5803d-5 |
---|
37 | a1 = +8.5258d-6 |
---|
38 | a2 = -6.836d-8 |
---|
39 | a3 = 6.6228d-10 |
---|
40 | |
---|
41 | b0 = +1.8932d-6 |
---|
42 | b1 = -4.2393d-8 |
---|
43 | |
---|
44 | c0 = +1.8741d-8 |
---|
45 | c1 = -6.7795d-10 |
---|
46 | c2 = +8.733d-12 |
---|
47 | c3 = -5.4481d-14 |
---|
48 | |
---|
49 | d0 = -1.1351d-10 |
---|
50 | d1 = 2.7759d-12 |
---|
51 | |
---|
52 | e0 = -4.6206d-13 |
---|
53 | e1 = +1.8676d-14 |
---|
54 | e2 = -2.1687d-16 |
---|
55 | |
---|
56 | sw_adtg = a0 + (a1 + (a2 + a3*T)*T)*T & |
---|
57 | + (b0 + b1*T)*(S-sref) & |
---|
58 | + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P & |
---|
59 | + ( e0 + (e1 + e2*T)*T )*P*P |
---|
60 | |
---|
61 | END FUNCTION sw_adtg |
---|
62 | |
---|
63 | ! ---------------------------------------------------------------------- |
---|
64 | ! SW_ADTG |
---|
65 | ! ---------------------------------------------------------------------- |
---|
66 | ! |
---|
67 | !> \file sw_ptmp.f90 |
---|
68 | !! \BRIEF |
---|
69 | !> Module with sw_ptmp function - compute potential T from in-situ T |
---|
70 | !> Function to calculate potential temperature [C] from in-situ temperature |
---|
71 | FUNCTION sw_ptmp (s,t,p,pr) |
---|
72 | |
---|
73 | ! ================================================================== |
---|
74 | ! Calculates potential temperature [C] from in-situ Temperature [C] |
---|
75 | ! From UNESCO 1983 report. |
---|
76 | ! Armin Koehl akoehl@ucsd.edu |
---|
77 | ! ================================================================== |
---|
78 | |
---|
79 | ! Input arguments: |
---|
80 | ! ------------------------------------- |
---|
81 | ! s = salinity [psu (PSS-78) ] |
---|
82 | ! t = temperature [degree C (IPTS-68)] |
---|
83 | ! p = pressure [db] |
---|
84 | ! pr = reference pressure [db] |
---|
85 | |
---|
86 | USE mocsy_singledouble |
---|
87 | IMPLICIT NONE |
---|
88 | |
---|
89 | ! Input arguments |
---|
90 | !> salinity [psu (PSS-78)] |
---|
91 | REAL(kind=wp) :: s |
---|
92 | !> temperature [degree C (IPTS-68)] |
---|
93 | REAL(kind=wp) :: t |
---|
94 | !> pressure [db] |
---|
95 | REAL(kind=wp) :: p |
---|
96 | !> reference pressure [db] |
---|
97 | REAL(kind=wp) :: pr |
---|
98 | |
---|
99 | ! local arguments |
---|
100 | REAL(kind=wp) :: del_P ,del_th, th, q |
---|
101 | REAL(kind=wp) :: onehalf, two, three |
---|
102 | PARAMETER (onehalf = 0.5d0, two = 2.d0, three = 3.d0 ) |
---|
103 | |
---|
104 | ! REAL(kind=wp) :: sw_adtg |
---|
105 | ! EXTERNAL sw_adtg |
---|
106 | |
---|
107 | ! Output |
---|
108 | REAL(kind=wp) :: sw_ptmp |
---|
109 | |
---|
110 | ! theta1 |
---|
111 | del_P = PR - P |
---|
112 | del_th = del_P*sw_adtg(S,T,P) |
---|
113 | th = T + onehalf*del_th |
---|
114 | q = del_th |
---|
115 | |
---|
116 | ! theta2 |
---|
117 | del_th = del_P*sw_adtg(S,th,P+onehalf*del_P) |
---|
118 | th = th + (1.d0 - 1.d0/SQRT(two))*(del_th - q) |
---|
119 | q = (two-SQRT(two))*del_th + (-two+three/SQRT(two))*q |
---|
120 | |
---|
121 | ! theta3 |
---|
122 | del_th = del_P*sw_adtg(S,th,P+onehalf*del_P) |
---|
123 | th = th + (1.d0 + 1.d0/SQRT(two))*(del_th - q) |
---|
124 | q = (two + SQRT(two))*del_th + (-two-three/SQRT(two))*q |
---|
125 | |
---|
126 | ! theta4 |
---|
127 | del_th = del_P*sw_adtg(S,th,P+del_P) |
---|
128 | sw_ptmp = th + (del_th - two*q)/(two*three) |
---|
129 | |
---|
130 | RETURN |
---|
131 | END FUNCTION sw_ptmp |
---|
132 | |
---|
133 | |
---|
134 | ! ---------------------------------------------------------------------- |
---|
135 | ! SW_TEMP |
---|
136 | ! ---------------------------------------------------------------------- |
---|
137 | ! |
---|
138 | !> \file sw_temp.f90 |
---|
139 | !! \BRIEF |
---|
140 | !> Module with sw_temp function - compute in-situ T from potential T |
---|
141 | !> Function to compute in-situ temperature [C] from potential temperature [C] |
---|
142 | FUNCTION sw_temp( s, t, p, pr ) |
---|
143 | ! ============================================================= |
---|
144 | ! SW_TEMP |
---|
145 | ! Computes in-situ temperature [C] from potential temperature [C] |
---|
146 | ! Routine available in seawater.f (used for MIT GCM) |
---|
147 | ! Downloaded seawater.f (on 17 April 2009) from |
---|
148 | ! http://ecco2.jpl.nasa.gov/data1/beaufort/MITgcm/bin/ |
---|
149 | ! ============================================================= |
---|
150 | |
---|
151 | ! REFERENCES: |
---|
152 | ! Fofonoff, P. and Millard, R.C. Jr |
---|
153 | ! Unesco 1983. Algorithms for computation of fundamental properties of |
---|
154 | ! seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp. |
---|
155 | ! Eqn.(31) p.39 |
---|
156 | |
---|
157 | ! Bryden, H. 1973. |
---|
158 | ! "New Polynomials for thermal expansion, adiabatic temperature gradient |
---|
159 | ! and potential temperature of sea water." |
---|
160 | ! DEEP-SEA RES., 1973, Vol20,401-408. |
---|
161 | ! ============================================================= |
---|
162 | |
---|
163 | ! Simple modifications: J. C. Orr, 16 April 2009 |
---|
164 | ! - combined fortran code from MITgcm site & simplification in |
---|
165 | ! CSIRO code (matlab equivalent) from Phil Morgan |
---|
166 | |
---|
167 | USE mocsy_singledouble |
---|
168 | IMPLICIT NONE |
---|
169 | |
---|
170 | ! Input arguments: |
---|
171 | ! ----------------------------------------------- |
---|
172 | ! s = salinity [psu (PSS-78) ] |
---|
173 | ! t = potential temperature [degree C (IPTS-68)] |
---|
174 | ! p = pressure [db] |
---|
175 | ! pr = reference pressure [db] |
---|
176 | |
---|
177 | !> salinity [psu (PSS-78)] |
---|
178 | REAL(kind=wp) :: s |
---|
179 | !> potential temperature [degree C (IPTS-68)] |
---|
180 | REAL(kind=wp) :: t |
---|
181 | !> pressure [db] |
---|
182 | REAL(kind=wp) :: p |
---|
183 | !> reference pressure [db] |
---|
184 | REAL(kind=wp) :: pr |
---|
185 | |
---|
186 | REAL(kind=wp) :: ds, dt, dp, dpr |
---|
187 | REAL(kind=wp) :: dsw_temp |
---|
188 | |
---|
189 | REAL(kind=wp) :: sw_temp |
---|
190 | ! EXTERNAL sw_ptmp |
---|
191 | ! REAL(kind=wp) :: sw_ptmp |
---|
192 | |
---|
193 | ds = DBLE(s) |
---|
194 | dt = DBLE(t) |
---|
195 | dp = DBLE(p) |
---|
196 | dpr = DBLE(pr) |
---|
197 | |
---|
198 | ! Simple solution |
---|
199 | ! (see https://svn.mpl.ird.fr/us191/oceano/tags/V0/lib/matlab/seawater/sw_temp.m) |
---|
200 | ! Carry out inverse calculation by swapping P_ref (pr) and Pressure (p) |
---|
201 | ! in routine that is normally used to compute potential temp from temp |
---|
202 | dsw_temp = sw_ptmp(ds, dt, dpr, dp) |
---|
203 | sw_temp = REAL(dsw_temp) |
---|
204 | |
---|
205 | ! The above simplification works extremely well (compared to Table in 1983 report) |
---|
206 | ! whereas the sw_temp routine from MIT GCM site does not seem to work right |
---|
207 | |
---|
208 | RETURN |
---|
209 | END FUNCTION sw_temp |
---|
210 | |
---|
211 | ! ---------------------------------------------------------------------- |
---|
212 | ! TPOT |
---|
213 | ! ---------------------------------------------------------------------- |
---|
214 | ! |
---|
215 | !> \file tpot.f90 |
---|
216 | !! \BRIEF |
---|
217 | !> Module with tpot subroutine - compute potential T from in situ T,S,P |
---|
218 | !> Compute potential temperature from arrays of in situ temp, salinity, and pressure. |
---|
219 | !! This subroutine is needed because sw_ptmp is a function (using scalars not arrays) |
---|
220 | SUBROUTINE tpot(salt, tempis, press, pressref, N, tempot) |
---|
221 | ! Purpose: |
---|
222 | ! Compute potential temperature from arrays of in situ temp, salinity, and pressure. |
---|
223 | ! Needed because sw_ptmp is a function |
---|
224 | |
---|
225 | USE mocsy_singledouble |
---|
226 | IMPLICIT NONE |
---|
227 | |
---|
228 | !> number of records |
---|
229 | INTEGER, intent(in) :: N |
---|
230 | |
---|
231 | ! INPUT variables |
---|
232 | !> salinity [psu] |
---|
233 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: salt |
---|
234 | !> in situ temperature [C] |
---|
235 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: tempis |
---|
236 | !> pressure [db] |
---|
237 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: press |
---|
238 | !f2py optional , depend(salt) :: n=len(salt) |
---|
239 | !> pressure reference level [db] |
---|
240 | REAL(kind=wp), INTENT(in) :: pressref |
---|
241 | |
---|
242 | ! OUTPUT variables: |
---|
243 | !> potential temperature [C] for pressref |
---|
244 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: tempot |
---|
245 | |
---|
246 | REAL(kind=wp) :: dsalt, dtempis, dpress, dpressref |
---|
247 | REAL(kind=wp) :: dtempot |
---|
248 | |
---|
249 | INTEGER :: i |
---|
250 | |
---|
251 | ! REAL(kind=wp) :: sw_ptmp |
---|
252 | ! EXTERNAL sw_ptmp |
---|
253 | |
---|
254 | DO i = 1,N |
---|
255 | dsalt = DBLE(salt(i)) |
---|
256 | dtempis = DBLE(tempis(i)) |
---|
257 | dpress = DBLE(press(i)) |
---|
258 | dpressref = DBLE(pressref) |
---|
259 | |
---|
260 | dtempot = sw_ptmp(dsalt, dtempis, dpress, dpressref) |
---|
261 | |
---|
262 | tempot(i) = REAL(dtempot) |
---|
263 | END DO |
---|
264 | |
---|
265 | RETURN |
---|
266 | END SUBROUTINE tpot |
---|
267 | |
---|
268 | ! ---------------------------------------------------------------------- |
---|
269 | ! TIS |
---|
270 | ! ---------------------------------------------------------------------- |
---|
271 | ! |
---|
272 | !> \file tis.f90 |
---|
273 | !! \BRIEF |
---|
274 | !> Module with tis subroutine - compute in situ T from S,T,P |
---|
275 | !> Compute in situ temperature from arrays of potential temp, salinity, and pressure. |
---|
276 | !! This subroutine is needed because sw_temp is a function (using scalars not arrays) |
---|
277 | SUBROUTINE tis(salt, tempot, press, pressref, N, tempis) |
---|
278 | ! Purpose: |
---|
279 | ! Compute in situ temperature from arrays of in situ temp, salinity, and pressure. |
---|
280 | ! Needed because sw_temp is a function |
---|
281 | |
---|
282 | USE mocsy_singledouble |
---|
283 | IMPLICIT NONE |
---|
284 | |
---|
285 | !> number of records |
---|
286 | INTEGER, intent(in) :: N |
---|
287 | |
---|
288 | ! INPUT variables |
---|
289 | !> salinity [psu] |
---|
290 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: salt |
---|
291 | !> potential temperature [C] |
---|
292 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: tempot |
---|
293 | !> pressure [db] |
---|
294 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: press |
---|
295 | !f2py optional , depend(salt) :: n=len(salt) |
---|
296 | !> pressure reference level [db] |
---|
297 | REAL(kind=wp), INTENT(in) :: pressref |
---|
298 | |
---|
299 | ! OUTPUT variables: |
---|
300 | !> in situ temperature [C] |
---|
301 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: tempis |
---|
302 | |
---|
303 | ! REAL(kind=wp) :: dsalt, dtempis, dpress, dpressref |
---|
304 | ! REAL(kind=wp) :: dtempot |
---|
305 | |
---|
306 | INTEGER :: i |
---|
307 | |
---|
308 | ! REAL(kind=wp) :: sw_temp |
---|
309 | ! REAL(kind=wp) :: sw_temp |
---|
310 | ! EXTERNAL sw_temp |
---|
311 | |
---|
312 | DO i = 1,N |
---|
313 | !dsalt = DBLE(salt(i)) |
---|
314 | !dtempot = DBLE(tempot(i)) |
---|
315 | !dpress = DBLE(press(i)) |
---|
316 | !dpressref = DBLE(pressref) |
---|
317 | !dtempis = sw_temp(dsalt, dtempot, dpress, dpressref) |
---|
318 | !tempis(i) = REAL(dtempis) |
---|
319 | |
---|
320 | tempis = sw_temp(salt(i), tempot(i), press(i), pressref) |
---|
321 | END DO |
---|
322 | |
---|
323 | RETURN |
---|
324 | END SUBROUTINE tis |
---|
325 | |
---|
326 | ! ---------------------------------------------------------------------- |
---|
327 | ! P80 |
---|
328 | ! ---------------------------------------------------------------------- |
---|
329 | ! |
---|
330 | !> \file p80.f90 |
---|
331 | !! \BRIEF |
---|
332 | !> Module with p80 function - compute pressure from depth |
---|
333 | !> Function to compute pressure from depth using Saunder's (1981) formula with eos80. |
---|
334 | FUNCTION p80(dpth,xlat) |
---|
335 | |
---|
336 | ! Compute Pressure from depth using Saunder's (1981) formula with eos80. |
---|
337 | |
---|
338 | ! Reference: |
---|
339 | ! Saunders, Peter M. (1981) Practical conversion of pressure |
---|
340 | ! to depth, J. Phys. Ooceanogr., 11, 573-574, (1981) |
---|
341 | |
---|
342 | ! Coded by: |
---|
343 | ! R. Millard |
---|
344 | ! March 9, 1983 |
---|
345 | ! check value: p80=7500.004 dbars at lat=30 deg., depth=7321.45 meters |
---|
346 | |
---|
347 | ! Modified (slight format changes + added ref. details): |
---|
348 | ! J. Orr, 16 April 2009 |
---|
349 | |
---|
350 | USE mocsy_singledouble |
---|
351 | IMPLICIT NONE |
---|
352 | |
---|
353 | ! Input variables: |
---|
354 | !> depth [m] |
---|
355 | REAL(kind=wp), INTENT(in) :: dpth |
---|
356 | !> latitude [degrees] |
---|
357 | REAL(kind=wp), INTENT(in) :: xlat |
---|
358 | |
---|
359 | ! Output variable: |
---|
360 | !> pressure [db] |
---|
361 | REAL(kind=wp) :: p80 |
---|
362 | |
---|
363 | ! Local variables: |
---|
364 | REAL(kind=wp) :: pi |
---|
365 | REAL(kind=wp) :: plat, d, c1 |
---|
366 | |
---|
367 | pi=3.141592654 |
---|
368 | |
---|
369 | plat = ABS(xlat*pi/180.) |
---|
370 | d = SIN(plat) |
---|
371 | c1 = 5.92e-3+d**2 * 5.25e-3 |
---|
372 | |
---|
373 | p80 = ((1-c1)-SQRT(((1-c1)**2)-(8.84e-6*dpth))) / 4.42e-6 |
---|
374 | |
---|
375 | RETURN |
---|
376 | END FUNCTION p80 |
---|
377 | |
---|
378 | ! ---------------------------------------------------------------------- |
---|
379 | ! RHO |
---|
380 | ! ---------------------------------------------------------------------- |
---|
381 | ! |
---|
382 | !> \file rho.f90 |
---|
383 | !! \BRIEF |
---|
384 | !> Module with rho function - computes in situ density from S, T, P |
---|
385 | !> Function to compute in situ density from salinity (psu), in situ temperature (C), & pressure (bar) |
---|
386 | FUNCTION rho(salt, temp, pbar) |
---|
387 | |
---|
388 | ! Compute in situ density from salinity (psu), in situ temperature (C), & pressure (bar) |
---|
389 | |
---|
390 | USE mocsy_singledouble |
---|
391 | IMPLICIT NONE |
---|
392 | |
---|
393 | !> salinity [psu] |
---|
394 | REAL(kind=wp) :: salt |
---|
395 | !> in situ temperature (C) |
---|
396 | REAL(kind=wp) :: temp |
---|
397 | !> pressure (bar) [Note units: this is NOT in decibars] |
---|
398 | REAL(kind=wp) :: pbar |
---|
399 | |
---|
400 | REAL(kind=wp) :: s, t, p |
---|
401 | ! REAL(kind=wp) :: t68 |
---|
402 | REAL(kind=wp) :: X |
---|
403 | REAL(kind=wp) :: rhow, rho0 |
---|
404 | REAL(kind=wp) :: a, b, c |
---|
405 | REAL(kind=wp) :: Ksbmw, Ksbm0, Ksbm |
---|
406 | REAL(kind=wp) :: drho |
---|
407 | |
---|
408 | REAL(kind=wp) :: rho |
---|
409 | |
---|
410 | ! Input arguments: |
---|
411 | ! ------------------------------------- |
---|
412 | ! s = salinity [psu (PSS-78) ] |
---|
413 | ! t = in situ temperature [degree C (IPTS-68)] |
---|
414 | ! p = pressure [bar] !!!! (not in [db] |
---|
415 | |
---|
416 | s = DBLE(salt) |
---|
417 | t = DBLE(temp) |
---|
418 | p = DBLE(pbar) |
---|
419 | |
---|
420 | ! Convert the temperature on today's "ITS 90" scale to older IPTS 68 scale |
---|
421 | ! (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote) |
---|
422 | ! According to Best-Practices guide, line above should be commented & 2 lines below should be uncommented |
---|
423 | ! Guide's answer of rho (s=35, t=25, p=0) = 1023.343 is for temperature given on ITPS-68 scale |
---|
424 | ! t68 = (T - 0.0002) / 0.99975 |
---|
425 | ! X = t68 |
---|
426 | ! Finally, don't do the ITS-90 to IPTS-68 conversion (T input var now already on IPTS-68 scale) |
---|
427 | X = T |
---|
428 | |
---|
429 | ! Density of pure water |
---|
430 | rhow = 999.842594d0 + 6.793952e-2_wp*X & |
---|
431 | -9.095290e-3_wp*X*X + 1.001685e-4_wp*X**3 & |
---|
432 | -1.120083e-6_wp*X**4 + 6.536332e-9_wp*X**5 |
---|
433 | |
---|
434 | ! Density of seawater at 1 atm, P=0 |
---|
435 | A = 8.24493e-1_wp - 4.0899e-3_wp*X & |
---|
436 | + 7.6438e-5_wp*X*X - 8.2467e-7_wp*X**3 + 5.3875e-9_wp*X**4 |
---|
437 | B = -5.72466e-3_wp + 1.0227e-4_wp*X - 1.6546e-6_wp*X*X |
---|
438 | C = 4.8314e-4_wp |
---|
439 | |
---|
440 | rho0 = rhow + A*S + B*S*SQRT(S) + C*S**2.0d0 |
---|
441 | |
---|
442 | ! Secant bulk modulus of pure water |
---|
443 | ! The secant bulk modulus is the average change in pressure |
---|
444 | ! divided by the total change in volume per unit of initial volume. |
---|
445 | Ksbmw = 19652.21d0 + 148.4206d0*X - 2.327105d0*X*X & |
---|
446 | + 1.360477e-2_wp*X**3 - 5.155288e-5_wp*X**4 |
---|
447 | |
---|
448 | ! Secant bulk modulus of seawater at 1 atm |
---|
449 | Ksbm0 = Ksbmw + S*( 54.6746d0 - 0.603459d0*X + 1.09987e-2_wp*X**2 & |
---|
450 | - 6.1670e-5_wp*X**3) & |
---|
451 | + S*SQRT(S)*( 7.944e-2_wp + 1.6483e-2_wp*X - 5.3009e-4_wp*X**2) |
---|
452 | |
---|
453 | ! Secant bulk modulus of seawater at S,T,P |
---|
454 | Ksbm = Ksbm0 & |
---|
455 | + P*(3.239908d0 + 1.43713e-3_wp*X + 1.16092e-4_wp*X**2 - 5.77905e-7_wp*X**3) & |
---|
456 | + P*S*(2.2838e-3_wp - 1.0981e-5_wp*X - 1.6078e-6_wp*X**2) & |
---|
457 | + P*S*SQRT(S)*1.91075e-4_wp & |
---|
458 | + P*P*(8.50935e-5_wp - 6.12293e-6_wp*X + 5.2787e-8_wp*X**2) & |
---|
459 | + P*P*S*(-9.9348e-7_wp + 2.0816e-8_wp*X + 9.1697e-10_wp*X**2) |
---|
460 | |
---|
461 | ! Density of seawater at S,T,P |
---|
462 | drho = rho0/(1.0d0 - P/Ksbm) |
---|
463 | rho = REAL(drho) |
---|
464 | |
---|
465 | RETURN |
---|
466 | END FUNCTION rho |
---|
467 | |
---|
468 | ! ---------------------------------------------------------------------- |
---|
469 | ! RHOINSITU |
---|
470 | ! ---------------------------------------------------------------------- |
---|
471 | ! |
---|
472 | !> \file rhoinsitu.f90 |
---|
473 | !! \BRIEF |
---|
474 | !> Module with rhoinsitu subroutine - compute in situ density from S, Tis, P |
---|
475 | !> Compute in situ density from salinity (psu), in situ temperature (C), & pressure (db). |
---|
476 | !! This subroutine is needed because rho is a function (using scalars not arrays) |
---|
477 | SUBROUTINE rhoinsitu(salt, tempis, pdbar, N, rhois) |
---|
478 | |
---|
479 | ! Purpose: |
---|
480 | ! Compute in situ density from salinity (psu), in situ temperature (C), & pressure (db) |
---|
481 | ! Needed because rho is a function |
---|
482 | |
---|
483 | USE mocsy_singledouble |
---|
484 | IMPLICIT NONE |
---|
485 | |
---|
486 | INTEGER :: N |
---|
487 | |
---|
488 | ! INPUT variables |
---|
489 | ! salt = salinity [psu] |
---|
490 | ! tempis = in situ temperature [C] |
---|
491 | ! pdbar = pressure [db] |
---|
492 | |
---|
493 | !> salinity [psu] |
---|
494 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: salt |
---|
495 | !> in situ temperature [C] |
---|
496 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: tempis |
---|
497 | !> pressure [db] |
---|
498 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: pdbar |
---|
499 | !f2py optional , depend(salt) :: n=len(salt) |
---|
500 | |
---|
501 | ! OUTPUT variables: |
---|
502 | ! rhois = in situ density |
---|
503 | |
---|
504 | !> in situ density [kg/m3] |
---|
505 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: rhois |
---|
506 | |
---|
507 | ! Local variables |
---|
508 | INTEGER :: i |
---|
509 | |
---|
510 | ! REAL(kind=wp) :: rho |
---|
511 | ! EXTERNAL rho |
---|
512 | |
---|
513 | DO i = 1,N |
---|
514 | rhois(i) = rho(salt(i), tempis(i), pdbar(i)/10.) |
---|
515 | END DO |
---|
516 | |
---|
517 | RETURN |
---|
518 | END SUBROUTINE rhoinsitu |
---|
519 | |
---|
520 | ! ---------------------------------------------------------------------- |
---|
521 | ! DEPTH2PRESS |
---|
522 | ! ---------------------------------------------------------------------- |
---|
523 | ! |
---|
524 | !> \file depth2press.f90 |
---|
525 | !! \BRIEF |
---|
526 | !> Module with depth2press subroutine - converts depth to pressure |
---|
527 | !! with Saunders (1981) formula |
---|
528 | !> Compute pressure [db] from depth [m] & latitude [degrees north]. |
---|
529 | !! This subroutine is needed because p80 is a function (using scalars not arrays) |
---|
530 | SUBROUTINE depth2press(depth, lat, pdbar, N) |
---|
531 | |
---|
532 | ! Purpose: |
---|
533 | ! Compute pressure [db] from depth [m] & latitude [degrees north]. |
---|
534 | ! Needed because p80 is a function |
---|
535 | |
---|
536 | USE mocsy_singledouble |
---|
537 | IMPLICIT NONE |
---|
538 | |
---|
539 | !> number of records |
---|
540 | INTEGER, intent(in) :: N |
---|
541 | |
---|
542 | ! INPUT variables |
---|
543 | !> depth [m] |
---|
544 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: depth |
---|
545 | !> latitude [degrees] |
---|
546 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: lat |
---|
547 | !f2py optional , depend(depth) :: n=len(depth) |
---|
548 | |
---|
549 | ! OUTPUT variables: |
---|
550 | !> pressure [db] |
---|
551 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: pdbar |
---|
552 | |
---|
553 | ! Local variables |
---|
554 | INTEGER :: i |
---|
555 | |
---|
556 | ! REAL(kind=wp) :: p80 |
---|
557 | ! EXTERNAL p80 |
---|
558 | |
---|
559 | DO i = 1,N |
---|
560 | pdbar(i) = p80(depth(i), lat(i)) |
---|
561 | END DO |
---|
562 | |
---|
563 | RETURN |
---|
564 | END SUBROUTINE depth2press |
---|
565 | |
---|
566 | ! ---------------------------------------------------------------------- |
---|
567 | ! CONSTANTS |
---|
568 | ! ---------------------------------------------------------------------- |
---|
569 | ! |
---|
570 | !> \file constants.f90 |
---|
571 | !! \BRIEF |
---|
572 | !> Module with contants subroutine - computes carbonate system constants |
---|
573 | !! from S,T,P |
---|
574 | !> Compute thermodynamic constants |
---|
575 | !! FROM temperature, salinity, and pressure (1D arrays) |
---|
576 | SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, & |
---|
577 | K1p, K2p, K3p, Ksi, & |
---|
578 | St, Ft, Bt, & |
---|
579 | temp, sal, Patm, & |
---|
580 | depth, lat, N, & |
---|
581 | optT, optP, optB, optK1K2, optKf, optGAS) |
---|
582 | |
---|
583 | ! Purpose: |
---|
584 | ! Compute thermodynamic constants |
---|
585 | ! FROM: temperature, salinity, and pressure (1D arrays) |
---|
586 | |
---|
587 | ! INPUT variables: |
---|
588 | ! ================ |
---|
589 | ! Patm = atmospheric pressure [atm] |
---|
590 | ! depth = depth [m] (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure) |
---|
591 | ! = pressure [db] (with optP='db') |
---|
592 | ! lat = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m') |
---|
593 | ! = dummy array (unused when optP='db') |
---|
594 | ! temp = potential temperature [degrees C] (with optT='Tpot', i.e., models carry tempot, not temp) |
---|
595 | ! = in situ temperature [degrees C] (with optT='Tinsitu', e.g., for data) |
---|
596 | ! sal = salinity in [psu] |
---|
597 | ! --------- |
---|
598 | ! optT: choose in situ vs. potential temperature as input |
---|
599 | ! --------- |
---|
600 | ! NOTE: Carbonate chem calculations require IN-SITU temperature (not potential Temperature) |
---|
601 | ! -> 'Tpot' means input is pot. Temperature (in situ Temp "tempis" is computed) |
---|
602 | ! -> 'Tinsitu' means input is already in-situ Temperature, not pot. Temp ("tempis" not computed) |
---|
603 | ! --------- |
---|
604 | ! optP: choose depth (m) vs pressure (db) as input |
---|
605 | ! --------- |
---|
606 | ! -> 'm' means "depth" input is in "m" (thus in situ Pressure "p" [db] is computed) |
---|
607 | ! -> 'db' means "depth" input is already in situ pressure [db], not m (p = depth) |
---|
608 | ! --------- |
---|
609 | ! optB: |
---|
610 | ! --------- |
---|
611 | ! -> 'u74' means use classic formulation of Uppström (1974) for total Boron |
---|
612 | ! -> 'l10' means use newer formulation of Lee et al. (2010) for total Boron |
---|
613 | ! --------- |
---|
614 | ! optK1K2: |
---|
615 | ! --------- |
---|
616 | ! -> 'l' means use Lueker et al. (2000) formulations for K1 & K2 (recommended by Dickson et al. 2007) |
---|
617 | ! **** BUT this should only be used when 2 < T < 35 and 19 < S < 43 |
---|
618 | ! -> 'm10' means use Millero (2010) formulation for K1 & K2 (see Dickson et al., 2007) |
---|
619 | ! **** Valid for 0 < T < 50°C and 1 < S < 50 psu |
---|
620 | ! -> 'w14' means use Waters (2014) formulation for K1 & K2 (see Dickson et al., 2007) |
---|
621 | ! **** Valid for 0 < T < 50°C and 1 < S < 50 psu |
---|
622 | ! ----------- |
---|
623 | ! optKf: |
---|
624 | ! ---------- |
---|
625 | ! -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007) |
---|
626 | ! **** BUT Valid for 9 < T < 33°C and 10 < S < 40. |
---|
627 | ! -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994) |
---|
628 | ! ----------- |
---|
629 | ! optGAS: choose in situ vs. potential fCO2 and pCO2 |
---|
630 | ! --------- |
---|
631 | ! PRESSURE corrections for K0 and the fugacity coefficient (Cf) |
---|
632 | ! -> 'Pzero' = 'zero order' fCO2 and pCO2 (typical approach, which is flawed) |
---|
633 | ! considers in situ T & only atm pressure (hydrostatic=0) |
---|
634 | ! -> 'Ppot' = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) |
---|
635 | ! considers potential T & only atm pressure (hydrostatic press = 0) |
---|
636 | ! -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure) |
---|
637 | ! considers in situ T & total pressure (atm + hydrostatic) |
---|
638 | ! --------- |
---|
639 | |
---|
640 | ! OUTPUT variables: |
---|
641 | ! ================= |
---|
642 | ! K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi |
---|
643 | ! St, Ft, Bt |
---|
644 | |
---|
645 | USE mocsy_singledouble |
---|
646 | IMPLICIT NONE |
---|
647 | |
---|
648 | ! Input variables |
---|
649 | !> number of records |
---|
650 | INTEGER, INTENT(in) :: N |
---|
651 | !> in <b>situ temperature</b> (when optT='Tinsitu', typical data) |
---|
652 | !! OR <b>potential temperature</b> (when optT='Tpot', typical models) [degree C] |
---|
653 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp |
---|
654 | !> depth in <b>meters</b> (when optP='m') or <b>decibars</b> (when optP='db') |
---|
655 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: depth |
---|
656 | !> latitude <b>[degrees north]</b> |
---|
657 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: lat |
---|
658 | !> salinity <b>[psu]</b> |
---|
659 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: sal |
---|
660 | !f2py optional , depend(sal) :: n=len(sal) |
---|
661 | |
---|
662 | !> atmospheric pressure <b>[atm]</b> |
---|
663 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm |
---|
664 | |
---|
665 | !> for temp input, choose \b 'Tinsitu' for in situ Temp or |
---|
666 | !! \b 'Tpot' for potential temperature (in situ Temp is computed, needed for models) |
---|
667 | CHARACTER(7), INTENT(in) :: optT |
---|
668 | !> for depth input, choose \b "db" for decibars (in situ pressure) or \b "m" for meters (pressure is computed, needed for models) |
---|
669 | CHARACTER(2), INTENT(in) :: optP |
---|
670 | !> for total boron, choose either \b 'u74' (Uppstrom, 1974) or \b 'l10' (Lee et al., 2010). |
---|
671 | !! The 'l10' formulation is based on 139 measurements (instead of 20), |
---|
672 | !! uses a more accurate method, and |
---|
673 | !! generally increases total boron in seawater by 4% |
---|
674 | !f2py character*3 optional, intent(in) :: optB='l10' |
---|
675 | CHARACTER(3), OPTIONAL, INTENT(in) :: optB |
---|
676 | !> for Kf, choose either \b 'pf' (Perez & Fraga, 1987) or \b 'dg' (Dickson & Riley, 1979) |
---|
677 | !f2py character*2 optional, intent(in) :: optKf='pf' |
---|
678 | CHARACTER(2), OPTIONAL, INTENT(in) :: optKf |
---|
679 | !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010) or \b 'w14' (Waters et al., 2014) |
---|
680 | !f2py character*3 optional, intent(in) :: optK1K2='l' |
---|
681 | CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2 |
---|
682 | !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction) |
---|
683 | !! 'Ppot' - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) |
---|
684 | !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf |
---|
685 | !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean |
---|
686 | !f2py character*7 optional, intent(in) :: optGAS='Pinsitu' |
---|
687 | CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS |
---|
688 | |
---|
689 | ! Ouput variables |
---|
690 | !> solubility of CO2 in seawater (Weiss, 1974), also known as K0 |
---|
691 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: K0 |
---|
692 | !> K1 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 |
---|
693 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: K1 |
---|
694 | !> K2 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 |
---|
695 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: K2 |
---|
696 | !> equilibrium constant for dissociation of boric acid |
---|
697 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kb |
---|
698 | !> equilibrium constant for the dissociation of water (Millero, 1995) |
---|
699 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kw |
---|
700 | !> equilibrium constant for the dissociation of bisulfate (Dickson, 1990) |
---|
701 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ks |
---|
702 | !> equilibrium constant for the dissociation of hydrogen fluoride |
---|
703 | !! either from Dickson and Riley (1979) or from Perez and Fraga (1987), depending on optKf |
---|
704 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kf |
---|
705 | !> solubility product for calcite (Mucci, 1983) |
---|
706 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kspc |
---|
707 | !> solubility product for aragonite (Mucci, 1983) |
---|
708 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kspa |
---|
709 | !> 1st dissociation constant for phosphoric acid (Millero, 1995) |
---|
710 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: K1p |
---|
711 | !> 2nd dissociation constant for phosphoric acid (Millero, 1995) |
---|
712 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: K2p |
---|
713 | !> 3rd dissociation constant for phosphoric acid (Millero, 1995) |
---|
714 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: K3p |
---|
715 | !> equilibrium constant for the dissociation of silicic acid (Millero, 1995) |
---|
716 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ksi |
---|
717 | !> total sulfate (Morris & Riley, 1966) |
---|
718 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: St |
---|
719 | !> total fluoride (Riley, 1965) |
---|
720 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ft |
---|
721 | !> total boron |
---|
722 | !! from either Uppstrom (1974) or Lee et al. (2010), depending on optB |
---|
723 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: Bt |
---|
724 | |
---|
725 | ! Local variables |
---|
726 | REAL(kind=wp) :: ssal |
---|
727 | REAL(kind=wp) :: p |
---|
728 | REAL(kind=wp) :: tempot, tempis68, tempot68 |
---|
729 | REAL(kind=wp) :: tempis |
---|
730 | REAL(kind=wp) :: is, invtk, dlogtk, is2, s2, sqrtis |
---|
731 | REAL(kind=wp) :: Ks_0p, Kf_0p |
---|
732 | REAL(kind=wp) :: total2free, free2SWS, total2SWS, SWS2total |
---|
733 | REAL(kind=wp) :: total2free_0p, free2SWS_0p, total2SWS_0p |
---|
734 | ! REAL(kind=wp) :: free2SWS, free2SWS_0p |
---|
735 | |
---|
736 | REAL(kind=wp) :: dtempot, dtempot68 |
---|
737 | REAL(kind=wp) :: R |
---|
738 | |
---|
739 | REAL(kind=wp) :: pK1o, ma1, mb1, mc1, pK1 |
---|
740 | REAL(kind=wp) :: pK2o, ma2, mb2, mc2, pK2 |
---|
741 | |
---|
742 | REAL(kind=wp), DIMENSION(12) :: a0, a1, a2, b0, b1, b2 |
---|
743 | REAL(kind=wp), DIMENSION(12) :: deltav, deltak, lnkpok0 |
---|
744 | REAL(kind=wp) :: tmp, nK0we74 |
---|
745 | |
---|
746 | INTEGER :: i, icount, ipc |
---|
747 | |
---|
748 | REAL(kind=wp) :: t, tk, tk0, prb |
---|
749 | REAL(kind=wp) :: s, sqrts, s15, scl |
---|
750 | |
---|
751 | REAL(kind=wp) :: Phydro_atm, Patmd, Ptot, Rgas_atm, vbarCO2 |
---|
752 | |
---|
753 | ! Arrays to pass optional arguments into or use defaults (Dickson et al., 2007) |
---|
754 | CHARACTER(3) :: opB |
---|
755 | CHARACTER(2) :: opKf |
---|
756 | CHARACTER(3) :: opK1K2 |
---|
757 | CHARACTER(7) :: opGAS |
---|
758 | |
---|
759 | ! CONSTANTS |
---|
760 | ! ========= |
---|
761 | ! Constants in formulation for Pressure effect on K's (Millero, 95) |
---|
762 | ! with corrected coefficients for Kb, Kw, Ksi, etc. |
---|
763 | |
---|
764 | ! index: 1) K1 , 2) K2, 3) Kb, 4) Kw, 5) Ks, 6) Kf, 7) Kspc, 8) Kspa, |
---|
765 | ! 9) K1P, 10) K2P, 11) K3P, 12) Ksi |
---|
766 | |
---|
767 | DATA a0 /-25.5_wp, -15.82_wp, -29.48_wp, -20.02_wp, & |
---|
768 | -18.03_wp, -9.78_wp, -48.76_wp, -45.96_wp, & |
---|
769 | -14.51_wp, -23.12_wp, -26.57_wp, -29.48_wp/ |
---|
770 | DATA a1 /0.1271_wp, -0.0219_wp, 0.1622_wp, 0.1119_wp, & |
---|
771 | 0.0466_wp, -0.0090_wp, 0.5304_wp, 0.5304_wp, & |
---|
772 | 0.1211_wp, 0.1758_wp, 0.2020_wp, 0.1622_wp/ |
---|
773 | DATA a2 / 0.0_wp, 0.0_wp, -2.608e-3_wp, -1.409e-3_wp, & |
---|
774 | 0.316e-3_wp, -0.942e-3_wp, 0.0_wp, 0.0_wp, & |
---|
775 | -0.321e-3_wp, -2.647e-3_wp, -3.042e-3_wp, -2.6080e-3_wp/ |
---|
776 | DATA b0 /-3.08e-3_wp, 1.13e-3_wp, -2.84e-3_wp, -5.13e-3_wp, & |
---|
777 | -4.53e-3_wp, -3.91e-3_wp, -11.76e-3_wp, -11.76e-3_wp, & |
---|
778 | -2.67e-3_wp, -5.15e-3_wp, -4.08e-3_wp, -2.84e-3_wp/ |
---|
779 | DATA b1 /0.0877e-3_wp, -0.1475e-3_wp, 0.0_wp, 0.0794e-3_wp, & |
---|
780 | 0.09e-3_wp, 0.054e-3_wp, 0.3692e-3_wp, 0.3692e-3_wp, & |
---|
781 | 0.0427e-3_wp, 0.09e-3_wp, 0.0714e-3_wp, 0.0_wp/ |
---|
782 | DATA b2 /12*0.0_wp/ |
---|
783 | |
---|
784 | ! Set defaults for optional arguments (in Fortran 90) |
---|
785 | ! Note: Optional arguments with f2py (python) are set above with |
---|
786 | ! the !f2py statements that precede each type declaraion |
---|
787 | IF (PRESENT(optB)) THEN |
---|
788 | opB = optB |
---|
789 | ELSE |
---|
790 | opB = 'l10' |
---|
791 | ENDIF |
---|
792 | IF (PRESENT(optKf)) THEN |
---|
793 | opKf = optKf |
---|
794 | ELSE |
---|
795 | opKf = 'pf' |
---|
796 | ENDIF |
---|
797 | IF (PRESENT(optK1K2)) THEN |
---|
798 | opK1K2 = optK1K2 |
---|
799 | ELSE |
---|
800 | opK1K2 = 'l' |
---|
801 | ENDIF |
---|
802 | IF (PRESENT(optGAS)) THEN |
---|
803 | opGAS = optGAS |
---|
804 | ELSE |
---|
805 | opGAS = 'Pinsitu' |
---|
806 | ENDIF |
---|
807 | |
---|
808 | R = 83.14472_wp |
---|
809 | |
---|
810 | icount = 0 |
---|
811 | DO i = 1, N |
---|
812 | icount = icount + 1 |
---|
813 | ! =============================================================== |
---|
814 | ! Convert model depth -> press; convert model Theta -> T in situ |
---|
815 | ! =============================================================== |
---|
816 | ! * Model temperature tracer is usually "potential temperature" |
---|
817 | ! * Model vertical grid is usually in meters |
---|
818 | ! BUT carbonate chem routines require pressure & in-situ T |
---|
819 | ! Thus before computing chemistry, if appropriate, |
---|
820 | ! convert these 2 model vars (input to this routine) |
---|
821 | ! - depth [m] => convert to pressure [db] |
---|
822 | ! - potential temperature (C) => convert to in-situ T (C) |
---|
823 | ! ------------------------------------------------------- |
---|
824 | ! 1) Compute pressure [db] from depth [m] and latitude [degrees] (if input is m, for models) |
---|
825 | IF (trim(optP) == 'm' ) THEN |
---|
826 | ! Compute pressure [db] from depth [m] and latitude [degrees] |
---|
827 | p = p80(depth(i), lat(i)) |
---|
828 | ELSEIF (trim(optP) == 'db' ) THEN |
---|
829 | ! In this case (where optP = 'db'), p is input & output (no depth->pressure conversion needed) |
---|
830 | p = depth(i) |
---|
831 | ELSE |
---|
832 | PRINT *,"optP must be 'm' or 'db'" |
---|
833 | STOP |
---|
834 | ENDIF |
---|
835 | |
---|
836 | ! 2) Convert potential T to in-situ T (if input is Tpot, i.e. case for models): |
---|
837 | IF (trim(optT) == 'Tpot' .OR. trim(optT) == 'tpot') THEN |
---|
838 | tempot = temp(i) |
---|
839 | ! This is the case for most models and some data |
---|
840 | ! a) Convert the pot. temp on today's "ITS 90" scale to older IPTS 68 scale |
---|
841 | ! (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote) |
---|
842 | tempot68 = (tempot - 0.0002) / 0.99975 |
---|
843 | ! b) Compute "in-situ Temperature" from "Potential Temperature" (both on IPTS 68) |
---|
844 | tempis68 = sw_temp(sal(i), tempot68, p, 0.0_wp ) |
---|
845 | ! c) Convert the in-situ temp on older IPTS 68 scale to modern scale (ITS 90) |
---|
846 | tempis = 0.99975*tempis68 + 0.0002 |
---|
847 | ! Note: parts (a) and (c) above are tiny corrections; |
---|
848 | ! part (b) is a big correction for deep waters (but zero at surface) |
---|
849 | ELSEIF (trim(optT) == 'Tinsitu' .OR. trim(optT) == 'tinsitu') THEN |
---|
850 | ! When optT = 'Tinsitu', tempis is input & output (no tempot needed) |
---|
851 | tempis = temp(i) |
---|
852 | tempis68 = (temp(i) - 0.0002) / 0.99975 |
---|
853 | ! dtempot68 = sw_ptmp(DBLE(sal(i)), DBLE(tempis68), DBLE(p), 0.0_wp) |
---|
854 | dtempot68 = sw_ptmp(sal(i), tempis68, p, 0.0_wp) |
---|
855 | dtempot = 0.99975*dtempot68 + 0.0002 |
---|
856 | ELSE |
---|
857 | PRINT *,"optT must be either 'Tpot' or 'Tinsitu'" |
---|
858 | PRINT *,"you specified optT =", trim(optT) |
---|
859 | STOP |
---|
860 | ENDIF |
---|
861 | |
---|
862 | ! Compute constants: |
---|
863 | IF (temp(i) >= -5. .AND. temp(i) < 1.0e+2) THEN |
---|
864 | ! Test to indicate if any of input variables are unreasonable |
---|
865 | IF ( sal(i) < 0. .OR. sal(i) > 1e+3) THEN |
---|
866 | PRINT *, 'i, icount, temp, sal =', i, icount, temp(i), sal(i) |
---|
867 | ENDIF |
---|
868 | ! Zero out negative salinity (prev case for OCMIP2 model w/ slightly negative S in some coastal cells) |
---|
869 | IF (sal(i) < 0.0) THEN |
---|
870 | ssal = 0.0 |
---|
871 | ELSE |
---|
872 | ssal = sal(i) |
---|
873 | ENDIF |
---|
874 | |
---|
875 | ! Absolute temperature (Kelvin) and related values |
---|
876 | t = DBLE(tempis) |
---|
877 | tk = 273.15d0 + t |
---|
878 | invtk=1.0d0/tk |
---|
879 | dlogtk=LOG(tk) |
---|
880 | |
---|
881 | ! Atmospheric pressure |
---|
882 | Patmd = DBLE(Patm(i)) |
---|
883 | |
---|
884 | ! Hydrostatic pressure (prb is in bars) |
---|
885 | prb = DBLE(p) / 10.0d0 |
---|
886 | |
---|
887 | ! Salinity and simply related values |
---|
888 | s = DBLE(ssal) |
---|
889 | s2=s*s |
---|
890 | sqrts=SQRT(s) |
---|
891 | s15=s**1.5d0 |
---|
892 | scl=s/1.80655d0 |
---|
893 | |
---|
894 | ! Ionic strength: |
---|
895 | is = 19.924d0*s/(1000.0d0 - 1.005d0*s) |
---|
896 | is2 = is*is |
---|
897 | sqrtis = SQRT(is) |
---|
898 | |
---|
899 | ! Total concentrations for sulfate, fluoride, and boron |
---|
900 | |
---|
901 | ! Sulfate: Morris & Riley (1966) |
---|
902 | St(i) = 0.14d0 * scl/96.062d0 |
---|
903 | |
---|
904 | ! Fluoride: Riley (1965) |
---|
905 | Ft(i) = 0.000067d0 * scl/18.9984d0 |
---|
906 | |
---|
907 | ! Boron: |
---|
908 | IF (trim(opB) == 'l10') THEN |
---|
909 | ! New formulation from Lee et al (2010) |
---|
910 | Bt(i) = 0.0002414d0 * scl/10.811d0 |
---|
911 | ELSEIF (trim(opB) == 'u74') THEN |
---|
912 | ! Classic formulation from Uppström (1974) |
---|
913 | Bt(i) = 0.000232d0 * scl/10.811d0 |
---|
914 | ELSE |
---|
915 | PRINT *,"optB must be 'l10' or 'u74'" |
---|
916 | STOP |
---|
917 | ENDIF |
---|
918 | |
---|
919 | ! K0 (K Henry) |
---|
920 | ! CO2(g) <-> CO2(aq.) |
---|
921 | ! K0 = [CO2]/ fCO2 |
---|
922 | ! Weiss (1974) [mol/kg/atm] |
---|
923 | IF (trim(opGAS) == 'Pzero' .OR. trim(opGAS) == 'pzero') THEN |
---|
924 | tk0 = tk !in situ temperature (K) for K0 calculation |
---|
925 | Ptot = Patmd !total pressure (in atm) = atmospheric pressure ONLY |
---|
926 | ELSEIF (trim(opGAS) == 'Ppot' .OR. trim(opGAS) == 'ppot') THEN |
---|
927 | tk0 = dtempot + 273.15d0 !potential temperature (K) for K0 calculation as needed for potential fCO2 & pCO2 |
---|
928 | Ptot = Patmd !total pressure (in atm) = atmospheric pressure ONLY |
---|
929 | ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN |
---|
930 | tk0 = tk !in situ temperature (K) for K0 calculation |
---|
931 | Phydro_atm = prb / 1.01325d0 !convert hydrostatic pressure from bar to atm (1.01325 bar / atm) |
---|
932 | Ptot = Patmd + Phydro_atm !total pressure (in atm) = atmospheric pressure + hydrostatic pressure |
---|
933 | ELSE |
---|
934 | PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'" |
---|
935 | STOP |
---|
936 | ENDIF |
---|
937 | tmp = 9345.17d0/tk0 - 60.2409d0 + 23.3585d0 * LOG(tk0/100.0d0) |
---|
938 | nK0we74 = tmp + s*(0.023517d0 - 0.00023656d0*tk0 + 0.0047036e-4_wp*tk0*tk0) |
---|
939 | K0(i) = EXP(nK0we74) |
---|
940 | |
---|
941 | ! K1 = [H][HCO3]/[H2CO3] |
---|
942 | ! K2 = [H][CO3]/[HCO3] |
---|
943 | IF (trim(opK1K2) == 'l') THEN |
---|
944 | ! Mehrbach et al. (1973) refit, by Lueker et al. (2000) (total pH scale) |
---|
945 | K1(i) = 10.0d0**(-1.0d0*(3633.86d0*invtk - 61.2172d0 + 9.6777d0*dlogtk & |
---|
946 | - 0.011555d0*s + 0.0001152d0*s2)) |
---|
947 | K2(i) = 10.0d0**(-1*(471.78d0*invtk + 25.9290d0 - 3.16967d0*dlogtk & |
---|
948 | - 0.01781d0*s + 0.0001122d0*s2)) |
---|
949 | ELSEIF (trim(opK1K2) == 'm10') THEN |
---|
950 | ! Millero (2010, Mar. Fresh Wat. Res.) (seawater pH scale) |
---|
951 | pK1o = 6320.813d0*invtk + 19.568224d0*dlogtk -126.34048d0 |
---|
952 | ma1 = 13.4038d0*sqrts + 0.03206d0*s - (5.242e-5)*s2 |
---|
953 | mb1 = -530.659d0*sqrts - 5.8210d0*s |
---|
954 | mc1 = -2.0664d0*sqrts |
---|
955 | pK1 = pK1o + ma1 + mb1*invtk + mc1*dlogtk |
---|
956 | K1(i) = 10.0d0**(-pK1) |
---|
957 | |
---|
958 | pK2o = 5143.692d0*invtk + 14.613358d0*dlogtk -90.18333d0 |
---|
959 | ma2 = 21.3728d0*sqrts + 0.1218d0*s - (3.688e-4)*s2 |
---|
960 | mb2 = -788.289d0*sqrts - 19.189d0*s |
---|
961 | mc2 = -3.374d0*sqrts |
---|
962 | pK2 = pK2o + ma2 + mb2*invtk + mc2*dlogtk |
---|
963 | K2(i) = 10.0d0**(-pK2) |
---|
964 | ELSEIF (trim(opK1K2) == 'w14') THEN |
---|
965 | ! Waters, Millero, Woosley (Mar. Chem., 165, 66-67, 2014) (seawater scale) |
---|
966 | pK1o = 6320.813d0*invtk + 19.568224d0*dlogtk -126.34048d0 |
---|
967 | ma1 = 13.409160d0*sqrts + 0.031646d0*s - (5.1895e-5)*s2 |
---|
968 | mb1 = -531.3642d0*sqrts - 5.713d0*s |
---|
969 | mc1 = -2.0669166d0*sqrts |
---|
970 | pK1 = pK1o + ma1 + mb1*invtk + mc1*dlogtk |
---|
971 | K1(i) = 10.0d0**(-pK1) |
---|
972 | |
---|
973 | pK2o = 5143.692d0*invtk + 14.613358d0*dlogtk -90.18333d0 |
---|
974 | ma2 = 21.225890d0*sqrts + 0.12450870d0*s - (3.7243e-4_r8)*s2 |
---|
975 | mb2 = -779.3444d0*sqrts - 19.91739d0*s |
---|
976 | mc2 = -3.3534679d0*sqrts |
---|
977 | pK2 = pK2o + ma2 + mb2*invtk + mc2*dlogtk |
---|
978 | K2(i) = 10.0d0**(-pK2) |
---|
979 | ELSE |
---|
980 | PRINT *, "optK1K2 must be either 'l' or 'm10', or 'w14'" |
---|
981 | STOP |
---|
982 | ENDIF |
---|
983 | |
---|
984 | ! Kb = [H][BO2]/[HBO2] |
---|
985 | ! (total scale) |
---|
986 | ! Millero p.669 (1995) using data from Dickson (1990) |
---|
987 | Kb(i) = EXP((-8966.90d0 - 2890.53d0*sqrts - 77.942d0*s + & |
---|
988 | 1.728d0*s15 - 0.0996d0*s2)*invtk + & |
---|
989 | (148.0248d0 + 137.1942d0*sqrts + 1.62142d0*s) + & |
---|
990 | (-24.4344d0 - 25.085d0*sqrts - 0.2474d0*s) * & |
---|
991 | dlogtk + 0.053105d0*sqrts*tk) |
---|
992 | |
---|
993 | ! K1p = [H][H2PO4]/[H3PO4] |
---|
994 | ! (seawater scale) |
---|
995 | ! DOE(1994) eq 7.2.20 with footnote using data from Millero (1974) |
---|
996 | ! Millero (1995), p.670, eq. 65 |
---|
997 | ! Use Millero equation's 115.540 constant instead of 115.525 (Dickson et al., 2007). |
---|
998 | ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) |
---|
999 | ! And we want to stay on the SWS scale anyway for the pressure correction later. |
---|
1000 | K1p(i) = EXP(-4576.752d0*invtk + 115.540d0 - 18.453d0*dlogtk + & |
---|
1001 | (-106.736d0*invtk + 0.69171d0) * sqrts + & |
---|
1002 | (-0.65643d0*invtk - 0.01844d0) * s) |
---|
1003 | |
---|
1004 | ! K2p = [H][HPO4]/[H2PO4] |
---|
1005 | ! (seawater scale) |
---|
1006 | ! DOE(1994) eq 7.2.23 with footnote using data from Millero (1974)) |
---|
1007 | ! Millero (1995), p.670, eq. 66 |
---|
1008 | ! Use Millero equation's 172.1033 constant instead of 172.0833 (Dickson et al., 2007). |
---|
1009 | ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) |
---|
1010 | ! And we want to stay on the SWS scale anyway for the pressure correction later. |
---|
1011 | K2p(i) = EXP(-8814.715d0*invtk + 172.1033d0 - 27.927d0*dlogtk + & |
---|
1012 | (-160.340d0*invtk + 1.3566d0)*sqrts + & |
---|
1013 | (0.37335d0*invtk - 0.05778d0)*s) |
---|
1014 | |
---|
1015 | ! K3p = [H][PO4]/[HPO4] |
---|
1016 | ! (seawater scale) |
---|
1017 | ! DOE(1994) eq 7.2.26 with footnote using data from Millero (1974) |
---|
1018 | ! Millero (1995), p.670, eq. 67 |
---|
1019 | ! Use Millero equation's 18.126 constant instead of 18.141 (Dickson et al., 2007). |
---|
1020 | ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) |
---|
1021 | ! And we want to stay on the SWS scale anyway for the pressure correction later. |
---|
1022 | K3p(i) = EXP(-3070.75d0*invtk - 18.126d0 + & |
---|
1023 | (17.27039d0*invtk + 2.81197d0) * & |
---|
1024 | sqrts + (-44.99486d0*invtk - 0.09984d0) * s) |
---|
1025 | |
---|
1026 | ! Ksi = [H][SiO(OH)3]/[Si(OH)4] |
---|
1027 | ! (seawater scale) |
---|
1028 | ! Millero (1995), p.671, eq. 72 |
---|
1029 | ! Use Millero equation's 117.400 constant instead of 117.385 (Dickson et al., 2007). |
---|
1030 | ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) |
---|
1031 | ! And we want to stay on the SWS scale anyway for the pressure correction later. |
---|
1032 | Ksi(i) = EXP(-8904.2d0*invtk + 117.400d0 - 19.334d0*dlogtk + & |
---|
1033 | (-458.79d0*invtk + 3.5913d0) * sqrtis + & |
---|
1034 | (188.74d0*invtk - 1.5998d0) * is + & |
---|
1035 | (-12.1652d0*invtk + 0.07871d0) * is2 + & |
---|
1036 | LOG(1.0 - 0.001005d0*s)) |
---|
1037 | |
---|
1038 | ! Kw = [H][OH] |
---|
1039 | ! (seawater scale) |
---|
1040 | ! Millero (1995) p.670, eq. 63 from composite data |
---|
1041 | ! Use Millero equation's 148.9802 constant instead of 148.9652 (Dickson et al., 2007). |
---|
1042 | ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) |
---|
1043 | ! And we want to stay on the SWS scale anyway for the pressure correction later. |
---|
1044 | Kw(i) = EXP(-13847.26d0*invtk + 148.9802d0 - 23.6521d0*dlogtk + & |
---|
1045 | (118.67d0*invtk - 5.977d0 + 1.0495d0 * dlogtk) * & |
---|
1046 | sqrts - 0.01615d0 * s) |
---|
1047 | |
---|
1048 | ! Ks = [H][SO4]/[HSO4] |
---|
1049 | ! (free scale) |
---|
1050 | ! Dickson (1990, J. chem. Thermodynamics 22, 113) |
---|
1051 | Ks_0p = EXP(-4276.1d0*invtk + 141.328d0 - 23.093d0*dlogtk & |
---|
1052 | + (-13856.d0*invtk + 324.57d0 - 47.986d0*dlogtk) * sqrtis & |
---|
1053 | + (35474.d0*invtk - 771.54 + 114.723d0*dlogtk) * is & |
---|
1054 | - 2698.d0*invtk*is**1.5 + 1776.d0*invtk*is2 & |
---|
1055 | + LOG(1.0d0 - 0.001005d0*s)) |
---|
1056 | |
---|
1057 | ! Kf = [H][F]/[HF] |
---|
1058 | ! (total scale) |
---|
1059 | IF (trim(opKf) == 'dg') THEN |
---|
1060 | ! Dickson and Riley (1979) -- change pH scale to total (following Dickson & Goyet, 1994) |
---|
1061 | Kf_0p = EXP(1590.2d0*invtk - 12.641d0 + 1.525d0*sqrtis + & |
---|
1062 | LOG(1.0d0 - 0.001005d0*s) + & |
---|
1063 | LOG(1.0d0 + St(i)/Ks_0p)) |
---|
1064 | ELSEIF (trim(opKf) == 'pf') THEN |
---|
1065 | ! Perez and Fraga (1987) - Already on Total scale (no need for last line above) |
---|
1066 | ! Formulation as given in Dickson et al. (2007) |
---|
1067 | Kf_0p = EXP(874.d0*invtk - 9.68d0 + 0.111d0*sqrts) |
---|
1068 | ELSE |
---|
1069 | PRINT *, "optKf must be either 'dg' or 'pf'" |
---|
1070 | STOP |
---|
1071 | ENDIF |
---|
1072 | |
---|
1073 | ! Kspc (calcite) - apparent solubility product of calcite |
---|
1074 | ! (no scale) |
---|
1075 | ! Kspc = [Ca2+] [CO32-] when soln is in equilibrium w/ calcite |
---|
1076 | ! Mucci 1983 mol/kg-soln |
---|
1077 | Kspc(i) = 10d0**(-171.9065d0 - 0.077993d0*tk + 2839.319d0/tk & |
---|
1078 | + 71.595d0*LOG10(tk) & |
---|
1079 | + (-0.77712d0 + 0.0028426d0*tk + 178.34d0/tk)*sqrts & |
---|
1080 | -0.07711d0*s + 0.0041249d0*s15 ) |
---|
1081 | |
---|
1082 | |
---|
1083 | ! Kspa (aragonite) - apparent solubility product of aragonite |
---|
1084 | ! (no scale) |
---|
1085 | ! Kspa = [Ca2+] [CO32-] when soln is in equilibrium w/ aragonite |
---|
1086 | ! Mucci 1983 mol/kg-soln |
---|
1087 | Kspa(i) = 10.d0**(-171.945d0 - 0.077993d0*tk + 2903.293d0/tk & |
---|
1088 | +71.595d0*LOG10(tk) & |
---|
1089 | +(-0.068393d0 + 0.0017276d0*tk + 88.135d0/tk)*sqrts & |
---|
1090 | -0.10018d0*s + 0.0059415d0*s15 ) |
---|
1091 | |
---|
1092 | ! Pressure effect on K0 based on Weiss (1974, equation 5) |
---|
1093 | Rgas_atm = 82.05736_wp ! (cm3 * atm) / (mol * K) CODATA (2006) |
---|
1094 | vbarCO2 = 32.3_wp ! partial molal volume (cm3 / mol) from Weiss (1974, Appendix, paragraph 3) |
---|
1095 | K0(i) = K0(i) * exp( ((1-Ptot)*vbarCO2)/(Rgas_atm*tk0) ) ! Weiss (1974, equation 5) |
---|
1096 | |
---|
1097 | ! Pressure effect on all other K's (based on Millero, (1995) |
---|
1098 | ! index: K1(1), K2(2), Kb(3), Kw(4), Ks(5), Kf(6), Kspc(7), Kspa(8), |
---|
1099 | ! K1p(9), K2p(10), K3p(11), Ksi(12) |
---|
1100 | DO ipc = 1, 12 |
---|
1101 | deltav(ipc) = a0(ipc) + a1(ipc) *t + a2(ipc) *t*t |
---|
1102 | deltak(ipc) = (b0(ipc) + b1(ipc) *t + b2(ipc) *t*t) |
---|
1103 | lnkpok0(ipc) = (-(deltav(ipc)) & |
---|
1104 | +(0.5d0*deltak(ipc) * prb) & |
---|
1105 | ) * prb/(R*tk) |
---|
1106 | END DO |
---|
1107 | |
---|
1108 | ! Pressure correction on Ks (Free scale) |
---|
1109 | Ks(i) = Ks_0p*EXP(lnkpok0(5)) |
---|
1110 | ! Conversion factor total -> free scale |
---|
1111 | total2free = 1.d0/(1.d0 + St(i)/Ks(i)) ! Kfree = Ktotal*total2free |
---|
1112 | ! Conversion factor total -> free scale at pressure zero |
---|
1113 | total2free_0p = 1.d0/(1.d0 + St(i)/Ks_0p) ! Kfree = Ktotal*total2free |
---|
1114 | |
---|
1115 | ! Pressure correction on Kf |
---|
1116 | ! Kf must be on FREE scale before correction |
---|
1117 | Kf_0p = Kf_0p * total2free_0p !Convert from Total to Free scale (pressure 0) |
---|
1118 | Kf(i) = Kf_0p * EXP(lnkpok0(6)) !Pressure correction (on Free scale) |
---|
1119 | Kf(i) = Kf(i)/total2free !Convert back from Free to Total scale |
---|
1120 | |
---|
1121 | ! Convert between seawater and total hydrogen (pH) scales |
---|
1122 | free2SWS = 1.d0 + St(i)/Ks(i) + Ft(i)/(Kf(i)*total2free) ! using Kf on free scale |
---|
1123 | total2SWS = total2free * free2SWS ! KSWS = Ktotal*total2SWS |
---|
1124 | SWS2total = 1.d0 / total2SWS |
---|
1125 | ! Conversion at pressure zero |
---|
1126 | free2SWS_0p = 1.d0 + St(i)/Ks_0p + Ft(i)/(Kf_0p) ! using Kf on free scale |
---|
1127 | total2SWS_0p = total2free_0p * free2SWS_0p ! KSWS = Ktotal*total2SWS |
---|
1128 | |
---|
1129 | ! Convert from Total to Seawater scale before pressure correction |
---|
1130 | ! Must change to SEAWATER scale: K1, K2, Kb |
---|
1131 | IF (trim(optK1K2) == 'l') THEN |
---|
1132 | K1(i) = K1(i)*total2SWS_0p |
---|
1133 | K2(i) = K2(i)*total2SWS_0p |
---|
1134 | !This conversion is unnecessary for the K1,K2 from Millero (2010), |
---|
1135 | !since we use here the formulation already on the seawater scale |
---|
1136 | ENDIF |
---|
1137 | Kb(i) = Kb(i)*total2SWS_0p |
---|
1138 | |
---|
1139 | ! Already on SEAWATER scale: K1p, K2p, K3p, Kb, Ksi, Kw |
---|
1140 | |
---|
1141 | ! Other contants (keep on another scale): |
---|
1142 | ! - K0 (independent of pH scale, already pressure corrected) |
---|
1143 | ! - Ks (already on Free scale; already pressure corrected) |
---|
1144 | ! - Kf (already on Total scale; already pressure corrected) |
---|
1145 | ! - Kspc, Kspa (independent of pH scale; pressure-corrected below) |
---|
1146 | |
---|
1147 | ! Perform actual pressure correction (on seawater scale) |
---|
1148 | K1(i) = K1(i)*EXP(lnkpok0(1)) |
---|
1149 | K2(i) = K2(i)*EXP(lnkpok0(2)) |
---|
1150 | Kb(i) = Kb(i)*EXP(lnkpok0(3)) |
---|
1151 | Kw(i) = Kw(i)*EXP(lnkpok0(4)) |
---|
1152 | Kspc(i) = Kspc(i)*EXP(lnkpok0(7)) |
---|
1153 | Kspa(i) = Kspa(i)*EXP(lnkpok0(8)) |
---|
1154 | K1p(i) = K1p(i)*EXP(lnkpok0(9)) |
---|
1155 | K2p(i) = K2p(i)*EXP(lnkpok0(10)) |
---|
1156 | K3p(i) = K3p(i)*EXP(lnkpok0(11)) |
---|
1157 | Ksi(i) = Ksi(i)*EXP(lnkpok0(12)) |
---|
1158 | |
---|
1159 | ! Convert back to original total scale: |
---|
1160 | K1(i) = K1(i) *SWS2total |
---|
1161 | K2(i) = K2(i) *SWS2total |
---|
1162 | K1p(i) = K1p(i)*SWS2total |
---|
1163 | K2p(i) = K2p(i)*SWS2total |
---|
1164 | K3p(i) = K3p(i)*SWS2total |
---|
1165 | Kb(i) = Kb(i) *SWS2total |
---|
1166 | Ksi(i) = Ksi(i)*SWS2total |
---|
1167 | Kw(i) = Kw(i) *SWS2total |
---|
1168 | |
---|
1169 | ELSE |
---|
1170 | |
---|
1171 | K0(i) = 1.e20_wp |
---|
1172 | K1(i) = 1.e20_wp |
---|
1173 | K2(i) = 1.e20_wp |
---|
1174 | Kb(i) = 1.e20_wp |
---|
1175 | Kw(i) = 1.e20_wp |
---|
1176 | Ks(i) = 1.e20_wp |
---|
1177 | Kf(i) = 1.e20_wp |
---|
1178 | Kspc(i) = 1.e20_wp |
---|
1179 | Kspa(i) = 1.e20_wp |
---|
1180 | K1p(i) = 1.e20_wp |
---|
1181 | K2p(i) = 1.e20_wp |
---|
1182 | K3p(i) = 1.e20_wp |
---|
1183 | Ksi(i) = 1.e20_wp |
---|
1184 | Bt(i) = 1.e20_wp |
---|
1185 | Ft(i) = 1.e20_wp |
---|
1186 | St(i) = 1.e20_wp |
---|
1187 | |
---|
1188 | ENDIF |
---|
1189 | |
---|
1190 | END DO |
---|
1191 | |
---|
1192 | RETURN |
---|
1193 | END SUBROUTINE constants |
---|
1194 | |
---|
1195 | ! ---------------------------------------------------------------------- |
---|
1196 | ! VARSOLVER |
---|
1197 | ! ---------------------------------------------------------------------- |
---|
1198 | ! |
---|
1199 | !> \file varsolver.f90 |
---|
1200 | !! \BRIEF |
---|
1201 | !> Module with varsolver subroutine - solve for pH and other carbonate system variables |
---|
1202 | !> Solve for pH and other carbonate system variables (with input from vars routine) |
---|
1203 | SUBROUTINE varsolver(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, & |
---|
1204 | temp, salt, ta, tc, pt, sit, & |
---|
1205 | Bt, St, Ft, & |
---|
1206 | K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi, & |
---|
1207 | Patm, Phydro_bar, rhodum, optGAS ) |
---|
1208 | |
---|
1209 | ! Purpose: Solve for pH and other carbonate system variables (with input from vars routine) |
---|
1210 | |
---|
1211 | ! INPUT variables: |
---|
1212 | ! ================ |
---|
1213 | ! temp = in situ temperature [degrees C] |
---|
1214 | ! ta = total alkalinity in [eq/m^3] or [eq/kg] based on optCON in calling routine (vars) |
---|
1215 | ! tc = dissolved inorganic carbon in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars) |
---|
1216 | ! pt = total dissolved inorganic phosphorus in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars) |
---|
1217 | ! sit = total dissolved inorganic silicon in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars) |
---|
1218 | ! Bt = total dissolved inorganic boron computed in calling routine (vars) |
---|
1219 | ! St = total dissolved inorganic sulfur computed in calling routine (vars) |
---|
1220 | ! Ft = total dissolved inorganic fluorine computed in calling routine (vars) |
---|
1221 | ! K's = K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi |
---|
1222 | ! Patm = atmospheric pressure [atm] |
---|
1223 | ! Phydro_bar = hydrostatic pressure [bar] |
---|
1224 | ! rhodum = density factor as computed in calling routine (vars) |
---|
1225 | ! ----------- |
---|
1226 | ! optGAS: choose in situ vs. potential fCO2 and pCO2 (default optGAS = 'Pinsitu') |
---|
1227 | ! --------- |
---|
1228 | ! PRESSURE & T corrections for K0 and the fugacity coefficient (Cf) |
---|
1229 | ! -> 'Pzero' = 'zero order' fCO2 and pCO2 (typical approach, which is flawed) |
---|
1230 | ! considers in situ T & only atm pressure (hydrostatic=0) |
---|
1231 | ! -> 'Ppot' = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) |
---|
1232 | ! considers potential T & only atm pressure (hydrostatic press = 0) |
---|
1233 | ! -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure) |
---|
1234 | ! considers in situ T & total pressure (atm + hydrostatic) |
---|
1235 | ! --------- |
---|
1236 | |
---|
1237 | ! OUTPUT variables: |
---|
1238 | ! ================= |
---|
1239 | ! ph = pH on total scale |
---|
1240 | ! pco2 = CO2 partial pressure (uatm) |
---|
1241 | ! fco2 = CO2 fugacity (uatm) |
---|
1242 | ! co2 = aqueous CO2 concentration in [mol/kg] or [mol/m^3] determined by rhodum (depends on optCON in calling routine) |
---|
1243 | ! hco3 = bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] determined by rhodum |
---|
1244 | ! co3 = carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] determined by rhodum |
---|
1245 | ! OmegaA = Omega for aragonite, i.e., the aragonite saturation state |
---|
1246 | ! OmegaC = Omega for calcite, i.e., the calcite saturation state |
---|
1247 | |
---|
1248 | USE mocsy_singledouble |
---|
1249 | USE mocsy_phsolvers |
---|
1250 | |
---|
1251 | IMPLICIT NONE |
---|
1252 | |
---|
1253 | ! Input variables |
---|
1254 | !> <b>in situ temperature</b> [degrees C] |
---|
1255 | REAL(kind=wp), INTENT(in) :: temp |
---|
1256 | !> <b>salinity</b> [on the practical salinity scale, dimensionless] |
---|
1257 | REAL(kind=wp), INTENT(in) :: salt |
---|
1258 | !> total alkalinity in <b>[eq/m^3]</b> OR in <b>[eq/kg]</b>, depending on optCON in calling routine |
---|
1259 | REAL(kind=wp), INTENT(in) :: ta |
---|
1260 | !> dissolved inorganic carbon in <b>[mol/m^3]</b> OR in <b>[mol/kg]</b>, depending on optCON in calling routine |
---|
1261 | REAL(kind=wp), INTENT(in) :: tc |
---|
1262 | !> phosphate concentration in <b>[mol/m^3]</b> OR in <b>[mol/kg]</b>, depending on optCON in calling routine |
---|
1263 | REAL(kind=wp), INTENT(in) :: pt |
---|
1264 | !> total dissolved inorganic silicon concentration in <b>[mol/m^3]</b> OR in <b>[mol/kg]</b>, depending on optCON in calling routine |
---|
1265 | REAL(kind=wp), INTENT(in) :: sit |
---|
1266 | !> total boron from either Uppstrom (1974) or Lee et al. (2010), depending on optB in calling routine |
---|
1267 | REAL(kind=wp), INTENT(in) :: Bt |
---|
1268 | !> total sulfate (Morris & Riley, 1966) |
---|
1269 | REAL(kind=wp), INTENT(in) :: St |
---|
1270 | !> total fluoride (Riley, 1965) |
---|
1271 | REAL(kind=wp), INTENT(in) :: Ft |
---|
1272 | !> solubility of CO2 in seawater (Weiss, 1974), also known as K0 |
---|
1273 | REAL(kind=wp), INTENT(in) :: K0 |
---|
1274 | !> K1 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 |
---|
1275 | REAL(kind=wp), INTENT(in) :: K1 |
---|
1276 | !> K2 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 |
---|
1277 | REAL(kind=wp), INTENT(in) :: K2 |
---|
1278 | !> equilibrium constant for dissociation of boric acid |
---|
1279 | REAL(kind=wp), INTENT(in) :: Kb |
---|
1280 | !> equilibrium constant for the dissociation of water (Millero, 1995) |
---|
1281 | REAL(kind=wp), INTENT(in) :: Kw |
---|
1282 | !> equilibrium constant for the dissociation of bisulfate (Dickson, 1990) |
---|
1283 | REAL(kind=wp), INTENT(in) :: Ks |
---|
1284 | !> equilibrium constant for the dissociation of hydrogen fluoride |
---|
1285 | !! from Dickson and Riley (1979) or Perez and Fraga (1987), depending on optKf |
---|
1286 | REAL(kind=wp), INTENT(in) :: Kf |
---|
1287 | !> solubility product for calcite (Mucci, 1983) |
---|
1288 | REAL(kind=wp), INTENT(in) :: Kspc |
---|
1289 | !> solubility product for aragonite (Mucci, 1983) |
---|
1290 | REAL(kind=wp), INTENT(in) :: Kspa |
---|
1291 | !> 1st dissociation constant for phosphoric acid (Millero, 1995) |
---|
1292 | REAL(kind=wp), INTENT(in) :: K1p |
---|
1293 | !> 2nd dissociation constant for phosphoric acid (Millero, 1995) |
---|
1294 | REAL(kind=wp), INTENT(in) :: K2p |
---|
1295 | !> 3rd dissociation constant for phosphoric acid (Millero, 1995) |
---|
1296 | REAL(kind=wp), INTENT(in) :: K3p |
---|
1297 | !> equilibrium constant for the dissociation of silicic acid (Millero, 1995) |
---|
1298 | REAL(kind=wp), INTENT(in) :: Ksi |
---|
1299 | !> total atmospheric pressure <b>[atm]</b> |
---|
1300 | REAL(kind=wp), INTENT(in) :: Patm |
---|
1301 | !> total hydrostatic pressure <b>[bar]</b> |
---|
1302 | REAL(kind=wp), INTENT(in) :: Phydro_bar |
---|
1303 | !> density factor as computed incalling routine (vars) |
---|
1304 | REAL(kind=wp), INTENT(in) :: rhodum |
---|
1305 | !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction) |
---|
1306 | !! 'Ppot' - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) |
---|
1307 | !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf |
---|
1308 | !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean |
---|
1309 | !f2py character*7 optional, intent(in) :: optGAS='Pinsitu' |
---|
1310 | CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS |
---|
1311 | |
---|
1312 | ! Output variables: |
---|
1313 | !> pH on the <b>total scale</b> |
---|
1314 | REAL(kind=wp), INTENT(out) :: ph |
---|
1315 | !> CO2 partial pressure <b>[uatm]</b> |
---|
1316 | REAL(kind=wp), INTENT(out) :: pco2 |
---|
1317 | !> CO2 fugacity <b>[uatm]</b> |
---|
1318 | REAL(kind=wp), INTENT(out) :: fco2 |
---|
1319 | !> aqueous CO2* concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg</b>] depending on choice for optCON |
---|
1320 | REAL(kind=wp), INTENT(out) :: co2 |
---|
1321 | !> bicarbonate ion (HCO3-) concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg]</b> depending on choice for optCON |
---|
1322 | REAL(kind=wp), INTENT(out) :: hco3 |
---|
1323 | !> carbonate ion (CO3--) concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg]</b> depending on choice for optCON |
---|
1324 | REAL(kind=wp), INTENT(out) :: co3 |
---|
1325 | !> Omega for aragonite, i.e., the aragonite saturation state |
---|
1326 | REAL(kind=wp), INTENT(out) :: OmegaA |
---|
1327 | !> Omega for calcite, i.e., the calcite saturation state |
---|
1328 | REAL(kind=wp), INTENT(out) :: OmegaC |
---|
1329 | |
---|
1330 | ! Local variables |
---|
1331 | REAL(kind=wp) :: Phydro_atm, Ptot |
---|
1332 | REAL(kind=wp) :: Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff |
---|
1333 | REAL(kind=wp) :: tk, tk0 |
---|
1334 | real(kind=wp) :: temp68, tempot, tempot68 |
---|
1335 | REAL(kind=wp) :: Hinit, H |
---|
1336 | REAL(kind=wp) :: HSO4, HF, HSI, HPO4 |
---|
1337 | REAL(kind=wp) :: ab, aw, ac |
---|
1338 | REAL(kind=wp) :: cu, cb, cc |
---|
1339 | REAL(kind=wp) :: Ca |
---|
1340 | ! Array to pass optional arguments |
---|
1341 | CHARACTER(7) :: opGAS |
---|
1342 | |
---|
1343 | IF (PRESENT(optGAS)) THEN |
---|
1344 | opGAS = optGAS |
---|
1345 | ELSE |
---|
1346 | opGAS = 'Pinsitu' |
---|
1347 | ENDIF |
---|
1348 | |
---|
1349 | ! Compute pH from constants and total concentrations |
---|
1350 | ! - use SolveSAPHE v1.0.1 routines from Munhoven (2013, GMD) modified to use mocsy's Ks instead of its own |
---|
1351 | ! 1) Compute best starting point for H+ calculation |
---|
1352 | call ahini_for_at(ta, tc, Bt, K1, K2, Kb, Hinit) |
---|
1353 | ! 2) Solve for H+ using above result as the initial H+ value |
---|
1354 | H = solve_at_general(ta, tc, Bt, & |
---|
1355 | pt, sit, & |
---|
1356 | St, Ft, & |
---|
1357 | K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & |
---|
1358 | Hinit) |
---|
1359 | ! 3) Calculate pH from from H+ concentration (mol/kg) |
---|
1360 | IF (H > 0.d0) THEN |
---|
1361 | pH = -1.*LOG10(H) |
---|
1362 | ELSE |
---|
1363 | pH = 1.e20_wp |
---|
1364 | ENDIF |
---|
1365 | |
---|
1366 | ! Compute carbonate Alk (Ac) by difference: from total Alk and other Alk components |
---|
1367 | HSO4 = St/(1.0d0 + Ks/(H/(1.0d0 + St/Ks))) |
---|
1368 | HF = 1.0d0/(1.0d0 + Kf/H) |
---|
1369 | HSI = 1.0d0/(1.0d0 + H/Ksi) |
---|
1370 | HPO4 = (K1p*K2p*(H + 2.*K3p) - H**3) / & |
---|
1371 | (H**3 + K1p*H**2 + K1p*K2p*H + K1p*K2p*K3p) |
---|
1372 | ab = Bt/(1.0d0 + H/Kb) |
---|
1373 | aw = Kw/H - H/(1.0d0 + St/Ks) |
---|
1374 | ac = ta + hso4 - sit*hsi - ab - aw + Ft*hf - pt*hpo4 |
---|
1375 | |
---|
1376 | ! Calculate CO2*, HCO3-, & CO32- (in mol/kg soln) from Ct, Ac, H+, K1, & K2 |
---|
1377 | cu = (2.0d0 * tc - ac) / (2.0d0 + K1 / H) |
---|
1378 | cb = K1 * cu / H |
---|
1379 | cc = K2 * cb / H |
---|
1380 | |
---|
1381 | ! When optCON = 'mol/m3' in calling routine (vars), then: |
---|
1382 | ! convert output var concentrations from mol/kg to mol/m^3 |
---|
1383 | ! e.g., for case when drho = 1028, multiply by [1.028 kg/L x 1000 L/m^3]) |
---|
1384 | co2 = cu * rhodum |
---|
1385 | hco3 = cb * rhodum |
---|
1386 | co3 = cc * rhodum |
---|
1387 | |
---|
1388 | ! Determine CO2 fugacity [uatm] |
---|
1389 | ! NOTE: equation just below requires CO2* in mol/kg |
---|
1390 | fCO2 = cu * 1.e6_wp/K0 |
---|
1391 | |
---|
1392 | ! Determine CO2 partial pressure from CO2 fugacity [uatm] |
---|
1393 | tk = 273.15d0 + temp |
---|
1394 | !Compute EITHER "potential pCO2" OR "in situ pCO2" (T and P used for calculations will differ) |
---|
1395 | IF (trim(opGAS) == 'Pzero' .OR. trim(opGAS) == 'pzero') THEN |
---|
1396 | tk0 = tk !in situ temperature (K) for K0 calculation |
---|
1397 | Ptot = Patm !total pressure (in atm) = atmospheric pressure ONLY |
---|
1398 | ELSEIF (trim(opGAS) == 'Ppot' .OR. trim(opGAS) == 'ppot') THEN |
---|
1399 | !Use potential temperature and atmospheric pressure (water parcel adiabatically brought back to surface) |
---|
1400 | !temp68 = (temp - 0.0002d0) / 0.99975d0 !temp = in situ T; temp68 is same converted to ITPS-68 scale |
---|
1401 | !tempot68 = sw_ptmp(salt, temp68, Phydro_bar*10d0, 0.0d0) !potential temperature (C) |
---|
1402 | !tempot = 0.99975*tempot68 + 0.0002 |
---|
1403 | !tk0 = tempot + 273.15d0 !potential temperature (K) for fugacity coeff. calc as needed for potential fCO2 & pCO2 |
---|
1404 | tempot = sw_ptmp(salt, temp, Phydro_bar*10._wp, 0.0_wp) !potential temperature (C) |
---|
1405 | tk0 = tempot + 273.15d0 !potential temperature (K) for fugacity coeff. calc as needed for potential fCO2 & pCO2 |
---|
1406 | Ptot = Patm !total pressure (in atm) = atmospheric pressure ONLY |
---|
1407 | ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN |
---|
1408 | !Use in situ temperature and total pressure |
---|
1409 | tk0 = tk !in situ temperature (K) for fugacity coefficient calculation |
---|
1410 | Phydro_atm = Phydro_bar / 1.01325d0 !convert hydrostatic pressure from bar to atm (1.01325 bar / atm) |
---|
1411 | Ptot = Patm + Phydro_atm !total pressure (in atm) = atmospheric pressure + hydrostatic pressure |
---|
1412 | ELSE |
---|
1413 | PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'" |
---|
1414 | STOP |
---|
1415 | ENDIF |
---|
1416 | |
---|
1417 | ! Now that we have T and P in the right form, continue with calculation of fugacity coefficient (and pCO2) |
---|
1418 | Rgas_atm = 82.05736_wp ! (cm3 * atm) / (mol * K) CODATA (2006) |
---|
1419 | ! To compute fugcoeff, we need 3 other terms (B, Del, xc2) in addition to 3 others above (tk, Ptot, Rgas_atm) |
---|
1420 | B = -1636.75d0 + 12.0408d0*tk0 - 0.0327957d0*(tk0*tk0) + 0.0000316528d0*(tk0*tk0*tk0) |
---|
1421 | Del = 57.7d0 - 0.118d0*tk0 |
---|
1422 | ! "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9 |
---|
1423 | ! x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite) |
---|
1424 | ! Let's assume that xCO2 = fCO2. Resulting fugcoeff is identical to 8th digit after the decimal. |
---|
1425 | xCO2approx = fCO2 * 1.e-6_wp |
---|
1426 | IF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN |
---|
1427 | ! xCO2approx = 400.0e-6_wp !a simple test (gives about same result as seacarb for pCO2insitu) |
---|
1428 | ! approximate surface xCO2 ~ surface fCO2 (i.e., in situ fCO2 d by exponential pressure correction) |
---|
1429 | xCO2approx = xCO2approx * exp( ((1-Ptot)*32.3_wp)/(82.05736_wp*tk0) ) ! of K0 press. correction, see Weiss (1974, equation 5) |
---|
1430 | ENDIF |
---|
1431 | xc2 = (1.0d0 - xCO2approx)**2 |
---|
1432 | fugcoeff = exp( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk0) ) |
---|
1433 | pCO2 = fCO2 / fugcoeff |
---|
1434 | |
---|
1435 | ! Determine Omega Calcite et Aragonite |
---|
1436 | ! OmegaA = ((0.01028d0*salt/35.0d0)*cc) / Kspa |
---|
1437 | ! OmegaC = ((0.01028d0*salt/35.0d0)*cc) / Kspc |
---|
1438 | ! - see comments from Munhoven on the best value "0.02128" which differs slightly from the best practices guide (0.02127) |
---|
1439 | Ca = (0.02128d0/40.078d0) * salt/1.80655d0 |
---|
1440 | OmegaA = (Ca*cc) / Kspa |
---|
1441 | OmegaC = (Ca*cc) / Kspc |
---|
1442 | |
---|
1443 | RETURN |
---|
1444 | END SUBROUTINE varsolver |
---|
1445 | |
---|
1446 | ! ---------------------------------------------------------------------- |
---|
1447 | ! VARS |
---|
1448 | ! ---------------------------------------------------------------------- |
---|
1449 | ! |
---|
1450 | !> \file vars.f90 |
---|
1451 | !! \BRIEF |
---|
1452 | !> Module with vars subroutine - compute carbonate system vars from DIC,Alk,T,S,P,nuts |
---|
1453 | !> Computes standard carbonate system variables (pH, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R) |
---|
1454 | !! as 1D arrays FROM |
---|
1455 | !! temperature, salinity, pressure, |
---|
1456 | !! total alkalinity (ALK), dissolved inorganic carbon (DIC), |
---|
1457 | !! silica and phosphate concentrations (all 1-D arrays) |
---|
1458 | SUBROUTINE vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis, & |
---|
1459 | temp, sal, alk, dic, sil, phos, Patm, depth, lat, N, & |
---|
1460 | optCON, optT, optP, optB, optK1K2, optKf, optGAS ) |
---|
1461 | |
---|
1462 | ! Purpose: |
---|
1463 | ! Computes other standard carbonate system variables (pH, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R) |
---|
1464 | ! as 1D arrays |
---|
1465 | ! FROM: |
---|
1466 | ! temperature, salinity, pressure, |
---|
1467 | ! total alkalinity (ALK), dissolved inorganic carbon (DIC), |
---|
1468 | ! silica and phosphate concentrations (all 1-D arrays) |
---|
1469 | |
---|
1470 | ! INPUT variables: |
---|
1471 | ! ================ |
---|
1472 | ! Patm = atmospheric pressure [atm] |
---|
1473 | ! depth = depth [m] (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure) |
---|
1474 | ! = pressure [db] (with optP='db') |
---|
1475 | ! lat = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m') |
---|
1476 | ! = dummy array (unused when optP='db') |
---|
1477 | ! temp = potential temperature [degrees C] (with optT='Tpot', i.e., models carry tempot, not in situ temp) |
---|
1478 | ! = in situ temperature [degrees C] (with optT='Tinsitu', e.g., for data) |
---|
1479 | ! sal = salinity in [psu] |
---|
1480 | ! alk = total alkalinity in [eq/m^3] with optCON = 'mol/m3' |
---|
1481 | ! = [eq/kg] with optCON = 'mol/kg' |
---|
1482 | ! dic = dissolved inorganic carbon [mol/m^3] with optCON = 'mol/m3' |
---|
1483 | ! = [mol/kg] with optCON = 'mol/kg' |
---|
1484 | ! sil = silica [mol/m^3] with optCON = 'mol/m3' |
---|
1485 | ! = [mol/kg] with optCON = 'mol/kg' |
---|
1486 | ! phos = phosphate [mol/m^3] with optCON = 'mol/m3' |
---|
1487 | ! = [mol/kg] with optCON = 'mol/kg' |
---|
1488 | ! INPUT options: |
---|
1489 | ! ============== |
---|
1490 | ! ----------- |
---|
1491 | ! optCON: choose input & output concentration units - mol/kg (data) vs. mol/m^3 (models) |
---|
1492 | ! ----------- |
---|
1493 | ! -> 'mol/kg' for DIC, ALK, sil, & phos given on mokal scale, i.e., in mol/kg (std DATA units) |
---|
1494 | ! -> 'mol/m3' for DIC, ALK, sil, & phos given in mol/m^3 (std MODEL units) |
---|
1495 | ! ----------- |
---|
1496 | ! optT: choose in situ vs. potential temperature as input |
---|
1497 | ! --------- |
---|
1498 | ! NOTE: Carbonate chem calculations require IN-SITU temperature (not potential Temperature) |
---|
1499 | ! -> 'Tpot' means input is pot. Temperature (in situ Temp "tempis" is computed) |
---|
1500 | ! -> 'Tinsitu' means input is already in-situ Temperature, not pot. Temp ("tempis" not computed) |
---|
1501 | ! --------- |
---|
1502 | ! optP: choose depth (m) vs pressure (db) as input |
---|
1503 | ! --------- |
---|
1504 | ! -> 'm' means "depth" input is in "m" (thus in situ Pressure "p" [db] is computed) |
---|
1505 | ! -> 'db' means "depth" input is already in situ pressure [db], not m (thus p = depth) |
---|
1506 | ! --------- |
---|
1507 | ! optB: choose total boron formulation - Uppström (1974) vs. Lee et al. (2010) |
---|
1508 | ! --------- |
---|
1509 | ! -> 'u74' means use classic formulation of Uppström (1974) for total Boron |
---|
1510 | ! -> 'l10' means use newer formulation of Lee et al. (2010) for total Boron |
---|
1511 | ! --------- |
---|
1512 | ! optK1K2: |
---|
1513 | ! --------- |
---|
1514 | ! -> 'l' means use Lueker et al. (2000) formulations for K1 & K2 (recommended by Dickson et al. 2007) |
---|
1515 | ! **** BUT this should only be used when 2 < T < 35 and 19 < S < 43 |
---|
1516 | ! -> 'm10' means use Millero (2010) formulation for K1 & K2 (see Dickson et al., 2007) |
---|
1517 | ! **** Valid for 0 < T < 50°C and 1 < S < 50 psu |
---|
1518 | ! ---------- |
---|
1519 | ! optKf: |
---|
1520 | ! ---------- |
---|
1521 | ! -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007) |
---|
1522 | ! **** BUT Valid for 9 < T < 33°C and 10 < S < 40. |
---|
1523 | ! -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994) |
---|
1524 | ! ----------- |
---|
1525 | ! optGAS: choose in situ vs. potential fCO2 and pCO2 |
---|
1526 | ! --------- |
---|
1527 | ! PRESSURE corrections for K0 and the fugacity coefficient (Cf) |
---|
1528 | ! -> 'Pzero' = 'zero order' fCO2 and pCO2 (typical approach, which is flawed) |
---|
1529 | ! considers in situ T & only atm pressure (hydrostatic=0) |
---|
1530 | ! -> 'Ppot' = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) |
---|
1531 | ! considers potential T & only atm pressure (hydrostatic press = 0) |
---|
1532 | ! -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure) |
---|
1533 | ! considers in situ T & total pressure (atm + hydrostatic) |
---|
1534 | ! --------- |
---|
1535 | |
---|
1536 | ! OUTPUT variables: |
---|
1537 | ! ================= |
---|
1538 | ! ph = pH on total scale |
---|
1539 | ! pco2 = CO2 partial pressure (uatm) |
---|
1540 | ! fco2 = CO2 fugacity (uatm) |
---|
1541 | ! co2 = aqueous CO2 concentration in [mol/kg] or [mol/m^3] depending on optCON |
---|
1542 | ! hco3 = bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] depending on optCON |
---|
1543 | ! co3 = carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] depending on optCON |
---|
1544 | ! OmegaA = Omega for aragonite, i.e., the aragonite saturation state |
---|
1545 | ! OmegaC = Omega for calcite, i.e., the calcite saturation state |
---|
1546 | ! BetaD = Revelle factor dpCO2/pCO2 / dDIC/DIC |
---|
1547 | ! rhoSW = in-situ density of seawater; rhoSW = f(s, t, p) |
---|
1548 | ! p = pressure [decibars]; p = f(depth, latitude) if computed from depth [m] OR p = depth if [db] |
---|
1549 | ! tempis = in-situ temperature [degrees C] |
---|
1550 | |
---|
1551 | USE mocsy_singledouble |
---|
1552 | |
---|
1553 | IMPLICIT NONE |
---|
1554 | |
---|
1555 | ! Input variables |
---|
1556 | !> number of records |
---|
1557 | INTEGER, INTENT(in) :: N |
---|
1558 | !> either <b>in situ temperature</b> (when optT='Tinsitu', typical data) |
---|
1559 | !! OR <b>potential temperature</b> (when optT='Tpot', typical models) <b>[degree C]</b> |
---|
1560 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp |
---|
1561 | !> salinity <b>[psu]</b> |
---|
1562 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: sal |
---|
1563 | !> total alkalinity in <b>[eq/m^3]</b> (when optCON = 'mol/m3') OR in <b>[eq/kg]</b> (when optCON = 'mol/kg') |
---|
1564 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: alk |
---|
1565 | !> dissolved inorganic carbon in <b>[mol/m^3]</b> (when optCON = 'mol/m3') OR in <b>[mol/kg]</b> (when optCON = 'mol/kg') |
---|
1566 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: dic |
---|
1567 | !> SiO2 concentration in <b>[mol/m^3]</b> (when optCON = 'mol/m3') OR in <b>[mol/kg]</b> (when optCON = 'mol/kg') |
---|
1568 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: sil |
---|
1569 | !> phosphate concentration in <b>[mol/m^3]</b> (when optCON = 'mol/m3') OR in <b>[mol/kg]</b> (when optCON = 'mol/kg') |
---|
1570 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: phos |
---|
1571 | !f2py optional , depend(sal) :: n=len(sal) |
---|
1572 | !> atmospheric pressure <b>[atm]</b> |
---|
1573 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm |
---|
1574 | !> depth in \b meters (when optP='m') or \b decibars (when optP='db') |
---|
1575 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: depth |
---|
1576 | !> latitude <b>[degrees north]</b> |
---|
1577 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: lat |
---|
1578 | |
---|
1579 | !> choose either \b 'mol/kg' (std DATA units) or \b 'mol/m3' (std MODEL units) to select |
---|
1580 | !! concentration units for input (for alk, dic, sil, phos) & output (co2, hco3, co3) |
---|
1581 | CHARACTER(6), INTENT(in) :: optCON |
---|
1582 | !> choose \b 'Tinsitu' for in situ temperature or \b 'Tpot' for potential temperature (in situ Temp is computed, needed for models) |
---|
1583 | CHARACTER(7), INTENT(in) :: optT |
---|
1584 | !> for depth input, choose \b "db" for decibars (in situ pressure) or \b "m" for meters (pressure is computed, needed for models) |
---|
1585 | CHARACTER(2), INTENT(in) :: optP |
---|
1586 | !> for total boron, choose either \b 'u74' (Uppstrom, 1974) or \b 'l10' (Lee et al., 2010). |
---|
1587 | !! The 'l10' formulation is based on 139 measurements (instead of 20), |
---|
1588 | !! uses a more accurate method, and |
---|
1589 | !! generally increases total boron in seawater by 4% |
---|
1590 | !f2py character*3 optional, intent(in) :: optB='l10' |
---|
1591 | CHARACTER(3), OPTIONAL, INTENT(in) :: optB |
---|
1592 | !> for Kf, choose either \b 'pf' (Perez & Fraga, 1987) or \b 'dg' (Dickson & Riley, 1979) |
---|
1593 | !f2py character*2 optional, intent(in) :: optKf='pf' |
---|
1594 | CHARACTER(2), OPTIONAL, INTENT(in) :: optKf |
---|
1595 | !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010) |
---|
1596 | !f2py character*3 optional, intent(in) :: optK1K2='l' |
---|
1597 | CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2 |
---|
1598 | !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction) |
---|
1599 | !! 'Ppot' - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) |
---|
1600 | !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf |
---|
1601 | !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean |
---|
1602 | !f2py character*7 optional, intent(in) :: optGAS='Pinsitu' |
---|
1603 | CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS |
---|
1604 | |
---|
1605 | ! Output variables: |
---|
1606 | !> pH on the <b>total scale</b> |
---|
1607 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: ph |
---|
1608 | !> CO2 partial pressure <b>[uatm]</b> |
---|
1609 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: pco2 |
---|
1610 | !> CO2 fugacity <b>[uatm]</b> |
---|
1611 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: fco2 |
---|
1612 | !> aqueous CO2* concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg</b>] depending on choice for optCON |
---|
1613 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: co2 |
---|
1614 | !> bicarbonate ion (HCO3-) concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg]</b> depending on choice for optCON |
---|
1615 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: hco3 |
---|
1616 | !> carbonate ion (CO3--) concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg]</b> depending on choice for optCON |
---|
1617 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: co3 |
---|
1618 | !> Omega for aragonite, i.e., the aragonite saturation state |
---|
1619 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: OmegaA |
---|
1620 | !> Omega for calcite, i.e., the calcite saturation state |
---|
1621 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: OmegaC |
---|
1622 | !> Revelle factor, i.e., dpCO2/pCO2 / dDIC/DIC |
---|
1623 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: BetaD |
---|
1624 | !> in-situ density of seawater; rhoSW = f(s, t, p) in <b>[kg/m3]</b> |
---|
1625 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: rhoSW |
---|
1626 | !> pressure <b>[decibars]</b>; p = f(depth, latitude) if computed from depth [m] (when optP='m') OR p = depth [db] (when optP='db') |
---|
1627 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: p |
---|
1628 | !> in-situ temperature \b <b>[degrees C]</b> |
---|
1629 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: tempis |
---|
1630 | |
---|
1631 | ! Local variables |
---|
1632 | REAL(kind=wp) :: ssal, salk, sdic, ssil, sphos |
---|
1633 | REAL(kind=wp) :: tempot, tempis68, tempot68 |
---|
1634 | REAL(kind=wp) :: drho |
---|
1635 | |
---|
1636 | REAL(kind=wp) :: K0, K1, K2, Kb, Kw, Ks, Kf, Kspc |
---|
1637 | REAL(kind=wp) :: Kspa, K1p, K2p, K3p, Ksi |
---|
1638 | REAL(kind=wp) :: St, Ft, Bt |
---|
1639 | |
---|
1640 | REAL(kind=wp), DIMENSION(1) :: aK0, aK1, aK2, aKb, aKw, aKs, aKf, aKspc |
---|
1641 | REAL(kind=wp), DIMENSION(1) :: aKspa, aK1p, aK2p, aK3p, aKsi |
---|
1642 | REAL(kind=wp), DIMENSION(1) :: aSt, aFt, aBt |
---|
1643 | |
---|
1644 | REAL(kind=wp) :: Patmd, Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff |
---|
1645 | REAL(kind=wp) :: Phydro_atm |
---|
1646 | |
---|
1647 | INTEGER :: i, icount |
---|
1648 | |
---|
1649 | REAL(kind=wp) :: t, tk, prb |
---|
1650 | REAL(kind=wp) :: s |
---|
1651 | REAL(kind=wp) :: tc, ta |
---|
1652 | REAL(kind=wp) :: sit, pt |
---|
1653 | REAL(kind=wp) :: Hinit |
---|
1654 | REAL(kind=wp) :: ah1 |
---|
1655 | |
---|
1656 | REAL(kind=wp) :: HSO4, HF, HSI, HPO4 |
---|
1657 | REAL(kind=wp) :: ab, aw, ac, ah2, erel |
---|
1658 | |
---|
1659 | REAL(kind=wp) :: cu, cb, cc |
---|
1660 | |
---|
1661 | REAL(kind=wp), DIMENSION(2) :: dicdel, pco2del |
---|
1662 | REAL(kind=wp) :: dx, Rf |
---|
1663 | REAL(kind=wp) :: dph, dpco2, dfco2, dco2, dhco3, dco3, dOmegaA, dOmegaC |
---|
1664 | |
---|
1665 | INTEGER :: kcomp |
---|
1666 | INTEGER :: j, minusplus |
---|
1667 | |
---|
1668 | ! Arrays to pass optional arguments into or use defaults (Dickson et al., 2007) |
---|
1669 | CHARACTER(3) :: opB |
---|
1670 | CHARACTER(2) :: opKf |
---|
1671 | CHARACTER(3) :: opK1K2 |
---|
1672 | CHARACTER(7) :: opGAS |
---|
1673 | |
---|
1674 | ! Set defaults for optional arguments (in Fortran 90) |
---|
1675 | ! Note: Optional arguments with f2py (python) are set above with |
---|
1676 | ! the !f2py statements that precede each type declaraion |
---|
1677 | IF (PRESENT(optB)) THEN |
---|
1678 | ! print *,"optB present:" |
---|
1679 | ! print *,"optB = ", optB |
---|
1680 | opB = optB |
---|
1681 | ELSE |
---|
1682 | ! Default is Lee et al (2010) for total boron |
---|
1683 | ! print *,"optB NOT present:" |
---|
1684 | opB = 'l10' |
---|
1685 | ! print *,"opB = ", opB |
---|
1686 | ENDIF |
---|
1687 | IF (PRESENT(optKf)) THEN |
---|
1688 | ! print *,"optKf = ", optKf |
---|
1689 | opKf = optKf |
---|
1690 | ELSE |
---|
1691 | ! print *,"optKf NOT present:" |
---|
1692 | ! Default is Perez & Fraga (1987) for Kf |
---|
1693 | opKf = 'pf' |
---|
1694 | ! print *,"opKf = ", opKf |
---|
1695 | ENDIF |
---|
1696 | IF (PRESENT(optK1K2)) THEN |
---|
1697 | ! print *,"optK1K2 = ", optK1K2 |
---|
1698 | opK1K2 = optK1K2 |
---|
1699 | ELSE |
---|
1700 | ! print *,"optK1K2 NOT present:" |
---|
1701 | ! Default is Lueker et al. 2000) for K1 & K2 |
---|
1702 | opK1K2 = 'l' |
---|
1703 | ! print *,"opK1K2 = ", opK1K2 |
---|
1704 | ENDIF |
---|
1705 | IF (PRESENT(optGAS)) THEN |
---|
1706 | opGAS = optGAS |
---|
1707 | ELSE |
---|
1708 | opGAS = 'Pinsitu' |
---|
1709 | ENDIF |
---|
1710 | |
---|
1711 | icount = 0 |
---|
1712 | DO i = 1, N |
---|
1713 | icount = icount + 1 |
---|
1714 | ! =============================================================== |
---|
1715 | ! Convert model depth -> press; convert model Theta -> T in situ |
---|
1716 | ! =============================================================== |
---|
1717 | ! * Model temperature tracer is usually "potential temperature" |
---|
1718 | ! * Model vertical grid is usually in meters |
---|
1719 | ! BUT carbonate chem routines require pressure & in-situ T |
---|
1720 | ! Thus before computing chemistry, if appropriate, |
---|
1721 | ! convert these 2 model vars (input to this routine) |
---|
1722 | ! - depth [m] => convert to pressure [db] |
---|
1723 | ! - potential temperature (C) => convert to in-situ T (C) |
---|
1724 | ! ------------------------------------------------------- |
---|
1725 | ! 1) Compute pressure [db] from depth [m] and latitude [degrees] (if input is m, for models) |
---|
1726 | !print *,"optP =", optP, "end" |
---|
1727 | IF (trim(optP) == 'm' ) THEN |
---|
1728 | ! Compute pressure [db] from depth [m] and latitude [degrees] |
---|
1729 | p(i) = p80(depth(i), lat(i)) |
---|
1730 | ELSEIF (trim(optP) == 'db') THEN |
---|
1731 | ! In this case (where optP = 'db'), p is input & output (no depth->pressure conversion needed) |
---|
1732 | p(i) = depth(i) |
---|
1733 | ELSE |
---|
1734 | !print *,"optP =", optP, "end" |
---|
1735 | PRINT *,"optP must be either 'm' or 'db'" |
---|
1736 | STOP |
---|
1737 | ENDIF |
---|
1738 | |
---|
1739 | ! 2) Convert potential T to in-situ T (if input is Tpot, i.e. case for models): |
---|
1740 | IF (trim(optT) == 'Tpot' .OR. trim(optT) == 'tpot') THEN |
---|
1741 | tempot = temp(i) |
---|
1742 | ! This is the case for most models and some data |
---|
1743 | ! a) Convert the pot. temp on today's "ITS 90" scale to older IPTS 68 scale |
---|
1744 | ! (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote) |
---|
1745 | tempot68 = (tempot - 0.0002) / 0.99975 |
---|
1746 | ! b) Compute "in-situ Temperature" from "Potential Temperature" (both on IPTS 68) |
---|
1747 | tempis68 = sw_temp(sal(i), tempot68, p(i), 0.0_wp ) |
---|
1748 | ! c) Convert the in-situ temp on older IPTS 68 scale to modern scale (ITS 90) |
---|
1749 | tempis(i) = 0.99975*tempis68 + 0.0002 |
---|
1750 | ! Note: parts (a) and (c) above are tiny corrections; |
---|
1751 | ! part (b) is a big correction for deep waters (but zero at surface) |
---|
1752 | ELSEIF (trim(optT) == 'Tinsitu' .OR. trim(optT) == 'tinsitu') THEN |
---|
1753 | ! When optT = 'Tinsitu', tempis is input & output (no tempot needed) |
---|
1754 | tempis(i) = temp(i) |
---|
1755 | tempis68 = (temp(i) - 0.0002) / 0.99975 |
---|
1756 | ! dtempot68 = sw_ptmp(DBLE(sal(i)), DBLE(tempis68), DBLE(p), 0.0d0) |
---|
1757 | ! dtempot = 0.99975*dtempot68 + 0.0002 |
---|
1758 | ELSE |
---|
1759 | PRINT *,"optT must be either 'Tpot' or 'Tinsitu'" |
---|
1760 | PRINT *,"you specified optT =", trim(optT) |
---|
1761 | STOP |
---|
1762 | ENDIF |
---|
1763 | |
---|
1764 | ! ================================================================ |
---|
1765 | ! Carbonate chemistry computations |
---|
1766 | ! ================================================================ |
---|
1767 | IF (dic(i) > 0. .AND. dic(i) < 1.0e+4) THEN |
---|
1768 | ! Test to indicate if any of input variables are unreasonable |
---|
1769 | IF ( sal(i) < 0. & |
---|
1770 | .OR. alk(i) < 0. & |
---|
1771 | .OR. dic(i) < 0. & |
---|
1772 | .OR. sil(i) < 0. & |
---|
1773 | .OR. phos(i) < 0. & |
---|
1774 | .OR. sal(i) > 1e+3 & |
---|
1775 | .OR. alk(i) > 1e+3 & |
---|
1776 | .OR. dic(i) > 1e+3 & |
---|
1777 | .OR. sil(i) > 1e+3 & |
---|
1778 | .OR. phos(i) > 1e+3) THEN |
---|
1779 | PRINT *, 'i, icount, tempot, sal, alk, dic, sil, phos =', & |
---|
1780 | i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i) |
---|
1781 | ENDIF |
---|
1782 | ! Zero out any negative salinity, phosphate, silica, dic, and alk |
---|
1783 | IF (sal(i) < 0.0) THEN |
---|
1784 | ssal = 0.0 |
---|
1785 | ELSE |
---|
1786 | ssal = sal(i) |
---|
1787 | ENDIF |
---|
1788 | IF (phos(i) < 0.0) THEN |
---|
1789 | sphos = 0.0 |
---|
1790 | ELSE |
---|
1791 | sphos = phos(i) |
---|
1792 | ENDIF |
---|
1793 | IF (sil(i) < 0.0) THEN |
---|
1794 | ssil = 0.0 |
---|
1795 | ELSE |
---|
1796 | ssil = sil(i) |
---|
1797 | ENDIF |
---|
1798 | IF (dic(i) < 0.0) THEN |
---|
1799 | sdic = 0.0 |
---|
1800 | ELSE |
---|
1801 | sdic = dic(i) |
---|
1802 | ENDIF |
---|
1803 | IF (alk(i) < 0.0) THEN |
---|
1804 | salk = 0.0 |
---|
1805 | ELSE |
---|
1806 | salk = alk(i) |
---|
1807 | ENDIF |
---|
1808 | |
---|
1809 | ! Absolute temperature (Kelvin) & related variables |
---|
1810 | t = DBLE(tempis(i)) |
---|
1811 | tk = 273.15d0 + t |
---|
1812 | |
---|
1813 | ! Atmospheric pressure |
---|
1814 | Patmd = DBLE(Patm(i)) |
---|
1815 | ! Hydrostatic pressure (prb is in bars) |
---|
1816 | prb = DBLE(p(i)) / 10.0d0 |
---|
1817 | Phydro_atm = prb / 1.01325d0 ! convert hydrostatic pressure from bar to atm (1.01325 bar / atm) |
---|
1818 | ! Total pressure [atm] |
---|
1819 | IF (trim(opGAS) == 'Pzero' .OR. trim(opGAS) == 'pzero') THEN |
---|
1820 | Ptot = Patmd ! total pressure (in atm) = atmospheric pressure ONLY |
---|
1821 | ELSEIF (trim(opGAS) == 'Ppot' .OR. trim(opGAS) == 'ppot') THEN |
---|
1822 | Ptot = Patmd ! total pressure (in atm) = atmospheric pressure ONLY |
---|
1823 | ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN |
---|
1824 | Ptot = Patmd + Phydro_atm ! total pressure (in atm) = atmospheric pressure + hydrostatic pressure |
---|
1825 | ELSE |
---|
1826 | PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'" |
---|
1827 | STOP |
---|
1828 | ENDIF |
---|
1829 | |
---|
1830 | ! Salinity (equivalent array in double precision) |
---|
1831 | s = DBLE(ssal) |
---|
1832 | |
---|
1833 | ! Get all equilibrium constants and total concentrations of SO4, F, B |
---|
1834 | CALL constants(aK0, aK1, aK2, aKb, aKw, aKs, aKf, aKspc, aKspa, & |
---|
1835 | aK1p, aK2p, aK3p, aKsi, & |
---|
1836 | aSt, aFt, aBt, & |
---|
1837 | temp(i), sal(i), Patm(i), & |
---|
1838 | depth(i), lat(i), 1, & |
---|
1839 | optT, optP, opB, opK1K2, opKf, opGAS ) |
---|
1840 | |
---|
1841 | ! Unlike f77, in F90 we can't assign an array (dimen=1) to a scalar in a routine argument |
---|
1842 | ! Thus, set scalar constants equal to array (dimension=1) values required as arguments |
---|
1843 | K0 = aK0(1) ; K1 = aK1(1) ; K2 = aK2(1) ; Kb = aKb(1) ; Kw = aKw(1) |
---|
1844 | Ks = aKs(1) ; Kf = aKs(1) ; Kspc = aKspc(1) ; Kspa = aKspa(1) |
---|
1845 | K1p = aK1p(1) ; K2p = aK2p(1) ; K3p = aK3p(1) ; Ksi = aKsi(1) |
---|
1846 | St = aSt(1) ; Ft = aFt(1) ; Bt = aBt(1) |
---|
1847 | |
---|
1848 | ! Compute in-situ density [kg/m^3] |
---|
1849 | rhoSW(i) = rho(ssal, tempis68, prb) |
---|
1850 | |
---|
1851 | ! Either convert units of DIC and ALK (MODEL case) or not (DATA case) |
---|
1852 | IF (trim(optCON) == 'mol/kg') THEN |
---|
1853 | ! No conversion: |
---|
1854 | ! print *,'DIC and ALK already given in mol/kg (std DATA units)' |
---|
1855 | drho = 1. |
---|
1856 | ELSEIF (trim(optCON) == 'mol/m3') THEN |
---|
1857 | ! Do conversion: |
---|
1858 | ! print *,"DIC and ALK given in mol/m^3 (std MODEL units)" |
---|
1859 | drho = DBLE(rhoSW(i)) |
---|
1860 | ELSE |
---|
1861 | PRINT *,"optCON must be either 'mol/kg' or 'mol/m3'" |
---|
1862 | STOP |
---|
1863 | ENDIF |
---|
1864 | |
---|
1865 | tc = DBLE(sdic)/drho |
---|
1866 | ta = DBLE(salk)/drho |
---|
1867 | sit = DBLE(ssil)/drho |
---|
1868 | pt = DBLE(sphos)/drho |
---|
1869 | |
---|
1870 | ! Solve for pH and all other variables |
---|
1871 | ! ------------------------------------ |
---|
1872 | CALL varsolver(dph, dpco2, dfco2, dco2, dhco3, dco3, dOmegaA, dOmegaC, & |
---|
1873 | t, s, ta, tc, pt, sit, & |
---|
1874 | Bt, St, Ft, & |
---|
1875 | K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi, & |
---|
1876 | Patmd, prb, drho, opGAS ) |
---|
1877 | |
---|
1878 | ! Convert all output variables from double to single precision |
---|
1879 | pH(i) = REAL(dph) |
---|
1880 | co2(i) = REAL(dco2) |
---|
1881 | hco3(i) = REAL(dhco3) |
---|
1882 | co3(i) = REAL(dco3) |
---|
1883 | fCO2(i) = REAL(dfCO2) |
---|
1884 | pCO2(i) = REAL(dpCO2) |
---|
1885 | OmegaA(i) = REAL(dOmegaA) |
---|
1886 | OmegaC(i) = REAL(dOmegaC) |
---|
1887 | |
---|
1888 | ! Compute Revelle factor numerically (derivative using centered-difference scheme) |
---|
1889 | DO j=1,2 |
---|
1890 | minusplus = (-1)**j |
---|
1891 | dx = 0.1 * 1e-6 ! Numerical tests found for DIC that optimal dx = 0.1 umol/kg (0.1e-6 mol/kg) |
---|
1892 | dicdel(j) = tc + DBLE(minusplus)*dx/2.0d0 |
---|
1893 | CALL varsolver(dph, dpco2, dfco2, dco2, dhco3, dco3, dOmegaA, dOmegaC, & |
---|
1894 | t, s, ta, dicdel(j), pt, sit, & |
---|
1895 | Bt, St, Ft, & |
---|
1896 | K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi, & |
---|
1897 | Patmd, prb, drho, optGAS ) |
---|
1898 | pco2del(j) = dpco2 |
---|
1899 | END DO |
---|
1900 | !Classic finite centered difference formula for derivative (2nd order accurate) |
---|
1901 | Rf = (pco2del(2) - pco2del(1)) / (dicdel(2) - dicdel(1)) ! dpCO2/dDIC |
---|
1902 | !Rf = (pco2del(2) - pco2del(1)) / (dx) ! dpCO2/dDIC (same as just above) |
---|
1903 | Rf = Rf * tc / dpco2 ! R = (dpCO2/dDIC) * (DIC/pCO2) |
---|
1904 | |
---|
1905 | BetaD(i) = REAL(Rf) |
---|
1906 | |
---|
1907 | ELSE |
---|
1908 | |
---|
1909 | ph(i) = 1.e20_wp |
---|
1910 | pco2(i) = 1.e20_wp |
---|
1911 | fco2(i) = 1.e20_wp |
---|
1912 | co2(i) = 1.e20_wp |
---|
1913 | hco3(i) = 1.e20_wp |
---|
1914 | co3(i) = 1.e20_wp |
---|
1915 | OmegaA(i) = 1.e20_wp |
---|
1916 | OmegaC(i) = 1.e20_wp |
---|
1917 | BetaD(i) = 1.e20_wp |
---|
1918 | rhoSW(i) = 1.e20_wp |
---|
1919 | p(i) = 1.e20_wp |
---|
1920 | tempis(i) = 1.e20_wp |
---|
1921 | |
---|
1922 | ENDIF |
---|
1923 | |
---|
1924 | END DO |
---|
1925 | |
---|
1926 | RETURN |
---|
1927 | END SUBROUTINE vars |
---|
1928 | |
---|
1929 | ! ---------------------------------------------------------------------- |
---|
1930 | ! P2FCO2 |
---|
1931 | ! ---------------------------------------------------------------------- |
---|
1932 | ! |
---|
1933 | !> \file p2fCO2.f90 |
---|
1934 | !! \BRIEF |
---|
1935 | !> Module with p2fCO2 subroutine - compute fCO2 from pCO2, in situ T, atm pressure, hydrostatic pressure |
---|
1936 | !> Compute fCO2 from arrays of pCO2, in situ temp, atm pressure, & hydrostatic pressure |
---|
1937 | SUBROUTINE p2fCO2(pCO2, temp, Patm, p, N, fCO2) |
---|
1938 | ! Purpose: |
---|
1939 | ! Compute fCO2 from arrays of pCO2, in situ temp, atm pressure, & hydrostatic pressure |
---|
1940 | |
---|
1941 | USE mocsy_singledouble |
---|
1942 | IMPLICIT NONE |
---|
1943 | |
---|
1944 | !> number of records |
---|
1945 | INTEGER, INTENT(in) :: N |
---|
1946 | |
---|
1947 | ! INPUT variables |
---|
1948 | !> oceanic partial pressure of CO2 [uatm] |
---|
1949 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: pCO2 |
---|
1950 | !> in situ temperature [C] |
---|
1951 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp |
---|
1952 | !> atmospheric pressure [atm] |
---|
1953 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm |
---|
1954 | !> hydrostatic pressure [db] |
---|
1955 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: p |
---|
1956 | |
---|
1957 | ! OUTPUT variables: |
---|
1958 | !> fugacity of CO2 [uatm] |
---|
1959 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: fCO2 |
---|
1960 | |
---|
1961 | ! LOCAL variables: |
---|
1962 | REAL(kind=wp) :: dpCO2, dtemp, tk, dPatm, prb |
---|
1963 | REAL(kind=wp) :: Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff |
---|
1964 | REAL(kind=wp) :: dfCO2 |
---|
1965 | |
---|
1966 | INTEGER :: i |
---|
1967 | |
---|
1968 | ! REAL(kind=wp) :: sw_ptmp |
---|
1969 | ! EXTERNAL sw_ptmp |
---|
1970 | |
---|
1971 | DO i = 1,N |
---|
1972 | dpCO2 = DBLE(pCO2(i)) |
---|
1973 | dtemp = DBLE(temp(i)) |
---|
1974 | dPatm = DBLE(Patm(i)) |
---|
1975 | tk = 273.15d0 + DBLE(temp(i)) !Absolute temperature (Kelvin) |
---|
1976 | prb = DBLE(p(i)) / 10.0d0 !Pressure effect (prb is in bars) |
---|
1977 | Ptot = dPatm + prb/1.01325d0 !Total pressure (atmospheric + hydrostatic) [atm] |
---|
1978 | Rgas_atm = 82.05736_wp !R in (cm3 * atm) / (mol * K) from CODATA (2006) |
---|
1979 | ! To compute fugcoeff, we need 3 other terms (B, Del, xc2) as well as 3 others above (tk, Ptot, Rgas_atm) |
---|
1980 | B = -1636.75d0 + 12.0408d0*tk - 0.0327957d0*(tk*tk) + 0.0000316528d0*(tk*tk*tk) |
---|
1981 | Del = 57.7d0 - 0.118d0*tk |
---|
1982 | ! "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9 |
---|
1983 | ! x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite) |
---|
1984 | ! Let's assume that xCO2 = pCO2. Resulting fugcoeff is identical to 8th digit after the decimal. |
---|
1985 | xCO2approx = dpCO2 * 1.e-6_wp |
---|
1986 | xc2 = (1.0d0 - xCO2approx)**2 |
---|
1987 | fugcoeff = EXP( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk) ) |
---|
1988 | dfCO2 = dpCO2 * fugcoeff |
---|
1989 | fCO2(i) = REAL(dfCO2) |
---|
1990 | END DO |
---|
1991 | |
---|
1992 | RETURN |
---|
1993 | END SUBROUTINE p2fCO2 |
---|
1994 | |
---|
1995 | ! ---------------------------------------------------------------------- |
---|
1996 | ! P2FCO2 |
---|
1997 | ! ---------------------------------------------------------------------- |
---|
1998 | ! |
---|
1999 | !> \file f2pCO2.f90 |
---|
2000 | !! \BRIEF |
---|
2001 | !> Module with f2pCO2 subroutine - compute pCO2 from fCO2, in situ T, atm pressure, hydrostatic pressure |
---|
2002 | !> Compute pCO2 from arrays of fCO2, in situ temp, atm pressure, & hydrostatic pressure |
---|
2003 | SUBROUTINE f2pCO2(fCO2, temp, Patm, p, N, pCO2) |
---|
2004 | ! Purpose: |
---|
2005 | ! Compute pCO2 from arrays of fCO2, in situ temp, atm pressure, & hydrostatic pressure |
---|
2006 | |
---|
2007 | USE mocsy_singledouble |
---|
2008 | IMPLICIT NONE |
---|
2009 | |
---|
2010 | !> number of records |
---|
2011 | INTEGER, intent(in) :: N |
---|
2012 | |
---|
2013 | ! INPUT variables |
---|
2014 | !> oceanic fugacity of CO2 [uatm] |
---|
2015 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: fCO2 |
---|
2016 | !> in situ temperature [C] |
---|
2017 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp |
---|
2018 | !> atmospheric pressure [atm] |
---|
2019 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm |
---|
2020 | !> hydrostatic pressure [db] |
---|
2021 | REAL(kind=wp), INTENT(in), DIMENSION(N) :: p |
---|
2022 | |
---|
2023 | ! OUTPUT variables: |
---|
2024 | !> oceanic partial pressure of CO2 [uatm] |
---|
2025 | REAL(kind=wp), INTENT(out), DIMENSION(N) :: pCO2 |
---|
2026 | |
---|
2027 | ! LOCAL variables: |
---|
2028 | REAL(kind=wp) :: dfCO2, dtemp, tk, dPatm, prb |
---|
2029 | REAL(kind=wp) :: Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff |
---|
2030 | REAL(kind=wp) :: dpCO2 |
---|
2031 | |
---|
2032 | INTEGER :: i |
---|
2033 | |
---|
2034 | ! REAL(kind=wp) :: sw_ptmp |
---|
2035 | ! EXTERNAL sw_ptmp |
---|
2036 | |
---|
2037 | DO i = 1,N |
---|
2038 | dfCO2 = DBLE(fCO2(i)) |
---|
2039 | dtemp = DBLE(temp(i)) |
---|
2040 | dPatm = DBLE(Patm(i)) |
---|
2041 | tk = 273.15d0 + DBLE(temp(i)) !Absolute temperature (Kelvin) |
---|
2042 | prb = DBLE(p(i)) / 10.0d0 !Pressure effect (prb is in bars) |
---|
2043 | Ptot = dPatm + prb/1.01325d0 !Total pressure (atmospheric + hydrostatic) [atm] |
---|
2044 | Rgas_atm = 82.05736_wp !R in (cm3 * atm) / (mol * K) from CODATA (2006) |
---|
2045 | ! To compute fugcoeff, we need 3 other terms (B, Del, xc2) as well as 3 others above (tk, Ptot, Rgas_atm) |
---|
2046 | B = -1636.75d0 + 12.0408d0*tk - 0.0327957d0*(tk*tk) + 0.0000316528d0*(tk*tk*tk) |
---|
2047 | Del = 57.7d0 - 0.118d0*tk |
---|
2048 | ! "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9 |
---|
2049 | ! x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite) |
---|
2050 | ! Let's assume that xCO2 = fCO2. Resulting fugcoeff is identical to 8th digit after the decimal. |
---|
2051 | xCO2approx = dfCO2 * 1.e-6_wp |
---|
2052 | xc2 = (1.0d0 - xCO2approx)**2 |
---|
2053 | fugcoeff = exp( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk) ) |
---|
2054 | dpCO2 = dfCO2 / fugcoeff |
---|
2055 | pCO2(i) = REAL(dpCO2) |
---|
2056 | END DO |
---|
2057 | |
---|
2058 | RETURN |
---|
2059 | END SUBROUTINE f2pCO2 |
---|
2060 | |
---|
2061 | END MODULE mocsy_mainmod |
---|