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

Last change on this file since 7894 was 7894, checked in by jpalmier, 4 years ago

JPALM — 11-04-2017 — MEDUSA spring tidy-up refreshning session

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