*IF DEF,OCEAN ORH1F305.4096
C ******************************COPYRIGHT****************************** GTS2F400.2485
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.2486
C GTS2F400.2487
C Use, duplication or disclosure of this code is subject to the GTS2F400.2488
C restrictions as set forth in the contract. GTS2F400.2489
C GTS2F400.2490
C Meteorological Office GTS2F400.2491
C London Road GTS2F400.2492
C BRACKNELL GTS2F400.2493
C Berkshire UK GTS2F400.2494
C RG12 2SZ GTS2F400.2495
C GTS2F400.2496
C If no contract has been raised with this copy of the code, the use, GTS2F400.2497
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.2498
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.2499
C Modelling at the above address. GTS2F400.2500
C ******************************COPYRIGHT****************************** GTS2F400.2501
C GTS2F400.2502
C*LL Subroutine EQUILIB_CONST EQ_CONST.3
CLL Can run on any FORTRAN 77 compiler with long lower case variables EQ_CONST.4
CLL EQ_CONST.5
CLL Author: N.K. TAYLOR EQ_CONST.6
CLL Date: 17 December 1993 EQ_CONST.7
CLL Version 3.3 EQ_CONST.8
CLL EQ_CONST.9
! Modification History: ORH1F305.4097
! Version Date Details ORH1F305.4098
! ------- ------- ------------------------------------------ ORH1F305.4099
! 3.5 16.01.95 Remove *IF dependency. R.Hill ORH1F305.4100
CLL Programming standards use Cox naming convention for Cox variables EQ_CONST.10
CLL with the addition that lower case variables are local to the EQ_CONST.11
CLL routine. EQ_CONST.12
CLL EQ_CONST.13
CLL This routine uses surface values of temperature and salinity along EQ_CONST.14
CLL a row to compute chemical equilibrium constants needed for the EQ_CONST.15
CLL ocean Carbon Cycle. The Henrys Law constant, REK0, is EQ_CONST.16
CLL calculated according to Weiss (1974). EQ_CONST.17
CLL The other equilibrium constants are computed from equations given EQ_CONST.18
CLL by Peng, Takahashi and Broecker (1987). EQ_CONST.19
CLL EQ_CONST.20
CLLEND ----------------------------------------------------------------- EQ_CONST.21
C* EQ_CONST.22
C*L-------------------------------- Arguments ------------------------ EQ_CONST.23
C EQ_CONST.24
SUBROUTINE EQUILIB_CONST (TA, REK0, REK, RG1, RG2, RG3, 1,5EQ_CONST.25
+ IMT, KM, NT) EQ_CONST.26
C EQ_CONST.27
IMPLICIT NONE EQ_CONST.28
C EQ_CONST.29
*CALL CNTLOCN
ORH1F305.4101
*CALL OARRYSIZ
ORH1F305.4102
C Define constants for array sizes EQ_CONST.30
C EQ_CONST.31
INTEGER EQ_CONST.32
+ IMT ! IN Number of points in horizontal EQ_CONST.33
+, KM ! IN Number of layers in model EQ_CONST.34
+, NT ! IN Number of tracers EQ_CONST.35
C EQ_CONST.36
C Physical arguments EQ_CONST.37
C EQ_CONST.38
REAL EQ_CONST.39
+ TA (IMT, KM, NT) ! IN Tracer values EQ_CONST.40
+, REK0 (IMT) ! OUT Solubility of CO2 (Henrys Law const) EQ_CONST.41
+, RG1 (IMT) ! OUT Products etc of the Equilibrium EQ_CONST.42
+, RG2 (IMT) ! OUT constants, which are used for EQ_CONST.43
+, RG3 (IMT) ! OUT the calculation of PCO2 EQ_CONST.44
+, REK (IMT) EQ_CONST.45
C* EQ_CONST.46
C EQ_CONST.47
C EQ_CONST.48
C Locally defined variables EQ_CONST.49
C EQ_CONST.50
INTEGER EQ_CONST.51
+ i ! Horizontal loop index EQ_CONST.52
C EQ_CONST.53
REAL rek1 ! Equilibrium constants for the EQ_CONST.54
+, rek2 ! reactions of CO2 with water EQ_CONST.55
+, rekb ! EQ_CONST.56
+, rekw ! EQ_CONST.57
+, tc ! Surface temperature in deg C OJP0F404.292
+, sal ! Surface salty in parts per mil EQ_CONST.59
EQ_CONST.61
C EQ_CONST.62
! Others OJP0F404.293
! OJP0F404.294
REAL CO2smolk,kone,ktwo,kboric,kwater ! FUNCTIONS to calc consts OJP0F404.295
OJP0F404.296
CL Compute equilibrium constants (used as intermediate constants only) EQ_CONST.63
C and the products etc of the constants (which are used later to EQ_CONST.64
C compute the partial pressure of CO2 EQ_CONST.65
! Acid dissoc. constants rek1,rek2,rekb are in units of moles per kg. OJP0F404.297
! Water dissociation constant is in units of (moles per kg)^2. OJP0F404.298
! Henry's law coeff. REK0 has units of moles per (kg.atm). OJP0F404.299
C EQ_CONST.66
OJP0F404.300
DO I = IFROM_CYC,ITO_CYC ORH1F305.4103
! Convert salinity into psu OJP0F404.301
sal = TA(I,1,2)*1000.0 + 35.0 EQ_CONST.75
tc = TA(I,1,1) OJP0F404.302
EQ_CONST.77
C Calculate dissociation constants for carbonate system EQ_CONST.78
EQ_CONST.79
REK0(I) = CO2smolk
(tc,sal) OJP0F404.303
rek1 = kone
(tc,sal) OJP0F404.304
rek2 = ktwo
(tc,sal) OJP0F404.305
rekb = kboric
(tc,sal) OJP0F404.306
rekw = kwater
(tc,sal) OJP0F404.307
EQ_CONST.83
C Calculate constants used for computation of pCO2 EQ_CONST.99
C EQ_CONST.100
REK(I) = sqrt(rek1/rek2) EQ_CONST.101
RG3(I) = sqrt(rek1*rek2) EQ_CONST.102
RG1(I) = RG3(I)/rekb EQ_CONST.103
RG2(I) = rekw/RG3(I) EQ_CONST.104
EQ_CONST.105
ENDDO EQ_CONST.106
C EQ_CONST.107
CL Return from EQUILIB_CONST EQ_CONST.108
C EQ_CONST.109
RETURN EQ_CONST.110
END EQ_CONST.111
! The following chemistry coefficient code is imported verbatim. OJP0F404.308
c OJP0F404.309
c Functions to calculate equilibrium constants and solubility of OJP0F404.310
c carbon dioxide in sea water. Taken from the U.S. Department of OJP0F404.311
c Energy's "Handbook of Methods for the Analysis of the OJP0F404.312
c Various Parameters of the Carbon Dioxide System in Sea Water", OJP0F404.313
c version 2, ed. by A.G. Dickson and C. Goyet, September, 1994. OJP0F404.314
c (Prepared for the U.S. Dept. of Energy, Special Research Grant OJP0F404.315
c Program 89-7A: Global Survey of Carbon dioxide in the Oceans). OJP0F404.316
OJP0F404.317
c N. Lefevre, OCMIP, LMCE/DSM/CE, Saclay FRANCE, April 1995 OJP0F404.318
OJP0F404.319
c REFERENCES: OJP0F404.320
c Roy,, R. N. L. N. Roy, M. Vogel, C. P. Moore, T. Pearson, OJP0F404.321
c C. E. Good, F. J. Millero, and D. J. Cambell, 1993. OJP0F404.322
c Determination of the ionization constants of carbonic acid in OJP0F404.323
c seawater, Mar. Chem., 44, 249--268. OJP0F404.324
OJP0F404.325
c Dickson, A.E., 1990. Thermodynamics of the dissociation of boric OJP0F404.326
c acid in syntheic sea water from 273.15 to 298.15 K. Deep-Sea OJP0F404.327
c Research, 37, 755-766. OJP0F404.328
OJP0F404.329
c Millero, F. J., 1994. The carbon dioxide system in the oceans, OJP0F404.330
c Geoch. Cosmch. Acta, (submitted) OJP0F404.331
OJP0F404.332
c Weiss, R. F., 1974. Carbon dioxide in water and seawater: the OJP0F404.333
c solubility of a non-ideal gas. Mar. Chem., 2, 203-215. OJP0F404.334
c OJP0F404.335
c ------------------------------------------------------------------ OJP0F404.336
real function kone(tc,sal) 1OJP0F404.337
c OJP0F404.338
c Calculates the 1st constant of the carbonic acid from Roy et al. OJP0F404.339
c (1993). It is expressed in mol/(kg.soln) on the total hydrogen OJP0F404.340
c ion scale. OJP0F404.341
c OJP0F404.342
c For this code Ln(kone(25.,35.)) = -13.4847 OJP0F404.343
c as it should according to Dickson and Goyet OJP0F404.344
c OJP0F404.345
c List of variables : OJP0F404.346
c tc temperature (celsius) OJP0F404.347
c tk temperature (Kelvin) OJP0F404.348
c sal salinity OJP0F404.349
c K1 first constant of the carbonic acid OJP0F404.350
c (mol/kg.soln) OJP0F404.351
c kone first constant of the carbonic acid OJP0F404.352
c (mol/kg.soln) OJP0F404.353
c A1 constant (celsius) OJP0F404.354
c A2 constant OJP0F404.355
c A3 constant OJP0F404.356
c B1 constant (celsius) OJP0F404.357
c B2 constant OJP0F404.358
c B3 constant OJP0F404.359
c B4 constant OJP0F404.360
c B5 constant OJP0F404.361
OJP0F404.362
implicit none OJP0F404.363
OJP0F404.364
real K1, tc, sal, tk, A1, A2, A3, B1, B2, B3 OJP0F404.365
real B4, B5 OJP0F404.366
data A1, A2, A3, B1, B2, B3, B4, B5 / -2307.1266, 2.83655, OJP0F404.367
. -1.5529413, -4.0484, -0.20760841, 0.08468345, OJP0F404.368
. -0.00654208, -0.001005/ OJP0F404.369
OJP0F404.370
tk = tc + 273.15 OJP0F404.371
K1 = exp(A1/tk + A2 + A3 * log(tk) OJP0F404.372
. + (B1/tk + B2) * sqrt(sal) + B3 * sal OJP0F404.373
. + B4 * sal**1.5 + log(1 + B5 * sal)) OJP0F404.374
kone = K1 OJP0F404.375
return OJP0F404.376
end OJP0F404.377
c OJP0F404.378
c ------------------------------------------------ OJP0F404.379
real function ktwo(tc,sal) 1OJP0F404.380
c OJP0F404.381
c Calculates the 2nd constant of the carbonic acid from Roy et al. OJP0F404.382
c (1993). It is expressed in mol/(kg.soln) on the total hydrogen OJP0F404.383
c ion scale. OJP0F404.384
c OJP0F404.385
c For this code Ln(ktwo(25.,35.)) = -20.5504 OJP0F404.386
c as it should according to Dickson and Goyet OJP0F404.387
c OJP0F404.388
c List of variables : OJP0F404.389
c tc temperature (celsius) OJP0F404.390
c tk temperature (Kelvin) OJP0F404.391
c sal salinity OJP0F404.392
c K2 second constant of the carbonic acid OJP0F404.393
c (mol/kg.soln) OJP0F404.394
c ktwo second constant of the carbonic acid OJP0F404.395
c (mol/kg.soln) OJP0F404.396
c A1 constant (celsius) OJP0F404.397
c A2 constant OJP0F404.398
c A3 constant OJP0F404.399
c B1 constant (celsius) OJP0F404.400
c B2 constant OJP0F404.401
c B3 constant OJP0F404.402
c B4 constant OJP0F404.403
c B5 constant OJP0F404.404
OJP0F404.405
implicit none OJP0F404.406
real K2, tc, sal, tk, A1, A2, A3, B1, B2, B3 OJP0F404.407
real B4, B5 OJP0F404.408
data A1, A2, A3, B1, B2, B3, B4, B5 / -3351.6106, -9.226508 OJP0F404.409
. ,-0.2005743, -23.9722, -0.106901773, 0.1130822, OJP0F404.410
. -0.00846934, -0.001005/ OJP0F404.411
OJP0F404.412
OJP0F404.413
tk = tc + 273.15 OJP0F404.414
K2 = exp(A1/tk + A2 + A3 * log(tk) OJP0F404.415
. + (B1/tk + B2) * sqrt(sal) + B3 * sal + OJP0F404.416
. B4 * sal**1.5 + log(1 + B5 * sal)) OJP0F404.417
ktwo = K2 OJP0F404.418
return OJP0F404.419
end OJP0F404.420
c OJP0F404.421
c ------------------------------------------------- OJP0F404.422
real function kboric(tc,sal) 1OJP0F404.423
c OJP0F404.424
c Calculates the equilibrium constant of the boric acid from OJP0F404.425
c Dickson (1990). It is expressed in mol/(kg.soln) on the total OJP0F404.426
c hydrogen ion scale. OJP0F404.427
c OJP0F404.428
c For this code, Ln(kboric(25.,35.)) = -19.7964 OJP0F404.429
c as it should according to Dickson and Goyet OJP0F404.430
c OJP0F404.431
c List of variables : OJP0F404.432
c tc temperature (celsius) OJP0F404.433
c tk temperature (Kelvin) OJP0F404.434
c sal salinity OJP0F404.435
c KB constant of the boric acid (mol/kg.soln) OJP0F404.436
c (mol/kg.soln) OJP0F404.437
c kboric constant of the carbonic acid (mol/kg.soln) OJP0F404.438
c A1 constant OJP0F404.439
c A2 constant OJP0F404.440
c A3 constant OJP0F404.441
c A4 constant OJP0F404.442
c A5 constant OJP0F404.443
c B1 constant OJP0F404.444
c B2 constant OJP0F404.445
c B3 constant OJP0F404.446
c C1 constant OJP0F404.447
c C2 constant OJP0F404.448
c C3 constant OJP0F404.449
c D constant OJP0F404.450
OJP0F404.451
implicit none OJP0F404.452
real KB, tc, sal, tk, A1, A2, A3, A4, A5, B1, OJP0F404.453
. B2, B3, C1, C2, C3, D OJP0F404.454
data A1, A2, A3, A4, A5, B1, B2, B3, C1, C2, C3, D / OJP0F404.455
. -8966.90, -2890.53, -77.942, 1.728, -0.0996, OJP0F404.456
. 148.0248, 137.1942, 1.62142, -24.4344, -25.085, OJP0F404.457
. -0.2474, 0.053105/ OJP0F404.458
OJP0F404.459
tk = tc + 273.15 OJP0F404.460
KB = exp((A1 + A2 * sqrt(sal) + A3 * sal + A4 * sal**1.5 OJP0F404.461
. + A5 * sal**2)/tk + (B1 + B2 * sqrt(sal) + B3 * sal) OJP0F404.462
. + (C1 + C2 * sqrt(sal) + C3 * sal) * log(tk) OJP0F404.463
. + D * sqrt(sal) * tk) OJP0F404.464
kboric = KB OJP0F404.465
return OJP0F404.466
end OJP0F404.467
c OJP0F404.468
c ------------------------------------------------- OJP0F404.469
real function kwater(tc,sal) 1OJP0F404.470
c OJP0F404.471
c Calculates the equilibrium constant of the water from Millero OJP0F404.472
c (1994). It is expressed in mol/(kg.soln) on the total hydrogen OJP0F404.473
c ion scale. OJP0F404.474
c OJP0F404.475
c For this code, Ln(kwater(25.,35.)) = -30.434 OJP0F404.476
c as it should according to Dickson and Goyet OJP0F404.477
c OJP0F404.478
c List of variables : OJP0F404.479
c tc temperature (celsius) OJP0F404.480
c tk temperature (Kelvin) OJP0F404.481
c sal salinity OJP0F404.482
c KW constant of the water (mol/kg.soln) OJP0F404.483
c (mol/kg.soln) OJP0F404.484
c kwater constant of the water (mol/kg.soln) OJP0F404.485
c A1 constant (celsius) OJP0F404.486
c A2 constant OJP0F404.487
c A3 constant OJP0F404.488
c B1 constant (celsius) OJP0F404.489
c B2 constant OJP0F404.490
c B3 constant OJP0F404.491
c C constant OJP0F404.492
OJP0F404.493
implicit none OJP0F404.494
real KW, tc, sal, tk, A1, A2, A3, B1, B2, B3, C OJP0F404.495
data A1, A2, A3, B1, B2, B3, C / -13847.26, 148.9652, OJP0F404.496
. -23.6521, 118.67, -5.977, 1.0495, -0.01615/ OJP0F404.497
OJP0F404.498
tk = tc + 273.15 OJP0F404.499
KW = exp(A1/tk + A2 + A3 * log(tk) OJP0F404.500
. + (B1/tk + B2 + B3 * log(tk)) * sqrt(sal) + C * sal) OJP0F404.501
kwater = KW OJP0F404.502
return OJP0F404.503
end OJP0F404.504
c OJP0F404.505
c -------------------------------------------------- OJP0F404.506
real function CO2smolk(tc,sal) 1OJP0F404.507
c OJP0F404.508
c Calculates the solubility of carbon dioxide in sea water from OJP0F404.509
c Weiss (1974). It is expressed in mol/(kg.atm). OJP0F404.510
c OJP0F404.511
c For this code, Ln(CO2smolk(25.,35.)) = -3.5617 OJP0F404.512
c as it should according to Dickson and Goyet OJP0F404.513
c OJP0F404.514
c List of variables : OJP0F404.515
c tc temperature (celsius) OJP0F404.516
c tk temperature (Kelvin) OJP0F404.517
c sal salinity OJP0F404.518
c K0 CO2 solubility coefficient (mol/kg.atm) OJP0F404.519
c CO2smolk CO2 solubility coefficient (mol/kg.atm) OJP0F404.520
c A1 constant (mol/kg.atm) OJP0F404.521
c A2 constant (mol/kg.atm) OJP0F404.522
c A3 constant (mol/kg.atm) OJP0F404.523
c B1 constant (mol/kg.atm) OJP0F404.524
c B2 constant (mol/kg.atm) OJP0F404.525
c B3 constant (mol/kg.atm) OJP0F404.526
OJP0F404.527
implicit none OJP0F404.528
real K0, tc, sal, tk, A1, A2, A3, B1, B2, B3 OJP0F404.529
data A1, A2, A3, B1, B2, B3 / -60.2409, 93.4517, OJP0F404.530
. 23.3585, 0.023517, -0.023656, 0.0047036/ OJP0F404.531
OJP0F404.532
tk = tc + 273.15 OJP0F404.533
K0 = exp(A1 + A2 * (100/tk) + A3 * log(tk/100) OJP0F404.534
. + sal * (B1 + B2 * (tk/100) + B3 * (tk/100)**2)) OJP0F404.535
CO2smolk = K0 OJP0F404.536
return OJP0F404.537
end OJP0F404.538
c OJP0F404.539
c ------------------------------------------------------ OJP0F404.540
real function CO2smoll(tc,sal) OJP0F404.541
c OJP0F404.542
c This function calculates the solubility of carbon OJP0F404.543
c dioxide in mol/(l.atm). OJP0F404.544
c OJP0F404.545
c For this code, CO2smoll(20.,35.) = 0.03322 OJP0F404.546
c as it given in Table II of Weiss (1974) OJP0F404.547
c OJP0F404.548
c List of variables : OJP0F404.549
c tc temperature (celsius) OJP0F404.550
c tk temperature (Kelvin) OJP0F404.551
c sal salinity OJP0F404.552
c K0 CO2 solubility coefficient (mol/l.atm) OJP0F404.553
c CO2smoll CO2 solubility coefficient (mol/l.atm) OJP0F404.554
c A1 constant (mol/l.atm) OJP0F404.555
c A2 constant (mol/l.atm) OJP0F404.556
c A3 constant (mol/l.atm) OJP0F404.557
c B1 constant (mol/l.atm) OJP0F404.558
c B2 constant (mol/l.atm) OJP0F404.559
c B3 constant (mol/l.atm) OJP0F404.560
OJP0F404.561
implicit none OJP0F404.562
real K0, tc, sal, tk, A1, A2, A3, B1, B2, B3 OJP0F404.563
data A1, A2, A3, B1, B2, B3 /-58.0931, 90.5069, OJP0F404.564
. 22.2940, 0.027766, -0.025888, 0.0050578/ OJP0F404.565
OJP0F404.566
tk = tc + 273.15 OJP0F404.567
K0 = exp(A1 + A2 * (100/tk) + A3 * log(tk/100) OJP0F404.568
. + sal * (B1 + B2 * (tk/100) + B3 * (tk/100)**2)) OJP0F404.569
CO2smoll = K0 OJP0F404.570
return OJP0F404.571
end OJP0F404.572
*ENDIF EQ_CONST.112