source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/mocsy_mainmod.F90 @ 5841

Last change on this file since 5841 was 5841, checked in by jpalmier, 5 years ago

JPALM —30-10-2015— Add MOCSY and DMS to MEDUSA-NEMO3.6

File size: 83.7 KB
Line 
1MODULE mocsy_mainmod
2
3      USE in_out_manager  ! I/O manager
4
5CONTAINS
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.
15FUNCTION 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
61END 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
71FUNCTION 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
131END 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]
142FUNCTION 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
209END 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)
220SUBROUTINE 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
266END 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)
277SUBROUTINE 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
324END 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.
334FUNCTION 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
376END 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)
386FUNCTION 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
466END 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)
477SUBROUTINE 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
518END 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)
530SUBROUTINE 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
564END 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)
576SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa,  &
577                     K1p, K2p, K3p, Ksi,                      &
578                     St, Ft, Bt,                              &
579                     temp, sal, Patm,                         &
580                     depth, lat, N,                           &
581                     optT, optP, optB, optK1K2, optKf, optGAS)
582
583  !   Purpose:
584  !     Compute thermodynamic constants
585  !     FROM: temperature, salinity, and pressure (1D arrays)
586
587  !     INPUT variables:
588  !     ================
589  !     Patm    = atmospheric pressure [atm]
590  !     depth   = depth [m]     (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure)
591  !             = pressure [db] (with optP='db')
592  !     lat     = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m')
593  !             = dummy array (unused when optP='db')
594  !     temp    = potential temperature [degrees C] (with optT='Tpot', i.e., models carry tempot, not temp)
595  !             = in situ   temperature [degrees C] (with optT='Tinsitu', e.g., for data)
596  !     sal     = salinity in [psu]
597  !     ---------
598  !     optT: choose in situ vs. potential temperature as input
599  !     ---------
600  !     NOTE: Carbonate chem calculations require IN-SITU temperature (not potential Temperature)
601  !       -> 'Tpot' means input is pot. Temperature (in situ Temp "tempis" is computed)
602  !       -> 'Tinsitu' means input is already in-situ Temperature, not pot. Temp ("tempis" not computed)
603  !     ---------
604  !     optP: choose depth (m) vs pressure (db) as input
605  !     ---------
606  !       -> 'm'  means "depth" input is in "m" (thus in situ Pressure "p" [db] is computed)
607  !       -> 'db' means "depth" input is already in situ pressure [db], not m (p = depth)
608  !     ---------
609  !     optB:
610  !     ---------
611  !       -> 'u74' means use classic formulation of Uppström (1974) for total Boron
612  !       -> 'l10' means use newer   formulation of Lee et al. (2010) for total Boron
613  !     ---------
614  !     optK1K2:
615  !     ---------
616  !       -> 'l'   means use Lueker et al. (2000) formulations for K1 & K2 (recommended by Dickson et al. 2007)
617  !                **** BUT this should only be used when 2 < T < 35 and 19 < S < 43
618  !       -> 'm10' means use Millero (2010) formulation for K1 & K2 (see Dickson et al., 2007)
619  !                **** Valid for 0 < T < 50°C and 1 < S < 50 psu
620  !     -----------
621  !     optKf:
622  !     ----------
623  !       -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007)
624  !               **** BUT Valid for  9 < T < 33°C and 10 < S < 40.
625  !       -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994)
626  !     -----------
627  !     optGAS: choose in situ vs. potential fCO2 and pCO2
628  !     ---------
629  !       PRESSURE corrections for K0 and the fugacity coefficient (Cf)
630  !       -> 'Pzero'   = 'zero order' fCO2 and pCO2 (typical approach, which is flawed)
631  !                      considers in situ T & only atm pressure (hydrostatic=0)
632  !       -> 'Ppot'    = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
633  !                      considers potential T & only atm pressure (hydrostatic press = 0)
634  !       -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure)
635  !                      considers in situ T & total pressure (atm + hydrostatic)
636  !     ---------
637
638  !     OUTPUT variables:
639  !     =================
640  !     K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi
641  !     St, Ft, Bt
642
643  USE mocsy_singledouble
644  IMPLICIT NONE
645
646! Input variables
647  !>     number of records
648  INTEGER, INTENT(in) :: N
649  !> in <b>situ temperature</b> (when optT='Tinsitu', typical data)
650  !! OR <b>potential temperature</b> (when optT='Tpot', typical models) [degree C]
651  REAL(kind=wp), INTENT(in),    DIMENSION(N) :: temp
652  !> depth in <b>meters</b> (when optP='m') or <b>decibars</b> (when optP='db')
653  REAL(kind=wp), INTENT(in),    DIMENSION(N) :: depth
654  !> latitude <b>[degrees north]</b>
655  REAL(kind=wp), INTENT(in),    DIMENSION(N) :: lat
656  !> salinity <b>[psu]</b>
657  REAL(kind=wp), INTENT(in), DIMENSION(N) :: sal
658!f2py optional , depend(sal) :: n=len(sal)
659
660  !> atmospheric pressure <b>[atm]</b>
661  REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm
662
663  !> for temp input, choose \b 'Tinsitu' for in situ Temp or
664  !! \b 'Tpot' for potential temperature (in situ Temp is computed, needed for models)
665  CHARACTER(7), INTENT(in) :: optT
666  !> for depth input, choose \b "db" for decibars (in situ pressure) or \b "m" for meters (pressure is computed, needed for models)
667  CHARACTER(2), INTENT(in) :: optP
668  !> for total boron, choose either \b 'u74' (Uppstrom, 1974) or \b 'l10' (Lee et al., 2010).
669  !! The 'l10' formulation is based on 139 measurements (instead of 20),
670  !! uses a more accurate method, and
671  !! generally increases total boron in seawater by 4%
672!f2py character*3 optional, intent(in) :: optB='l10'
673  CHARACTER(3), OPTIONAL, INTENT(in) :: optB
674  !> for Kf, choose either \b 'pf' (Perez & Fraga, 1987) or \b 'dg' (Dickson & Riley, 1979)
675!f2py character*2 optional, intent(in) :: optKf='pf'
676  CHARACTER(2), OPTIONAL, INTENT(in) :: optKf
677  !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010)
678!f2py character*3 optional, intent(in) :: optK1K2='l'
679  CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2
680  !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction)
681  !! 'Ppot'    - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
682  !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf
683  !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean
684!f2py character*7 optional, intent(in) :: optGAS='Pinsitu'
685  CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS
686
687! Ouput variables
688  !> solubility of CO2 in seawater (Weiss, 1974), also known as K0
689  REAL(kind=wp), INTENT(out), DIMENSION(N) :: K0
690  !> K1 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2
691  REAL(kind=wp), INTENT(out), DIMENSION(N) :: K1
692  !> K2 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) :: K2
694  !> equilibrium constant for dissociation of boric acid
695  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kb
696  !> equilibrium constant for the dissociation of water (Millero, 1995)
697  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kw
698  !> equilibrium constant for the dissociation of bisulfate (Dickson, 1990)
699  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ks
700  !> equilibrium constant for the dissociation of hydrogen fluoride
701  !! either from Dickson and Riley (1979) or from Perez and Fraga (1987), depending on optKf
702  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kf
703  !> solubility product for calcite (Mucci, 1983)
704  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kspc
705  !> solubility product for aragonite (Mucci, 1983)
706  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kspa
707  !> 1st dissociation constant for phosphoric acid (Millero, 1995)
708  REAL(kind=wp), INTENT(out), DIMENSION(N) :: K1p
709  !> 2nd dissociation constant for phosphoric acid (Millero, 1995)
710  REAL(kind=wp), INTENT(out), DIMENSION(N) :: K2p
711  !> 3rd dissociation constant for phosphoric acid (Millero, 1995)
712  REAL(kind=wp), INTENT(out), DIMENSION(N) :: K3p
713  !> equilibrium constant for the dissociation of silicic acid (Millero, 1995)
714  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ksi
715  !> total sulfate (Morris & Riley, 1966)
716  REAL(kind=wp), INTENT(out), DIMENSION(N) :: St
717  !> total fluoride  (Riley, 1965)
718  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ft
719  !> total boron
720  !! from either Uppstrom (1974) or Lee et al. (2010), depending on optB
721  REAL(kind=wp), INTENT(out), DIMENSION(N) :: Bt
722
723! Local variables
724  REAL(kind=wp) :: ssal
725  REAL(kind=wp) :: p
726  REAL(kind=wp) :: tempot, tempis68, tempot68
727  REAL(kind=wp) :: tempis
728  REAL(kind=wp) :: is, invtk, dlogtk, is2, s2, sqrtis
729  REAL(kind=wp) :: Ks_0p, Kf_0p
730  REAL(kind=wp) :: total2free, free2SWS, total2SWS, SWS2total
731  REAL(kind=wp) :: total2free_0p, free2SWS_0p, total2SWS_0p
732! REAL(kind=wp) :: free2SWS, free2SWS_0p
733
734  REAL(kind=wp) :: dtempot, dtempot68
735  REAL(kind=wp) :: R
736
737  REAL(kind=wp) :: pK1o, ma1, mb1, mc1, pK1
738  REAL(kind=wp) :: pK2o, ma2, mb2, mc2, pK2
739
740  REAL(kind=wp), DIMENSION(12) :: a0, a1, a2, b0, b1, b2
741  REAL(kind=wp), DIMENSION(12) :: deltav, deltak, lnkpok0
742  REAL(kind=wp) :: tmp, nK0we74
743
744  INTEGER :: i, icount, ipc
745
746  REAL(kind=wp) :: t, tk, tk0, prb
747  REAL(kind=wp) :: s, sqrts, s15, scl
748
749  REAL(kind=wp) :: Phydro_atm, Patmd, Ptot, Rgas_atm, vbarCO2
750
751! Arrays to pass optional arguments into or use defaults (Dickson et al., 2007)
752  CHARACTER(3) :: opB
753  CHARACTER(2) :: opKf
754  CHARACTER(3) :: opK1K2
755  CHARACTER(7) :: opGAS
756
757  ! CONSTANTS
758  ! =========
759  ! Constants in formulation for Pressure effect on K's (Millero, 95)
760  ! with corrected coefficients for Kb, Kw, Ksi, etc.
761
762  ! index: 1) K1 , 2) K2, 3) Kb, 4) Kw, 5) Ks, 6) Kf, 7) Kspc, 8) Kspa,
763  !            9) K1P, 10) K2P, 11) K3P, 12) Ksi
764
765  DATA a0 /-25.5_wp, -15.82_wp, -29.48_wp, -20.02_wp, &
766          -18.03_wp,  -9.78_wp, -48.76_wp, -45.96_wp, &
767          -14.51_wp, -23.12_wp, -26.57_wp, -29.48_wp/
768  DATA a1 /0.1271_wp, -0.0219_wp, 0.1622_wp, 0.1119_wp, &
769           0.0466_wp, -0.0090_wp, 0.5304_wp, 0.5304_wp, &
770           0.1211_wp, 0.1758_wp, 0.2020_wp, 0.1622_wp/
771  DATA a2 /     0.0_wp,       0.0_wp, -2.608e-3_wp, -1.409e-3_wp, &
772           0.316e-3_wp, -0.942e-3_wp,  0.0_wp,       0.0_wp, &
773          -0.321e-3_wp, -2.647e-3_wp, -3.042e-3_wp, -2.6080e-3_wp/
774  DATA b0 /-3.08e-3_wp, 1.13e-3_wp,  -2.84e-3_wp,   -5.13e-3_wp, &
775           -4.53e-3_wp, -3.91e-3_wp, -11.76e-3_wp, -11.76e-3_wp, &
776           -2.67e-3_wp, -5.15e-3_wp,  -4.08e-3_wp,  -2.84e-3_wp/
777  DATA b1 /0.0877e-3_wp, -0.1475e-3_wp, 0.0_wp,       0.0794e-3_wp, &
778           0.09e-3_wp,    0.054e-3_wp,  0.3692e-3_wp, 0.3692e-3_wp, &
779           0.0427e-3_wp,  0.09e-3_wp,   0.0714e-3_wp, 0.0_wp/
780  DATA b2 /12*0.0_wp/
781
782! Set defaults for optional arguments (in Fortran 90)
783! Note:  Optional arguments with f2py (python) are set above with
784!        the !f2py statements that precede each type declaraion
785  IF (PRESENT(optB)) THEN
786    opB = optB
787  ELSE
788    opB = 'l10'
789  ENDIF
790  IF (PRESENT(optKf)) THEN
791    opKf = optKf
792  ELSE
793    opKf = 'pf'
794  ENDIF
795  IF (PRESENT(optK1K2)) THEN
796    opK1K2 = optK1K2
797  ELSE
798    opK1K2 = 'l'
799  ENDIF
800  IF (PRESENT(optGAS)) THEN
801    opGAS = optGAS
802  ELSE
803    opGAS = 'Pinsitu'
804  ENDIF
805
806  R = 83.14472_wp
807
808  icount = 0
809  DO i = 1, N
810     icount = icount + 1
811!    ===============================================================
812!    Convert model depth -> press; convert model Theta -> T in situ
813!    ===============================================================
814!    * Model temperature tracer is usually "potential temperature"
815!    * Model vertical grid is usually in meters
816!    BUT carbonate chem routines require pressure & in-situ T
817!    Thus before computing chemistry, if appropriate,
818!    convert these 2 model vars (input to this routine)
819!     - depth [m] => convert to pressure [db]
820!     - potential temperature (C) => convert to in-situ T (C)
821!    -------------------------------------------------------
822!    1)  Compute pressure [db] from depth [m] and latitude [degrees] (if input is m, for models)
823     IF (trim(optP) == 'm' ) THEN
824!       Compute pressure [db] from depth [m] and latitude [degrees]
825        p = p80(depth(i), lat(i))
826     ELSEIF (trim(optP) == 'db' ) THEN
827!       In this case (where optP = 'db'), p is input & output (no depth->pressure conversion needed)
828        p = depth(i)
829     ELSE
830        PRINT *,"optP must be 'm' or 'db'"
831        STOP
832     ENDIF
833
834!    2) Convert potential T to in-situ T (if input is Tpot, i.e. case for models):
835     IF (trim(optT) == 'Tpot' .OR. trim(optT) == 'tpot') THEN
836        tempot = temp(i)
837!       This is the case for most models and some data
838!       a) Convert the pot. temp on today's "ITS 90" scale to older IPTS 68 scale
839!          (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote)
840        tempot68 = (tempot - 0.0002) / 0.99975
841!       b) Compute "in-situ Temperature" from "Potential Temperature" (both on IPTS 68)
842        tempis68 = sw_temp(sal(i), tempot68, p, 0.0_wp )
843!       c) Convert the in-situ temp on older IPTS 68 scale to modern scale (ITS 90)
844        tempis = 0.99975*tempis68 + 0.0002
845!       Note: parts (a) and (c) above are tiny corrections;
846!             part  (b) is a big correction for deep waters (but zero at surface)
847     ELSEIF (trim(optT) == 'Tinsitu' .OR. trim(optT) == 'tinsitu') THEN
848!       When optT = 'Tinsitu', tempis is input & output (no tempot needed)
849        tempis    = temp(i)
850        tempis68  = (temp(i) - 0.0002) / 0.99975
851!       dtempot68 = sw_ptmp(DBLE(sal(i)), DBLE(tempis68), DBLE(p), 0.0_wp)
852        dtempot68 = sw_ptmp(sal(i), tempis68, p, 0.0_wp)
853        dtempot   = 0.99975*dtempot68 + 0.0002
854     ELSE
855        PRINT *,"optT must be either 'Tpot' or 'Tinsitu'"
856        PRINT *,"you specified optT =", trim(optT) 
857        STOP
858     ENDIF
859
860!    Compute constants:
861     IF (temp(i) >= -5. .AND. temp(i) < 1.0e+2) THEN
862!       Test to indicate if any of input variables are unreasonable
863        IF (      sal(i) < 0.  .OR.  sal(i) > 1e+3) THEN
864           PRINT *, 'i, icount, temp, sal =', i, icount, temp(i), sal(i)
865        ENDIF
866!       Zero out negative salinity (prev case for OCMIP2 model w/ slightly negative S in some coastal cells)
867        IF (sal(i) < 0.0) THEN
868           ssal = 0.0
869        ELSE
870           ssal = sal(i)
871        ENDIF
872
873!       Absolute temperature (Kelvin) and related values
874        t = DBLE(tempis)
875        tk = 273.15d0 + t
876        invtk=1.0d0/tk
877        dlogtk=LOG(tk)
878
879!       Atmospheric pressure
880        Patmd = DBLE(Patm(i))
881
882!       Hydrostatic pressure (prb is in bars)
883        prb = DBLE(p) / 10.0d0
884
885!       Salinity and simply related values
886        s = DBLE(ssal)
887        s2=s*s
888        sqrts=SQRT(s)
889        s15=s**1.5d0
890        scl=s/1.80655d0
891
892!       Ionic strength:
893        is = 19.924d0*s/(1000.0d0 - 1.005d0*s)
894        is2 = is*is
895        sqrtis = SQRT(is)
896
897!       Total concentrations for sulfate, fluoride, and boron
898
899!       Sulfate: Morris & Riley (1966)
900        St(i) = 0.14d0 * scl/96.062d0
901
902!       Fluoride:  Riley (1965)
903        Ft(i) = 0.000067d0 * scl/18.9984d0
904
905!       Boron:
906        IF (trim(opB) == 'l10') THEN
907!          New formulation from Lee et al (2010)
908           Bt(i) = 0.0002414d0 * scl/10.811d0
909        ELSEIF (trim(opB) == 'u74') THEN
910!          Classic formulation from Uppström (1974)
911           Bt(i) = 0.000232d0  * scl/10.811d0
912        ELSE
913           PRINT *,"optB must be 'l10' or 'u74'"
914           STOP
915        ENDIF
916
917!       K0 (K Henry)
918!       CO2(g) <-> CO2(aq.)
919!       K0  = [CO2]/ fCO2
920!       Weiss (1974)   [mol/kg/atm]
921        IF     (trim(opGAS) == 'Pzero'   .OR. trim(opGAS) == 'pzero') THEN
922           tk0 = tk                   !in situ temperature (K) for K0 calculation
923           Ptot = Patmd               !total pressure (in atm) = atmospheric pressure ONLY
924        ELSEIF (trim(opGAS) == 'Ppot'    .OR. trim(opGAS) == 'ppot') THEN
925           tk0 = dtempot + 273.15d0   !potential temperature (K) for K0 calculation as needed for potential fCO2 & pCO2
926           Ptot = Patmd               !total pressure (in atm) = atmospheric pressure ONLY
927        ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN
928           tk0 = tk                     !in situ temperature (K) for K0 calculation
929           Phydro_atm = prb / 1.01325d0 !convert hydrostatic pressure from bar to atm (1.01325 bar / atm)
930           Ptot = Patmd + Phydro_atm    !total pressure (in atm) = atmospheric pressure + hydrostatic pressure
931        ELSE
932           PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'"
933           STOP
934        ENDIF
935        tmp = 9345.17d0/tk0 - 60.2409d0 + 23.3585d0 * LOG(tk0/100.0d0)
936        nK0we74 = tmp + s*(0.023517d0 - 0.00023656d0*tk0 + 0.0047036e-4_wp*tk0*tk0)
937        K0(i) = EXP(nK0we74)
938
939!       K1 = [H][HCO3]/[H2CO3]
940!       K2 = [H][CO3]/[HCO3]
941        IF (trim(opK1K2) == 'l') THEN
942!         Mehrbach et al. (1973) refit, by Lueker et al. (2000) (total pH scale)
943          K1(i) = 10.0d0**(-1.0d0*(3633.86d0*invtk - 61.2172d0 + 9.6777d0*dlogtk  &
944                  - 0.011555d0*s + 0.0001152d0*s2))
945          K2(i) = 10.0d0**(-1*(471.78d0*invtk + 25.9290d0 - 3.16967d0*dlogtk      &
946                  - 0.01781d0*s + 0.0001122d0*s2))
947        ELSEIF (trim(opK1K2) == 'm10') THEN
948!         Millero (2010, Mar. Fresh Wat. Res.) (total pH scale)
949!         pK1o = 6320.813d0*invtk + 19.568224d0*dlogtk -126.34048d0
950!         ma1 = 13.4051d0*sqrts + 0.03185d0*s - (5.218e-5)*s2
951!         mb1 = -531.095d0*sqrts - 5.7789d0*s
952!         mc1 = -2.0663d0*sqrts
953!         pK1 = pK1o + ma1 + mb1*invtk + mc1*dlogtk
954!         K1(i) = 10.0d0**(-pK1)
955
956!         pK2o = 5143.692d0*invtk + 14.613358d0*dlogtk -90.18333d0
957!         ma2 = 21.5724d0*sqrts + 0.1212d0*s - (3.714e-4)*s2
958!         mb2 = -798.292d0*sqrts - 18.951d0*s
959!         mc2 = -3.403d0*sqrts
960!         pK2 = pK2o + ma2 + mb2*invtk + mc2*dlogtk
961!         K2(i) = 10.0d0**(-pK2)
962
963!         Millero (2010, Mar. Fresh Wat. Res.) (seawater pH scale)
964          pK1o = 6320.813d0*invtk + 19.568224d0*dlogtk -126.34048d0
965          ma1 = 13.4038d0*sqrts + 0.03206d0*s - (5.242e-5)*s2
966          mb1 = -530.659d0*sqrts - 5.8210d0*s
967          mc1 = -2.0664d0*sqrts
968          pK1 = pK1o + ma1 + mb1*invtk + mc1*dlogtk
969          K1(i) = 10.0d0**(-pK1) 
970
971          pK2o = 5143.692d0*invtk + 14.613358d0*dlogtk -90.18333d0
972          ma2 = 21.3728d0*sqrts + 0.1218d0*s - (3.688e-4)*s2
973          mb2 = -788.289d0*sqrts - 19.189d0*s
974          mc2 = -3.374d0*sqrts
975          pK2 = pK2o + ma2 + mb2*invtk + mc2*dlogtk
976          K2(i) = 10.0d0**(-pK2)
977        ELSE
978           PRINT *, "optK1K2 must be either 'l' or 'm10'"
979           STOP
980        ENDIF
981
982!       Kb = [H][BO2]/[HBO2]
983!       (total scale)
984!       Millero p.669 (1995) using data from Dickson (1990)
985        Kb(i) = EXP((-8966.90d0 - 2890.53d0*sqrts - 77.942d0*s +  &
986                1.728d0*s15 - 0.0996d0*s2)*invtk +              &
987                (148.0248d0 + 137.1942d0*sqrts + 1.62142d0*s) +   &
988                (-24.4344d0 - 25.085d0*sqrts - 0.2474d0*s) *      &
989                dlogtk + 0.053105d0*sqrts*tk)
990
991!       K1p = [H][H2PO4]/[H3PO4]
992!       (seawater scale)
993!       DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
994!       Millero (1995), p.670, eq. 65
995!       Use Millero equation's 115.540 constant instead of 115.525 (Dickson et al., 2007).
996!       The latter is only an crude approximation to convert to Total scale (by subtracting 0.015)
997!       And we want to stay on the SWS scale anyway for the pressure correction later.
998        K1p(i) = EXP(-4576.752d0*invtk + 115.540d0 - 18.453d0*dlogtk +  &
999                 (-106.736d0*invtk + 0.69171d0) * sqrts +             &
1000                 (-0.65643d0*invtk - 0.01844d0) * s)
1001
1002!       K2p = [H][HPO4]/[H2PO4]
1003!       (seawater scale)
1004!       DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
1005!       Millero (1995), p.670, eq. 66
1006!       Use Millero equation's 172.1033 constant instead of 172.0833 (Dickson et al., 2007).
1007!       The latter is only an crude approximation to convert to Total scale (by subtracting 0.015)
1008!       And we want to stay on the SWS scale anyway for the pressure correction later.
1009        K2p(i) = EXP(-8814.715d0*invtk + 172.1033d0 - 27.927d0*dlogtk +  &
1010                 (-160.340d0*invtk + 1.3566d0)*sqrts +                 &
1011                 (0.37335d0*invtk - 0.05778d0)*s)
1012
1013!       K3p = [H][PO4]/[HPO4]
1014!       (seawater scale)
1015!       DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
1016!       Millero (1995), p.670, eq. 67
1017!       Use Millero equation's 18.126 constant instead of 18.141 (Dickson et al., 2007).
1018!       The latter is only an crude approximation to convert to Total scale (by subtracting 0.015)
1019!       And we want to stay on the SWS scale anyway for the pressure correction later.
1020        K3p(i) = EXP(-3070.75d0*invtk - 18.126d0 +            &
1021                 (17.27039d0*invtk + 2.81197d0) *             &
1022                 sqrts + (-44.99486d0*invtk - 0.09984d0) * s)
1023
1024!       Ksi = [H][SiO(OH)3]/[Si(OH)4]
1025!       (seawater scale)
1026!       Millero (1995), p.671, eq. 72
1027!       Use Millero equation's 117.400 constant instead of 117.385 (Dickson et al., 2007).
1028!       The latter is only an crude approximation to convert to Total scale (by subtracting 0.015)
1029!       And we want to stay on the SWS scale anyway for the pressure correction later.
1030        Ksi(i) = EXP(-8904.2d0*invtk  + 117.400d0 - 19.334d0*dlogtk +  &
1031                 (-458.79d0*invtk + 3.5913d0) * sqrtis +             &
1032                 (188.74d0*invtk - 1.5998d0) * is +                  &
1033                 (-12.1652d0*invtk + 0.07871d0) * is2 +              &
1034                 LOG(1.0 - 0.001005d0*s))
1035
1036!       Kw = [H][OH]
1037!       (seawater scale)
1038!       Millero (1995) p.670, eq. 63 from composite data
1039!       Use Millero equation's 148.9802 constant instead of 148.9652 (Dickson et al., 2007).
1040!       The latter is only an crude approximation to convert to Total scale (by subtracting 0.015)
1041!       And we want to stay on the SWS scale anyway for the pressure correction later.
1042        Kw(i) = EXP(-13847.26d0*invtk + 148.9802d0 - 23.6521d0*dlogtk +  &
1043               (118.67d0*invtk - 5.977d0 + 1.0495d0 * dlogtk) *          &
1044               sqrts - 0.01615d0 * s)
1045
1046!       Ks = [H][SO4]/[HSO4]
1047!       (free scale)
1048!       Dickson (1990, J. chem. Thermodynamics 22, 113)
1049        Ks_0p = EXP(-4276.1d0*invtk + 141.328d0 - 23.093d0*dlogtk          &
1050                + (-13856.d0*invtk + 324.57d0 - 47.986d0*dlogtk) * sqrtis  &
1051                + (35474.d0*invtk - 771.54 + 114.723d0*dlogtk) * is      &
1052                - 2698.d0*invtk*is**1.5 + 1776.d0*invtk*is2              &
1053                + LOG(1.0d0 - 0.001005d0*s))
1054
1055!       Kf = [H][F]/[HF]
1056!       (total scale)
1057        IF (trim(opKf) == 'dg') THEN
1058!          Dickson and Riley (1979) -- change pH scale to total (following Dickson & Goyet, 1994)
1059           Kf_0p = EXP(1590.2d0*invtk - 12.641d0 + 1.525d0*sqrtis +  &
1060                   LOG(1.0d0 - 0.001005d0*s) +                     &
1061                   LOG(1.0d0 + St(i)/Ks_0p))
1062        ELSEIF (trim(opKf) == 'pf') THEN
1063!          Perez and Fraga (1987) - Already on Total scale (no need for last line above)
1064!          Formulation as given in Dickson et al. (2007)
1065           Kf_0p = EXP(874.d0*invtk - 9.68d0 + 0.111d0*sqrts)
1066        ELSE
1067           PRINT *, "optKf must be either 'dg' or 'pf'"
1068           STOP
1069        ENDIF
1070
1071!       Kspc (calcite) - apparent solubility product of calcite
1072!       (no scale)
1073!       Kspc = [Ca2+] [CO32-] when soln is in equilibrium w/ calcite
1074!       Mucci 1983 mol/kg-soln
1075        Kspc(i) = 10d0**(-171.9065d0 - 0.077993d0*tk + 2839.319d0/tk    &
1076                 + 71.595d0*LOG10(tk)                             &
1077                 + (-0.77712d0 + 0.0028426d0*tk + 178.34d0/tk)*sqrts  &
1078                 -0.07711d0*s + 0.0041249d0*s15 )
1079
1080
1081!       Kspa (aragonite) - apparent solubility product of aragonite
1082!       (no scale)
1083!       Kspa = [Ca2+] [CO32-] when soln is in equilibrium w/ aragonite
1084!       Mucci 1983 mol/kg-soln
1085        Kspa(i) = 10.d0**(-171.945d0 - 0.077993d0*tk + 2903.293d0/tk &
1086             +71.595d0*LOG10(tk) &
1087             +(-0.068393d0 + 0.0017276d0*tk + 88.135d0/tk)*sqrts &
1088             -0.10018d0*s + 0.0059415d0*s15 )
1089
1090!       Pressure effect on K0 based on Weiss (1974, equation 5)
1091        Rgas_atm = 82.05736_wp      ! (cm3 * atm) / (mol * K)  CODATA (2006)
1092        vbarCO2 = 32.3_wp           ! partial molal volume (cm3 / mol) from Weiss (1974, Appendix, paragraph 3)
1093        K0(i) = K0(i) * exp( ((1-Ptot)*vbarCO2)/(Rgas_atm*tk0) )   ! Weiss (1974, equation 5)
1094
1095!       Pressure effect on all other K's (based on Millero, (1995)
1096!           index: K1(1), K2(2), Kb(3), Kw(4), Ks(5), Kf(6), Kspc(7), Kspa(8),
1097!                  K1p(9), K2p(10), K3p(11), Ksi(12)
1098        DO ipc = 1, 12
1099           deltav(ipc)  =  a0(ipc) + a1(ipc) *t + a2(ipc) *t*t
1100           deltak(ipc)   = (b0(ipc)  + b1(ipc) *t + b2(ipc) *t*t)
1101           lnkpok0(ipc)  = (-(deltav(ipc)) &
1102                +(0.5d0*deltak(ipc) * prb) &
1103                )                         * prb/(R*tk)
1104        END DO
1105
1106!       Pressure correction on Ks (Free scale)
1107        Ks(i) = Ks_0p*EXP(lnkpok0(5))
1108!       Conversion factor total -> free scale
1109        total2free     = 1.d0/(1.d0 + St(i)/Ks(i))   ! Kfree = Ktotal*total2free
1110!       Conversion factor total -> free scale at pressure zero
1111        total2free_0p  = 1.d0/(1.d0 + St(i)/Ks_0p)   ! Kfree = Ktotal*total2free
1112
1113!       Pressure correction on Kf
1114!       Kf must be on FREE scale before correction
1115        Kf_0p = Kf_0p * total2free_0p   !Convert from Total to Free scale (pressure 0)
1116        Kf(i) = Kf_0p * EXP(lnkpok0(6)) !Pressure correction (on Free scale)
1117        Kf(i) = Kf(i)/total2free        !Convert back from Free to Total scale
1118
1119!       Convert between seawater and total hydrogen (pH) scales
1120        free2SWS  = 1.d0 + St(i)/Ks(i) + Ft(i)/(Kf(i)*total2free)  ! using Kf on free scale
1121        total2SWS = total2free * free2SWS                          ! KSWS = Ktotal*total2SWS
1122        SWS2total = 1.d0 / total2SWS
1123!       Conversion at pressure zero
1124        free2SWS_0p  = 1.d0 + St(i)/Ks_0p + Ft(i)/(Kf_0p)  ! using Kf on free scale
1125        total2SWS_0p = total2free_0p * free2SWS_0p         ! KSWS = Ktotal*total2SWS
1126
1127!       Convert from Total to Seawater scale before pressure correction
1128!       Must change to SEAWATER scale: K1, K2, Kb
1129        IF (trim(optK1K2) == 'l') THEN
1130          K1(i)  = K1(i)*total2SWS_0p
1131          K2(i)  = K2(i)*total2SWS_0p
1132          !This conversion is unnecessary for the K1,K2 from Millero (2010),
1133          !since we use here the formulation already on the seawater scale
1134        ENDIF
1135        Kb(i)  = Kb(i)*total2SWS_0p
1136
1137!       Already on SEAWATER scale: K1p, K2p, K3p, Kb, Ksi, Kw
1138
1139!       Other contants (keep on another scale):
1140!          - K0         (independent of pH scale, already pressure corrected)
1141!          - Ks         (already on Free scale;   already pressure corrected)
1142!          - Kf         (already on Total scale;  already pressure corrected)
1143!          - Kspc, Kspa (independent of pH scale; pressure-corrected below)
1144
1145!       Perform actual pressure correction (on seawater scale)
1146        K1(i)   = K1(i)*EXP(lnkpok0(1))
1147        K2(i)   = K2(i)*EXP(lnkpok0(2))
1148        Kb(i)   = Kb(i)*EXP(lnkpok0(3))
1149        Kw(i)   = Kw(i)*EXP(lnkpok0(4))
1150        Kspc(i) = Kspc(i)*EXP(lnkpok0(7))
1151        Kspa(i) = Kspa(i)*EXP(lnkpok0(8))
1152        K1p(i)  = K1p(i)*EXP(lnkpok0(9))
1153        K2p(i)  = K2p(i)*EXP(lnkpok0(10))
1154        K3p(i)  = K3p(i)*EXP(lnkpok0(11))
1155        Ksi(i)  = Ksi(i)*EXP(lnkpok0(12))
1156
1157!       Convert back to original total scale:
1158        K1(i)  = K1(i) *SWS2total
1159        K2(i)  = K2(i) *SWS2total
1160        K1p(i) = K1p(i)*SWS2total
1161        K2p(i) = K2p(i)*SWS2total
1162        K3p(i) = K3p(i)*SWS2total
1163        Kb(i)  = Kb(i) *SWS2total
1164        Ksi(i) = Ksi(i)*SWS2total
1165        Kw(i)  = Kw(i) *SWS2total
1166
1167     ELSE
1168
1169        K0(i)   = 1.e20_wp
1170        K1(i)   = 1.e20_wp
1171        K2(i)   = 1.e20_wp
1172        Kb(i)   = 1.e20_wp
1173        Kw(i)   = 1.e20_wp
1174        Ks(i)   = 1.e20_wp
1175        Kf(i)   = 1.e20_wp
1176        Kspc(i) = 1.e20_wp
1177        Kspa(i) = 1.e20_wp
1178        K1p(i)  = 1.e20_wp
1179        K2p(i)  = 1.e20_wp
1180        K3p(i)  = 1.e20_wp
1181        Ksi(i)  = 1.e20_wp
1182        Bt(i)   = 1.e20_wp
1183        Ft(i)   = 1.e20_wp
1184        St(i)   = 1.e20_wp
1185
1186     ENDIF
1187
1188  END DO
1189
1190  RETURN
1191END SUBROUTINE constants
1192
1193! ----------------------------------------------------------------------
1194!  VARSOLVER
1195! ----------------------------------------------------------------------
1196!
1197!> \file varsolver.f90
1198!! \BRIEF
1199!> Module with varsolver subroutine - solve for pH and other carbonate system variables
1200!>    Solve for pH and other carbonate system variables (with input from vars routine)
1201SUBROUTINE varsolver(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC,             &
1202                    temp, salt, ta, tc, pt, sit,                                 &
1203                    Bt, St, Ft,                                                  &
1204                    K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi,  & 
1205                    Patm, Phydro_bar, rhodum, optGAS                             )
1206
1207  !   Purpose: Solve for pH and other carbonate system variables (with input from vars routine)
1208
1209  !     INPUT variables:
1210  !     ================
1211  !     temp    = in situ temperature [degrees C]
1212  !     ta      = total alkalinity                     in [eq/m^3] or [eq/kg]   based on optCON in calling routine (vars)
1213  !     tc      = dissolved inorganic carbon           in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars)
1214  !     pt      = total dissolved inorganic phosphorus in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars)
1215  !     sit     = total dissolved inorganic silicon    in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars)
1216  !     Bt      = total dissolved inorganic boron      computed in calling routine (vars)
1217  !     St      = total dissolved inorganic sulfur     computed in calling routine (vars)
1218  !     Ft      = total dissolved inorganic fluorine   computed in calling routine (vars)
1219  !     K's     = K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi
1220  !     Patm    = atmospheric pressure [atm]
1221  !     Phydro_bar = hydrostatic pressure [bar]
1222  !     rhodum  = density factor as computed in calling routine  (vars)
1223  !     -----------
1224  !     optGAS: choose in situ vs. potential fCO2 and pCO2 (default optGAS = 'Pinsitu')
1225  !     ---------
1226  !       PRESSURE & T corrections for K0 and the fugacity coefficient (Cf)
1227  !       -> 'Pzero'   = 'zero order' fCO2 and pCO2 (typical approach, which is flawed)
1228  !                      considers in situ T & only atm pressure (hydrostatic=0)
1229  !       -> 'Ppot'    = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
1230  !                      considers potential T & only atm pressure (hydrostatic press = 0)
1231  !       -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure)
1232  !                      considers in situ T & total pressure (atm + hydrostatic)
1233  !     ---------
1234
1235  !     OUTPUT variables:
1236  !     =================
1237  !     ph   = pH on total scale
1238  !     pco2 = CO2 partial pressure (uatm)
1239  !     fco2 = CO2 fugacity (uatm)
1240  !     co2  = aqueous CO2 concentration in [mol/kg] or [mol/m^3] determined by rhodum (depends on optCON in calling routine)
1241  !     hco3 = bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] determined by rhodum
1242  !     co3  = carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] determined by rhodum
1243  !     OmegaA = Omega for aragonite, i.e., the aragonite saturation state
1244  !     OmegaC = Omega for calcite, i.e., the   calcite saturation state
1245
1246  USE mocsy_singledouble
1247  USE mocsy_phsolvers
1248
1249  IMPLICIT NONE
1250
1251! Input variables
1252  !> <b>in situ temperature</b> [degrees C]
1253  REAL(kind=wp), INTENT(in) :: temp
1254  !> <b>salinity</b> [on the practical salinity scale, dimensionless]
1255  REAL(kind=wp), INTENT(in) :: salt
1256  !> total alkalinity in <b>[eq/m^3]</b> OR in <b>[eq/kg]</b>, depending on optCON in calling routine
1257  REAL(kind=wp), INTENT(in) :: ta
1258  !> dissolved inorganic carbon in <b>[mol/m^3]</b> OR in <b>[mol/kg]</b>, depending on optCON in calling routine
1259  REAL(kind=wp), INTENT(in) :: tc
1260  !> phosphate concentration in <b>[mol/m^3]</b> OR in <b>[mol/kg]</b>, depending on optCON in calling routine
1261  REAL(kind=wp), INTENT(in) :: pt
1262  !> total dissolved inorganic silicon 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) :: sit
1264  !> total boron from either Uppstrom (1974) or Lee et al. (2010), depending on optB in calling routine
1265  REAL(kind=wp), INTENT(in) :: Bt
1266  !> total sulfate (Morris & Riley, 1966)
1267  REAL(kind=wp), INTENT(in) :: St
1268  !> total fluoride  (Riley, 1965)
1269  REAL(kind=wp), INTENT(in) :: Ft
1270  !> solubility of CO2 in seawater (Weiss, 1974), also known as K0
1271  REAL(kind=wp), INTENT(in) :: K0
1272  !> K1 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2
1273  REAL(kind=wp), INTENT(in) :: K1
1274  !> K2 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2
1275  REAL(kind=wp), INTENT(in) :: K2
1276  !> equilibrium constant for dissociation of boric acid
1277  REAL(kind=wp), INTENT(in) :: Kb
1278  !> equilibrium constant for the dissociation of water (Millero, 1995)
1279  REAL(kind=wp), INTENT(in) :: Kw
1280  !> equilibrium constant for the dissociation of bisulfate (Dickson, 1990)
1281  REAL(kind=wp), INTENT(in) :: Ks
1282  !> equilibrium constant for the dissociation of hydrogen fluoride
1283  !! from Dickson and Riley (1979) or Perez and Fraga (1987), depending on optKf
1284  REAL(kind=wp), INTENT(in) :: Kf
1285  !> solubility product for calcite (Mucci, 1983)
1286  REAL(kind=wp), INTENT(in) :: Kspc
1287  !> solubility product for aragonite (Mucci, 1983)
1288  REAL(kind=wp), INTENT(in) :: Kspa
1289  !> 1st dissociation constant for phosphoric acid (Millero, 1995)
1290  REAL(kind=wp), INTENT(in) :: K1p
1291  !> 2nd dissociation constant for phosphoric acid (Millero, 1995)
1292  REAL(kind=wp), INTENT(in) :: K2p
1293  !> 3rd dissociation constant for phosphoric acid (Millero, 1995)
1294  REAL(kind=wp), INTENT(in) :: K3p
1295  !> equilibrium constant for the dissociation of silicic acid (Millero, 1995)
1296  REAL(kind=wp), INTENT(in) :: Ksi
1297  !> total atmospheric pressure <b>[atm]</b>
1298  REAL(kind=wp), INTENT(in) :: Patm
1299  !> total hydrostatic pressure <b>[bar]</b>
1300  REAL(kind=wp), INTENT(in) :: Phydro_bar
1301  !> density factor as computed incalling routine  (vars)
1302  REAL(kind=wp), INTENT(in) :: rhodum
1303  !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction)
1304  !! 'Ppot'    - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
1305  !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf
1306  !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean
1307!f2py character*7 optional, intent(in) :: optGAS='Pinsitu'
1308  CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS
1309
1310! Output variables:
1311  !> pH on the <b>total scale</b>
1312  REAL(kind=wp), INTENT(out) :: ph
1313  !> CO2 partial pressure <b>[uatm]</b>
1314  REAL(kind=wp), INTENT(out) :: pco2
1315  !> CO2 fugacity <b>[uatm]</b>
1316  REAL(kind=wp), INTENT(out) :: fco2
1317  !> aqueous CO2* concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg</b>] depending on choice for optCON
1318  REAL(kind=wp), INTENT(out) :: co2
1319  !> bicarbonate ion (HCO3-) concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg]</b> depending on choice for optCON
1320  REAL(kind=wp), INTENT(out) :: hco3
1321  !> carbonate ion (CO3--) concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg]</b> depending on choice for optCON
1322  REAL(kind=wp), INTENT(out) :: co3
1323  !> Omega for aragonite, i.e., the aragonite saturation state
1324  REAL(kind=wp), INTENT(out) :: OmegaA
1325  !> Omega for calcite, i.e., the calcite saturation state
1326  REAL(kind=wp), INTENT(out) :: OmegaC
1327
1328! Local variables
1329  REAL(kind=wp) :: Phydro_atm, Ptot
1330  REAL(kind=wp) :: Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff
1331  REAL(kind=wp) :: tk, tk0
1332  real(kind=wp) :: temp68, tempot, tempot68
1333  REAL(kind=wp) :: Hinit, H
1334  REAL(kind=wp) :: HSO4, HF, HSI, HPO4
1335  REAL(kind=wp) :: ab, aw, ac
1336  REAL(kind=wp) :: cu, cb, cc
1337  REAL(kind=wp) :: Ca
1338! Array to pass optional arguments
1339  CHARACTER(7) :: opGAS
1340
1341  IF (PRESENT(optGAS)) THEN
1342    opGAS = optGAS
1343  ELSE
1344    opGAS = 'Pinsitu'
1345  ENDIF
1346
1347! Compute pH from constants and total concentrations
1348! - use SolveSAPHE v1.0.1 routines from Munhoven (2013, GMD) modified to use mocsy's Ks instead of its own
1349! 1) Compute best starting point for H+ calculation
1350  call ahini_for_at(ta, tc, Bt, K1, K2, Kb, Hinit)
1351! 2) Solve for H+ using above result as the initial H+ value
1352  H = solve_at_general(ta, tc, Bt,                                         & 
1353                       pt,     sit,                                        &
1354                       St, Ft,                                             &
1355                       K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi,     &
1356                       Hinit)
1357! 3) Calculate pH from from H+ concentration (mol/kg)
1358  IF (H > 0.d0) THEN
1359     pH = -1.*LOG10(H)
1360  ELSE
1361     pH = 1.e20_wp
1362  ENDIF
1363
1364! Compute carbonate Alk (Ac) by difference: from total Alk and other Alk components
1365  HSO4 = St/(1.0d0 + Ks/(H/(1.0d0 + St/Ks)))
1366  HF = 1.0d0/(1.0d0 + Kf/H)
1367  HSI = 1.0d0/(1.0d0 + H/Ksi)
1368  HPO4 = (K1p*K2p*(H + 2.*K3p) - H**3) /                &
1369  (H**3 + K1p*H**2 + K1p*K2p*H + K1p*K2p*K3p)
1370  ab = Bt/(1.0d0 + H/Kb)
1371  aw = Kw/H - H/(1.0d0 + St/Ks)
1372  ac = ta + hso4 - sit*hsi - ab - aw + Ft*hf - pt*hpo4
1373
1374! Calculate CO2*, HCO3-, & CO32- (in mol/kg soln) from Ct, Ac, H+, K1, & K2
1375  cu = (2.0d0 * tc - ac) / (2.0d0 + K1 / H)
1376  cb = K1 * cu / H
1377  cc = K2 * cb / H
1378
1379! When optCON = 'mol/m3' in calling routine (vars), then:
1380! convert output var concentrations from mol/kg to mol/m^3
1381! e.g., for case when drho = 1028, multiply by [1.028 kg/L  x  1000 L/m^3])
1382  co2  = cu * rhodum
1383  hco3 = cb * rhodum
1384  co3  = cc * rhodum
1385
1386! Determine CO2 fugacity [uatm]
1387! NOTE: equation just below requires CO2* in mol/kg
1388  fCO2 = cu * 1.e6_wp/K0
1389
1390! Determine CO2 partial pressure from CO2 fugacity [uatm]
1391  tk = 273.15d0 + temp
1392  !Compute EITHER "potential pCO2" OR "in situ pCO2" (T and P used for calculations will differ)
1393  IF     (trim(opGAS) == 'Pzero'   .OR. trim(opGAS) == 'pzero') THEN
1394     tk0 = tk                 !in situ temperature (K) for K0 calculation
1395     Ptot = Patm              !total pressure (in atm) = atmospheric pressure ONLY
1396  ELSEIF (trim(opGAS) == 'Ppot' .OR. trim(opGAS) == 'ppot') THEN
1397     !Use potential temperature and atmospheric pressure (water parcel adiabatically brought back to surface)
1398     !temp68 = (temp - 0.0002d0) / 0.99975d0          !temp = in situ T; temp68 is same converted to ITPS-68 scale
1399     !tempot68 = sw_ptmp(salt, temp68, Phydro_bar*10d0, 0.0d0) !potential temperature (C)
1400     !tempot   = 0.99975*tempot68 + 0.0002
1401     !tk0 = tempot + 273.15d0  !potential temperature (K) for fugacity coeff. calc as needed for potential fCO2 & pCO2
1402     tempot = sw_ptmp(salt, temp, Phydro_bar*10._wp, 0.0_wp) !potential temperature (C)
1403     tk0 = tempot + 273.15d0  !potential temperature (K) for fugacity coeff. calc as needed for potential fCO2 & pCO2
1404     Ptot = Patm              !total pressure (in atm) = atmospheric pressure ONLY
1405  ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN
1406     !Use in situ temperature and total pressure
1407     tk0 = tk                             !in situ temperature (K) for fugacity coefficient calculation
1408     Phydro_atm = Phydro_bar / 1.01325d0  !convert hydrostatic pressure from bar to atm (1.01325 bar / atm)
1409     Ptot = Patm + Phydro_atm            !total pressure (in atm) = atmospheric pressure + hydrostatic pressure
1410  ELSE
1411     PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'"
1412     STOP
1413  ENDIF
1414
1415! Now that we have T and P in the right form, continue with calculation of fugacity coefficient (and pCO2)
1416  Rgas_atm = 82.05736_wp      ! (cm3 * atm) / (mol * K)  CODATA (2006)
1417! To compute fugcoeff, we need 3 other terms (B, Del, xc2) in addition to 3 others above (tk, Ptot, Rgas_atm)
1418  B = -1636.75d0 + 12.0408d0*tk0 - 0.0327957d0*(tk0*tk0) + 0.0000316528d0*(tk0*tk0*tk0)
1419  Del = 57.7d0 - 0.118d0*tk0
1420! "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9
1421! x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite)
1422! Let's assume that xCO2 = fCO2. Resulting fugcoeff is identical to 8th digit after the decimal.
1423  xCO2approx = fCO2 * 1.e-6_wp
1424  IF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN
1425!    xCO2approx = 400.0e-6_wp      !a simple test (gives about same result as seacarb for pCO2insitu)
1426!    approximate surface xCO2 ~ surface fCO2 (i.e., in situ fCO2 d by exponential pressure correction)
1427     xCO2approx = xCO2approx * exp( ((1-Ptot)*32.3_wp)/(82.05736_wp*tk0) )   ! of K0 press. correction, see Weiss (1974, equation 5)
1428  ENDIF
1429  xc2 = (1.0d0 - xCO2approx)**2 
1430  fugcoeff = exp( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk0) )
1431  pCO2 = fCO2 / fugcoeff
1432
1433! Determine Omega Calcite et Aragonite
1434! OmegaA = ((0.01028d0*salt/35.0d0)*cc) / Kspa
1435! OmegaC = ((0.01028d0*salt/35.0d0)*cc) / Kspc
1436! - see comments from Munhoven on the best value "0.02128" which differs slightly from the best practices guide (0.02127)
1437  Ca = (0.02128d0/40.078d0) * salt/1.80655d0
1438  OmegaA = (Ca*cc) / Kspa
1439  OmegaC = (Ca*cc) / Kspc
1440
1441  RETURN
1442END SUBROUTINE varsolver
1443
1444! ----------------------------------------------------------------------
1445!  VARS
1446! ----------------------------------------------------------------------
1447!
1448!> \file vars.f90
1449!! \BRIEF
1450!> Module with vars subroutine - compute carbonate system vars from DIC,Alk,T,S,P,nuts
1451!>    Computes standard carbonate system variables (pH, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R)
1452!!    as 1D arrays FROM
1453!!    temperature, salinity, pressure,
1454!!    total alkalinity (ALK), dissolved inorganic carbon (DIC),
1455!!    silica and phosphate concentrations (all 1-D arrays)
1456SUBROUTINE vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis,  &
1457                temp, sal, alk, dic, sil, phos, Patm, depth, lat, N,                      &
1458                optCON, optT, optP, optB, optK1K2, optKf, optGAS                          )
1459
1460  !   Purpose:
1461  !     Computes other standard carbonate system variables (pH, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R)
1462  !     as 1D arrays
1463  !     FROM:
1464  !     temperature, salinity, pressure,
1465  !     total alkalinity (ALK), dissolved inorganic carbon (DIC),
1466  !     silica and phosphate concentrations (all 1-D arrays)
1467
1468  !     INPUT variables:
1469  !     ================
1470  !     Patm    = atmospheric pressure [atm]
1471  !     depth   = depth [m]     (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure)
1472  !             = pressure [db] (with optP='db')
1473  !     lat     = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m')
1474  !             = dummy array (unused when optP='db')
1475  !     temp    = potential temperature [degrees C] (with optT='Tpot', i.e., models carry tempot, not in situ temp)
1476  !             = in situ   temperature [degrees C] (with optT='Tinsitu', e.g., for data)
1477  !     sal     = salinity in [psu]
1478  !     alk     = total alkalinity in [eq/m^3] with optCON = 'mol/m3'
1479  !             =               [eq/kg]  with optCON = 'mol/kg'
1480  !     dic     = dissolved inorganic carbon [mol/m^3] with optCON = 'mol/m3'
1481  !             =                            [mol/kg]  with optCON = 'mol/kg'
1482  !     sil     = silica    [mol/m^3] with optCON = 'mol/m3'
1483  !             =           [mol/kg]  with optCON = 'mol/kg'
1484  !     phos    = phosphate [mol/m^3] with optCON = 'mol/m3'
1485  !             =           [mol/kg]  with optCON = 'mol/kg'
1486  !     INPUT options:
1487  !     ==============
1488  !     -----------
1489  !     optCON: choose input & output concentration units - mol/kg (data) vs. mol/m^3 (models)
1490  !     -----------
1491  !       -> 'mol/kg' for DIC, ALK, sil, & phos given on mokal scale, i.e., in mol/kg  (std DATA units)
1492  !       -> 'mol/m3' for DIC, ALK, sil, & phos given in mol/m^3 (std MODEL units)
1493  !     -----------
1494  !     optT: choose in situ vs. potential temperature as input
1495  !     ---------
1496  !     NOTE: Carbonate chem calculations require IN-SITU temperature (not potential Temperature)
1497  !       -> 'Tpot' means input is pot. Temperature (in situ Temp "tempis" is computed)
1498  !       -> 'Tinsitu' means input is already in-situ Temperature, not pot. Temp ("tempis" not computed)
1499  !     ---------
1500  !     optP: choose depth (m) vs pressure (db) as input
1501  !     ---------
1502  !       -> 'm'  means "depth" input is in "m" (thus in situ Pressure "p" [db] is computed)
1503  !       -> 'db' means "depth" input is already in situ pressure [db], not m (thus p = depth)
1504  !     ---------
1505  !     optB: choose total boron formulation - Uppström (1974) vs. Lee et al. (2010)
1506  !     ---------
1507  !       -> 'u74' means use classic formulation of Uppström (1974) for total Boron
1508  !       -> 'l10' means use newer   formulation of Lee et al. (2010) for total Boron
1509  !     ---------
1510  !     optK1K2:
1511  !     ---------
1512  !       -> 'l'   means use Lueker et al. (2000) formulations for K1 & K2 (recommended by Dickson et al. 2007)
1513  !                **** BUT this should only be used when 2 < T < 35 and 19 < S < 43
1514  !       -> 'm10' means use Millero (2010) formulation for K1 & K2 (see Dickson et al., 2007)
1515  !                **** Valid for 0 < T < 50°C and 1 < S < 50 psu
1516  !     ----------
1517  !     optKf:
1518  !     ----------
1519  !       -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007)
1520  !               **** BUT Valid for  9 < T < 33°C and 10 < S < 40.
1521  !       -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994)
1522  !     -----------
1523  !     optGAS: choose in situ vs. potential fCO2 and pCO2
1524  !     ---------
1525  !       PRESSURE corrections for K0 and the fugacity coefficient (Cf)
1526  !       -> 'Pzero'   = 'zero order' fCO2 and pCO2 (typical approach, which is flawed)
1527  !                      considers in situ T & only atm pressure (hydrostatic=0)
1528  !       -> 'Ppot'    = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
1529  !                      considers potential T & only atm pressure (hydrostatic press = 0)
1530  !       -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure)
1531  !                      considers in situ T & total pressure (atm + hydrostatic)
1532  !     ---------
1533
1534  !     OUTPUT variables:
1535  !     =================
1536  !     ph   = pH on total scale
1537  !     pco2 = CO2 partial pressure (uatm)
1538  !     fco2 = CO2 fugacity (uatm)
1539  !     co2  = aqueous CO2 concentration in [mol/kg] or [mol/m^3] depending on optCON
1540  !     hco3 = bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] depending on optCON
1541  !     co3  = carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] depending on optCON
1542  !     OmegaA = Omega for aragonite, i.e., the aragonite saturation state
1543  !     OmegaC = Omega for calcite, i.e., the   calcite saturation state
1544  !     BetaD = Revelle factor   dpCO2/pCO2 / dDIC/DIC
1545  !     rhoSW  = in-situ density of seawater; rhoSW = f(s, t, p)
1546  !     p = pressure [decibars]; p = f(depth, latitude) if computed from depth [m] OR p = depth if [db]
1547  !     tempis  = in-situ temperature [degrees C]
1548
1549  USE mocsy_singledouble
1550
1551  IMPLICIT NONE
1552
1553! Input variables
1554  !>     number of records
1555  INTEGER, INTENT(in) :: N
1556  !> either <b>in situ temperature</b> (when optT='Tinsitu', typical data)
1557  !! OR <b>potential temperature</b> (when optT='Tpot', typical models) <b>[degree C]</b>
1558  REAL(kind=wp), INTENT(in),    DIMENSION(N) :: temp
1559  !> salinity <b>[psu]</b>
1560  REAL(kind=wp), INTENT(in), DIMENSION(N) :: sal
1561  !> total alkalinity in <b>[eq/m^3]</b> (when optCON = 'mol/m3') OR in <b>[eq/kg]</b>  (when optCON = 'mol/kg')
1562  REAL(kind=wp), INTENT(in), DIMENSION(N) :: alk
1563  !> dissolved inorganic carbon in <b>[mol/m^3]</b> (when optCON = 'mol/m3') OR in <b>[mol/kg]</b> (when optCON = 'mol/kg')
1564  REAL(kind=wp), INTENT(in), DIMENSION(N) :: dic
1565  !> SiO2 concentration 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) :: sil
1567  !> phosphate 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) :: phos
1569!f2py optional , depend(sal) :: n=len(sal)
1570  !> atmospheric pressure <b>[atm]</b>
1571  REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm
1572  !> depth in \b meters (when optP='m') or \b decibars (when optP='db')
1573  REAL(kind=wp), INTENT(in),    DIMENSION(N) :: depth
1574  !> latitude <b>[degrees north]</b>
1575  REAL(kind=wp), INTENT(in),    DIMENSION(N) :: lat
1576
1577  !> choose either \b 'mol/kg' (std DATA units) or \b 'mol/m3' (std MODEL units) to select
1578  !! concentration units for input (for alk, dic, sil, phos) & output (co2, hco3, co3)
1579  CHARACTER(6), INTENT(in) :: optCON
1580  !> choose \b 'Tinsitu' for in situ temperature or \b 'Tpot' for potential temperature (in situ Temp is computed, needed for models)
1581  CHARACTER(7), INTENT(in) :: optT
1582  !> for depth input, choose \b "db" for decibars (in situ pressure) or \b "m" for meters (pressure is computed, needed for models)
1583  CHARACTER(2), INTENT(in) :: optP
1584  !> for total boron, choose either \b 'u74' (Uppstrom, 1974) or \b 'l10' (Lee et al., 2010).
1585  !! The 'l10' formulation is based on 139 measurements (instead of 20),
1586  !! uses a more accurate method, and
1587  !! generally increases total boron in seawater by 4%
1588!f2py character*3 optional, intent(in) :: optB='l10'
1589  CHARACTER(3), OPTIONAL, INTENT(in) :: optB
1590  !> for Kf, choose either \b 'pf' (Perez & Fraga, 1987) or \b 'dg' (Dickson & Riley, 1979)
1591!f2py character*2 optional, intent(in) :: optKf='pf'
1592  CHARACTER(2), OPTIONAL, INTENT(in) :: optKf
1593  !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010)
1594!f2py character*3 optional, intent(in) :: optK1K2='l'
1595  CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2
1596  !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction)
1597  !! 'Ppot'    - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface)
1598  !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf
1599  !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean
1600!f2py character*7 optional, intent(in) :: optGAS='Pinsitu'
1601  CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS
1602
1603! Output variables:
1604  !> pH on the <b>total scale</b>
1605  REAL(kind=wp), INTENT(out), DIMENSION(N) :: ph
1606  !> CO2 partial pressure <b>[uatm]</b>
1607  REAL(kind=wp), INTENT(out), DIMENSION(N) :: pco2
1608  !> CO2 fugacity <b>[uatm]</b>
1609  REAL(kind=wp), INTENT(out), DIMENSION(N) :: fco2
1610  !> aqueous CO2* concentration, either in <b>[mol/m^3]</b> or <b>[mol/kg</b>] depending on choice for optCON
1611  REAL(kind=wp), INTENT(out), DIMENSION(N) :: co2
1612  !> bicarbonate ion (HCO3-) 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) :: hco3
1614  !> carbonate ion (CO3--) 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) :: co3
1616  !> Omega for aragonite, i.e., the aragonite saturation state
1617  REAL(kind=wp), INTENT(out), DIMENSION(N) :: OmegaA
1618  !> Omega for calcite, i.e., the calcite saturation state
1619  REAL(kind=wp), INTENT(out), DIMENSION(N) :: OmegaC
1620  !> Revelle factor, i.e., dpCO2/pCO2 / dDIC/DIC
1621  REAL(kind=wp), INTENT(out), DIMENSION(N) :: BetaD
1622  !> in-situ density of seawater; rhoSW = f(s, t, p) in <b>[kg/m3]</b>
1623  REAL(kind=wp), INTENT(out), DIMENSION(N) :: rhoSW
1624  !> pressure <b>[decibars]</b>; p = f(depth, latitude) if computed from depth [m] (when optP='m') OR p = depth [db] (when optP='db')
1625  REAL(kind=wp), INTENT(out), DIMENSION(N) :: p
1626  !> in-situ temperature \b <b>[degrees C]</b>
1627  REAL(kind=wp), INTENT(out), DIMENSION(N) :: tempis
1628
1629! Local variables
1630  REAL(kind=wp) :: ssal, salk, sdic, ssil, sphos
1631
1632  REAL(kind=wp) :: tempot, tempis68, tempot68
1633! REAL(kind=wp) :: dtempot, dtempot68
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
1927END 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
1937SUBROUTINE 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!f2py optional , depend(pCO2) :: n=len(pCO2)
1957
1958! OUTPUT variables:
1959  !> fugacity of CO2 [uatm]
1960  REAL(kind=wp), INTENT(out), DIMENSION(N) :: fCO2
1961
1962! LOCAL variables:
1963  REAL(kind=wp) :: dpCO2, dtemp, tk, dPatm, prb
1964  REAL(kind=wp) :: Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff
1965  REAL(kind=wp) :: dfCO2
1966
1967  INTEGER :: i
1968
1969! REAL(kind=wp) :: sw_ptmp
1970! EXTERNAL sw_ptmp
1971
1972  DO i = 1,N
1973     dpCO2     = DBLE(pCO2(i))
1974     dtemp     = DBLE(temp(i))
1975     dPatm     = DBLE(Patm(i))
1976     tk = 273.15d0 + DBLE(temp(i))     !Absolute temperature (Kelvin)
1977     prb = DBLE(p(i)) / 10.0d0         !Pressure effect (prb is in bars)
1978     Ptot = dPatm + prb/1.01325d0      !Total pressure (atmospheric + hydrostatic) [atm]
1979     Rgas_atm = 82.05736_wp            !R in (cm3 * atm) / (mol * K)  from CODATA (2006)
1980!    To compute fugcoeff, we need 3 other terms (B, Del, xc2) as well as 3 others above (tk, Ptot, Rgas_atm)
1981     B = -1636.75d0 + 12.0408d0*tk - 0.0327957d0*(tk*tk) + 0.0000316528d0*(tk*tk*tk)
1982     Del = 57.7d0 - 0.118d0*tk
1983!    "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9
1984!    x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite)
1985!    Let's assume that xCO2 = pCO2. Resulting fugcoeff is identical to 8th digit after the decimal.
1986     xCO2approx = dpCO2 * 1.e-6_wp
1987     xc2 = (1.0d0 - xCO2approx)**2 
1988     fugcoeff = EXP( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk) )
1989     dfCO2 = dpCO2 * fugcoeff
1990     fCO2(i) = REAL(dfCO2)
1991  END DO
1992
1993  RETURN
1994END SUBROUTINE p2fCO2
1995
1996! ----------------------------------------------------------------------
1997!  P2FCO2
1998! ----------------------------------------------------------------------
1999!
2000!> \file f2pCO2.f90
2001!! \BRIEF
2002!>    Module with f2pCO2 subroutine - compute pCO2 from fCO2, in situ T, atm pressure, hydrostatic pressure
2003!>    Compute pCO2 from arrays of fCO2, in situ temp, atm pressure, & hydrostatic pressure
2004SUBROUTINE f2pCO2(fCO2, temp, Patm, p, N, pCO2)
2005  !    Purpose:
2006  !    Compute pCO2 from arrays of fCO2, in situ temp, atm pressure, & hydrostatic pressure
2007
2008  USE mocsy_singledouble
2009  IMPLICIT NONE
2010
2011  !> number of records
2012  INTEGER, intent(in) :: N
2013
2014! INPUT variables
2015  !> oceanic fugacity of CO2 [uatm]
2016  REAL(kind=wp), INTENT(in), DIMENSION(N) :: fCO2
2017  !> in situ temperature [C]
2018  REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp
2019  !> atmospheric pressure [atm]
2020  REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm
2021  !> hydrostatic pressure [db]
2022  REAL(kind=wp), INTENT(in), DIMENSION(N) :: p
2023!f2py optional , depend(pCO2) :: n=len(pCO2)
2024
2025! OUTPUT variables:
2026  !> oceanic partial pressure of CO2 [uatm]
2027  REAL(kind=wp), INTENT(out), DIMENSION(N) :: pCO2
2028
2029! LOCAL variables:
2030  REAL(kind=wp) :: dfCO2, dtemp, tk, dPatm, prb
2031  REAL(kind=wp) :: Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff
2032  REAL(kind=wp) :: dpCO2
2033
2034  INTEGER :: i
2035
2036! REAL(kind=wp) :: sw_ptmp
2037! EXTERNAL sw_ptmp
2038
2039  DO i = 1,N
2040     dfCO2     = DBLE(fCO2(i))
2041     dtemp     = DBLE(temp(i))
2042     dPatm     = DBLE(Patm(i))
2043     tk = 273.15d0 + DBLE(temp(i))     !Absolute temperature (Kelvin)
2044     prb = DBLE(p(i)) / 10.0d0         !Pressure effect (prb is in bars)
2045     Ptot = dPatm + prb/1.01325d0      !Total pressure (atmospheric + hydrostatic) [atm]
2046     Rgas_atm = 82.05736_wp            !R in (cm3 * atm) / (mol * K)  from CODATA (2006)
2047!    To compute fugcoeff, we need 3 other terms (B, Del, xc2) as well as 3 others above (tk, Ptot, Rgas_atm)
2048     B = -1636.75d0 + 12.0408d0*tk - 0.0327957d0*(tk*tk) + 0.0000316528d0*(tk*tk*tk)
2049     Del = 57.7d0 - 0.118d0*tk
2050!    "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9
2051!    x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite)
2052!    Let's assume that xCO2 = fCO2. Resulting fugcoeff is identical to 8th digit after the decimal.
2053     xCO2approx = dfCO2 * 1.e-6_wp
2054     xc2 = (1.0d0 - xCO2approx)**2 
2055     fugcoeff = exp( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk) )
2056     dpCO2 = dfCO2 / fugcoeff
2057     pCO2(i) = REAL(dpCO2)
2058  END DO
2059
2060  RETURN
2061END SUBROUTINE f2pCO2
2062
2063END MODULE mocsy_mainmod
Note: See TracBrowser for help on using the repository browser.