MODULE mocsy_mainmod USE in_out_manager ! I/O manager CONTAINS ! ---------------------------------------------------------------------- ! SW_ADTG ! ---------------------------------------------------------------------- ! !> \file sw_adtg.f90 !! \BRIEF !> Module with sw_adtg function - compute adiabatic temp. gradient from S,T,P !> Function to calculate adiabatic temperature gradient as per UNESCO 1983 routines. FUNCTION sw_adtg (s,t,p) ! ================================================================== ! Calculates adiabatic temperature gradient as per UNESCO 1983 routines. ! Armin Koehl akoehl@ucsd.edu ! ================================================================== USE mocsy_singledouble IMPLICIT NONE !> salinity [psu (PSU-78)] REAL(kind=wp) :: s !> temperature [degree C (IPTS-68)] REAL(kind=wp) :: t !> pressure [db] REAL(kind=wp) :: p REAL(kind=wp) :: a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2 REAL(kind=wp) :: sref REAL(kind=wp) :: sw_adtg sref = 35.d0 a0 = 3.5803d-5 a1 = +8.5258d-6 a2 = -6.836d-8 a3 = 6.6228d-10 b0 = +1.8932d-6 b1 = -4.2393d-8 c0 = +1.8741d-8 c1 = -6.7795d-10 c2 = +8.733d-12 c3 = -5.4481d-14 d0 = -1.1351d-10 d1 = 2.7759d-12 e0 = -4.6206d-13 e1 = +1.8676d-14 e2 = -2.1687d-16 sw_adtg = a0 + (a1 + (a2 + a3*T)*T)*T & + (b0 + b1*T)*(S-sref) & + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P & + ( e0 + (e1 + e2*T)*T )*P*P END FUNCTION sw_adtg ! ---------------------------------------------------------------------- ! SW_ADTG ! ---------------------------------------------------------------------- ! !> \file sw_ptmp.f90 !! \BRIEF !> Module with sw_ptmp function - compute potential T from in-situ T !> Function to calculate potential temperature [C] from in-situ temperature FUNCTION sw_ptmp (s,t,p,pr) ! ================================================================== ! Calculates potential temperature [C] from in-situ Temperature [C] ! From UNESCO 1983 report. ! Armin Koehl akoehl@ucsd.edu ! ================================================================== ! Input arguments: ! ------------------------------------- ! s = salinity [psu (PSS-78) ] ! t = temperature [degree C (IPTS-68)] ! p = pressure [db] ! pr = reference pressure [db] USE mocsy_singledouble IMPLICIT NONE ! Input arguments !> salinity [psu (PSS-78)] REAL(kind=wp) :: s !> temperature [degree C (IPTS-68)] REAL(kind=wp) :: t !> pressure [db] REAL(kind=wp) :: p !> reference pressure [db] REAL(kind=wp) :: pr ! local arguments REAL(kind=wp) :: del_P ,del_th, th, q REAL(kind=wp) :: onehalf, two, three PARAMETER (onehalf = 0.5d0, two = 2.d0, three = 3.d0 ) ! REAL(kind=wp) :: sw_adtg ! EXTERNAL sw_adtg ! Output REAL(kind=wp) :: sw_ptmp ! theta1 del_P = PR - P del_th = del_P*sw_adtg(S,T,P) th = T + onehalf*del_th q = del_th ! theta2 del_th = del_P*sw_adtg(S,th,P+onehalf*del_P) th = th + (1.d0 - 1.d0/SQRT(two))*(del_th - q) q = (two-SQRT(two))*del_th + (-two+three/SQRT(two))*q ! theta3 del_th = del_P*sw_adtg(S,th,P+onehalf*del_P) th = th + (1.d0 + 1.d0/SQRT(two))*(del_th - q) q = (two + SQRT(two))*del_th + (-two-three/SQRT(two))*q ! theta4 del_th = del_P*sw_adtg(S,th,P+del_P) sw_ptmp = th + (del_th - two*q)/(two*three) RETURN END FUNCTION sw_ptmp ! ---------------------------------------------------------------------- ! SW_TEMP ! ---------------------------------------------------------------------- ! !> \file sw_temp.f90 !! \BRIEF !> Module with sw_temp function - compute in-situ T from potential T !> Function to compute in-situ temperature [C] from potential temperature [C] FUNCTION sw_temp( s, t, p, pr ) ! ============================================================= ! SW_TEMP ! Computes in-situ temperature [C] from potential temperature [C] ! Routine available in seawater.f (used for MIT GCM) ! Downloaded seawater.f (on 17 April 2009) from ! http://ecco2.jpl.nasa.gov/data1/beaufort/MITgcm/bin/ ! ============================================================= ! REFERENCES: ! Fofonoff, P. and Millard, R.C. Jr ! Unesco 1983. Algorithms for computation of fundamental properties of ! seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp. ! Eqn.(31) p.39 ! Bryden, H. 1973. ! "New Polynomials for thermal expansion, adiabatic temperature gradient ! and potential temperature of sea water." ! DEEP-SEA RES., 1973, Vol20,401-408. ! ============================================================= ! Simple modifications: J. C. Orr, 16 April 2009 ! - combined fortran code from MITgcm site & simplification in ! CSIRO code (matlab equivalent) from Phil Morgan USE mocsy_singledouble IMPLICIT NONE ! Input arguments: ! ----------------------------------------------- ! s = salinity [psu (PSS-78) ] ! t = potential temperature [degree C (IPTS-68)] ! p = pressure [db] ! pr = reference pressure [db] !> salinity [psu (PSS-78)] REAL(kind=wp) :: s !> potential temperature [degree C (IPTS-68)] REAL(kind=wp) :: t !> pressure [db] REAL(kind=wp) :: p !> reference pressure [db] REAL(kind=wp) :: pr REAL(kind=wp) :: ds, dt, dp, dpr REAL(kind=wp) :: dsw_temp REAL(kind=wp) :: sw_temp ! EXTERNAL sw_ptmp ! REAL(kind=wp) :: sw_ptmp ds = DBLE(s) dt = DBLE(t) dp = DBLE(p) dpr = DBLE(pr) ! Simple solution ! (see https://svn.mpl.ird.fr/us191/oceano/tags/V0/lib/matlab/seawater/sw_temp.m) ! Carry out inverse calculation by swapping P_ref (pr) and Pressure (p) ! in routine that is normally used to compute potential temp from temp dsw_temp = sw_ptmp(ds, dt, dpr, dp) sw_temp = REAL(dsw_temp) ! The above simplification works extremely well (compared to Table in 1983 report) ! whereas the sw_temp routine from MIT GCM site does not seem to work right RETURN END FUNCTION sw_temp ! ---------------------------------------------------------------------- ! TPOT ! ---------------------------------------------------------------------- ! !> \file tpot.f90 !! \BRIEF !> Module with tpot subroutine - compute potential T from in situ T,S,P !> Compute potential temperature from arrays of in situ temp, salinity, and pressure. !! This subroutine is needed because sw_ptmp is a function (using scalars not arrays) SUBROUTINE tpot(salt, tempis, press, pressref, N, tempot) ! Purpose: ! Compute potential temperature from arrays of in situ temp, salinity, and pressure. ! Needed because sw_ptmp is a function USE mocsy_singledouble IMPLICIT NONE !> number of records INTEGER, intent(in) :: N ! INPUT variables !> salinity [psu] REAL(kind=wp), INTENT(in), DIMENSION(N) :: salt !> in situ temperature [C] REAL(kind=wp), INTENT(in), DIMENSION(N) :: tempis !> pressure [db] REAL(kind=wp), INTENT(in), DIMENSION(N) :: press !f2py optional , depend(salt) :: n=len(salt) !> pressure reference level [db] REAL(kind=wp), INTENT(in) :: pressref ! OUTPUT variables: !> potential temperature [C] for pressref REAL(kind=wp), INTENT(out), DIMENSION(N) :: tempot REAL(kind=wp) :: dsalt, dtempis, dpress, dpressref REAL(kind=wp) :: dtempot INTEGER :: i ! REAL(kind=wp) :: sw_ptmp ! EXTERNAL sw_ptmp DO i = 1,N dsalt = DBLE(salt(i)) dtempis = DBLE(tempis(i)) dpress = DBLE(press(i)) dpressref = DBLE(pressref) dtempot = sw_ptmp(dsalt, dtempis, dpress, dpressref) tempot(i) = REAL(dtempot) END DO RETURN END SUBROUTINE tpot ! ---------------------------------------------------------------------- ! TIS ! ---------------------------------------------------------------------- ! !> \file tis.f90 !! \BRIEF !> Module with tis subroutine - compute in situ T from S,T,P !> Compute in situ temperature from arrays of potential temp, salinity, and pressure. !! This subroutine is needed because sw_temp is a function (using scalars not arrays) SUBROUTINE tis(salt, tempot, press, pressref, N, tempis) ! Purpose: ! Compute in situ temperature from arrays of in situ temp, salinity, and pressure. ! Needed because sw_temp is a function USE mocsy_singledouble IMPLICIT NONE !> number of records INTEGER, intent(in) :: N ! INPUT variables !> salinity [psu] REAL(kind=wp), INTENT(in), DIMENSION(N) :: salt !> potential temperature [C] REAL(kind=wp), INTENT(in), DIMENSION(N) :: tempot !> pressure [db] REAL(kind=wp), INTENT(in), DIMENSION(N) :: press !f2py optional , depend(salt) :: n=len(salt) !> pressure reference level [db] REAL(kind=wp), INTENT(in) :: pressref ! OUTPUT variables: !> in situ temperature [C] REAL(kind=wp), INTENT(out), DIMENSION(N) :: tempis ! REAL(kind=wp) :: dsalt, dtempis, dpress, dpressref ! REAL(kind=wp) :: dtempot INTEGER :: i ! REAL(kind=wp) :: sw_temp ! REAL(kind=wp) :: sw_temp ! EXTERNAL sw_temp DO i = 1,N !dsalt = DBLE(salt(i)) !dtempot = DBLE(tempot(i)) !dpress = DBLE(press(i)) !dpressref = DBLE(pressref) !dtempis = sw_temp(dsalt, dtempot, dpress, dpressref) !tempis(i) = REAL(dtempis) tempis = sw_temp(salt(i), tempot(i), press(i), pressref) END DO RETURN END SUBROUTINE tis ! ---------------------------------------------------------------------- ! P80 ! ---------------------------------------------------------------------- ! !> \file p80.f90 !! \BRIEF !> Module with p80 function - compute pressure from depth !> Function to compute pressure from depth using Saunder's (1981) formula with eos80. FUNCTION p80(dpth,xlat) ! Compute Pressure from depth using Saunder's (1981) formula with eos80. ! Reference: ! Saunders, Peter M. (1981) Practical conversion of pressure ! to depth, J. Phys. Ooceanogr., 11, 573-574, (1981) ! Coded by: ! R. Millard ! March 9, 1983 ! check value: p80=7500.004 dbars at lat=30 deg., depth=7321.45 meters ! Modified (slight format changes + added ref. details): ! J. Orr, 16 April 2009 USE mocsy_singledouble IMPLICIT NONE ! Input variables: !> depth [m] REAL(kind=wp), INTENT(in) :: dpth !> latitude [degrees] REAL(kind=wp), INTENT(in) :: xlat ! Output variable: !> pressure [db] REAL(kind=wp) :: p80 ! Local variables: REAL(kind=wp) :: pi REAL(kind=wp) :: plat, d, c1 pi=3.141592654 plat = ABS(xlat*pi/180.) d = SIN(plat) c1 = 5.92e-3+d**2 * 5.25e-3 p80 = ((1-c1)-SQRT(((1-c1)**2)-(8.84e-6*dpth))) / 4.42e-6 RETURN END FUNCTION p80 ! ---------------------------------------------------------------------- ! RHO ! ---------------------------------------------------------------------- ! !> \file rho.f90 !! \BRIEF !> Module with rho function - computes in situ density from S, T, P !> Function to compute in situ density from salinity (psu), in situ temperature (C), & pressure (bar) FUNCTION rho(salt, temp, pbar) ! Compute in situ density from salinity (psu), in situ temperature (C), & pressure (bar) USE mocsy_singledouble IMPLICIT NONE !> salinity [psu] REAL(kind=wp) :: salt !> in situ temperature (C) REAL(kind=wp) :: temp !> pressure (bar) [Note units: this is NOT in decibars] REAL(kind=wp) :: pbar REAL(kind=wp) :: s, t, p ! REAL(kind=wp) :: t68 REAL(kind=wp) :: X REAL(kind=wp) :: rhow, rho0 REAL(kind=wp) :: a, b, c REAL(kind=wp) :: Ksbmw, Ksbm0, Ksbm REAL(kind=wp) :: drho REAL(kind=wp) :: rho ! Input arguments: ! ------------------------------------- ! s = salinity [psu (PSS-78) ] ! t = in situ temperature [degree C (IPTS-68)] ! p = pressure [bar] !!!! (not in [db] s = DBLE(salt) t = DBLE(temp) p = DBLE(pbar) ! Convert the temperature on today's "ITS 90" scale to older IPTS 68 scale ! (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote) ! According to Best-Practices guide, line above should be commented & 2 lines below should be uncommented ! Guide's answer of rho (s=35, t=25, p=0) = 1023.343 is for temperature given on ITPS-68 scale ! t68 = (T - 0.0002) / 0.99975 ! X = t68 ! Finally, don't do the ITS-90 to IPTS-68 conversion (T input var now already on IPTS-68 scale) X = T ! Density of pure water rhow = 999.842594d0 + 6.793952e-2_wp*X & -9.095290e-3_wp*X*X + 1.001685e-4_wp*X**3 & -1.120083e-6_wp*X**4 + 6.536332e-9_wp*X**5 ! Density of seawater at 1 atm, P=0 A = 8.24493e-1_wp - 4.0899e-3_wp*X & + 7.6438e-5_wp*X*X - 8.2467e-7_wp*X**3 + 5.3875e-9_wp*X**4 B = -5.72466e-3_wp + 1.0227e-4_wp*X - 1.6546e-6_wp*X*X C = 4.8314e-4_wp rho0 = rhow + A*S + B*S*SQRT(S) + C*S**2.0d0 ! Secant bulk modulus of pure water ! The secant bulk modulus is the average change in pressure ! divided by the total change in volume per unit of initial volume. Ksbmw = 19652.21d0 + 148.4206d0*X - 2.327105d0*X*X & + 1.360477e-2_wp*X**3 - 5.155288e-5_wp*X**4 ! Secant bulk modulus of seawater at 1 atm Ksbm0 = Ksbmw + S*( 54.6746d0 - 0.603459d0*X + 1.09987e-2_wp*X**2 & - 6.1670e-5_wp*X**3) & + S*SQRT(S)*( 7.944e-2_wp + 1.6483e-2_wp*X - 5.3009e-4_wp*X**2) ! Secant bulk modulus of seawater at S,T,P Ksbm = Ksbm0 & + P*(3.239908d0 + 1.43713e-3_wp*X + 1.16092e-4_wp*X**2 - 5.77905e-7_wp*X**3) & + P*S*(2.2838e-3_wp - 1.0981e-5_wp*X - 1.6078e-6_wp*X**2) & + P*S*SQRT(S)*1.91075e-4_wp & + P*P*(8.50935e-5_wp - 6.12293e-6_wp*X + 5.2787e-8_wp*X**2) & + P*P*S*(-9.9348e-7_wp + 2.0816e-8_wp*X + 9.1697e-10_wp*X**2) ! Density of seawater at S,T,P drho = rho0/(1.0d0 - P/Ksbm) rho = REAL(drho) RETURN END FUNCTION rho ! ---------------------------------------------------------------------- ! RHOINSITU ! ---------------------------------------------------------------------- ! !> \file rhoinsitu.f90 !! \BRIEF !> Module with rhoinsitu subroutine - compute in situ density from S, Tis, P !> Compute in situ density from salinity (psu), in situ temperature (C), & pressure (db). !! This subroutine is needed because rho is a function (using scalars not arrays) SUBROUTINE rhoinsitu(salt, tempis, pdbar, N, rhois) ! Purpose: ! Compute in situ density from salinity (psu), in situ temperature (C), & pressure (db) ! Needed because rho is a function USE mocsy_singledouble IMPLICIT NONE INTEGER :: N ! INPUT variables ! salt = salinity [psu] ! tempis = in situ temperature [C] ! pdbar = pressure [db] !> salinity [psu] REAL(kind=wp), INTENT(in), DIMENSION(N) :: salt !> in situ temperature [C] REAL(kind=wp), INTENT(in), DIMENSION(N) :: tempis !> pressure [db] REAL(kind=wp), INTENT(in), DIMENSION(N) :: pdbar !f2py optional , depend(salt) :: n=len(salt) ! OUTPUT variables: ! rhois = in situ density !> in situ density [kg/m3] REAL(kind=wp), INTENT(out), DIMENSION(N) :: rhois ! Local variables INTEGER :: i ! REAL(kind=wp) :: rho ! EXTERNAL rho DO i = 1,N rhois(i) = rho(salt(i), tempis(i), pdbar(i)/10.) END DO RETURN END SUBROUTINE rhoinsitu ! ---------------------------------------------------------------------- ! DEPTH2PRESS ! ---------------------------------------------------------------------- ! !> \file depth2press.f90 !! \BRIEF !> Module with depth2press subroutine - converts depth to pressure !! with Saunders (1981) formula !> Compute pressure [db] from depth [m] & latitude [degrees north]. !! This subroutine is needed because p80 is a function (using scalars not arrays) SUBROUTINE depth2press(depth, lat, pdbar, N) ! Purpose: ! Compute pressure [db] from depth [m] & latitude [degrees north]. ! Needed because p80 is a function USE mocsy_singledouble IMPLICIT NONE !> number of records INTEGER, intent(in) :: N ! INPUT variables !> depth [m] REAL(kind=wp), INTENT(in), DIMENSION(N) :: depth !> latitude [degrees] REAL(kind=wp), INTENT(in), DIMENSION(N) :: lat !f2py optional , depend(depth) :: n=len(depth) ! OUTPUT variables: !> pressure [db] REAL(kind=wp), INTENT(out), DIMENSION(N) :: pdbar ! Local variables INTEGER :: i ! REAL(kind=wp) :: p80 ! EXTERNAL p80 DO i = 1,N pdbar(i) = p80(depth(i), lat(i)) END DO RETURN END SUBROUTINE depth2press ! ---------------------------------------------------------------------- ! CONSTANTS ! ---------------------------------------------------------------------- ! !> \file constants.f90 !! \BRIEF !> Module with contants subroutine - computes carbonate system constants !! from S,T,P !> Compute thermodynamic constants !! FROM temperature, salinity, and pressure (1D arrays) SUBROUTINE constants(K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, & K1p, K2p, K3p, Ksi, & St, Ft, Bt, & temp, sal, Patm, & depth, lat, N, & optT, optP, optB, optK1K2, optKf, optGAS) ! Purpose: ! Compute thermodynamic constants ! FROM: temperature, salinity, and pressure (1D arrays) ! INPUT variables: ! ================ ! Patm = atmospheric pressure [atm] ! depth = depth [m] (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure) ! = pressure [db] (with optP='db') ! lat = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m') ! = dummy array (unused when optP='db') ! temp = potential temperature [degrees C] (with optT='Tpot', i.e., models carry tempot, not temp) ! = in situ temperature [degrees C] (with optT='Tinsitu', e.g., for data) ! sal = salinity in [psu] ! --------- ! optT: choose in situ vs. potential temperature as input ! --------- ! NOTE: Carbonate chem calculations require IN-SITU temperature (not potential Temperature) ! -> 'Tpot' means input is pot. Temperature (in situ Temp "tempis" is computed) ! -> 'Tinsitu' means input is already in-situ Temperature, not pot. Temp ("tempis" not computed) ! --------- ! optP: choose depth (m) vs pressure (db) as input ! --------- ! -> 'm' means "depth" input is in "m" (thus in situ Pressure "p" [db] is computed) ! -> 'db' means "depth" input is already in situ pressure [db], not m (p = depth) ! --------- ! optB: ! --------- ! -> 'u74' means use classic formulation of Uppström (1974) for total Boron ! -> 'l10' means use newer formulation of Lee et al. (2010) for total Boron ! --------- ! optK1K2: ! --------- ! -> 'l' means use Lueker et al. (2000) formulations for K1 & K2 (recommended by Dickson et al. 2007) ! **** BUT this should only be used when 2 < T < 35 and 19 < S < 43 ! -> 'm10' means use Millero (2010) formulation for K1 & K2 (see Dickson et al., 2007) ! **** Valid for 0 < T < 50°C and 1 < S < 50 psu ! -> 'w14' means use Waters (2014) formulation for K1 & K2 (see Dickson et al., 2007) ! **** Valid for 0 < T < 50°C and 1 < S < 50 psu ! ----------- ! optKf: ! ---------- ! -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007) ! **** BUT Valid for 9 < T < 33°C and 10 < S < 40. ! -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994) ! ----------- ! optGAS: choose in situ vs. potential fCO2 and pCO2 ! --------- ! PRESSURE corrections for K0 and the fugacity coefficient (Cf) ! -> 'Pzero' = 'zero order' fCO2 and pCO2 (typical approach, which is flawed) ! considers in situ T & only atm pressure (hydrostatic=0) ! -> 'Ppot' = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) ! considers potential T & only atm pressure (hydrostatic press = 0) ! -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure) ! considers in situ T & total pressure (atm + hydrostatic) ! --------- ! OUTPUT variables: ! ================= ! K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi ! St, Ft, Bt USE mocsy_singledouble IMPLICIT NONE ! Input variables !> number of records INTEGER, INTENT(in) :: N !> in situ temperature (when optT='Tinsitu', typical data) !! OR potential temperature (when optT='Tpot', typical models) [degree C] REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp !> depth in meters (when optP='m') or decibars (when optP='db') REAL(kind=wp), INTENT(in), DIMENSION(N) :: depth !> latitude [degrees north] REAL(kind=wp), INTENT(in), DIMENSION(N) :: lat !> salinity [psu] REAL(kind=wp), INTENT(in), DIMENSION(N) :: sal !f2py optional , depend(sal) :: n=len(sal) !> atmospheric pressure [atm] REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm !> for temp input, choose \b 'Tinsitu' for in situ Temp or !! \b 'Tpot' for potential temperature (in situ Temp is computed, needed for models) CHARACTER(7), INTENT(in) :: optT !> for depth input, choose \b "db" for decibars (in situ pressure) or \b "m" for meters (pressure is computed, needed for models) CHARACTER(2), INTENT(in) :: optP !> for total boron, choose either \b 'u74' (Uppstrom, 1974) or \b 'l10' (Lee et al., 2010). !! The 'l10' formulation is based on 139 measurements (instead of 20), !! uses a more accurate method, and !! generally increases total boron in seawater by 4% !f2py character*3 optional, intent(in) :: optB='l10' CHARACTER(3), OPTIONAL, INTENT(in) :: optB !> for Kf, choose either \b 'pf' (Perez & Fraga, 1987) or \b 'dg' (Dickson & Riley, 1979) !f2py character*2 optional, intent(in) :: optKf='pf' CHARACTER(2), OPTIONAL, INTENT(in) :: optKf !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010) or \b 'w14' (Waters et al., 2014) !f2py character*3 optional, intent(in) :: optK1K2='l' CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2 !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction) !! 'Ppot' - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean !f2py character*7 optional, intent(in) :: optGAS='Pinsitu' CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS ! Ouput variables !> solubility of CO2 in seawater (Weiss, 1974), also known as K0 REAL(kind=wp), INTENT(out), DIMENSION(N) :: K0 !> K1 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 REAL(kind=wp), INTENT(out), DIMENSION(N) :: K1 !> K2 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 REAL(kind=wp), INTENT(out), DIMENSION(N) :: K2 !> equilibrium constant for dissociation of boric acid REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kb !> equilibrium constant for the dissociation of water (Millero, 1995) REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kw !> equilibrium constant for the dissociation of bisulfate (Dickson, 1990) REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ks !> equilibrium constant for the dissociation of hydrogen fluoride !! either from Dickson and Riley (1979) or from Perez and Fraga (1987), depending on optKf REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kf !> solubility product for calcite (Mucci, 1983) REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kspc !> solubility product for aragonite (Mucci, 1983) REAL(kind=wp), INTENT(out), DIMENSION(N) :: Kspa !> 1st dissociation constant for phosphoric acid (Millero, 1995) REAL(kind=wp), INTENT(out), DIMENSION(N) :: K1p !> 2nd dissociation constant for phosphoric acid (Millero, 1995) REAL(kind=wp), INTENT(out), DIMENSION(N) :: K2p !> 3rd dissociation constant for phosphoric acid (Millero, 1995) REAL(kind=wp), INTENT(out), DIMENSION(N) :: K3p !> equilibrium constant for the dissociation of silicic acid (Millero, 1995) REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ksi !> total sulfate (Morris & Riley, 1966) REAL(kind=wp), INTENT(out), DIMENSION(N) :: St !> total fluoride (Riley, 1965) REAL(kind=wp), INTENT(out), DIMENSION(N) :: Ft !> total boron !! from either Uppstrom (1974) or Lee et al. (2010), depending on optB REAL(kind=wp), INTENT(out), DIMENSION(N) :: Bt ! Local variables REAL(kind=wp) :: ssal REAL(kind=wp) :: p REAL(kind=wp) :: tempot, tempis68, tempot68 REAL(kind=wp) :: tempis REAL(kind=wp) :: is, invtk, dlogtk, is2, s2, sqrtis REAL(kind=wp) :: Ks_0p, Kf_0p REAL(kind=wp) :: total2free, free2SWS, total2SWS, SWS2total REAL(kind=wp) :: total2free_0p, free2SWS_0p, total2SWS_0p ! REAL(kind=wp) :: free2SWS, free2SWS_0p REAL(kind=wp) :: dtempot, dtempot68 REAL(kind=wp) :: R REAL(kind=wp) :: pK1o, ma1, mb1, mc1, pK1 REAL(kind=wp) :: pK2o, ma2, mb2, mc2, pK2 REAL(kind=wp), DIMENSION(12) :: a0, a1, a2, b0, b1, b2 REAL(kind=wp), DIMENSION(12) :: deltav, deltak, lnkpok0 REAL(kind=wp) :: tmp, nK0we74 INTEGER :: i, icount, ipc REAL(kind=wp) :: t, tk, tk0, prb REAL(kind=wp) :: s, sqrts, s15, scl REAL(kind=wp) :: Phydro_atm, Patmd, Ptot, Rgas_atm, vbarCO2 ! Arrays to pass optional arguments into or use defaults (Dickson et al., 2007) CHARACTER(3) :: opB CHARACTER(2) :: opKf CHARACTER(3) :: opK1K2 CHARACTER(7) :: opGAS ! CONSTANTS ! ========= ! Constants in formulation for Pressure effect on K's (Millero, 95) ! with corrected coefficients for Kb, Kw, Ksi, etc. ! index: 1) K1 , 2) K2, 3) Kb, 4) Kw, 5) Ks, 6) Kf, 7) Kspc, 8) Kspa, ! 9) K1P, 10) K2P, 11) K3P, 12) Ksi DATA a0 /-25.5_wp, -15.82_wp, -29.48_wp, -20.02_wp, & -18.03_wp, -9.78_wp, -48.76_wp, -45.96_wp, & -14.51_wp, -23.12_wp, -26.57_wp, -29.48_wp/ DATA a1 /0.1271_wp, -0.0219_wp, 0.1622_wp, 0.1119_wp, & 0.0466_wp, -0.0090_wp, 0.5304_wp, 0.5304_wp, & 0.1211_wp, 0.1758_wp, 0.2020_wp, 0.1622_wp/ DATA a2 / 0.0_wp, 0.0_wp, -2.608e-3_wp, -1.409e-3_wp, & 0.316e-3_wp, -0.942e-3_wp, 0.0_wp, 0.0_wp, & -0.321e-3_wp, -2.647e-3_wp, -3.042e-3_wp, -2.6080e-3_wp/ DATA b0 /-3.08e-3_wp, 1.13e-3_wp, -2.84e-3_wp, -5.13e-3_wp, & -4.53e-3_wp, -3.91e-3_wp, -11.76e-3_wp, -11.76e-3_wp, & -2.67e-3_wp, -5.15e-3_wp, -4.08e-3_wp, -2.84e-3_wp/ DATA b1 /0.0877e-3_wp, -0.1475e-3_wp, 0.0_wp, 0.0794e-3_wp, & 0.09e-3_wp, 0.054e-3_wp, 0.3692e-3_wp, 0.3692e-3_wp, & 0.0427e-3_wp, 0.09e-3_wp, 0.0714e-3_wp, 0.0_wp/ DATA b2 /12*0.0_wp/ ! Set defaults for optional arguments (in Fortran 90) ! Note: Optional arguments with f2py (python) are set above with ! the !f2py statements that precede each type declaraion IF (PRESENT(optB)) THEN opB = optB ELSE opB = 'l10' ENDIF IF (PRESENT(optKf)) THEN opKf = optKf ELSE opKf = 'pf' ENDIF IF (PRESENT(optK1K2)) THEN opK1K2 = optK1K2 ELSE opK1K2 = 'l' ENDIF IF (PRESENT(optGAS)) THEN opGAS = optGAS ELSE opGAS = 'Pinsitu' ENDIF R = 83.14472_wp icount = 0 DO i = 1, N icount = icount + 1 ! =============================================================== ! Convert model depth -> press; convert model Theta -> T in situ ! =============================================================== ! * Model temperature tracer is usually "potential temperature" ! * Model vertical grid is usually in meters ! BUT carbonate chem routines require pressure & in-situ T ! Thus before computing chemistry, if appropriate, ! convert these 2 model vars (input to this routine) ! - depth [m] => convert to pressure [db] ! - potential temperature (C) => convert to in-situ T (C) ! ------------------------------------------------------- ! 1) Compute pressure [db] from depth [m] and latitude [degrees] (if input is m, for models) IF (trim(optP) == 'm' ) THEN ! Compute pressure [db] from depth [m] and latitude [degrees] p = p80(depth(i), lat(i)) ELSEIF (trim(optP) == 'db' ) THEN ! In this case (where optP = 'db'), p is input & output (no depth->pressure conversion needed) p = depth(i) ELSE PRINT *,"optP must be 'm' or 'db'" STOP ENDIF ! 2) Convert potential T to in-situ T (if input is Tpot, i.e. case for models): IF (trim(optT) == 'Tpot' .OR. trim(optT) == 'tpot') THEN tempot = temp(i) ! This is the case for most models and some data ! a) Convert the pot. temp on today's "ITS 90" scale to older IPTS 68 scale ! (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote) tempot68 = (tempot - 0.0002) / 0.99975 ! b) Compute "in-situ Temperature" from "Potential Temperature" (both on IPTS 68) tempis68 = sw_temp(sal(i), tempot68, p, 0.0_wp ) ! c) Convert the in-situ temp on older IPTS 68 scale to modern scale (ITS 90) tempis = 0.99975*tempis68 + 0.0002 ! Note: parts (a) and (c) above are tiny corrections; ! part (b) is a big correction for deep waters (but zero at surface) ELSEIF (trim(optT) == 'Tinsitu' .OR. trim(optT) == 'tinsitu') THEN ! When optT = 'Tinsitu', tempis is input & output (no tempot needed) tempis = temp(i) tempis68 = (temp(i) - 0.0002) / 0.99975 ! dtempot68 = sw_ptmp(DBLE(sal(i)), DBLE(tempis68), DBLE(p), 0.0_wp) dtempot68 = sw_ptmp(sal(i), tempis68, p, 0.0_wp) dtempot = 0.99975*dtempot68 + 0.0002 ELSE PRINT *,"optT must be either 'Tpot' or 'Tinsitu'" PRINT *,"you specified optT =", trim(optT) STOP ENDIF ! Compute constants: IF (temp(i) >= -5. .AND. temp(i) < 1.0e+2) THEN ! Test to indicate if any of input variables are unreasonable IF ( sal(i) < 0. .OR. sal(i) > 1e+3) THEN PRINT *, 'i, icount, temp, sal =', i, icount, temp(i), sal(i) ENDIF ! Zero out negative salinity (prev case for OCMIP2 model w/ slightly negative S in some coastal cells) IF (sal(i) < 0.0) THEN ssal = 0.0 ELSE ssal = sal(i) ENDIF ! Absolute temperature (Kelvin) and related values t = DBLE(tempis) tk = 273.15d0 + t invtk=1.0d0/tk dlogtk=LOG(tk) ! Atmospheric pressure Patmd = DBLE(Patm(i)) ! Hydrostatic pressure (prb is in bars) prb = DBLE(p) / 10.0d0 ! Salinity and simply related values s = DBLE(ssal) s2=s*s sqrts=SQRT(s) s15=s**1.5d0 scl=s/1.80655d0 ! Ionic strength: is = 19.924d0*s/(1000.0d0 - 1.005d0*s) is2 = is*is sqrtis = SQRT(is) ! Total concentrations for sulfate, fluoride, and boron ! Sulfate: Morris & Riley (1966) St(i) = 0.14d0 * scl/96.062d0 ! Fluoride: Riley (1965) Ft(i) = 0.000067d0 * scl/18.9984d0 ! Boron: IF (trim(opB) == 'l10') THEN ! New formulation from Lee et al (2010) Bt(i) = 0.0002414d0 * scl/10.811d0 ELSEIF (trim(opB) == 'u74') THEN ! Classic formulation from Uppström (1974) Bt(i) = 0.000232d0 * scl/10.811d0 ELSE PRINT *,"optB must be 'l10' or 'u74'" STOP ENDIF ! K0 (K Henry) ! CO2(g) <-> CO2(aq.) ! K0 = [CO2]/ fCO2 ! Weiss (1974) [mol/kg/atm] IF (trim(opGAS) == 'Pzero' .OR. trim(opGAS) == 'pzero') THEN tk0 = tk !in situ temperature (K) for K0 calculation Ptot = Patmd !total pressure (in atm) = atmospheric pressure ONLY ELSEIF (trim(opGAS) == 'Ppot' .OR. trim(opGAS) == 'ppot') THEN tk0 = dtempot + 273.15d0 !potential temperature (K) for K0 calculation as needed for potential fCO2 & pCO2 Ptot = Patmd !total pressure (in atm) = atmospheric pressure ONLY ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN tk0 = tk !in situ temperature (K) for K0 calculation Phydro_atm = prb / 1.01325d0 !convert hydrostatic pressure from bar to atm (1.01325 bar / atm) Ptot = Patmd + Phydro_atm !total pressure (in atm) = atmospheric pressure + hydrostatic pressure ELSE PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'" STOP ENDIF tmp = 9345.17d0/tk0 - 60.2409d0 + 23.3585d0 * LOG(tk0/100.0d0) nK0we74 = tmp + s*(0.023517d0 - 0.00023656d0*tk0 + 0.0047036e-4_wp*tk0*tk0) K0(i) = EXP(nK0we74) ! K1 = [H][HCO3]/[H2CO3] ! K2 = [H][CO3]/[HCO3] IF (trim(opK1K2) == 'l') THEN ! Mehrbach et al. (1973) refit, by Lueker et al. (2000) (total pH scale) K1(i) = 10.0d0**(-1.0d0*(3633.86d0*invtk - 61.2172d0 + 9.6777d0*dlogtk & - 0.011555d0*s + 0.0001152d0*s2)) K2(i) = 10.0d0**(-1*(471.78d0*invtk + 25.9290d0 - 3.16967d0*dlogtk & - 0.01781d0*s + 0.0001122d0*s2)) ELSEIF (trim(opK1K2) == 'm10') THEN ! Millero (2010, Mar. Fresh Wat. Res.) (seawater pH scale) pK1o = 6320.813d0*invtk + 19.568224d0*dlogtk -126.34048d0 ma1 = 13.4038d0*sqrts + 0.03206d0*s - (5.242e-5)*s2 mb1 = -530.659d0*sqrts - 5.8210d0*s mc1 = -2.0664d0*sqrts pK1 = pK1o + ma1 + mb1*invtk + mc1*dlogtk K1(i) = 10.0d0**(-pK1) pK2o = 5143.692d0*invtk + 14.613358d0*dlogtk -90.18333d0 ma2 = 21.3728d0*sqrts + 0.1218d0*s - (3.688e-4)*s2 mb2 = -788.289d0*sqrts - 19.189d0*s mc2 = -3.374d0*sqrts pK2 = pK2o + ma2 + mb2*invtk + mc2*dlogtk K2(i) = 10.0d0**(-pK2) ELSEIF (trim(opK1K2) == 'w14') THEN ! Waters, Millero, Woosley (Mar. Chem., 165, 66-67, 2014) (seawater scale) pK1o = 6320.813d0*invtk + 19.568224d0*dlogtk -126.34048d0 ma1 = 13.409160d0*sqrts + 0.031646d0*s - (5.1895e-5)*s2 mb1 = -531.3642d0*sqrts - 5.713d0*s mc1 = -2.0669166d0*sqrts pK1 = pK1o + ma1 + mb1*invtk + mc1*dlogtk K1(i) = 10.0d0**(-pK1) pK2o = 5143.692d0*invtk + 14.613358d0*dlogtk -90.18333d0 ma2 = 21.225890d0*sqrts + 0.12450870d0*s - (3.7243e-4_r8)*s2 mb2 = -779.3444d0*sqrts - 19.91739d0*s mc2 = -3.3534679d0*sqrts pK2 = pK2o + ma2 + mb2*invtk + mc2*dlogtk K2(i) = 10.0d0**(-pK2) ELSE PRINT *, "optK1K2 must be either 'l' or 'm10', or 'w14'" STOP ENDIF ! Kb = [H][BO2]/[HBO2] ! (total scale) ! Millero p.669 (1995) using data from Dickson (1990) Kb(i) = EXP((-8966.90d0 - 2890.53d0*sqrts - 77.942d0*s + & 1.728d0*s15 - 0.0996d0*s2)*invtk + & (148.0248d0 + 137.1942d0*sqrts + 1.62142d0*s) + & (-24.4344d0 - 25.085d0*sqrts - 0.2474d0*s) * & dlogtk + 0.053105d0*sqrts*tk) ! K1p = [H][H2PO4]/[H3PO4] ! (seawater scale) ! DOE(1994) eq 7.2.20 with footnote using data from Millero (1974) ! Millero (1995), p.670, eq. 65 ! Use Millero equation's 115.540 constant instead of 115.525 (Dickson et al., 2007). ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) ! And we want to stay on the SWS scale anyway for the pressure correction later. K1p(i) = EXP(-4576.752d0*invtk + 115.540d0 - 18.453d0*dlogtk + & (-106.736d0*invtk + 0.69171d0) * sqrts + & (-0.65643d0*invtk - 0.01844d0) * s) ! K2p = [H][HPO4]/[H2PO4] ! (seawater scale) ! DOE(1994) eq 7.2.23 with footnote using data from Millero (1974)) ! Millero (1995), p.670, eq. 66 ! Use Millero equation's 172.1033 constant instead of 172.0833 (Dickson et al., 2007). ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) ! And we want to stay on the SWS scale anyway for the pressure correction later. K2p(i) = EXP(-8814.715d0*invtk + 172.1033d0 - 27.927d0*dlogtk + & (-160.340d0*invtk + 1.3566d0)*sqrts + & (0.37335d0*invtk - 0.05778d0)*s) ! K3p = [H][PO4]/[HPO4] ! (seawater scale) ! DOE(1994) eq 7.2.26 with footnote using data from Millero (1974) ! Millero (1995), p.670, eq. 67 ! Use Millero equation's 18.126 constant instead of 18.141 (Dickson et al., 2007). ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) ! And we want to stay on the SWS scale anyway for the pressure correction later. K3p(i) = EXP(-3070.75d0*invtk - 18.126d0 + & (17.27039d0*invtk + 2.81197d0) * & sqrts + (-44.99486d0*invtk - 0.09984d0) * s) ! Ksi = [H][SiO(OH)3]/[Si(OH)4] ! (seawater scale) ! Millero (1995), p.671, eq. 72 ! Use Millero equation's 117.400 constant instead of 117.385 (Dickson et al., 2007). ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) ! And we want to stay on the SWS scale anyway for the pressure correction later. Ksi(i) = EXP(-8904.2d0*invtk + 117.400d0 - 19.334d0*dlogtk + & (-458.79d0*invtk + 3.5913d0) * sqrtis + & (188.74d0*invtk - 1.5998d0) * is + & (-12.1652d0*invtk + 0.07871d0) * is2 + & LOG(1.0 - 0.001005d0*s)) ! Kw = [H][OH] ! (seawater scale) ! Millero (1995) p.670, eq. 63 from composite data ! Use Millero equation's 148.9802 constant instead of 148.9652 (Dickson et al., 2007). ! The latter is only an crude approximation to convert to Total scale (by subtracting 0.015) ! And we want to stay on the SWS scale anyway for the pressure correction later. Kw(i) = EXP(-13847.26d0*invtk + 148.9802d0 - 23.6521d0*dlogtk + & (118.67d0*invtk - 5.977d0 + 1.0495d0 * dlogtk) * & sqrts - 0.01615d0 * s) ! Ks = [H][SO4]/[HSO4] ! (free scale) ! Dickson (1990, J. chem. Thermodynamics 22, 113) Ks_0p = EXP(-4276.1d0*invtk + 141.328d0 - 23.093d0*dlogtk & + (-13856.d0*invtk + 324.57d0 - 47.986d0*dlogtk) * sqrtis & + (35474.d0*invtk - 771.54 + 114.723d0*dlogtk) * is & - 2698.d0*invtk*is**1.5 + 1776.d0*invtk*is2 & + LOG(1.0d0 - 0.001005d0*s)) ! Kf = [H][F]/[HF] ! (total scale) IF (trim(opKf) == 'dg') THEN ! Dickson and Riley (1979) -- change pH scale to total (following Dickson & Goyet, 1994) Kf_0p = EXP(1590.2d0*invtk - 12.641d0 + 1.525d0*sqrtis + & LOG(1.0d0 - 0.001005d0*s) + & LOG(1.0d0 + St(i)/Ks_0p)) ELSEIF (trim(opKf) == 'pf') THEN ! Perez and Fraga (1987) - Already on Total scale (no need for last line above) ! Formulation as given in Dickson et al. (2007) Kf_0p = EXP(874.d0*invtk - 9.68d0 + 0.111d0*sqrts) ELSE PRINT *, "optKf must be either 'dg' or 'pf'" STOP ENDIF ! Kspc (calcite) - apparent solubility product of calcite ! (no scale) ! Kspc = [Ca2+] [CO32-] when soln is in equilibrium w/ calcite ! Mucci 1983 mol/kg-soln Kspc(i) = 10d0**(-171.9065d0 - 0.077993d0*tk + 2839.319d0/tk & + 71.595d0*LOG10(tk) & + (-0.77712d0 + 0.0028426d0*tk + 178.34d0/tk)*sqrts & -0.07711d0*s + 0.0041249d0*s15 ) ! Kspa (aragonite) - apparent solubility product of aragonite ! (no scale) ! Kspa = [Ca2+] [CO32-] when soln is in equilibrium w/ aragonite ! Mucci 1983 mol/kg-soln Kspa(i) = 10.d0**(-171.945d0 - 0.077993d0*tk + 2903.293d0/tk & +71.595d0*LOG10(tk) & +(-0.068393d0 + 0.0017276d0*tk + 88.135d0/tk)*sqrts & -0.10018d0*s + 0.0059415d0*s15 ) ! Pressure effect on K0 based on Weiss (1974, equation 5) Rgas_atm = 82.05736_wp ! (cm3 * atm) / (mol * K) CODATA (2006) vbarCO2 = 32.3_wp ! partial molal volume (cm3 / mol) from Weiss (1974, Appendix, paragraph 3) K0(i) = K0(i) * exp( ((1-Ptot)*vbarCO2)/(Rgas_atm*tk0) ) ! Weiss (1974, equation 5) ! Pressure effect on all other K's (based on Millero, (1995) ! index: K1(1), K2(2), Kb(3), Kw(4), Ks(5), Kf(6), Kspc(7), Kspa(8), ! K1p(9), K2p(10), K3p(11), Ksi(12) DO ipc = 1, 12 deltav(ipc) = a0(ipc) + a1(ipc) *t + a2(ipc) *t*t deltak(ipc) = (b0(ipc) + b1(ipc) *t + b2(ipc) *t*t) lnkpok0(ipc) = (-(deltav(ipc)) & +(0.5d0*deltak(ipc) * prb) & ) * prb/(R*tk) END DO ! Pressure correction on Ks (Free scale) Ks(i) = Ks_0p*EXP(lnkpok0(5)) ! Conversion factor total -> free scale total2free = 1.d0/(1.d0 + St(i)/Ks(i)) ! Kfree = Ktotal*total2free ! Conversion factor total -> free scale at pressure zero total2free_0p = 1.d0/(1.d0 + St(i)/Ks_0p) ! Kfree = Ktotal*total2free ! Pressure correction on Kf ! Kf must be on FREE scale before correction Kf_0p = Kf_0p * total2free_0p !Convert from Total to Free scale (pressure 0) Kf(i) = Kf_0p * EXP(lnkpok0(6)) !Pressure correction (on Free scale) Kf(i) = Kf(i)/total2free !Convert back from Free to Total scale ! Convert between seawater and total hydrogen (pH) scales free2SWS = 1.d0 + St(i)/Ks(i) + Ft(i)/(Kf(i)*total2free) ! using Kf on free scale total2SWS = total2free * free2SWS ! KSWS = Ktotal*total2SWS SWS2total = 1.d0 / total2SWS ! Conversion at pressure zero free2SWS_0p = 1.d0 + St(i)/Ks_0p + Ft(i)/(Kf_0p) ! using Kf on free scale total2SWS_0p = total2free_0p * free2SWS_0p ! KSWS = Ktotal*total2SWS ! Convert from Total to Seawater scale before pressure correction ! Must change to SEAWATER scale: K1, K2, Kb IF (trim(optK1K2) == 'l') THEN K1(i) = K1(i)*total2SWS_0p K2(i) = K2(i)*total2SWS_0p !This conversion is unnecessary for the K1,K2 from Millero (2010), !since we use here the formulation already on the seawater scale ENDIF Kb(i) = Kb(i)*total2SWS_0p ! Already on SEAWATER scale: K1p, K2p, K3p, Kb, Ksi, Kw ! Other contants (keep on another scale): ! - K0 (independent of pH scale, already pressure corrected) ! - Ks (already on Free scale; already pressure corrected) ! - Kf (already on Total scale; already pressure corrected) ! - Kspc, Kspa (independent of pH scale; pressure-corrected below) ! Perform actual pressure correction (on seawater scale) K1(i) = K1(i)*EXP(lnkpok0(1)) K2(i) = K2(i)*EXP(lnkpok0(2)) Kb(i) = Kb(i)*EXP(lnkpok0(3)) Kw(i) = Kw(i)*EXP(lnkpok0(4)) Kspc(i) = Kspc(i)*EXP(lnkpok0(7)) Kspa(i) = Kspa(i)*EXP(lnkpok0(8)) K1p(i) = K1p(i)*EXP(lnkpok0(9)) K2p(i) = K2p(i)*EXP(lnkpok0(10)) K3p(i) = K3p(i)*EXP(lnkpok0(11)) Ksi(i) = Ksi(i)*EXP(lnkpok0(12)) ! Convert back to original total scale: K1(i) = K1(i) *SWS2total K2(i) = K2(i) *SWS2total K1p(i) = K1p(i)*SWS2total K2p(i) = K2p(i)*SWS2total K3p(i) = K3p(i)*SWS2total Kb(i) = Kb(i) *SWS2total Ksi(i) = Ksi(i)*SWS2total Kw(i) = Kw(i) *SWS2total ELSE K0(i) = 1.e20_wp K1(i) = 1.e20_wp K2(i) = 1.e20_wp Kb(i) = 1.e20_wp Kw(i) = 1.e20_wp Ks(i) = 1.e20_wp Kf(i) = 1.e20_wp Kspc(i) = 1.e20_wp Kspa(i) = 1.e20_wp K1p(i) = 1.e20_wp K2p(i) = 1.e20_wp K3p(i) = 1.e20_wp Ksi(i) = 1.e20_wp Bt(i) = 1.e20_wp Ft(i) = 1.e20_wp St(i) = 1.e20_wp ENDIF END DO RETURN END SUBROUTINE constants ! ---------------------------------------------------------------------- ! VARSOLVER ! ---------------------------------------------------------------------- ! !> \file varsolver.f90 !! \BRIEF !> Module with varsolver subroutine - solve for pH and other carbonate system variables !> Solve for pH and other carbonate system variables (with input from vars routine) SUBROUTINE varsolver(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, & temp, salt, ta, tc, pt, sit, & Bt, St, Ft, & K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi, & Patm, Phydro_bar, rhodum, optGAS ) ! Purpose: Solve for pH and other carbonate system variables (with input from vars routine) ! INPUT variables: ! ================ ! temp = in situ temperature [degrees C] ! ta = total alkalinity in [eq/m^3] or [eq/kg] based on optCON in calling routine (vars) ! tc = dissolved inorganic carbon in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars) ! pt = total dissolved inorganic phosphorus in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars) ! sit = total dissolved inorganic silicon in [mol/m^3] or [mol/kg] based on optCON in calling routine (vars) ! Bt = total dissolved inorganic boron computed in calling routine (vars) ! St = total dissolved inorganic sulfur computed in calling routine (vars) ! Ft = total dissolved inorganic fluorine computed in calling routine (vars) ! K's = K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi ! Patm = atmospheric pressure [atm] ! Phydro_bar = hydrostatic pressure [bar] ! rhodum = density factor as computed in calling routine (vars) ! ----------- ! optGAS: choose in situ vs. potential fCO2 and pCO2 (default optGAS = 'Pinsitu') ! --------- ! PRESSURE & T corrections for K0 and the fugacity coefficient (Cf) ! -> 'Pzero' = 'zero order' fCO2 and pCO2 (typical approach, which is flawed) ! considers in situ T & only atm pressure (hydrostatic=0) ! -> 'Ppot' = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) ! considers potential T & only atm pressure (hydrostatic press = 0) ! -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure) ! considers in situ T & total pressure (atm + hydrostatic) ! --------- ! OUTPUT variables: ! ================= ! ph = pH on total scale ! pco2 = CO2 partial pressure (uatm) ! fco2 = CO2 fugacity (uatm) ! co2 = aqueous CO2 concentration in [mol/kg] or [mol/m^3] determined by rhodum (depends on optCON in calling routine) ! hco3 = bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] determined by rhodum ! co3 = carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] determined by rhodum ! OmegaA = Omega for aragonite, i.e., the aragonite saturation state ! OmegaC = Omega for calcite, i.e., the calcite saturation state USE mocsy_singledouble USE mocsy_phsolvers IMPLICIT NONE ! Input variables !> in situ temperature [degrees C] REAL(kind=wp), INTENT(in) :: temp !> salinity [on the practical salinity scale, dimensionless] REAL(kind=wp), INTENT(in) :: salt !> total alkalinity in [eq/m^3] OR in [eq/kg], depending on optCON in calling routine REAL(kind=wp), INTENT(in) :: ta !> dissolved inorganic carbon in [mol/m^3] OR in [mol/kg], depending on optCON in calling routine REAL(kind=wp), INTENT(in) :: tc !> phosphate concentration in [mol/m^3] OR in [mol/kg], depending on optCON in calling routine REAL(kind=wp), INTENT(in) :: pt !> total dissolved inorganic silicon concentration in [mol/m^3] OR in [mol/kg], depending on optCON in calling routine REAL(kind=wp), INTENT(in) :: sit !> total boron from either Uppstrom (1974) or Lee et al. (2010), depending on optB in calling routine REAL(kind=wp), INTENT(in) :: Bt !> total sulfate (Morris & Riley, 1966) REAL(kind=wp), INTENT(in) :: St !> total fluoride (Riley, 1965) REAL(kind=wp), INTENT(in) :: Ft !> solubility of CO2 in seawater (Weiss, 1974), also known as K0 REAL(kind=wp), INTENT(in) :: K0 !> K1 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 REAL(kind=wp), INTENT(in) :: K1 !> K2 for the dissociation of carbonic acid from Lueker et al. (2000) or Millero (2010), depending on optK1K2 REAL(kind=wp), INTENT(in) :: K2 !> equilibrium constant for dissociation of boric acid REAL(kind=wp), INTENT(in) :: Kb !> equilibrium constant for the dissociation of water (Millero, 1995) REAL(kind=wp), INTENT(in) :: Kw !> equilibrium constant for the dissociation of bisulfate (Dickson, 1990) REAL(kind=wp), INTENT(in) :: Ks !> equilibrium constant for the dissociation of hydrogen fluoride !! from Dickson and Riley (1979) or Perez and Fraga (1987), depending on optKf REAL(kind=wp), INTENT(in) :: Kf !> solubility product for calcite (Mucci, 1983) REAL(kind=wp), INTENT(in) :: Kspc !> solubility product for aragonite (Mucci, 1983) REAL(kind=wp), INTENT(in) :: Kspa !> 1st dissociation constant for phosphoric acid (Millero, 1995) REAL(kind=wp), INTENT(in) :: K1p !> 2nd dissociation constant for phosphoric acid (Millero, 1995) REAL(kind=wp), INTENT(in) :: K2p !> 3rd dissociation constant for phosphoric acid (Millero, 1995) REAL(kind=wp), INTENT(in) :: K3p !> equilibrium constant for the dissociation of silicic acid (Millero, 1995) REAL(kind=wp), INTENT(in) :: Ksi !> total atmospheric pressure [atm] REAL(kind=wp), INTENT(in) :: Patm !> total hydrostatic pressure [bar] REAL(kind=wp), INTENT(in) :: Phydro_bar !> density factor as computed incalling routine (vars) REAL(kind=wp), INTENT(in) :: rhodum !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction) !! 'Ppot' - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean !f2py character*7 optional, intent(in) :: optGAS='Pinsitu' CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS ! Output variables: !> pH on the total scale REAL(kind=wp), INTENT(out) :: ph !> CO2 partial pressure [uatm] REAL(kind=wp), INTENT(out) :: pco2 !> CO2 fugacity [uatm] REAL(kind=wp), INTENT(out) :: fco2 !> aqueous CO2* concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON REAL(kind=wp), INTENT(out) :: co2 !> bicarbonate ion (HCO3-) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON REAL(kind=wp), INTENT(out) :: hco3 !> carbonate ion (CO3--) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON REAL(kind=wp), INTENT(out) :: co3 !> Omega for aragonite, i.e., the aragonite saturation state REAL(kind=wp), INTENT(out) :: OmegaA !> Omega for calcite, i.e., the calcite saturation state REAL(kind=wp), INTENT(out) :: OmegaC ! Local variables REAL(kind=wp) :: Phydro_atm, Ptot REAL(kind=wp) :: Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff REAL(kind=wp) :: tk, tk0 real(kind=wp) :: temp68, tempot, tempot68 REAL(kind=wp) :: Hinit, H REAL(kind=wp) :: HSO4, HF, HSI, HPO4 REAL(kind=wp) :: ab, aw, ac REAL(kind=wp) :: cu, cb, cc REAL(kind=wp) :: Ca ! Array to pass optional arguments CHARACTER(7) :: opGAS IF (PRESENT(optGAS)) THEN opGAS = optGAS ELSE opGAS = 'Pinsitu' ENDIF ! Compute pH from constants and total concentrations ! - use SolveSAPHE v1.0.1 routines from Munhoven (2013, GMD) modified to use mocsy's Ks instead of its own ! 1) Compute best starting point for H+ calculation call ahini_for_at(ta, tc, Bt, K1, K2, Kb, Hinit) ! 2) Solve for H+ using above result as the initial H+ value H = solve_at_general(ta, tc, Bt, & pt, sit, & St, Ft, & K0, K1, K2, Kb, Kw, Ks, Kf, K1p, K2p, K3p, Ksi, & Hinit) ! 3) Calculate pH from from H+ concentration (mol/kg) IF (H > 0.d0) THEN pH = -1.*LOG10(H) ELSE pH = 1.e20_wp ENDIF ! Compute carbonate Alk (Ac) by difference: from total Alk and other Alk components HSO4 = St/(1.0d0 + Ks/(H/(1.0d0 + St/Ks))) HF = 1.0d0/(1.0d0 + Kf/H) HSI = 1.0d0/(1.0d0 + H/Ksi) HPO4 = (K1p*K2p*(H + 2.*K3p) - H**3) / & (H**3 + K1p*H**2 + K1p*K2p*H + K1p*K2p*K3p) ab = Bt/(1.0d0 + H/Kb) aw = Kw/H - H/(1.0d0 + St/Ks) ac = ta + hso4 - sit*hsi - ab - aw + Ft*hf - pt*hpo4 ! Calculate CO2*, HCO3-, & CO32- (in mol/kg soln) from Ct, Ac, H+, K1, & K2 cu = (2.0d0 * tc - ac) / (2.0d0 + K1 / H) cb = K1 * cu / H cc = K2 * cb / H ! When optCON = 'mol/m3' in calling routine (vars), then: ! convert output var concentrations from mol/kg to mol/m^3 ! e.g., for case when drho = 1028, multiply by [1.028 kg/L x 1000 L/m^3]) co2 = cu * rhodum hco3 = cb * rhodum co3 = cc * rhodum ! Determine CO2 fugacity [uatm] ! NOTE: equation just below requires CO2* in mol/kg fCO2 = cu * 1.e6_wp/K0 ! Determine CO2 partial pressure from CO2 fugacity [uatm] tk = 273.15d0 + temp !Compute EITHER "potential pCO2" OR "in situ pCO2" (T and P used for calculations will differ) IF (trim(opGAS) == 'Pzero' .OR. trim(opGAS) == 'pzero') THEN tk0 = tk !in situ temperature (K) for K0 calculation Ptot = Patm !total pressure (in atm) = atmospheric pressure ONLY ELSEIF (trim(opGAS) == 'Ppot' .OR. trim(opGAS) == 'ppot') THEN !Use potential temperature and atmospheric pressure (water parcel adiabatically brought back to surface) !temp68 = (temp - 0.0002d0) / 0.99975d0 !temp = in situ T; temp68 is same converted to ITPS-68 scale !tempot68 = sw_ptmp(salt, temp68, Phydro_bar*10d0, 0.0d0) !potential temperature (C) !tempot = 0.99975*tempot68 + 0.0002 !tk0 = tempot + 273.15d0 !potential temperature (K) for fugacity coeff. calc as needed for potential fCO2 & pCO2 tempot = sw_ptmp(salt, temp, Phydro_bar*10._wp, 0.0_wp) !potential temperature (C) tk0 = tempot + 273.15d0 !potential temperature (K) for fugacity coeff. calc as needed for potential fCO2 & pCO2 Ptot = Patm !total pressure (in atm) = atmospheric pressure ONLY ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN !Use in situ temperature and total pressure tk0 = tk !in situ temperature (K) for fugacity coefficient calculation Phydro_atm = Phydro_bar / 1.01325d0 !convert hydrostatic pressure from bar to atm (1.01325 bar / atm) Ptot = Patm + Phydro_atm !total pressure (in atm) = atmospheric pressure + hydrostatic pressure ELSE PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'" STOP ENDIF ! Now that we have T and P in the right form, continue with calculation of fugacity coefficient (and pCO2) Rgas_atm = 82.05736_wp ! (cm3 * atm) / (mol * K) CODATA (2006) ! To compute fugcoeff, we need 3 other terms (B, Del, xc2) in addition to 3 others above (tk, Ptot, Rgas_atm) B = -1636.75d0 + 12.0408d0*tk0 - 0.0327957d0*(tk0*tk0) + 0.0000316528d0*(tk0*tk0*tk0) Del = 57.7d0 - 0.118d0*tk0 ! "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9 ! x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite) ! Let's assume that xCO2 = fCO2. Resulting fugcoeff is identical to 8th digit after the decimal. xCO2approx = fCO2 * 1.e-6_wp IF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN ! xCO2approx = 400.0e-6_wp !a simple test (gives about same result as seacarb for pCO2insitu) ! approximate surface xCO2 ~ surface fCO2 (i.e., in situ fCO2 d by exponential pressure correction) xCO2approx = xCO2approx * exp( ((1-Ptot)*32.3_wp)/(82.05736_wp*tk0) ) ! of K0 press. correction, see Weiss (1974, equation 5) ENDIF xc2 = (1.0d0 - xCO2approx)**2 fugcoeff = exp( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk0) ) pCO2 = fCO2 / fugcoeff ! Determine Omega Calcite et Aragonite ! OmegaA = ((0.01028d0*salt/35.0d0)*cc) / Kspa ! OmegaC = ((0.01028d0*salt/35.0d0)*cc) / Kspc ! - see comments from Munhoven on the best value "0.02128" which differs slightly from the best practices guide (0.02127) Ca = (0.02128d0/40.078d0) * salt/1.80655d0 OmegaA = (Ca*cc) / Kspa OmegaC = (Ca*cc) / Kspc RETURN END SUBROUTINE varsolver ! ---------------------------------------------------------------------- ! VARS ! ---------------------------------------------------------------------- ! !> \file vars.f90 !! \BRIEF !> Module with vars subroutine - compute carbonate system vars from DIC,Alk,T,S,P,nuts !> Computes standard carbonate system variables (pH, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R) !! as 1D arrays FROM !! temperature, salinity, pressure, !! total alkalinity (ALK), dissolved inorganic carbon (DIC), !! silica and phosphate concentrations (all 1-D arrays) SUBROUTINE vars(ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis, & temp, sal, alk, dic, sil, phos, Patm, depth, lat, N, & optCON, optT, optP, optB, optK1K2, optKf, optGAS ) ! Purpose: ! Computes other standard carbonate system variables (pH, CO2*, HCO3- and CO32-, OmegaA, OmegaC, R) ! as 1D arrays ! FROM: ! temperature, salinity, pressure, ! total alkalinity (ALK), dissolved inorganic carbon (DIC), ! silica and phosphate concentrations (all 1-D arrays) ! INPUT variables: ! ================ ! Patm = atmospheric pressure [atm] ! depth = depth [m] (with optP='m', i.e., for a z-coordinate model vertical grid is depth, not pressure) ! = pressure [db] (with optP='db') ! lat = latitude [degrees] (needed to convert depth to pressure, i.e., when optP='m') ! = dummy array (unused when optP='db') ! temp = potential temperature [degrees C] (with optT='Tpot', i.e., models carry tempot, not in situ temp) ! = in situ temperature [degrees C] (with optT='Tinsitu', e.g., for data) ! sal = salinity in [psu] ! alk = total alkalinity in [eq/m^3] with optCON = 'mol/m3' ! = [eq/kg] with optCON = 'mol/kg' ! dic = dissolved inorganic carbon [mol/m^3] with optCON = 'mol/m3' ! = [mol/kg] with optCON = 'mol/kg' ! sil = silica [mol/m^3] with optCON = 'mol/m3' ! = [mol/kg] with optCON = 'mol/kg' ! phos = phosphate [mol/m^3] with optCON = 'mol/m3' ! = [mol/kg] with optCON = 'mol/kg' ! INPUT options: ! ============== ! ----------- ! optCON: choose input & output concentration units - mol/kg (data) vs. mol/m^3 (models) ! ----------- ! -> 'mol/kg' for DIC, ALK, sil, & phos given on mokal scale, i.e., in mol/kg (std DATA units) ! -> 'mol/m3' for DIC, ALK, sil, & phos given in mol/m^3 (std MODEL units) ! ----------- ! optT: choose in situ vs. potential temperature as input ! --------- ! NOTE: Carbonate chem calculations require IN-SITU temperature (not potential Temperature) ! -> 'Tpot' means input is pot. Temperature (in situ Temp "tempis" is computed) ! -> 'Tinsitu' means input is already in-situ Temperature, not pot. Temp ("tempis" not computed) ! --------- ! optP: choose depth (m) vs pressure (db) as input ! --------- ! -> 'm' means "depth" input is in "m" (thus in situ Pressure "p" [db] is computed) ! -> 'db' means "depth" input is already in situ pressure [db], not m (thus p = depth) ! --------- ! optB: choose total boron formulation - Uppström (1974) vs. Lee et al. (2010) ! --------- ! -> 'u74' means use classic formulation of Uppström (1974) for total Boron ! -> 'l10' means use newer formulation of Lee et al. (2010) for total Boron ! --------- ! optK1K2: ! --------- ! -> 'l' means use Lueker et al. (2000) formulations for K1 & K2 (recommended by Dickson et al. 2007) ! **** BUT this should only be used when 2 < T < 35 and 19 < S < 43 ! -> 'm10' means use Millero (2010) formulation for K1 & K2 (see Dickson et al., 2007) ! **** Valid for 0 < T < 50°C and 1 < S < 50 psu ! ---------- ! optKf: ! ---------- ! -> 'pf' means use Perez & Fraga (1987) formulation for Kf (recommended by Dickson et al., 2007) ! **** BUT Valid for 9 < T < 33°C and 10 < S < 40. ! -> 'dg' means use Dickson & Riley (1979) formulation for Kf (recommended by Dickson & Goyet, 1994) ! ----------- ! optGAS: choose in situ vs. potential fCO2 and pCO2 ! --------- ! PRESSURE corrections for K0 and the fugacity coefficient (Cf) ! -> 'Pzero' = 'zero order' fCO2 and pCO2 (typical approach, which is flawed) ! considers in situ T & only atm pressure (hydrostatic=0) ! -> 'Ppot' = 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) ! considers potential T & only atm pressure (hydrostatic press = 0) ! -> 'Pinsitu' = 'in situ' fCO2 and pCO2 (accounts for huge effects of pressure) ! considers in situ T & total pressure (atm + hydrostatic) ! --------- ! OUTPUT variables: ! ================= ! ph = pH on total scale ! pco2 = CO2 partial pressure (uatm) ! fco2 = CO2 fugacity (uatm) ! co2 = aqueous CO2 concentration in [mol/kg] or [mol/m^3] depending on optCON ! hco3 = bicarbonate (HCO3-) concentration in [mol/kg] or [mol/m^3] depending on optCON ! co3 = carbonate (CO3--) concentration in [mol/kg] or [mol/m^3] depending on optCON ! OmegaA = Omega for aragonite, i.e., the aragonite saturation state ! OmegaC = Omega for calcite, i.e., the calcite saturation state ! BetaD = Revelle factor dpCO2/pCO2 / dDIC/DIC ! rhoSW = in-situ density of seawater; rhoSW = f(s, t, p) ! p = pressure [decibars]; p = f(depth, latitude) if computed from depth [m] OR p = depth if [db] ! tempis = in-situ temperature [degrees C] USE mocsy_singledouble IMPLICIT NONE ! Input variables !> number of records INTEGER, INTENT(in) :: N !> either in situ temperature (when optT='Tinsitu', typical data) !! OR potential temperature (when optT='Tpot', typical models) [degree C] REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp !> salinity [psu] REAL(kind=wp), INTENT(in), DIMENSION(N) :: sal !> total alkalinity in [eq/m^3] (when optCON = 'mol/m3') OR in [eq/kg] (when optCON = 'mol/kg') REAL(kind=wp), INTENT(in), DIMENSION(N) :: alk !> dissolved inorganic carbon in [mol/m^3] (when optCON = 'mol/m3') OR in [mol/kg] (when optCON = 'mol/kg') REAL(kind=wp), INTENT(in), DIMENSION(N) :: dic !> SiO2 concentration in [mol/m^3] (when optCON = 'mol/m3') OR in [mol/kg] (when optCON = 'mol/kg') REAL(kind=wp), INTENT(in), DIMENSION(N) :: sil !> phosphate concentration in [mol/m^3] (when optCON = 'mol/m3') OR in [mol/kg] (when optCON = 'mol/kg') REAL(kind=wp), INTENT(in), DIMENSION(N) :: phos !f2py optional , depend(sal) :: n=len(sal) !> atmospheric pressure [atm] REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm !> depth in \b meters (when optP='m') or \b decibars (when optP='db') REAL(kind=wp), INTENT(in), DIMENSION(N) :: depth !> latitude [degrees north] REAL(kind=wp), INTENT(in), DIMENSION(N) :: lat !> choose either \b 'mol/kg' (std DATA units) or \b 'mol/m3' (std MODEL units) to select !! concentration units for input (for alk, dic, sil, phos) & output (co2, hco3, co3) CHARACTER(6), INTENT(in) :: optCON !> choose \b 'Tinsitu' for in situ temperature or \b 'Tpot' for potential temperature (in situ Temp is computed, needed for models) CHARACTER(7), INTENT(in) :: optT !> for depth input, choose \b "db" for decibars (in situ pressure) or \b "m" for meters (pressure is computed, needed for models) CHARACTER(2), INTENT(in) :: optP !> for total boron, choose either \b 'u74' (Uppstrom, 1974) or \b 'l10' (Lee et al., 2010). !! The 'l10' formulation is based on 139 measurements (instead of 20), !! uses a more accurate method, and !! generally increases total boron in seawater by 4% !f2py character*3 optional, intent(in) :: optB='l10' CHARACTER(3), OPTIONAL, INTENT(in) :: optB !> for Kf, choose either \b 'pf' (Perez & Fraga, 1987) or \b 'dg' (Dickson & Riley, 1979) !f2py character*2 optional, intent(in) :: optKf='pf' CHARACTER(2), OPTIONAL, INTENT(in) :: optKf !> for K1,K2 choose either \b 'l' (Lueker et al., 2000) or \b 'm10' (Millero, 2010) !f2py character*3 optional, intent(in) :: optK1K2='l' CHARACTER(3), OPTIONAL, INTENT(in) :: optK1K2 !> for K0,fugacity coefficient choose either \b 'Ppot' (no pressure correction) or \b 'Pinsitu' (with pressure correction) !! 'Ppot' - for 'potential' fCO2 and pCO2 (water parcel brought adiabatically to the surface) !! 'Pinsitu' - for 'in situ' values of fCO2 and pCO2, accounting for pressure on K0 and Cf !! with 'Pinsitu' the fCO2 and pCO2 will be many times higher in the deep ocean !f2py character*7 optional, intent(in) :: optGAS='Pinsitu' CHARACTER(7), OPTIONAL, INTENT(in) :: optGAS ! Output variables: !> pH on the total scale REAL(kind=wp), INTENT(out), DIMENSION(N) :: ph !> CO2 partial pressure [uatm] REAL(kind=wp), INTENT(out), DIMENSION(N) :: pco2 !> CO2 fugacity [uatm] REAL(kind=wp), INTENT(out), DIMENSION(N) :: fco2 !> aqueous CO2* concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON REAL(kind=wp), INTENT(out), DIMENSION(N) :: co2 !> bicarbonate ion (HCO3-) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON REAL(kind=wp), INTENT(out), DIMENSION(N) :: hco3 !> carbonate ion (CO3--) concentration, either in [mol/m^3] or [mol/kg] depending on choice for optCON REAL(kind=wp), INTENT(out), DIMENSION(N) :: co3 !> Omega for aragonite, i.e., the aragonite saturation state REAL(kind=wp), INTENT(out), DIMENSION(N) :: OmegaA !> Omega for calcite, i.e., the calcite saturation state REAL(kind=wp), INTENT(out), DIMENSION(N) :: OmegaC !> Revelle factor, i.e., dpCO2/pCO2 / dDIC/DIC REAL(kind=wp), INTENT(out), DIMENSION(N) :: BetaD !> in-situ density of seawater; rhoSW = f(s, t, p) in [kg/m3] REAL(kind=wp), INTENT(out), DIMENSION(N) :: rhoSW !> pressure [decibars]; p = f(depth, latitude) if computed from depth [m] (when optP='m') OR p = depth [db] (when optP='db') REAL(kind=wp), INTENT(out), DIMENSION(N) :: p !> in-situ temperature \b [degrees C] REAL(kind=wp), INTENT(out), DIMENSION(N) :: tempis ! Local variables REAL(kind=wp) :: ssal, salk, sdic, ssil, sphos REAL(kind=wp) :: tempot, tempis68, tempot68 REAL(kind=wp) :: drho REAL(kind=wp) :: K0, K1, K2, Kb, Kw, Ks, Kf, Kspc REAL(kind=wp) :: Kspa, K1p, K2p, K3p, Ksi REAL(kind=wp) :: St, Ft, Bt REAL(kind=wp), DIMENSION(1) :: aK0, aK1, aK2, aKb, aKw, aKs, aKf, aKspc REAL(kind=wp), DIMENSION(1) :: aKspa, aK1p, aK2p, aK3p, aKsi REAL(kind=wp), DIMENSION(1) :: aSt, aFt, aBt REAL(kind=wp) :: Patmd, Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff REAL(kind=wp) :: Phydro_atm INTEGER :: i, icount REAL(kind=wp) :: t, tk, prb REAL(kind=wp) :: s REAL(kind=wp) :: tc, ta REAL(kind=wp) :: sit, pt REAL(kind=wp) :: Hinit REAL(kind=wp) :: ah1 REAL(kind=wp) :: HSO4, HF, HSI, HPO4 REAL(kind=wp) :: ab, aw, ac, ah2, erel REAL(kind=wp) :: cu, cb, cc REAL(kind=wp), DIMENSION(2) :: dicdel, pco2del REAL(kind=wp) :: dx, Rf REAL(kind=wp) :: dph, dpco2, dfco2, dco2, dhco3, dco3, dOmegaA, dOmegaC INTEGER :: kcomp INTEGER :: j, minusplus ! Arrays to pass optional arguments into or use defaults (Dickson et al., 2007) CHARACTER(3) :: opB CHARACTER(2) :: opKf CHARACTER(3) :: opK1K2 CHARACTER(7) :: opGAS ! Set defaults for optional arguments (in Fortran 90) ! Note: Optional arguments with f2py (python) are set above with ! the !f2py statements that precede each type declaraion IF (PRESENT(optB)) THEN ! print *,"optB present:" ! print *,"optB = ", optB opB = optB ELSE ! Default is Lee et al (2010) for total boron ! print *,"optB NOT present:" opB = 'l10' ! print *,"opB = ", opB ENDIF IF (PRESENT(optKf)) THEN ! print *,"optKf = ", optKf opKf = optKf ELSE ! print *,"optKf NOT present:" ! Default is Perez & Fraga (1987) for Kf opKf = 'pf' ! print *,"opKf = ", opKf ENDIF IF (PRESENT(optK1K2)) THEN ! print *,"optK1K2 = ", optK1K2 opK1K2 = optK1K2 ELSE ! print *,"optK1K2 NOT present:" ! Default is Lueker et al. 2000) for K1 & K2 opK1K2 = 'l' ! print *,"opK1K2 = ", opK1K2 ENDIF IF (PRESENT(optGAS)) THEN opGAS = optGAS ELSE opGAS = 'Pinsitu' ENDIF icount = 0 DO i = 1, N icount = icount + 1 ! =============================================================== ! Convert model depth -> press; convert model Theta -> T in situ ! =============================================================== ! * Model temperature tracer is usually "potential temperature" ! * Model vertical grid is usually in meters ! BUT carbonate chem routines require pressure & in-situ T ! Thus before computing chemistry, if appropriate, ! convert these 2 model vars (input to this routine) ! - depth [m] => convert to pressure [db] ! - potential temperature (C) => convert to in-situ T (C) ! ------------------------------------------------------- ! 1) Compute pressure [db] from depth [m] and latitude [degrees] (if input is m, for models) !print *,"optP =", optP, "end" IF (trim(optP) == 'm' ) THEN ! Compute pressure [db] from depth [m] and latitude [degrees] p(i) = p80(depth(i), lat(i)) ELSEIF (trim(optP) == 'db') THEN ! In this case (where optP = 'db'), p is input & output (no depth->pressure conversion needed) p(i) = depth(i) ELSE !print *,"optP =", optP, "end" PRINT *,"optP must be either 'm' or 'db'" STOP ENDIF ! 2) Convert potential T to in-situ T (if input is Tpot, i.e. case for models): IF (trim(optT) == 'Tpot' .OR. trim(optT) == 'tpot') THEN tempot = temp(i) ! This is the case for most models and some data ! a) Convert the pot. temp on today's "ITS 90" scale to older IPTS 68 scale ! (see Dickson et al., Best Practices Guide, 2007, Chap. 5, p. 7, including footnote) tempot68 = (tempot - 0.0002) / 0.99975 ! b) Compute "in-situ Temperature" from "Potential Temperature" (both on IPTS 68) tempis68 = sw_temp(sal(i), tempot68, p(i), 0.0_wp ) ! c) Convert the in-situ temp on older IPTS 68 scale to modern scale (ITS 90) tempis(i) = 0.99975*tempis68 + 0.0002 ! Note: parts (a) and (c) above are tiny corrections; ! part (b) is a big correction for deep waters (but zero at surface) ELSEIF (trim(optT) == 'Tinsitu' .OR. trim(optT) == 'tinsitu') THEN ! When optT = 'Tinsitu', tempis is input & output (no tempot needed) tempis(i) = temp(i) tempis68 = (temp(i) - 0.0002) / 0.99975 ! dtempot68 = sw_ptmp(DBLE(sal(i)), DBLE(tempis68), DBLE(p), 0.0d0) ! dtempot = 0.99975*dtempot68 + 0.0002 ELSE PRINT *,"optT must be either 'Tpot' or 'Tinsitu'" PRINT *,"you specified optT =", trim(optT) STOP ENDIF ! ================================================================ ! Carbonate chemistry computations ! ================================================================ IF (dic(i) > 0. .AND. dic(i) < 1.0e+4) THEN ! Test to indicate if any of input variables are unreasonable IF ( sal(i) < 0. & .OR. alk(i) < 0. & .OR. dic(i) < 0. & .OR. sil(i) < 0. & .OR. phos(i) < 0. & .OR. sal(i) > 1e+3 & .OR. alk(i) > 1e+3 & .OR. dic(i) > 1e+3 & .OR. sil(i) > 1e+3 & .OR. phos(i) > 1e+3) THEN PRINT *, 'i, icount, tempot, sal, alk, dic, sil, phos =', & i, icount, tempot, sal(i), alk(i), dic(i), sil(i), phos(i) ENDIF ! Zero out any negative salinity, phosphate, silica, dic, and alk IF (sal(i) < 0.0) THEN ssal = 0.0 ELSE ssal = sal(i) ENDIF IF (phos(i) < 0.0) THEN sphos = 0.0 ELSE sphos = phos(i) ENDIF IF (sil(i) < 0.0) THEN ssil = 0.0 ELSE ssil = sil(i) ENDIF IF (dic(i) < 0.0) THEN sdic = 0.0 ELSE sdic = dic(i) ENDIF IF (alk(i) < 0.0) THEN salk = 0.0 ELSE salk = alk(i) ENDIF ! Absolute temperature (Kelvin) & related variables t = DBLE(tempis(i)) tk = 273.15d0 + t ! Atmospheric pressure Patmd = DBLE(Patm(i)) ! Hydrostatic pressure (prb is in bars) prb = DBLE(p(i)) / 10.0d0 Phydro_atm = prb / 1.01325d0 ! convert hydrostatic pressure from bar to atm (1.01325 bar / atm) ! Total pressure [atm] IF (trim(opGAS) == 'Pzero' .OR. trim(opGAS) == 'pzero') THEN Ptot = Patmd ! total pressure (in atm) = atmospheric pressure ONLY ELSEIF (trim(opGAS) == 'Ppot' .OR. trim(opGAS) == 'ppot') THEN Ptot = Patmd ! total pressure (in atm) = atmospheric pressure ONLY ELSEIF (trim(opGAS) == 'Pinsitu' .OR. trim(opGAS) == 'pinsitu') THEN Ptot = Patmd + Phydro_atm ! total pressure (in atm) = atmospheric pressure + hydrostatic pressure ELSE PRINT *, "optGAS must be 'Pzero', 'Ppot', or 'Pinsitu'" STOP ENDIF ! Salinity (equivalent array in double precision) s = DBLE(ssal) ! Get all equilibrium constants and total concentrations of SO4, F, B CALL constants(aK0, aK1, aK2, aKb, aKw, aKs, aKf, aKspc, aKspa, & aK1p, aK2p, aK3p, aKsi, & aSt, aFt, aBt, & temp(i), sal(i), Patm(i), & depth(i), lat(i), 1, & optT, optP, opB, opK1K2, opKf, opGAS ) ! Unlike f77, in F90 we can't assign an array (dimen=1) to a scalar in a routine argument ! Thus, set scalar constants equal to array (dimension=1) values required as arguments K0 = aK0(1) ; K1 = aK1(1) ; K2 = aK2(1) ; Kb = aKb(1) ; Kw = aKw(1) Ks = aKs(1) ; Kf = aKs(1) ; Kspc = aKspc(1) ; Kspa = aKspa(1) K1p = aK1p(1) ; K2p = aK2p(1) ; K3p = aK3p(1) ; Ksi = aKsi(1) St = aSt(1) ; Ft = aFt(1) ; Bt = aBt(1) ! Compute in-situ density [kg/m^3] rhoSW(i) = rho(ssal, tempis68, prb) ! Either convert units of DIC and ALK (MODEL case) or not (DATA case) IF (trim(optCON) == 'mol/kg') THEN ! No conversion: ! print *,'DIC and ALK already given in mol/kg (std DATA units)' drho = 1. ELSEIF (trim(optCON) == 'mol/m3') THEN ! Do conversion: ! print *,"DIC and ALK given in mol/m^3 (std MODEL units)" drho = DBLE(rhoSW(i)) ELSE PRINT *,"optCON must be either 'mol/kg' or 'mol/m3'" STOP ENDIF tc = DBLE(sdic)/drho ta = DBLE(salk)/drho sit = DBLE(ssil)/drho pt = DBLE(sphos)/drho ! Solve for pH and all other variables ! ------------------------------------ CALL varsolver(dph, dpco2, dfco2, dco2, dhco3, dco3, dOmegaA, dOmegaC, & t, s, ta, tc, pt, sit, & Bt, St, Ft, & K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi, & Patmd, prb, drho, opGAS ) ! Convert all output variables from double to single precision pH(i) = REAL(dph) co2(i) = REAL(dco2) hco3(i) = REAL(dhco3) co3(i) = REAL(dco3) fCO2(i) = REAL(dfCO2) pCO2(i) = REAL(dpCO2) OmegaA(i) = REAL(dOmegaA) OmegaC(i) = REAL(dOmegaC) ! Compute Revelle factor numerically (derivative using centered-difference scheme) DO j=1,2 minusplus = (-1)**j dx = 0.1 * 1e-6 ! Numerical tests found for DIC that optimal dx = 0.1 umol/kg (0.1e-6 mol/kg) dicdel(j) = tc + DBLE(minusplus)*dx/2.0d0 CALL varsolver(dph, dpco2, dfco2, dco2, dhco3, dco3, dOmegaA, dOmegaC, & t, s, ta, dicdel(j), pt, sit, & Bt, St, Ft, & K0, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi, & Patmd, prb, drho, optGAS ) pco2del(j) = dpco2 END DO !Classic finite centered difference formula for derivative (2nd order accurate) Rf = (pco2del(2) - pco2del(1)) / (dicdel(2) - dicdel(1)) ! dpCO2/dDIC !Rf = (pco2del(2) - pco2del(1)) / (dx) ! dpCO2/dDIC (same as just above) Rf = Rf * tc / dpco2 ! R = (dpCO2/dDIC) * (DIC/pCO2) BetaD(i) = REAL(Rf) ELSE ph(i) = 1.e20_wp pco2(i) = 1.e20_wp fco2(i) = 1.e20_wp co2(i) = 1.e20_wp hco3(i) = 1.e20_wp co3(i) = 1.e20_wp OmegaA(i) = 1.e20_wp OmegaC(i) = 1.e20_wp BetaD(i) = 1.e20_wp rhoSW(i) = 1.e20_wp p(i) = 1.e20_wp tempis(i) = 1.e20_wp ENDIF END DO RETURN END SUBROUTINE vars ! ---------------------------------------------------------------------- ! P2FCO2 ! ---------------------------------------------------------------------- ! !> \file p2fCO2.f90 !! \BRIEF !> Module with p2fCO2 subroutine - compute fCO2 from pCO2, in situ T, atm pressure, hydrostatic pressure !> Compute fCO2 from arrays of pCO2, in situ temp, atm pressure, & hydrostatic pressure SUBROUTINE p2fCO2(pCO2, temp, Patm, p, N, fCO2) ! Purpose: ! Compute fCO2 from arrays of pCO2, in situ temp, atm pressure, & hydrostatic pressure USE mocsy_singledouble IMPLICIT NONE !> number of records INTEGER, INTENT(in) :: N ! INPUT variables !> oceanic partial pressure of CO2 [uatm] REAL(kind=wp), INTENT(in), DIMENSION(N) :: pCO2 !> in situ temperature [C] REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp !> atmospheric pressure [atm] REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm !> hydrostatic pressure [db] REAL(kind=wp), INTENT(in), DIMENSION(N) :: p ! OUTPUT variables: !> fugacity of CO2 [uatm] REAL(kind=wp), INTENT(out), DIMENSION(N) :: fCO2 ! LOCAL variables: REAL(kind=wp) :: dpCO2, dtemp, tk, dPatm, prb REAL(kind=wp) :: Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff REAL(kind=wp) :: dfCO2 INTEGER :: i ! REAL(kind=wp) :: sw_ptmp ! EXTERNAL sw_ptmp DO i = 1,N dpCO2 = DBLE(pCO2(i)) dtemp = DBLE(temp(i)) dPatm = DBLE(Patm(i)) tk = 273.15d0 + DBLE(temp(i)) !Absolute temperature (Kelvin) prb = DBLE(p(i)) / 10.0d0 !Pressure effect (prb is in bars) Ptot = dPatm + prb/1.01325d0 !Total pressure (atmospheric + hydrostatic) [atm] Rgas_atm = 82.05736_wp !R in (cm3 * atm) / (mol * K) from CODATA (2006) ! To compute fugcoeff, we need 3 other terms (B, Del, xc2) as well as 3 others above (tk, Ptot, Rgas_atm) B = -1636.75d0 + 12.0408d0*tk - 0.0327957d0*(tk*tk) + 0.0000316528d0*(tk*tk*tk) Del = 57.7d0 - 0.118d0*tk ! "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9 ! x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite) ! Let's assume that xCO2 = pCO2. Resulting fugcoeff is identical to 8th digit after the decimal. xCO2approx = dpCO2 * 1.e-6_wp xc2 = (1.0d0 - xCO2approx)**2 fugcoeff = EXP( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk) ) dfCO2 = dpCO2 * fugcoeff fCO2(i) = REAL(dfCO2) END DO RETURN END SUBROUTINE p2fCO2 ! ---------------------------------------------------------------------- ! P2FCO2 ! ---------------------------------------------------------------------- ! !> \file f2pCO2.f90 !! \BRIEF !> Module with f2pCO2 subroutine - compute pCO2 from fCO2, in situ T, atm pressure, hydrostatic pressure !> Compute pCO2 from arrays of fCO2, in situ temp, atm pressure, & hydrostatic pressure SUBROUTINE f2pCO2(fCO2, temp, Patm, p, N, pCO2) ! Purpose: ! Compute pCO2 from arrays of fCO2, in situ temp, atm pressure, & hydrostatic pressure USE mocsy_singledouble IMPLICIT NONE !> number of records INTEGER, intent(in) :: N ! INPUT variables !> oceanic fugacity of CO2 [uatm] REAL(kind=wp), INTENT(in), DIMENSION(N) :: fCO2 !> in situ temperature [C] REAL(kind=wp), INTENT(in), DIMENSION(N) :: temp !> atmospheric pressure [atm] REAL(kind=wp), INTENT(in), DIMENSION(N) :: Patm !> hydrostatic pressure [db] REAL(kind=wp), INTENT(in), DIMENSION(N) :: p ! OUTPUT variables: !> oceanic partial pressure of CO2 [uatm] REAL(kind=wp), INTENT(out), DIMENSION(N) :: pCO2 ! LOCAL variables: REAL(kind=wp) :: dfCO2, dtemp, tk, dPatm, prb REAL(kind=wp) :: Ptot, Rgas_atm, B, Del, xCO2approx, xc2, fugcoeff REAL(kind=wp) :: dpCO2 INTEGER :: i ! REAL(kind=wp) :: sw_ptmp ! EXTERNAL sw_ptmp DO i = 1,N dfCO2 = DBLE(fCO2(i)) dtemp = DBLE(temp(i)) dPatm = DBLE(Patm(i)) tk = 273.15d0 + DBLE(temp(i)) !Absolute temperature (Kelvin) prb = DBLE(p(i)) / 10.0d0 !Pressure effect (prb is in bars) Ptot = dPatm + prb/1.01325d0 !Total pressure (atmospheric + hydrostatic) [atm] Rgas_atm = 82.05736_wp !R in (cm3 * atm) / (mol * K) from CODATA (2006) ! To compute fugcoeff, we need 3 other terms (B, Del, xc2) as well as 3 others above (tk, Ptot, Rgas_atm) B = -1636.75d0 + 12.0408d0*tk - 0.0327957d0*(tk*tk) + 0.0000316528d0*(tk*tk*tk) Del = 57.7d0 - 0.118d0*tk ! "x2" term often neglected (assumed = 1) in applications of Weiss's (1974) equation 9 ! x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite) ! Let's assume that xCO2 = fCO2. Resulting fugcoeff is identical to 8th digit after the decimal. xCO2approx = dfCO2 * 1.e-6_wp xc2 = (1.0d0 - xCO2approx)**2 fugcoeff = exp( Ptot*(B + 2.0d0*xc2*Del)/(Rgas_atm*tk) ) dpCO2 = dfCO2 / fugcoeff pCO2(i) = REAL(dpCO2) END DO RETURN END SUBROUTINE f2pCO2 END MODULE mocsy_mainmod