[5841] | 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 |
---|
[7894] | 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 |
---|
[5841] | 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 |
---|
[7894] | 679 | !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010) or \b 'w14' (Waters et al., 2014) |
---|
[5841] | 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) |
---|
[7894] | 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) |
---|
[5841] | 979 | ELSE |
---|
[7894] | 980 | PRINT *, "optK1K2 must be either 'l' or 'm10', or 'w14'" |
---|
[5841] | 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 |
---|