Back to home page

MITgcm
 
 

    

File indexing completed on 2023-10-26 05:10:19 UTC

view on githubraw file Latest commit a2c35952 on 2023-10-25 14:50:13 UTC
6acab690ae Jona*0001 #include "DIC_OPTIONS.h"
                0002 
                0003 C--  File dic_solvesaphe.F:
                0004 C--   Contents:
                0005 C--   o AHINI_FOR_AT
                0006 C--   o ANW_INFSUP
                0007 C--   o EQUATION_AT
                0008 C--   o CALC_PCO2_SOLVESAPHE
                0009 C--   o DIC_COEFFS_SURF
                0010 C--   o DIC_COEFFS_DEEP
                0011 C--   o SOLVE_AT_FAST
                0012 C--   o SOLVE_AT_GENERAL
                0013 C--   o SOLVE_AT_GENERAL_SEC
                0014 
                0015 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0016 C    SolveSAPHE is free software: you can redistribute it and/or modify
                0017 C    it under the terms of the GNU Lesser General Public License as published by
                0018 C    the Free Software Foundation, either version 3 of the License, or
                0019 C    (at your option) any later version.
                0020 C
                0021 C    SolveSAPHE is distributed in the hope that it will be useful,
                0022 C    but WITHOUT ANY WARRANTY; without even the implied warranty of
                0023 C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
                0024 C    GNU Lesser General Public License for more details.
                0025 C
                0026 C    You should have received a copy of the GNU Lesser General Public License
                0027 C    along with SolveSAPHE.  If not, see <http://www.gnu.org/licenses/>.
                0028 C
                0029 C    Copyright 2013 Guy Munhoven
                0030 C
                0031 C    From the paper: Munhoven, G (2013),
                0032 C        Mathematics of the total alkalinity–pH equation – pathway to robust
                0033 C        and universal solution algorithms: the SolveSAPHE package v1.0.1.
                0034 C      Geosci. Model Dev., 6, 1367–1388,  doi:10.5194/gmd-6-1367-2013
                0035 C
                0036 C    JMLauderdale Summer 2018:
                0037 C    - Modified for MITGCM style, just keeping the "general" equations from
                0038 C        MOD_PHSOLVERS SOLVE_AT_GENERAL SOLVE_AT_GENERAL_SEC and SOLVE_AT_FAST.
                0039 C        Use runtime flags in data.dic to use GENERAL (selectPHsolver=1,
                0040 C        default), SEC (selectPHsolver=2) or FAST (selectPHsolver=3).
                0041 C
                0042 C    - MOD_CHEMCONST is included here as DIC_COEFFS_SURF, with the
                0043 C        style brought in line with S/R CARBON_COEFF (i.e all ak values are
                0044 C        filled in one call, rather than all as separate functions).
                0045 C
                0046 C    - MOD_CHEMCONST pressure corrections have been split off into the new
                0047 C        S/R DIC_COEFFS_DEEP, call both to get pressure
                0048 C        adjusted dissociation constants
                0049 C
                0050 C    - Different coefficients can be accessed using runtime flags in data.dic:
                0051 C
                0052 C    Borate concentration from salinity:
                0053 C     selectBTconst=1 for the default formulation of Uppström (1974)
                0054 C                     i.e. the same as S/R CARBON_COEFFS
                0055 C     selectBTconst=2 for the new formulation from Lee et al (2010)
                0056 C
                0057 C    Fluoride concentration from salinity:
                0058 C     selectFTconst=1 for the default formulation of Riley (1965)
                0059 C                     i.e. the same as S/R CARBON_COEFFS
                0060 C     selectFTconst=2 for the new formulation from Culkin (1965)
                0061 C
                0062 C    First dissociation constant for hydrogen fluoride:
                0063 C     selectHFconst=1 for the default  Dickson and Riley (1979)
                0064 C                     i.e. the same as S/R CARBON_COEFFS
                0065 C     selectHFconst=2 for the formulation of Perez and Fraga (1987)
                0066 C
                0067 C    First and second dissociation constants of carbonic acid:
                0068 C     selectK1K2const=1 for the default formulation of Millero (1995) with data
                0069 C              from Mehrbach et al. (1973),  i.e. the same as S/R CARBON_COEFFS
                0070 C     selectK1K2const=2 for the formulation of Roy et al. (1993)
                0071 C     selectK1K2const=3 for the "combination" formulation of Millero (1995)
                0072 C     selectK1K2const=4 for the formulation of Luecker et al. (2000)
                0073 C     selectK1K2const=5 for the formulation of Millero
                0074 C                                              (2010, Mar. Fresh Wat. Res.)
                0075 C     selectK1K2const=6 for the formulation of Waters, Millero, Woosley
                0076 C                                              (2014, Mar. Chem.)
                0077 C      =========================================================================
                0078       SUBROUTINE AHINI_FOR_AT(p_alkcb, p_dictot, p_bortot, p_hini,
                0079      &                    i, j, k, bi, bj, myIter, myThid )
                0080 
                0081 C      Subroutine returns the root for the 2nd order approximation of the
                0082 C      DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic polynomial)
                0083 C      around the local minimum, if it exists.
                0084 
                0085 C      Returns * 1E-03 if p_alkcb <= 0
                0086 C              * 1E-10 if p_alkcb >= 2*p_dictot + p_bortot
                0087 C              * 1E-07 if 0 < p_alkcb < 2*p_dictot + p_bortot
                0088 C                       and the 2nd order approximation does not have a solution
                0089 
                0090       IMPLICIT NONE
                0091 
                0092 C     MITgcm GLobal variables
                0093 #include "SIZE.h"
                0094 #include "EEPARAMS.h"
                0095 #include "DIC_VARS.h"
                0096 
                0097 C -----------------------------------------------------------------------
                0098 C     Argument variables
                0099 C -----------------------------------------------------------------------
                0100 
                0101       _RL p_alkcb
                0102       _RL p_dictot
                0103       _RL p_bortot
                0104       _RL p_hini
                0105       INTEGER i,j,k,bi,bj,myIter,myThid
                0106 
                0107 #ifdef CARBONCHEM_SOLVESAPHE
                0108 C -----------------------------------------------------------------------
                0109 C     Local variables
                0110 C -----------------------------------------------------------------------
                0111 
                0112       _RL zcar, zbor
                0113       _RL zd, zsqrtd, zhmin
                0114       _RL za2, za1, za0
                0115 
                0116       IF (p_alkcb .LE. 0. _d 0) THEN
                0117        p_hini = 1. _d -3
                0118       ELSEIF (p_alkcb .GE. (2. _d 0*p_dictot + p_bortot)) THEN
                0119        p_hini = 1. _d -10
                0120       ELSE
                0121        zcar = p_dictot/p_alkcb
                0122        zbor = p_bortot/p_alkcb
                0123 
                0124 C      Coefficients of the cubic polynomial
                0125        za2 = akb(i,j,bi,bj)*(1. _d 0 - zbor) + ak1(i,j,bi,bj)
                0126      &       *(1. _d 0-zcar)
                0127        za1 = ak1(i,j,bi,bj)*akb(i,j,bi,bj)*(1. _d 0 - zbor - zcar)
                0128      &     + ak1(i,j,bi,bj)*ak2(i,j,bi,bj)*(1. _d 0 - (zcar+zcar))
                0129        za0 = ak1(i,j,bi,bj)*ak2(i,j,bi,bj)*akb(i,j,bi,bj)
                0130      &       *(1. _d 0 - zbor - (zcar+zcar))
                0131 
                0132 C      Taylor expansion around the minimum
                0133        zd = za2*za2 - 3. _d 0*za1
                0134 C            Discriminant of the quadratic equation
                0135 C                 for the minimum close to the root
                0136 
                0137        IF(zd .GT. 0. _d 0) THEN
                0138 C      If the discriminant is positive
                0139          zsqrtd = SQRT(zd)
                0140          IF(za2 .LT. 0) THEN
                0141            zhmin = (-za2 + zsqrtd)/3. _d 0
                0142          ELSE
                0143            zhmin = -za1/(za2 + zsqrtd)
                0144         ENDIF
                0145         p_hini = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))
                0146      &            /zsqrtd)
                0147        ELSE
                0148          p_hini = 1. _d -7
                0149        ENDIF
                0150       ENDIF
                0151 
                0152 #endif /* CARBONCHEM_SOLVESAPHE */
                0153       RETURN
                0154       END
                0155 C      END SUBROUTINE AHINI_FOR_AT
                0156 C      =========================================================================
                0157 
                0158 C      =========================================================================
                0159       SUBROUTINE ANW_INFSUP(p_dictot, p_bortot,
                0160      &                      p_po4tot, p_siltot,
                0161      &                      p_nh4tot, p_h2stot,
                0162      &                      p_so4tot, p_flutot,
                0163      &                      p_alknw_inf, p_alknw_sup,
                0164      &                      i, j, k, bi, bj, myIter,
                0165      &                      myThid )
                0166 
                0167 C      Subroutine returns the lower and upper bounds of "non-water-selfionization"
                0168 C      Contributions to total alkalinity (the infimum and the supremum), i.e
                0169 C      inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
                0170 
                0171       IMPLICIT NONE
                0172 
                0173 C     MITgcm GLobal variables
                0174 #include "SIZE.h"
                0175 #include "EEPARAMS.h"
                0176 #include "DIC_VARS.h"
                0177 
                0178 C      --------------------
                0179 C      Argument variables
                0180 C      --------------------
                0181 
                0182       _RL p_dictot
                0183       _RL p_bortot
                0184       _RL p_po4tot
                0185       _RL p_siltot
                0186       _RL p_nh4tot
                0187       _RL p_h2stot
                0188       _RL p_so4tot
                0189       _RL p_flutot
                0190       _RL p_alknw_inf
                0191       _RL p_alknw_sup
                0192       INTEGER i,j,k,bi,bj,myIter,myThid
                0193 
                0194 #ifdef CARBONCHEM_SOLVESAPHE
                0195 C      p_alknw_inf = -\Sum_i m_i Xtot_i
                0196 C
                0197 C      p_alknw_inf =-p_dictot*0. _d 0 & ! n = 2, m = 0
                0198 C                  -p_bortot*0. _d 0 &  ! n = 1, m = 0
                0199 C                  -p_po4tot*1. _d 0 &  ! n = 3, m = 1
                0200 C                  -p_siltot*0. _d 0 &  ! n = 1, m = 0
                0201 C                  -p_nh4tot*0. _d 0 &  ! n = 1, m = 0
                0202 C                  -p_h2stot*0. _d 0 &  ! n = 1, m = 0
                0203 C                  -p_so4tot*1. _d 0 &  ! n = 1, m = 1
                0204 C                  -p_flutot*1. _d 0    ! n = 1, m = 1
                0205 
                0206       p_alknw_inf =    -p_po4tot - p_so4tot - p_flutot
                0207 
                0208 C      p_alknw_sup = \Sum_i (n_i - m_i) Xtot_i
                0209 C
                0210 C      p_alknw_sup = p_dictot*(2. _d 0-0. _d 0) &  ! n = 2, m = 0
                0211 C                  p_bortot*(1. _d 0-0. _d 0) &    ! n = 1, m = 0
                0212 C                  p_po4tot*(3. _d 0-1. _d 0) &    ! n = 3, m = 1
                0213 C                  p_siltot*(1. _d 0-0. _d 0) &    ! n = 1, m = 0
                0214 C                  p_nh4tot*(1. _d 0-0. _d 0) &    ! n = 1, m = 0
                0215 C                  p_h2stot*(1. _d 0-0. _d 0) &    ! n = 1, m = 0
                0216 C                  p_so4tot*(1. _d 0-1. _d 0) &    ! n = 1, m = 1
                0217 C                  p_flutot*(1. _d 0-1. _d 0)      ! n = 1, m = 1
                0218 
                0219       p_alknw_sup =  p_dictot + p_dictot + p_bortot
                0220      &             + p_po4tot + p_po4tot + p_siltot
                0221      &             + p_nh4tot + p_h2stot
                0222 
                0223 #endif /* CARBONCHEM_SOLVESAPHE */
                0224       RETURN
                0225       END
                0226 C      END SUBROUTINE ANW_INFSUP
                0227 C      =========================================================================
                0228 
b5953f02df Jona*0229 C      =========================================================================
6acab690ae Jona*0230       SUBROUTINE CALC_PCO2_SOLVESAPHE(
                0231      I                    t,s,z_dictot,z_po4tot,z_siltot,z_alktot,
                0232      U                    pHlocal,z_pco2,z_co3,
                0233      I                    i,j,k,bi,bj,debugPrt,myIter,myThid )
                0234 
                0235       IMPLICIT NONE
                0236 
                0237 C     MITgcm GLobal variables
                0238 #include "SIZE.h"
                0239 #include "EEPARAMS.h"
                0240 #include "DIC_VARS.h"
                0241 
                0242 C -----------------------------------------------------------------------
                0243 C General parameters
                0244 C -----------------------------------------------------------------------
                0245 C       diclocal = total inorganic carbon (mol/m^3)
                0246 C             where 1 T = 1 metric ton = 1000 kg
                0247 C       ta  = total alkalinity (eq/m^3)
                0248 C       pt  = inorganic phosphate (mol/^3)
                0249 C       sit = inorganic silicate (mol/^3)
                0250 C       t   = temperature (degrees C)
ba0b047096 Mart*0251 C       s   = salinity (g/kg)
6acab690ae Jona*0252       _RL  t, s, z_po4tot, z_siltot, z_alktot
                0253       _RL  z_pco2, z_dictot, pHlocal
                0254       _RL  z_co3
                0255       INTEGER i,j,k,bi,bj
                0256       LOGICAL debugPrt
                0257       INTEGER myIter,myThid
                0258 
                0259 #ifdef CARBONCHEM_SOLVESAPHE
                0260 
                0261 C     == Local variables ==
                0262       _RL  z_h
                0263       _RL  z_hini
                0264       _RL  z_bortot
                0265       _RL  z_nh4tot
                0266       _RL  z_h2stot
                0267       _RL  z_so4tot
                0268       _RL  z_flutot
                0269       _RL  z_val
                0270       _RL  z_co2s
                0271       _RL  z_fco2
                0272       _RL  z_hco3
                0273       _RL  z_co2aq
                0274       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0275 C -----------------------------------------------------------------------
                0276 C Change units from the input of mol/m^3 -> mol/kg:
                0277 c (1 mol/m^3)  x (1 m^3/1024.5 kg)
                0278 c where the ocean mean surface density is 1024.5 kg/m^3
                0279 c Note: mol/kg are actually what the body of this routine uses
                0280 c for calculations.  Units are reconverted back to mol/m^3 at the
                0281 c end of this routine.
                0282 c To convert input in mol/m^3 -> mol/kg
                0283       z_po4tot=z_po4tot*permil
                0284       z_siltot=z_siltot*permil
                0285       z_alktot=z_alktot*permil
                0286       z_dictot=z_dictot*permil
                0287 
                0288 C Load from the carbon_chem common block
                0289       z_so4tot = st(i,j,bi,bj)
                0290       z_flutot = ft(i,j,bi,bj)
                0291       z_bortot = bt(i,j,bi,bj)
                0292 
                0293       z_nh4tot = 0. _d 0
                0294       z_h2stot = 0. _d 0
                0295 
                0296       z_hini = 10. _d 0**(-1. _d 0 * pHlocal)
                0297 
                0298       at_maxniter  = 50
                0299 
                0300       IF ( selectPHsolver.EQ.1 ) THEN
                0301 #ifdef ALLOW_DEBUG
                0302         IF (debugPrt) CALL DEBUG_CALL('SOLVE_AT_GENERAL',myThid)
                0303 #endif
                0304         CALL SOLVE_AT_GENERAL(z_alktot, z_dictot, z_bortot,
                0305      &                z_po4tot, z_siltot, z_nh4tot, z_h2stot,
                0306      &                z_so4tot, z_flutot, z_hini,   z_val,
                0307      &                z_h,
                0308      &                i, j, k, bi, bj, debugPrt, myIter, myThid )
                0309       ELSEIF ( selectPHsolver.EQ.2 ) THEN
                0310 #ifdef ALLOW_DEBUG
                0311         IF (debugPrt) CALL DEBUG_CALL('SOLVE_AT_GENERAL_SEC',myThid)
                0312 #endif
                0313         CALL SOLVE_AT_GENERAL_SEC(z_alktot, z_dictot, z_bortot,
                0314      &                z_po4tot, z_siltot, z_nh4tot, z_h2stot,
                0315      &                z_so4tot, z_flutot, z_hini,   z_val,
                0316      &                z_h,
                0317      &                i, j, k, bi, bj, debugPrt, myIter, myThid )
                0318       ELSEIF ( selectPHsolver.EQ.3 ) THEN
                0319 #ifdef ALLOW_DEBUG
                0320         IF (debugPrt) CALL DEBUG_CALL('SOLVE_AT_FAST',myThid)
                0321 #endif
                0322         CALL SOLVE_AT_FAST(z_alktot, z_dictot, z_bortot,
                0323      &                z_po4tot, z_siltot, z_nh4tot, z_h2stot,
                0324      &                z_so4tot, z_flutot, z_hini,   z_val,
                0325      &                z_h,
                0326      &                i, j, k, bi, bj, debugPrt, myIter, myThid )
                0327       ENDIF
                0328 
                0329 C Unlikely, but it might happen...
                0330       IF ( z_h .LT. 0. _d 0 ) THEN
                0331         WRITE(msgBuf,'(A,A,A)')
                0332      &    'S/R CALC_PCO2_SOLVESAPHE:',
                0333      &    ' H+ ion concentration less than 0',
                0334      &    ' Divergence or too many iterations'
                0335         CALL PRINT_ERROR( msgBuf, myThid )
                0336         STOP 'ABNORMAL END: S/R CALC_PCO2_SOLVESAPHE'
                0337       ENDIF
                0338 
                0339 C Return update pH to main routine
                0340       phlocal = -log10(z_h)
                0341 
                0342 C now determine [CO2*]
                0343       z_co2s  = z_dictot/
                0344      &   (1.0 _d 0 + (ak1(i,j,bi,bj)/z_h)
                0345      & + (ak1(i,j,bi,bj)*ak2(i,j,bi,bj)/(z_h*z_h)))
                0346 C Evaluate HCO3- and CO32- , carbonate ion concentration
                0347 C used in determination of calcite compensation depth
                0348       z_hco3 = ak1(i,j,bi,bj)*z_dictot /
                0349      &         (z_h*z_h + ak1(i,j,bi,bj)*z_h
                0350      &        + ak1(i,j,bi,bj)*ak2(i,j,bi,bj))
                0351       z_co3 = ak1(i,j,bi,bj)*ak2(i,j,bi,bj)*z_dictot /
                0352      &         (z_h*z_h + ak1(i,j,bi,bj)*z_h
                0353      &        + ak1(i,j,bi,bj)*ak2(i,j,bi,bj))
                0354 
                0355 c ---------------------------------------------------------------
                0356 c surface pCO2 (following Dickson and Goyet, DOE...)
                0357 #ifdef WATERVAP_BUG
                0358       z_pco2 = z_co2s/ff(i,j,bi,bj)
                0359 #else
                0360 c bug fix by Bennington
                0361       z_fco2 = z_co2s/ak0(i,j,bi,bj)
                0362       z_pco2 = z_fco2/fugf(i,j,bi,bj)
                0363 #endif
                0364 
                0365 C ----------------------------------------------------------------
                0366 C Reconvert from mol/kg -> mol/m^3
                0367       z_po4tot = z_po4tot/permil
                0368       z_siltot = z_siltot/permil
                0369       z_alktot = z_alktot/permil
                0370       z_dictot = z_dictot/permil
                0371       z_hco3   = z_hco3/permil
                0372       z_co3    = z_co3/permil
                0373       z_co2aq  = z_co2s/permil
                0374 
                0375 #endif /* CARBONCHEM_SOLVESAPHE */
                0376       RETURN
                0377       END
                0378 C      END SUBROUTINE CALC_PCO2_SOLVESAPHE
                0379 C      =========================================================================
                0380 
                0381 C      =========================================================================
                0382       SUBROUTINE DIC_COEFFS_SURF(
                0383      &                   ttemp,stemp,
                0384      &                   bi,bj,iMin,iMax,jMin,jMax,myThid)
                0385 
                0386 C     Determine coefficients for surface carbon chemistry,
                0387 C     loaded into common block
                0388 
                0389       IMPLICIT NONE
                0390 
                0391 C     MITgcm GLobal variables
                0392 #include "SIZE.h"
                0393 #include "EEPARAMS.h"
                0394 #include "PARAMS.h"
                0395 #include "GRID.h"
                0396 #include "DIC_VARS.h"
                0397 
                0398       _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0399       _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0400       INTEGER bi,bj,iMin,iMax,jMin,jMax,myThid
                0401 
                0402 #ifdef CARBONCHEM_SOLVESAPHE
                0403 
                0404 C LOCAL VARIABLES
                0405       INTEGER i, j
                0406       _RL t
                0407       _RL s
                0408       _RL t_k
                0409       _RL t_k_o_100
                0410       _RL t_k_o_100_2
                0411       _RL dlog_t_k
a2c35952f2 Jona*0412       _RL dlog10_t_k
6acab690ae Jona*0413       _RL inv_t_k
                0414       _RL ion_st, is_2, sqrtis
                0415       _RL s_2, sqrts, s_15, scl, s35
                0416       _RL log_fw2sw
                0417       _RL B, B1
                0418       _RL delta
                0419       _RL P1atm
                0420       _RL RT
                0421       _RL Rgas
                0422 C conversions for different pH scales
                0423       _RL total2free
8d230ec051 Jona*0424       _RL free2total
6acab690ae Jona*0425       _RL free2sw
8d230ec051 Jona*0426       _RL sw2free
6acab690ae Jona*0427       _RL total2sw
8d230ec051 Jona*0428       _RL sw2total
92031a2908 Jean*0429 
6acab690ae Jona*0430         do i=imin,imax
                0431          do j=jmin,jmax
                0432           if (hFacC(i,j,1,bi,bj).gt.0. _d 0) then
                0433            t = ttemp(i,j)
                0434            s = stemp(i,j)
                0435 C terms used more than once for:
                0436 C temperature
                0437            t_k         = 273.15 _d 0 + t
                0438            t_k_o_100   = t_k/100. _d 0
                0439            t_k_o_100_2 = t_k_o_100*t_k_o_100
                0440            inv_t_k=1.0 _d 0/t_k
a2c35952f2 Jona*0441            dlog_t_k    = LOG(t_k)
                0442            dlog10_t_k  = LOG10(t_k)
6acab690ae Jona*0443 C ionic strength (converted to kgsw)
                0444            ion_st=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)
                0445            is_2=ion_st*ion_st
                0446            sqrtis=sqrt(ion_st)
                0447 C salinity
                0448            s_2  = s*s
                0449            sqrts=sqrt(s)
                0450            s_15 = s*sqrts
                0451            scl  = s/1.80655 _d 0
                0452            s35  = s/35. _d 0
                0453 
                0454            log_fw2sw = LOG( 1. _d 0 - 0.001005 _d 0*s )
                0455 
                0456       IF ( selectBTconst.EQ.1 ) THEN
                0457 C -----------------------------------------------------------------------
                0458 C       Calculate the total borate concentration in mol/kg-SW
                0459 C       given the salinity of a sample
                0460 C       Ref: Uppström (1974), cited by  Dickson et al. (2007, chapter 5, p 10)
                0461 C            Millero (1982) cited in Millero (1995)
                0462 C       pH scale  : N/A
                0463            bt(i,j,bi,bj) = 0.000232 _d 0* scl/10.811 _d 0
                0464       ELSEIF ( selectBTconst.EQ.2 ) THEN
                0465 C -----------------------------------------------------------------------
                0466 C       Calculate the total borate concentration in mol/kg-SW
                0467 C       given the salinity of a sample
                0468 C       References: New formulation from Lee et al (2010)
                0469 C       pH scale  : N/A
                0470            bt(i,j,bi,bj) = 0.0002414 _d 0* scl/10.811 _d 0
                0471       ENDIF
                0472 
                0473       IF ( selectFTconst.EQ.1 ) THEN
                0474 C -----------------------------------------------------------------------
                0475 C       Calculate the total fluoride concentration in mol/kg-SW
                0476 C       given the salinity of a sample
                0477 C       References: Riley (1965)
                0478 C       pH scale  : N/A
                0479            ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0
                0480       ELSEIF ( selectFTconst.EQ.2 ) THEN
                0481 C -----------------------------------------------------------------------
                0482 C       Calculate the total fluoride concentration in mol/kg-SW
                0483 C       given the salinity of a sample
                0484 C       References: Culkin (1965) (???)
                0485 C       pH scale  : N/A
                0486            ft(i,j,bi,bj) = 0.000068 _d 0*(s35)
                0487       ENDIF
                0488 
                0489 C -----------------------------------------------------------------------
                0490 C       Calculate the total sulfate concentration in mol/kg-SW
                0491 C       given the salinity of a sample
                0492 C       References: Morris & Riley (1966) quoted in Handbook (2007)
                0493 C       pH scale  : N/A
                0494            st(i,j,bi,bj) = 0.14 _d 0 * scl/96.062 _d 0
                0495 
                0496 C -----------------------------------------------------------------------
                0497 C       Calculate the total calcium concentration in mol/kg-SW
                0498 C       given the salinity of a sample
                0499 C       References: Culkin (1965) (???)
                0500 C       pH scale  : N/A
                0501            cat(i,j,bi,bj) = 0.010282 _d 0*(s35)
                0502 
                0503 C -----------------------------------------------------------------------
                0504 C       Calculate K0 in (mol/kg-SW)/atmosphere
                0505 C       References: Weiss (1979) [(mol/kg-SW)/atm]
                0506 C       pH scale  : N/A
                0507 C       Note      : currently no pressure correction
                0508            ak0(i,j,bi,bj)  = EXP( 93.4517 _d 0/t_k_o_100 - 60.2409 _d 0
                0509      &                 + 23.3585 _d 0*LOG(t_k_o_100)
                0510      &                 + s * (0.023517 _d 0 - 0.023656 _d 0*t_k_o_100
                0511      &                 + 0.0047036 _d 0*t_k_o_100_2))
                0512 
                0513 C------------------------------------------------------------------------
                0514 C       Calculate f = k0(1-pH2O)*correction term for non-ideality
                0515 C       References: Weiss & Price (1980, Mar. Chem., 8, 347-359
                0516 C                   Eq 13 with table 6 values)
                0517 C       pH scale  : N/A
                0518            ff(i,j,bi,bj) = exp(-162.8301 _d 0 + 218.2968 _d 0/t_k_o_100
                0519      &          + 90.9241 _d 0*log(t_k_o_100) - 1.47696 _d 0*t_k_o_100_2
                0520      &          + s * (.025695 _d 0 - .025225 _d 0*t_k_o_100
                0521      &          + 0.0049867 _d 0*t_k_o_100_2))
                0522 
                0523 C------------------------------------------------------------------------
                0524 C       Calculate Fugacity Factor needed for non-ideality in ocean
                0525 C       References: Weiss (1974) Marine Chemistry
                0526 C       pH scale  : N/A
                0527            P1atm = 1.01325 _d 0 ! bars
                0528            Rgas = 83.1451 _d 0  ! bar*cm3/(mol*K)
                0529            RT = Rgas*t_k
                0530            delta = (57.7 _d 0 - 0.118 _d 0*t_k)
92031a2908 Jean*0531            B1 = -1636.75 _d 0 + 12.0408 _d 0*t_k
791313b485 Jona*0532      &                  - 0.0327957 _d 0*t_k*t_k
6acab690ae Jona*0533            B  = B1 + 3.16528 _d 0*t_k*t_k*t_k*(0.00001 _d 0)
                0534 
                0535 C   "x2" term often neglected (assumed=1) in applications of Weiss's (1974) eq.9
                0536 C    x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite)
                0537            fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * P1atm / RT)
                0538 
                0539       IF ( selectK1K2const.EQ.1 ) THEN
                0540 C -----------------------------------------------------------------------
                0541 C       Calculate first dissociation constant of carbonic acid
                0542 C       in mol/kg-SW on the Seawater pH-scale.
                0543 C       References: Millero (1995, eq 35 -- pK1(MEHR));
                0544 C                   Mehrbach et al. (1973) data
                0545 C       pH scale:   Seawater
                0546            ak1(i,j,bi,bj)=10.**(-1. _d 0*(3670.7 _d 0*inv_t_k
                0547      &          - 62.008 _d 0 + 9.7944 _d 0*dlog_t_k
                0548      &          - 0.0118 _d 0 * s + 0.000116 _d 0*s_2))
                0549 
                0550 C       Calculate second dissociation constant of carbonic acid
                0551 C       in mol/kg-SW on the Seawater pH-scale.
                0552 C       References: Millero (1995, eq 36 -- pK2(MEHR))
                0553 C                   Mehrbach et al. (1973) data
                0554 C       pH scale:   Seawater
                0555            ak2(i,j,bi,bj)=10.**(-1. _d 0*(1394.7 _d 0*inv_t_k
                0556      &           + 4.777 _d 0 - 0.0184 _d 0*s + 0.000118 _d 0*s_2))
                0557 
                0558       ELSEIF ( selectK1K2const.EQ.2 ) THEN
                0559 C -----------------------------------------------------------------------
                0560 C       Calculate first dissociation constant of carbonic acid
                0561 C       in mol/kg-SW on the Total pH-scale.
                0562 C       References: Roy et al. (1993) -- also Handbook (1994)
                0563 C       pH scale  : Total
                0564 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0565            ak1(i,j,bi,bj)  =  EXP(-2307.1255 _d 0*inv_t_k + 2.83655 _d 0
                0566      &              - 1.5529413 _d 0*dlog_t_k
                0567      &              + (-4.0484 _d 0*inv_t_k - 0.20760841)*sqrts
                0568      &              + 0.08468345*s
                0569      &              - 0.00654208*s_15
                0570      &              + log_fw2sw )
                0571 
                0572 C       Calculate second dissociation constant of carbonic acid
                0573 C       in mol/kg-SW on the Total pH-scale.
                0574 C       References: Roy et al. (1993) -- also Handbook (1994)
                0575 C       pH scale  : Total
                0576 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0577            ak2(i,j,bi,bj) = EXP(-3351.6106 _d 0*inv_t_k - 9.226508 _d 0
                0578      &              - 0.2005743 _d 0*dlog_t_k
                0579      &              + ( -23.9722 _d 0*inv_t_k - 0.106901773 _d 0)*sqrts
                0580      &              + 0.1130822*s - 0.00846934 _d 0*s_15
                0581      &              + log_fw2sw )
                0582 
                0583       ELSEIF ( selectK1K2const.EQ.3 ) THEN
                0584 C -----------------------------------------------------------------------
                0585 C       Calculate first dissociation constant of carbonic acid
                0586 C       in mol/kg-SW on the Seawater pH-scale.
                0587 C       References: Millero (1995, eq 50 -- ln K1(COM))
                0588 C       pH scale:   Seawater
                0589            ak1(i,j,bi,bj)  = EXP(2.18867 _d 0 - 2275.0360 _d 0*inv_t_k
                0590      &              - 1.468591 _d 0*dlog_t_k
                0591      &              + ( -0.138681 _d 0 - 9.33291 _d 0*inv_t_k)*sqrts
                0592      &              + 0.0726483 _d 0*s - 0.00574938 _d 0*s_15)
                0593 
                0594 C       Calculate second dissociation constant of carbonic acid
                0595 C       in mol/kg-SW on the Seawater pH-scale.
                0596 C       References: Millero (1995, eq 51 -- ln K2(COM))
                0597 C       pH scale:   Seawater
                0598            ak2(i,j,bi,bj)  = EXP(-0.84226 _d 0 - 3741.1288 _d 0*inv_t_k
                0599      &              -  1.437139 _d 0*dlog_t_k
                0600      &              + (-0.128417 _d 0 - 24.41239 _d 0*inv_t_k)*sqrts
                0601      &              +  0.1195308 _d 0*s - 0.00912840 _d 0*s_15)
                0602 
                0603       ELSEIF ( selectK1K2const.EQ.4 ) THEN
                0604 C -----------------------------------------------------------------------
                0605 C       Calculate first dissociation constant of carbonic acid
                0606 C       in mol/kg-SW on the Total pH-scale.
                0607 C       Suitable when 2 < T < 35 and 19 < S < 43
                0608 C       References: Luecker et al. (2000) -- also Handbook (2007)
                0609 C       pH scale:   Total
                0610            ak1(i,j,bi,bj)  = 10. _d 0**( 61.2172 _d 0
                0611      &                 - 3633.86 _d 0*inv_t_k - 9.67770 _d 0*dlog_t_k
                0612      &                 + s*(0.011555 - s*0.0001152 _d 0))
                0613 
                0614 C       Calculate second dissociation constant of carbonic acid
                0615 C       in mol/kg-SW on the Total pH-scale.
                0616 C       Suitable when 2 < T < 35 and 19 < S < 43
                0617 C       References: Luecker et al. (2000) -- also Handbook (2007)
                0618 C       pH scale:   Total
                0619            ak2(i,j,bi,bj)  = 10. _d 0**(-25.9290 _d 0
                0620      &                 - 471.78 _d 0*inv_t_k + 3.16967 _d 0*dlog_t_k
                0621      &                 + s*(0.01781 _d 0 - s*0.0001122 _d 0))
                0622 
                0623       ELSEIF ( selectK1K2const.EQ.5 ) THEN
                0624 C -----------------------------------------------------------------------
                0625 C       Calculate first dissociation constant of carbonic acid
                0626 C       in mol/kg-SW on the Seawater pH-scale.
                0627 C       Suitable when 0 < T < 50 and 1 < S < 50
                0628 C       References: Millero (2010, Mar. Fresh Wat. Res.)
                0629 C                   Millero (1979) pressure correction
                0630 C       pH scale:   Seawater
                0631            ak1(i,j,bi,bj) = 10.0 _d 0**(-1*(6320.813 _d 0*inv_t_k
                0632      &      + 19.568224 _d 0*dlog_t_k -126.34048 _d 0
                0633      &      + 13.4038 _d 0*sqrts + 0.03206 _d 0*s - (5.242 _d -5)*s_2
                0634      &      + (-530.659 _d 0*sqrts - 5.8210 _d 0*s)*inv_t_k
                0635      &      -2.0664 _d 0*sqrts*dlog_t_k))
                0636 
                0637 C       Calculate second dissociation constant of carbonic acid
                0638 C       in mol/kg-SW on the Seawater pH-scale.
                0639 C       Suitable when 0 < T < 50 and 1 < S < 50
                0640 C       References: Millero (2010, Mar. Fresh Wat. Res.)
                0641 C                   Millero (1979) pressure correction
                0642 C       pH scale:   Seawater
                0643            ak2(i,j,bi,bj) = 10.0 _d 0**(-1*(5143.692 _d 0*inv_t_k
                0644      &     + 14.613358 _d 0*dlog_t_k -90.18333 _d 0
                0645      &     + 21.3728 _d 0*sqrts + 0.1218 _d 0*s - (3.688 _d -4)*s_2
                0646      &     + (-788.289 _d 0*sqrts - 19.189 _d 0*s)*inv_t_k
                0647      &     -3.374 _d 0*sqrts*dlog_t_k))
                0648 
                0649       ELSEIF ( selectK1K2const.EQ.6 ) THEN
                0650 C -----------------------------------------------------------------------
                0651 C       Calculate first dissociation constant of carbonic acid
                0652 C       in mol/kg-SW on the Seawater pH-scale.
                0653 C       Suitable when 0 < T < 50 and 1 < S < 50
                0654 C       References: Waters, Millero, Woosley (Mar. Chem., 165, 66-67, 2014)
                0655 C                   Millero (1979) pressure correction
                0656 C       pH scale:   Seawater
                0657            ak1(i,j,bi,bj) = 10.0 _d 0**(-1.*(6320.813 _d 0*inv_t_k
                0658      &     + 19.568224 _d 0*dlog_t_k -126.34048 _d 0
                0659      &     + 13.409160 _d 0*sqrts + 0.031646 _d 0*s - (5.1895 _d -5)*s_2
                0660      &     + (-531.3642 _d 0*sqrts - 5.713 _d 0*s)*inv_t_k
                0661      &     -2.0669166 _d 0*sqrts*dlog_t_k))
                0662 
                0663 C       Calculate second dissociation constant of carbonic acid
                0664 C       in mol/kg-SW on the Seawater pH-scale.
                0665 C       Suitable when 0 < T < 50 and 1 < S < 50
                0666 C       References: Waters, Millero, Woosley (Mar. Chem., 165, 66-67, 2014)
                0667 C                   Millero (1979) pressure correction
                0668 C       pH scale:   Seawater
                0669            ak2(i,j,bi,bj) = 10.0 _d 0**(-1.*
                0670      &     ( 5143.692 _d 0*inv_t_k + 14.613358 _d 0*dlog_t_k
                0671      &      - 90.18333 _d 0 + 21.225890 _d 0*sqrts + 0.12450870 _d 0*s
                0672      &      - (3.7243 _d -4)*s_2
                0673      &      + (-779.3444 _d 0*sqrts - 19.91739 _d 0*s)*inv_t_k
                0674      &      - 3.3534679 _d 0*sqrts*dlog_t_k ) )
                0675 
                0676       ENDIF /* selectK1K2const */
                0677 
                0678 C -----------------------------------------------------------------------
                0679 C       Calculate boric acid dissociation constant KB
                0680 C       in mol/kg-SW on the total pH-scale.
                0681 C       References: Dickson (1990, eq. 23) -- also Handbook (2007, eq. 37)
                0682 C       pH scale  : total
                0683            akb(i,j,bi,bj)  = EXP(( -8966.90 _d 0 - 2890.53 _d 0*sqrts
                0684      &      -77.942 _d 0*s + 1.728 _d 0*s_15 - 0.0996 _d 0*s_2 )*inv_t_k
                0685      &      + (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s)
                0686      &      + (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s)*
                0687      &      dlog_t_k + 0.053105 _d 0*sqrts*t_k )
                0688 
                0689 C -----------------------------------------------------------------------
                0690 C       Calculate the first dissociation constant
                0691 C       of phosphoric acid (H3PO4) in seawater
                0692 C       References: Yao and Millero (1995)
                0693 C       pH scale  : Seawater
                0694            ak1p(i,j,bi,bj) = EXP(115.54 _d 0
                0695      &              - 4576.752 _d 0*inv_t_k - 18.453 _d 0*dlog_t_k
                0696      &              + ( 0.69171 _d 0 -  106.736 _d 0*inv_t_k)*sqrts
                0697      &              + (-0.01844 _d 0 -  0.65643 _d 0*inv_t_k)*s)
                0698 
                0699 C -----------------------------------------------------------------------
                0700 C       Calculate the second dissociation constant
                0701 C       of phosphoric acid (H3PO4) in seawater
                0702 C       References: Yao and Millero (1995)
                0703 C       pH scale  : Seawater
                0704            ak2p(i,j,bi,bj) = EXP( 172.1033 _d 0
                0705      &              - 8814.715 _d 0*inv_t_k
                0706      &              -   27.927 _d 0*dlog_t_k
                0707      &              + (  1.3566 _d 0 -  160.340 _d 0*inv_t_k)*sqrts
                0708      &              + (-0.05778 _d 0 +  0.37335 _d 0*inv_t_k)*s)
                0709 
                0710 C -----------------------------------------------------------------------
                0711 C       Calculate the third dissociation constant
                0712 C       of phosphoric acid (H3PO4) in seawater
                0713 C       References: Yao and Millero (1995)
                0714 C       pH scale  : Seawater
                0715            ak3p(i,j,bi,bj) = EXP(-18.126 _d 0  -  3070.75 _d 0*inv_t_k
                0716      &                + ( 2.81197 _d 0 + 17.27039 _d 0*inv_t_k)*sqrts
                0717      &                + (-0.09984 _d 0 - 44.99486 _d 0*inv_t_k)*s)
                0718 
                0719 C -----------------------------------------------------------------------
                0720 C       Calculate the first dissociation constant
                0721 C       of silicic acid (H4SiO4) in seawater
                0722 C       References: Yao and Millero (1995) cited by Millero (1995)
                0723 C       pH scale  : Seawater (according to Dickson et al, 2007)
                0724 C       Note      : converted here from mol/kg-H2O to mol/kg-sw
                0725            aksi(i,j,bi,bj) = EXP(
                0726      &           117.40 _d 0 - 8904.2 _d 0*inv_t_k
                0727      &         - 19.334 _d 0 * dlog_t_k
                0728      &         + ( 3.5913 _d 0 -  458.79 _d 0*inv_t_k) * sqrtis
                0729      &         + (-1.5998 _d 0 +  188.74 _d 0*inv_t_k) * ion_st
                0730      &         + (0.07871 _d 0 - 12.1652 _d 0*inv_t_k) * ion_st*ion_st
                0731      &         + log_fw2sw )
                0732 
                0733 C -----------------------------------------------------------------------
                0734 C       Calculate the dissociation constant
                0735 C       of ammonium in sea-water [mol/kg-SW]
                0736 C       References: Yao and Millero (1995)
                0737 C       pH scale  : Seawater
                0738            akn(i,j,bi,bj)  = EXP(-0.25444 _d 0 -  6285.33 _d 0*inv_t_k
                0739      &                + 0.0001635 _d 0 * t_k
                0740      &                + ( 0.46532 _d 0 - 123.7184 _d 0*inv_t_k)*sqrts
                0741      &                + (-0.01992 _d 0 +  3.17556 _d 0*inv_t_k)*s)
                0742 
                0743 C -----------------------------------------------------------------------
                0744 C       Calculate the dissociation constant of hydrogen sulfide in sea-water
                0745 C       References: Millero et al. (1988) (cited by Millero (1995)
                0746 C       pH scale  : - Seawater (according to Yao and Millero, 1995,
                0747 C                               p. 82: "refitted if necessary")
                0748 C                   - Total (according to Lewis and Wallace, 1998)
                0749 C       Note      : we stick to Seawater here for the time being
                0750 C       Note      : the fits from Millero (1995) and Yao and Millero (1995)
                0751 C                   derive from Millero et al. (1998), with all the coefficients
                0752 C                   multiplied by -ln(10)
                0753            akhs(i,j,bi,bj) = EXP( 225.838 _d 0 - 13275.3 _d 0*inv_t_k
                0754      &               - 34.6435 _d 0 * dlog_t_k
                0755      &               +  0.3449 _d 0*sqrts -  0.0274 _d 0*s)
                0756 
                0757 C -----------------------------------------------------------------------
                0758 C       Calculate the dissociation constant of hydrogen sulfate (bisulfate)
                0759 C       References: Dickson (1990) -- also Handbook (2007)
                0760 C       pH scale  : free
                0761 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0762            aks(i,j,bi,bj) = EXP(141.328 _d 0
                0763      &                -   4276.1 _d 0*inv_t_k -  23.093 _d 0*dlog_t_k
                0764      &                + ( 324.57 _d 0 - 13856. _d 0*inv_t_k
                0765      &                -   47.986 _d 0*dlog_t_k) * sqrtis
                0766      &                + (-771.54 _d 0 + 35474. _d 0*inv_t_k
                0767      &                +  114.723 _d 0*dlog_t_k) * ion_st
                0768      &                - 2698. _d 0*inv_t_k*ion_st**1.5 _d 0
                0769      &                + 1776. _d 0*inv_t_k*ion_st*ion_st
                0770      &                + log_fw2sw )
                0771 
                0772       IF ( selectHFconst.EQ.1 ) THEN
                0773 C -----------------------------------------------------------------------
                0774 C       Calculate the dissociation constant \beta_{HF} [(mol/kg-SW)^{-1}]
                0775 C       in (mol/kg-SW)^{-1}, where
                0776 C         \beta_{HF} = \frac{ [HF] }{ [H^{+}] [F^{-}] }
                0777 C       References: Dickson and Riley (1979)
                0778 C       pH scale  : free
                0779 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0780            akf(i,j,bi,bj) = EXP(1590.2 _d 0*inv_t_k - 12.641 _d 0
                0781      &                 + 1.525 _d 0*sqrtis
                0782      &                 + log_fw2sw )
                0783       ELSEIF ( selectHFconst.EQ.2 ) THEN
                0784 C -----------------------------------------------------------------------
                0785 C       Calculate the dissociation constant for hydrogen fluoride
                0786 C       in mol/kg-SW
                0787 C       References: Perez and Fraga (1987)
                0788 C       pH scale  : Total (according to Handbook, 2007)
                0789            akf(i,j,bi,bj) = EXP( 874. _d 0*inv_t_k - 9.68 _d 0
                0790      &                           + 0.111 _d 0*sqrts )
                0791       ENDIF
                0792 
                0793 C -----------------------------------------------------------------------
                0794 C       Calculate the water dissociation constant Kw in (mol/kg-SW)^2
                0795 C       References: Millero (1995) for value at pressc = 0
                0796 C       pH scale  : Seawater
                0797            akw(i,j,bi,bj) =  EXP(148.9802 _d 0 - 13847.26 _d 0*inv_t_k
                0798      &              -   23.6521 _d 0*dlog_t_k
                0799      &              + (-5.977 _d 0 + 118.67 _d 0*inv_t_k
                0800      &              + 1.0495 _d 0*dlog_t_k)*sqrts
                0801      &              -   0.01615 _d 0*s )
                0802 
                0803 C -----------------------------------------------------------------------
                0804 C pH scale conversion factors and conversions. Everything on total pH scale
                0805            total2free = 1. _d 0/
                0806      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                0807 
92031a2908 Jean*0808            free2total = (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
8d230ec051 Jona*0809 
6acab690ae Jona*0810            free2sw = 1. _d 0
                0811      &                 + (st(i,j,bi,bj)/ aks(i,j,bi,bj))
                0812      &                 + (ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free))
                0813 
8d230ec051 Jona*0814            sw2free = 1. _d 0 / free2sw
92031a2908 Jean*0815 
6acab690ae Jona*0816            total2sw = total2free * free2sw
                0817 
8d230ec051 Jona*0818            sw2total = 1. _d 0 / total2sw
                0819 
6acab690ae Jona*0820            aphscale(i,j,bi,bj) = 1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj)
                0821 
                0822       IF ( selectK1K2const.NE.2 .AND. selectK1K2const.NE.4 ) THEN
                0823 C Convert to the total pH scale
8d230ec051 Jona*0824            ak1(i,j,bi,bj)  = ak1(i,j,bi,bj)*sw2total
                0825            ak2(i,j,bi,bj)  = ak2(i,j,bi,bj)*sw2total
6acab690ae Jona*0826       ENDIF
8d230ec051 Jona*0827            ak1p(i,j,bi,bj) = ak1p(i,j,bi,bj)*sw2total
                0828            ak2p(i,j,bi,bj) = ak2p(i,j,bi,bj)*sw2total
                0829            ak3p(i,j,bi,bj) = ak3p(i,j,bi,bj)*sw2total
                0830            aksi(i,j,bi,bj) = aksi(i,j,bi,bj)*sw2total
                0831            akn (i,j,bi,bj) = akn (i,j,bi,bj)*sw2total
                0832            akhs(i,j,bi,bj) = akhs(i,j,bi,bj)*sw2total
                0833            aks (i,j,bi,bj) = aks (i,j,bi,bj)*free2total
6acab690ae Jona*0834       IF ( selectHFconst.EQ.1 ) THEN
8d230ec051 Jona*0835            akf (i,j,bi,bj) = akf (i,j,bi,bj)*free2total
6acab690ae Jona*0836       ENDIF
8d230ec051 Jona*0837            akw (i,j,bi,bj) = akw (i,j,bi,bj)*sw2total
6acab690ae Jona*0838 
                0839 C -----------------------------------------------------------------------
                0840 C       Calculate the stoichiometric solubility product
                0841 C       of calcite in seawater
                0842 C       References: Mucci (1983)
                0843 C       pH scale  : N/A
                0844 C       Units     : (mol/kg-SW)^2
                0845 
                0846            Ksp_TP_Calc(i,j,bi,bj) = 10. _d 0**(-171.9065 _d 0
                0847      &             - 0.077993 _d 0*t_k
a2c35952f2 Jona*0848      &             + 2839.319 _d 0*inv_t_k + 71.595 _d 0*dlog10_t_k
6acab690ae Jona*0849      &             + ( -0.77712 _d 0 + 0.0028426 _d 0*t_k
                0850      &             + 178.34 _d 0*inv_t_k)*sqrts
                0851      &             - 0.07711 _d 0*s + 0.0041249 _d 0*s_15)
                0852 
                0853 C -----------------------------------------------------------------------
                0854 C       Calculate the stoichiometric solubility product
                0855 C       of aragonite in seawater
                0856 C       References: Mucci (1983)
                0857 C       pH scale  : N/A
                0858 C       Units     : (mol/kg-SW)^2
                0859 
                0860            Ksp_TP_Arag(i,j,bi,bj) = 10. _d 0**(-171.945 _d 0
                0861      &             - 0.077993 _d 0*t_k
a2c35952f2 Jona*0862      &             + 2903.293 _d 0*inv_t_k + 71.595 _d 0*dlog10_t_k
6acab690ae Jona*0863      &             + ( -0.068393 _d 0 + 0.0017276 _d 0*t_k
                0864      &             + 88.135 _d 0*inv_t_k)*sqrts
                0865      &             - 0.10018 _d 0*s + 0.0059415 _d 0*s_15)
                0866          else
                0867            bt(i,j,bi,bj)  = 0. _d 0
                0868            st(i,j,bi,bj)  = 0. _d 0
                0869            ft(i,j,bi,bj)  = 0. _d 0
                0870            cat(i,j,bi,bj) = 0. _d 0
                0871            fugf(i,j,bi,bj)= 0. _d 0
                0872            ff(i,j,bi,bj)  = 0. _d 0
                0873            ak0(i,j,bi,bj) = 0. _d 0
                0874            ak1(i,j,bi,bj) = 0. _d 0
                0875            ak2(i,j,bi,bj) = 0. _d 0
                0876            akb(i,j,bi,bj) = 0. _d 0
                0877            ak1p(i,j,bi,bj)= 0. _d 0
                0878            ak2p(i,j,bi,bj)= 0. _d 0
                0879            ak3p(i,j,bi,bj)= 0. _d 0
                0880            aksi(i,j,bi,bj)= 0. _d 0
                0881            akw(i,j,bi,bj) = 0. _d 0
                0882            aks(i,j,bi,bj) = 0. _d 0
                0883            akf(i,j,bi,bj) = 0. _d 0
                0884            akn(i,j,bi,bj) = 0. _d 0
                0885            akhs(i,j,bi,bj)= 0. _d 0
                0886            Ksp_TP_Calc(i,j,bi,bj) = 0. _d 0
                0887            Ksp_TP_Arag(i,j,bi,bj) = 0. _d 0
                0888            aphscale(i,j,bi,bj)    = 0. _d 0
                0889          endif
                0890          end do
                0891         end do
                0892 #endif /* CARBONCHEM_SOLVESAPHE */
                0893       RETURN
                0894       END
                0895 C     END SUBROUTINE DIC_COEFFS_SURF
                0896 C -----------------------------------------------------------------------
                0897 
                0898 C -----------------------------------------------------------------------
                0899       SUBROUTINE DIC_COEFFS_DEEP(
                0900      &                   ttemp,stemp,
                0901      &                   bi,bj,iMin,iMax,jMin,jMax,
                0902      &                   Klevel,myThid)
                0903 
                0904 C     Add depth dependence to carbon chemistry coefficients loaded into
                0905 C     common block. Corrections occur on the Seawater pH scale and are
                0906 C     converted back to the total scale
                0907       IMPLICIT NONE
                0908 
                0909 C     MITgcm GLobal variables
                0910 #include "SIZE.h"
                0911 #include "EEPARAMS.h"
                0912 #include "PARAMS.h"
                0913 #include "GRID.h"
                0914 #include "DIC_VARS.h"
                0915 
                0916       _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0917       _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0918       INTEGER bi,bj,iMin,iMax,jMin,jMax,Klevel,myThid
                0919 
                0920 #ifdef CARBONCHEM_SOLVESAPHE
                0921 C LOCAL VARIABLES
                0922       INTEGER i, j, k
                0923       _RL  bdepth
                0924       _RL  cdepth
                0925       _RL  pressc
                0926       _RL t
                0927       _RL s
                0928       _RL zds
                0929       _RL t_k
                0930       _RL t_k_o_100
                0931       _RL t_k_o_100_2
                0932       _RL dlog_t_k
                0933       _RL sqrtis
                0934       _RL sqrts
                0935       _RL s_2
                0936       _RL s_15
                0937       _RL inv_t_k
                0938       _RL ion_st
                0939       _RL is_2
                0940       _RL scl
                0941       _RL zrt
                0942       _RL B1
                0943       _RL B
                0944       _RL delta
                0945       _RL Pzatm
                0946       _RL zdvi
                0947       _RL zdki
                0948       _RL pfac
                0949 C pH scale converstions
                0950       _RL total2free_surf
                0951       _RL free2sw_surf
                0952       _RL total2sw_surf
                0953       _RL total2free
8d230ec051 Jona*0954       _RL free2total
6acab690ae Jona*0955       _RL free2sw
8d230ec051 Jona*0956       _RL sw2free
6acab690ae Jona*0957       _RL total2sw
8d230ec051 Jona*0958       _RL sw2total
92031a2908 Jean*0959 
6acab690ae Jona*0960 c determine pressure (bar) from depth
                0961 c 1 BAR at z=0m (atmos pressure)
                0962 c use UPPER surface of cell so top layer pressure = 0 bar
                0963 c for surface exchange coeffs
                0964 
                0965 cmick..............................
                0966 c        write(6,*)'Klevel ',klevel
                0967 
                0968         bdepth = 0. _d 0
                0969         cdepth = 0. _d 0
                0970         pressc = 1. _d 0
                0971         do k = 1,Klevel
                0972             cdepth = bdepth + 0.5 _d 0*drF(k)
                0973             bdepth = bdepth + drF(k)
                0974             pressc = 1. _d 0 + 0.1 _d 0*cdepth
                0975         end do
                0976 cmick...................................................
                0977 c        write(6,*)'depth,pressc ',cdepth,pressc
                0978 cmick....................................................
                0979 
                0980         do i=imin,imax
                0981          do j=jmin,jmax
                0982           if (hFacC(i,j,Klevel,bi,bj).gt.0. _d 0) then
                0983            t = ttemp(i,j)
                0984            s = stemp(i,j)
                0985 C terms used more than once for:
                0986 C temperature
                0987            t_k = 273.15 _d 0 + t
                0988            zrt= 83.14472 _d 0 * t_k
                0989            t_k_o_100 = t_k/100. _d 0
                0990            t_k_o_100_2=t_k_o_100*t_k_o_100
                0991            inv_t_k=1.0 _d 0/t_k
                0992            dlog_t_k=log(t_k)
                0993 C ionic strength
                0994            ion_st=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)
                0995            is_2=ion_st*ion_st
                0996            sqrtis=sqrt(ion_st)
                0997 C salinity
                0998            s_2=s*s
                0999            sqrts=sqrt(s)
                1000            s_15=s**1.5 _d 0
                1001            scl=s/1.80655 _d 0
                1002            zds = s - 34.8 _d 0
                1003 
                1004            total2free_surf = 1. _d 0/
                1005      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                1006 
                1007            free2sw_surf = 1. _d 0
                1008      &                 + st(i,j,bi,bj)/ aks(i,j,bi,bj)
                1009      &                 + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free_surf)
                1010 
                1011            total2sw_surf = total2free_surf * free2sw_surf
                1012 
                1013 C------------------------------------------------------------------------
                1014 C       Recalculate Fugacity Factor needed for non-ideality in ocean
                1015 C       with pressure dependence.
                1016 C       Reference : Weiss (1974) Marine Chemistry
                1017 C       pH scale  : N/A
                1018            Pzatm = 1.01325 _d 0+pressc ! bars
                1019            delta = (57.7 _d 0 - 0.118 _d 0*t_k)
                1020            B1 = -1636.75 _d 0 + 12.0408 _d 0*t_k -0.0327957 _d 0*t_k*t_k
                1021            B  = B1 + 3.16528 _d 0*t_k*t_k*t_k*(0.00001 _d 0)
                1022 C   "x2" term often neglected (assumed=1) in applications of Weiss's (1974) eq.9
                1023 C    x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite)
                1024            fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * Pzatm / zrt)
                1025 
                1026 C -----------------------------------------------------------------------
                1027 C       Apply pressure dependence to the dissociation constant of hydrogen
                1028 C       sulfate (bisulfate).  Ref: Millero (1995) for pressure correction
                1029            zdvi         =  -18.03 _d 0 + t*(0.0466 _d 0 + t*0.316 _d -3)
                1030            zdki         = ( -4.53 _d 0 + t*0.0900 _d 0)*1. _d -3
                1031            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1032            aks(i,j,bi,bj) = total2free_surf*aks(i,j,bi,bj) * exp(pfac)
                1033 
                1034            total2free = 1. _d 0/
                1035      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
92031a2908 Jean*1036 
8d230ec051 Jona*1037            free2total = (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
6acab690ae Jona*1038 
                1039 C free2sw has an additional component from fluoride
                1040            free2sw    = 1. _d 0
                1041      &                 + st(i,j,bi,bj)/ aks(i,j,bi,bj)
92031a2908 Jean*1042 
8d230ec051 Jona*1043            aks(i,j,bi,bj) = aks(i,j,bi,bj)*free2total
6acab690ae Jona*1044 
                1045 C -----------------------------------------------------------------------
                1046 C       Apply pressure dependence to dissociation constant for hydrogen fluoride
                1047 C       References: Millero (1995) for pressure correction
                1048            zdvi   =   -9.78 _d 0 + t*(-0.0090 _d 0 - t*0.942 _d -3)
                1049            zdki   = ( -3.91 _d 0 + t*0.054 _d 0)*1. _d -3
                1050            pfac   = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1051            akf(i,j,bi,bj) = total2free_surf*akf(i,j,bi,bj) * exp(pfac)
                1052 
b5953f02df Jona*1053 C free2sw has an additional component from fluoride add it here
6acab690ae Jona*1054            free2sw = free2sw
                1055      &               + ft(i,j,bi,bj)/akf(i,j,bi,bj)
92031a2908 Jean*1056 
8d230ec051 Jona*1057            sw2free = 1. _d 0 / free2sw
92031a2908 Jean*1058 
6acab690ae Jona*1059            total2sw = total2free * free2sw
                1060 
92031a2908 Jean*1061            sw2total = 1. _d 0 / total2sw
8d230ec051 Jona*1062 
                1063            akf(i,j,bi,bj) = akf(i,j,bi,bj)*free2total
6acab690ae Jona*1064 
                1065 C -----------------------------------------------------------------------
                1066 C       Apply pressure dependence to 1rst dissociation constant of carbonic acid
                1067 C       References: Millero (1982) pressure correction
                1068            zdvi =  -25.50 _d 0 -0.151 _d 0*zds + 0.1271 _d 0*t
                1069            zdki = ( -3.08 _d 0 -0.578 _d 0*zds + 0.0877 _d 0*t)*1. _d -3
                1070            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1071            ak1(i,j,bi,bj) = (total2sw_surf*ak1(i,j,bi,bj)
8d230ec051 Jona*1072      &                      * exp(pfac))*sw2total
92031a2908 Jean*1073 
6acab690ae Jona*1074 C -----------------------------------------------------------------------
                1075 C       Apply pressure dependence to 2nd dissociation constant of carbonic acid
                1076 C       References: Millero (1979) pressure correction
                1077            zdvi = -15.82 _d 0 + 0.321 _d 0*zds - 0.0219 _d 0*t
                1078            zdki = ( 1.13 _d 0 - 0.314 _d 0*zds - 0.1475 _d 0*t)*1. _d -3
                1079            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1080            ak2(i,j,bi,bj) = (total2sw_surf*ak2(i,j,bi,bj)
8d230ec051 Jona*1081      &                      * exp(pfac))*sw2total
6acab690ae Jona*1082 
                1083 C -----------------------------------------------------------------------
                1084 C       Apply pressure dependence to the boric acid dissociation constant KB
                1085 C       References: Millero (1979) pressure correction
                1086            zdvi      = -29.48 _d 0 + 0.295 _d 0*zds + 0.1622 _d 0*t
                1087      &                 - 0.002608 _d 0*t*t
                1088            zdki      = (-2.84 _d 0 + 0.354 _d 0*zds)*1. _d -3
                1089            pfac =  (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1090            akb(i,j,bi,bj) = (total2sw_surf*akb(i,j,bi,bj)
8d230ec051 Jona*1091      &                      * exp(pfac))*sw2total
6acab690ae Jona*1092 
                1093 C -----------------------------------------------------------------------
                1094 C       Apply pressure dependence to water dissociation constant Kw in
                1095 C       (mol/kg-SW)^2. Ref.: Millero (pers. comm. 1996) for pressure correction
                1096            zdvi    =  -20.02 _d 0 + 0.1119 _d 0*t - 0.1409 _d -2*t*t
                1097            zdki    = ( -5.13 _d 0 + 0.0794 _d 0*t)*1. _d -3
                1098            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1099            akw(i,j,bi,bj) = (total2sw_surf*akw(i,j,bi,bj)
8d230ec051 Jona*1100      &                      * exp(pfac))*sw2total
6acab690ae Jona*1101 
                1102 C -----------------------------------------------------------------------
                1103 C       Apply pressure dependence to the first dissociation constant
                1104 C       of phosphoric acid (H3PO4) in seawater
                1105 C       References: Millero (1995) for pressure correction
                1106            zdvi      =  -14.51 _d 0 + 0.1211 _d 0*t - 0.321 _d -3*t*t
                1107            zdki      = ( -2.67 _d 0 + 0.0427 _d 0*t)*1. _d -3
                1108            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1109            ak1p(i,j,bi,bj) = (total2sw_surf*ak1p(i,j,bi,bj)
8d230ec051 Jona*1110      &                       * exp(pfac))*sw2total
6acab690ae Jona*1111 
                1112 C -----------------------------------------------------------------------
                1113 C       Apply pressure dependence to the second dissociation constant
                1114 C       of phosphoric acid (H3PO4) in seawater
                1115 C       References: Millero (1995) for pressure correction
                1116            zdvi       =  -23.12 _d 0 + 0.1758 _d 0*t -2.647 _d -3*t*t
                1117            zdki       = ( -5.15 _d 0 +   0.09 _d 0*t)*1. _d -3
                1118            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1119            ak2p(i,j,bi,bj) = (total2sw_surf*ak2p(i,j,bi,bj)
8d230ec051 Jona*1120      &                       * exp(pfac))*sw2total
6acab690ae Jona*1121 
                1122 C -----------------------------------------------------------------------
                1123 C       Apply pressure dependence to the third dissociation constant
                1124 C       of phosphoric acid (H3PO4) in seawater
                1125 C       References: Millero (1995) for pressure correction
                1126            zdvi      =  -26.57 _d 0 + 0.2020 _d 0*t -3.042 _d -3*t*t
                1127            zdki      = ( -4.08 _d 0 + 0.0714 _d 0*t)*1. _d -3
                1128            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1129            ak3p(i,j,bi,bj) = (total2sw_surf*ak3p(i,j,bi,bj)
8d230ec051 Jona*1130      &                       * exp(pfac))*sw2total
6acab690ae Jona*1131 
                1132 C -----------------------------------------------------------------------
                1133 C       Apply pressure dependence to the first dissociation constant
                1134 C       of silicic acid (H4SiO4) in seawater
                1135 C       References: Millero (1979) pressure correction. Note: Pressure
                1136 C        correction estimated to be the same as borate (Millero, 1995)
                1137            zdvi      = -29.48 _d 0 + 0.295 _d 0*zds + 0.1622 _d 0*t
                1138      &                 - 0.002608 _d 0*t*t
                1139            zdki      = (-2.84 _d 0 + 0.354 _d 0*zds)*1. _d -3
                1140            pfac =  (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1141            aksi(i,j,bi,bj) = (total2sw_surf*aksi(i,j,bi,bj)
8d230ec051 Jona*1142      &                       * exp(pfac))*sw2total
6acab690ae Jona*1143 
                1144 C -----------------------------------------------------------------------
                1145 C       Apply pressure dependence to the dissociation constant of hydrogen
                1146 C       sulfide in sea-water
                1147 C       References: Millero (1995) for pressure correction
                1148            zdvi         =  -14.80 _d 0 + t*(0.0020 _d 0 - t*0.400 _d -3)
                1149            zdki         = (  2.89 _d 0 + t*0.054 _d 0)*1. _d -3
                1150            pfac  = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1151            akhs(i,j,bi,bj) = (total2sw_surf*akhs(i,j,bi,bj)
8d230ec051 Jona*1152      &                       * exp(pfac))*sw2total
6acab690ae Jona*1153 
                1154 C -----------------------------------------------------------------------
                1155 C       Apply pressure dependence to the dissociation constant
                1156 C       of ammonium in sea-water [mol/kg-SW]
                1157 C       References: Millero (1995) for pressure correction
                1158            zdvi         =  -26.43 _d 0 + t*(0.0889 _d 0 - t*0.905 _d -3)
                1159            zdki         = ( -5.03 _d 0 + t*0.0814 _d 0)*1. _d -3
                1160            pfac  = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
079499602a Jona*1161            akn(i,j,bi,bj) = (total2sw_surf*akn(i,j,bi,bj)
8d230ec051 Jona*1162      &                      * exp(pfac))*sw2total
6acab690ae Jona*1163 
                1164 C -----------------------------------------------------------------------
                1165 C       Apply pressure dependence to the stoichiometric solubility product
                1166 C       of calcite in seawater
                1167 C       References: Millero (1995) for pressure correction
                1168            zdvi      =  -48.76 _d 0 + 0.5304 _d 0*t
                1169            zdki      = (-11.76 _d 0 + 0.3692 _d 0*t)*1. _d -3
                1170            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1171            Ksp_TP_Calc(i,j,bi,bj) = Ksp_TP_Calc(i,j,bi,bj) * exp(pfac)
                1172 
                1173 C -----------------------------------------------------------------------
                1174 C       Apply pressure dependence to the stoichiometric solubility product
                1175 C       of aragonite in seawater
                1176 C       References: Millero (1979) for pressure correction
                1177            zdvi      =  -48.76 _d 0 + 0.5304 _d 0*t  + 2.8 _d 0
                1178            zdki      = (-11.76 _d 0 + 0.3692 _d 0*t)*1. _d -3
                1179            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1180            Ksp_TP_Arag(i,j,bi,bj) = Ksp_TP_Arag(i,j,bi,bj) * exp(pfac)
                1181          else
                1182            bt(i,j,bi,bj)  = 0. _d 0
                1183            st(i,j,bi,bj)  = 0. _d 0
                1184            ft(i,j,bi,bj)  = 0. _d 0
                1185            cat(i,j,bi,bj) = 0. _d 0
                1186            fugf(i,j,bi,bj)= 0. _d 0
                1187            ff(i,j,bi,bj)  = 0. _d 0
                1188            ak0(i,j,bi,bj) = 0. _d 0
                1189            ak1(i,j,bi,bj) = 0. _d 0
                1190            ak2(i,j,bi,bj) = 0. _d 0
                1191            akb(i,j,bi,bj) = 0. _d 0
                1192            ak1p(i,j,bi,bj)= 0. _d 0
                1193            ak2p(i,j,bi,bj)= 0. _d 0
                1194            ak3p(i,j,bi,bj)= 0. _d 0
                1195            aksi(i,j,bi,bj)= 0. _d 0
                1196            akw(i,j,bi,bj) = 0. _d 0
                1197            aks(i,j,bi,bj) = 0. _d 0
                1198            akf(i,j,bi,bj) = 0. _d 0
                1199            akn(i,j,bi,bj) = 0. _d 0
                1200            akhs(i,j,bi,bj)= 0. _d 0
                1201            Ksp_TP_Calc(i,j,bi,bj) = 0. _d 0
                1202            Ksp_TP_Arag(i,j,bi,bj) = 0. _d 0
                1203            aphscale(i,j,bi,bj)    = 0. _d 0
                1204          endif
                1205        enddo
                1206       enddo
                1207 #endif /* CARBONCHEM_SOLVESAPHE */
                1208 
                1209       RETURN
                1210       END
                1211 C     END SUBROUTINE DIC_COEFFS_DEEP
                1212 C      =========================================================================
                1213 
                1214 C      =========================================================================
                1215       SUBROUTINE EQUATION_AT(p_alktot, p_h,      p_dictot,
                1216      &                       p_bortot, p_po4tot, p_siltot,
                1217      &                       p_nh4tot, p_h2stot, p_so4tot,
                1218      &                       p_flutot, p_eqn   , p_deriveqn,
                1219      &                       i, j, k, bi, bj, myIter,
                1220      &                       myThid )
                1221 
                1222       IMPLICIT NONE
                1223 
                1224 C     MITgcm GLobal variables
                1225 #include "SIZE.h"
                1226 #include "EEPARAMS.h"
                1227 #include "DIC_VARS.h"
                1228 
                1229 C      --------------------
                1230 C      Argument variables
                1231 C      --------------------
                1232 
                1233       _RL p_alktot
                1234       _RL p_h
                1235       _RL p_dictot
                1236       _RL p_bortot
                1237       _RL p_po4tot
                1238       _RL p_siltot
                1239       _RL p_nh4tot
                1240       _RL p_h2stot
                1241       _RL p_so4tot
                1242       _RL p_flutot
                1243       _RL p_deriveqn
                1244       _RL p_eqn
                1245       INTEGER i,j,k,bi,bj,myIter,myThid
                1246 
                1247 #ifdef CARBONCHEM_SOLVESAPHE
                1248 
                1249 C      -----------------
                1250 C      Local variables
                1251 C      -----------------
                1252 
                1253       _RL znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
                1254       _RL znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
                1255       _RL znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
                1256       _RL znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
                1257       _RL znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4
                1258       _RL znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s
                1259       _RL znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
                1260       _RL znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
                1261       _RL                                      zalk_wat
                1262 
                1263 C      H2CO3 - HCO3 - CO3 : n=2, m=0
                1264       znumer_dic = 2. _d 0*ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1265      &             + p_h*  ak1(i,j,bi,bj)
                1266       zdenom_dic =      ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1267      &             + p_h*(ak1(i,j,bi,bj) + p_h)
                1268       zalk_dic   = p_dictot * (znumer_dic/zdenom_dic)
                1269 
                1270 C      B(OH)3 - B(OH)4 : n=1, m=0
                1271       znumer_bor =      akb(i,j,bi,bj)
                1272       zdenom_bor =      akb(i,j,bi,bj) + p_h
                1273       zalk_bor   = p_bortot * (znumer_bor/zdenom_bor)
                1274 
                1275 C      H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
                1276       znumer_po4 =
                1277      &        3. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)*ak3p(i,j,bi,bj)
                1278      &      + p_h*(2. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1279      &               + p_h*ak1p(i,j,bi,bj))
                1280       zdenom_po4 =    ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)*ak3p(i,j,bi,bj)
                1281      &           + p_h*(ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1282      &           + p_h*(ak1p(i,j,bi,bj) + p_h))
                1283       zalk_po4   = p_po4tot * (znumer_po4/zdenom_po4 - 1. _d 0)
                1284 C       Zero level of H3PO4 = 1
                1285 
                1286 C      H4SiO4 - H3SiO4 : n=1, m=0
                1287       znumer_sil =      aksi(i,j,bi,bj)
                1288       zdenom_sil =      aksi(i,j,bi,bj) + p_h
                1289       zalk_sil   = p_siltot * (znumer_sil/zdenom_sil)
                1290 
                1291 C      NH4 - NH3 : n=1, m=0
                1292       znumer_nh4 =      akn(i,j,bi,bj)
                1293       zdenom_nh4 =      akn(i,j,bi,bj) + p_h
                1294       zalk_nh4   = p_nh4tot * (znumer_nh4/zdenom_nh4)
                1295 
                1296 C      H2S - HS : n=1, m=0
                1297       znumer_h2s =      akhs(i,j,bi,bj)
                1298       zdenom_h2s =      akhs(i,j,bi,bj) + p_h
                1299       zalk_h2s   = p_h2stot * (znumer_h2s/zdenom_h2s)
                1300 
                1301 C      HSO4 - SO4 : n=1, m=1
                1302       znumer_so4 =      aks(i,j,bi,bj)
                1303       zdenom_so4 =      aks(i,j,bi,bj) + p_h
                1304       zalk_so4   = p_so4tot * (znumer_so4/zdenom_so4 - 1. _d 0)
                1305 
                1306 C      HF - F : n=1, m=1
                1307       znumer_flu =      akf(i,j,bi,bj)
                1308       zdenom_flu =      akf(i,j,bi,bj) + p_h
                1309       zalk_flu   = p_flutot * (znumer_flu/zdenom_flu - 1. _d 0)
                1310 
                1311 C      H2O - OH
                1312       zalk_wat   = akw(i,j,bi,bj)/p_h - p_h/aphscale(i,j,bi,bj)
                1313 
                1314 C      EQUATION_AT =    zalk_dic + zalk_bor + zalk_po4 + zalk_sil &
                1315 C              + zalk_nh4 + zalk_h2s + zalk_so4 + zalk_flu &
                1316 C                  + zalk_wat - p_alktot
                1317       p_eqn =   zalk_dic + zalk_bor + zalk_po4 + zalk_sil
                1318      &         + zalk_nh4 + zalk_h2s + zalk_so4 + zalk_flu
                1319      &         + zalk_wat - p_alktot
                1320 
                1321 C      IF(PRESENT(p_deriveqn)) THEN
                1322 
                1323 C      H2CO3 - HCO3 - CO3 : n=2
                1324         zdnumer_dic = ak1(i,j,bi,bj)*ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1325      &                + p_h*(4. _d 0*ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1326      &                               + p_h*      ak1(i,j,bi,bj))
                1327         zdalk_dic   = -p_dictot*(zdnumer_dic/zdenom_dic**2)
                1328 
                1329 C      B(OH)3 - B(OH)4 : n=1
                1330         zdnumer_bor = akb(i,j,bi,bj)
                1331         zdalk_bor   = -p_bortot*(zdnumer_bor/zdenom_bor**2)
                1332 
                1333 C      H3PO4 - H2PO4 - HPO4 - PO4 : n=3
                1334         zdnumer_po4 = ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1335      &               *ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1336      &               *ak3p(i,j,bi,bj)
                1337      & + p_h*(4. _d 0*ak1p(i,j,bi,bj)*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1338      &               *ak3p(i,j,bi,bj)
                1339      & + p_h*(9. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)*ak3p(i,j,bi,bj)
                1340      &              + ak1p(i,j,bi,bj)*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1341      & + p_h*(4. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1342      & + p_h*         ak1p(i,j,bi,bj))))
                1343         zdalk_po4   = -p_po4tot * (zdnumer_po4/zdenom_po4**2)
                1344 
                1345 C      H4SiO4 - H3SiO4 : n=1
                1346         zdnumer_sil = aksi(i,j,bi,bj)
                1347         zdalk_sil   = -p_siltot * (zdnumer_sil/zdenom_sil**2)
                1348 
                1349 C      NH4 - NH3 : n=1
                1350         zdnumer_nh4 = akn(i,j,bi,bj)
                1351         zdalk_nh4   = -p_nh4tot * (zdnumer_nh4/zdenom_nh4**2)
                1352 
                1353 C      H2S - HS : n=1
                1354         zdnumer_h2s = akhs(i,j,bi,bj)
                1355         zdalk_h2s   = -p_h2stot * (zdnumer_h2s/zdenom_h2s**2)
                1356 
                1357 C      HSO4 - SO4 : n=1
                1358         zdnumer_so4 = aks(i,j,bi,bj)
                1359         zdalk_so4   = -p_so4tot * (zdnumer_so4/zdenom_so4**2)
                1360 
                1361 C      HF - F : n=1
                1362         zdnumer_flu = akf(i,j,bi,bj)
                1363         zdalk_flu   = -p_flutot * (zdnumer_flu/zdenom_flu**2)
                1364 
                1365         p_deriveqn = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil
                1366      &             + zdalk_nh4 + zdalk_h2s + zdalk_so4 + zdalk_flu
                1367      &             - akw(i,j,bi,bj)/p_h**2 - 1. _d 0/aphscale(i,j,bi,bj)
                1368 C      ENDIF
                1369 
                1370 #endif /* CARBONCHEM_SOLVESAPHE */
                1371       RETURN
                1372       END
                1373 C     END SUBROUTINE EQUATION_AT
                1374 C      =========================================================================
                1375 
                1376 C      =========================================================================
                1377       SUBROUTINE SOLVE_AT_FAST(p_alktot, p_dictot, p_bortot,
                1378      &                         p_po4tot, p_siltot, p_nh4tot,
                1379      &                         p_h2stot, p_so4tot, p_flutot,
                1380      &                         p_hini,   p_val,    p_hnew,
                1381      &                         i, j, k, bi, bj,
                1382      &                         debugPrt, myIter, myThid )
                1383 C      =========================================================================
                1384 C      Fast version of SOLVE_AT_GENERAL, without any bounds checking.
                1385 
                1386       IMPLICIT NONE
                1387 
                1388 C     MITgcm GLobal variables
                1389 #include "SIZE.h"
                1390 #include "EEPARAMS.h"
                1391 #include "DIC_VARS.h"
                1392 
                1393 C     _RL SOLVE_AT_FAST
                1394 
                1395 C     --------------------
                1396 C     Argument variables
                1397 C     --------------------
                1398 
                1399       _RL p_alktot
                1400       _RL p_dictot
                1401       _RL p_bortot
                1402       _RL p_po4tot
                1403       _RL p_siltot
                1404       _RL p_nh4tot
                1405       _RL p_h2stot
                1406       _RL p_so4tot
                1407       _RL p_flutot
                1408       _RL p_hini
                1409       _RL p_val
                1410       _RL p_hnew
                1411       INTEGER i, j, k, bi, bj
                1412       LOGICAL debugPrt
                1413       INTEGER myIter, myThid
                1414 
                1415 #ifdef CARBONCHEM_SOLVESAPHE
                1416 
                1417 C     --------------------
                1418 C     Local variables
                1419 C     --------------------
                1420 
                1421       _RL zh_ini, zh, zh_prev, zh_lnfactor
                1422       _RL zhdelta
                1423       _RL zeqn, zdeqndh
                1424 
                1425 C     LOGICAL l_exitnow
                1426       LOGICAL iterate4pH
                1427       _RL pz_exp_threshold
                1428       _RL pp_rdel_ah_target
                1429       INTEGER niter_atfast
                1430 
                1431       pp_rdel_ah_target = 1. _d -8
                1432       pz_exp_threshold = 1.0 _d 0
                1433 C      =========================================================================
                1434 
                1435 C     IF(PRESENT(p_hini)) THEN
                1436       zh_ini = p_hini
                1437 C     ELSE
                1438 C     #if defined(DEBUG_PHSOLVERS)
                1439 C        PRINT*, '[SOLVE_AT_FAST] Calling AHINI_FOR_AT for h_ini'
                1440 C     #endif
                1441 C
                1442 C        CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, zh_ini)
                1443 C
                1444 C     #if defined(DEBUG_PHSOLVERS)
                1445 C        PRINT*, '[SOLVE_AT_FAST] h_ini :', zh_ini
                1446 C     #endif
                1447 C     ENDIF
                1448 
                1449       zh = zh_ini
                1450 C Reset counters of iterations
                1451       niter_atfast    = 0
                1452 
                1453       iterate4pH = .TRUE.
                1454       DO WHILE ( iterate4pH )
                1455 
                1456        niter_atfast = niter_atfast + 1
                1457        IF ( niter_atfast .GT. at_maxniter ) THEN
                1458          zh = -1. _d 0
                1459          iterate4pH = .FALSE.
                1460        ELSE
                1461 
                1462          zh_prev = zh
                1463 
                1464 C        zeqn = EQUATION_AT(p_alktot, zh,       p_dictot, p_bortot,
                1465 C     &                      p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1466 C     &                      p_so4tot, p_flutot, P_DERIVEQN = zdeqndh)
                1467 
                1468         CALL EQUATION_AT( p_alktot, zh      , p_dictot, p_bortot,
                1469      &                    p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1470      &                    p_so4tot, p_flutot,
                1471      &                    zeqn    , zdeqndh,
                1472      &                    i, j, k, bi, bj, myIter, myThid )
                1473 C zh is the root
                1474         iterate4pH = ( zeqn .NE. 0. _d 0 )
                1475        ENDIF
                1476 
                1477        IF ( iterate4pH ) THEN
                1478         zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
                1479         IF(ABS(zh_lnfactor) .GT. pz_exp_threshold) THEN
                1480            zh          = zh_prev*EXP(zh_lnfactor)
                1481         ELSE
                1482            zhdelta     = zh_lnfactor*zh_prev
                1483            zh          = zh_prev + zhdelta
                1484         ENDIF
                1485 
                1486 C     #if defined(DEBUG_PHSOLVERS)
                1487 C        PRINT*, '[SOLVE_AT_FAST] testing zh :', zh, zeqn, zh_lnfactor
                1488 C     #endif
                1489 C
                1490 C      ! Stop iterations once |\delta{[H]}/[H]| < rdel
                1491 C      ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
                1492 C      ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
                1493 C
                1494 C      ! Alternatively:
                1495 C      ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
                1496 C      !             ~ 1/LOG(10) * |\Delta [H]|/[H]
                1497 C      !             < 1/LOG(10) * rdel
                1498 C
                1499 C      ! Hence |zeqn/(zdeqndh*zh)| < rdel
                1500 C
                1501 C      ! rdel <- pp_rdel_ah_target
                1502 
                1503 C        l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
                1504 C        IF(l_exitnow) EXIT
                1505         iterate4pH = ( ABS(zh_lnfactor) .GE. pp_rdel_ah_target )
                1506 
                1507        ENDIF
                1508       ENDDO
                1509 
                1510 C     SOLVE_AT_FAST = zh
                1511       p_hnew = zh
                1512 
                1513 C     IF(PRESENT(p_val)) THEN
                1514         IF(zh .GT. 0. _d 0) THEN
                1515 C           p_val = EQUATION_AT(p_alktot, zh,       p_dictot, p_bortot,
                1516 C     &                          p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1517 C     &                          p_so4tot, p_flutot)
                1518            CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                1519      &                       p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1520      &                       p_so4tot, p_flutot,
                1521      &                       p_val   , zdeqndh,
                1522      &                       i, j, k, bi, bj, myIter, myThid )
                1523         ELSE
                1524 c          p_val = HUGE(1. _d 0)
                1525            p_val = 1. _d +65
                1526         ENDIF
                1527 C     ENDIF
                1528 
                1529 #endif /* CARBONCHEM_SOLVESAPHE */
                1530       RETURN
                1531       END
                1532 C END FUNCTION SOLVE_AT_FAST
                1533 C      =========================================================================
                1534 
                1535 C      =========================================================================
                1536       SUBROUTINE SOLVE_AT_GENERAL(p_alktot, p_dictot, p_bortot,
                1537      &                         p_po4tot, p_siltot, p_nh4tot,
                1538      &                         p_h2stot, p_so4tot, p_flutot,
                1539      &                         p_hini,   p_val,    p_hnew,
                1540      &                         i, j, k, bi, bj,
                1541      &                         debugPrt, myIter, myThid )
                1542 C      Universal pH solver that converges from any given initial value,
                1543 C      determines upper an lower bounds for the solution if required
                1544 
                1545       IMPLICIT NONE
                1546 
                1547 C     MITgcm GLobal variables
                1548 #include "SIZE.h"
                1549 #include "EEPARAMS.h"
                1550 #include "DIC_VARS.h"
                1551 
                1552 C      _RL SOLVE_AT_GENERAL
                1553 
                1554 C     --------------------
                1555 C     Argument variables
                1556 C     --------------------
                1557 
                1558        _RL p_alktot
                1559        _RL p_dictot
                1560        _RL p_bortot
                1561        _RL p_po4tot
                1562        _RL p_siltot
                1563        _RL p_nh4tot
                1564        _RL p_h2stot
                1565        _RL p_so4tot
                1566        _RL p_flutot
                1567        _RL p_hini
                1568        _RL p_hnew
                1569        _RL p_val
                1570       INTEGER i, j, k, bi, bj
                1571       LOGICAL debugPrt
                1572       INTEGER myIter, myThid
                1573 
                1574 #ifdef CARBONCHEM_SOLVESAPHE
                1575 
                1576 C     -----------------
                1577 C     Local variables
                1578 C     -----------------
                1579 
                1580       _RL zh_ini, zh, zh_prev, zh_lnfactor
                1581       _RL p_alknw_inf, p_alknw_sup
                1582       _RL zh_min, zh_max
                1583       _RL zdelta, zh_delta
                1584       _RL zeqn, zdeqndh, zeqn_absmin
                1585       INTEGER niter_atgen
                1586 
                1587 C      LOGICAL       :: l_exitnow
                1588       LOGICAL iterate4pH
                1589       _RL pz_exp_threshold
                1590       _RL pp_rdel_ah_target
                1591 
                1592       pp_rdel_ah_target = 1. _d -8
                1593       pz_exp_threshold = 1. _d 0
                1594 
                1595 C      IF(PRESENT(p_hini)) THEN
                1596         zh_ini = p_hini
                1597 C      ELSE
                1598 C#if defined(DEBUG_PHSOLVERS)
                1599 C        PRINT*, '[SOLVE_AT_GENERAL] Calling AHINI_FOR_AT for h_ini'
                1600 C#endif
                1601 C        CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, zh_ini)
                1602 C#if defined(DEBUG_PHSOLVERS)
                1603 C        PRINT*, '[SOLVE_AT_GENERAL] h_ini :', zh_ini
                1604 C#endif
                1605 C      ENDIF
                1606 #ifdef ALLOW_DEBUG
                1607       IF (debugPrt) CALL DEBUG_CALL('ANW_INFSUP',myThid)
                1608 #endif
                1609       CALL ANW_INFSUP(p_dictot, p_bortot,
                1610      &                p_po4tot, p_siltot,  p_nh4tot, p_h2stot,
                1611      &                p_so4tot, p_flutot,
                1612      &                p_alknw_inf, p_alknw_sup,
                1613      &                i, j, k, bi, bj, myIter, myThid )
                1614 
                1615       zdelta = (p_alktot-p_alknw_inf)**2 + 4. _d 0*akw(i,j,bi,bj)
                1616      &         /aphscale(i,j,bi,bj)
                1617 
                1618       IF(p_alktot .GE. p_alknw_inf) THEN
                1619         zh_min = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_inf
                1620      &           + SQRT(zdelta) )
                1621       ELSE
                1622         zh_min = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_inf)
                1623      &           + SQRT(zdelta) ) / 2. _d 0
                1624       ENDIF
                1625 
                1626       zdelta = (p_alktot-p_alknw_sup)**2 + 4. _d 0*akw(i,j,bi,bj)
                1627      &         /aphscale(i,j,bi,bj)
                1628 
                1629       IF(p_alktot .LE. p_alknw_sup) THEN
                1630         zh_max = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_sup)
                1631      &           + SQRT(zdelta) ) / 2. _d 0
                1632       ELSE
                1633         zh_max = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_sup
                1634      &           + SQRT(zdelta) )
                1635       ENDIF
                1636 
                1637 C#if defined(DEBUG_PHSOLVERS)
                1638 C           PRINT*, '[SOLVE_AT_GENERAL] h_min :', zh_min
                1639 C           PRINT*, '[SOLVE_AT_GENERAL] h_max :', zh_max
                1640 C#endif
                1641 
                1642       zh = MAX(MIN(zh_max, zh_ini), zh_min)
                1643 C     Uncomment this line for the "safe" initialisation test
                1644 C      zh = SQRT(zh_max*zh_min)
                1645 
                1646 C     Reset counters of iterations
                1647       niter_atgen       = 0
                1648 c     zeqn_absmin       = HUGE(1. _d 0)
                1649       zeqn_absmin       = 1. _d +65
                1650 
                1651       iterate4pH = .TRUE.
                1652       DO WHILE ( iterate4pH )
                1653 
                1654        IF ( niter_atgen .GE. at_maxniter ) THEN
                1655            zh = -1. _d 0
                1656            iterate4pH = .FALSE.
                1657        ELSE
                1658 
                1659         zh_prev = zh
                1660 #ifdef ALLOW_DEBUG
                1661         IF (debugPrt) CALL DEBUG_CALL('EQUATION_AT',myThid)
                1662 #endif
                1663 C        zeqn = EQUATION_AT(p_alktot, zh,      p_dictot, p_bortot,
                1664 C     &                   p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1665 C     &                   p_so4tot, p_flutot, P_DERIVEQN = zdeqndh)
                1666         CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                1667      &                    p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1668      &                    p_so4tot, p_flutot,
                1669      &                    zeqn   , zdeqndh,
                1670      &                    i, j, k, bi, bj, myIter, myThid )
                1671 
                1672 C         Adapt bracketing interval
                1673         IF(zeqn .GT. 0. _d 0) THEN
                1674            zh_min = zh_prev
                1675         ELSEIF(zeqn .LT. 0. _d 0) THEN
                1676            zh_max = zh_prev
                1677         ELSE
                1678 C          zh is the root; unlikely but, one never knows
                1679            iterate4pH = .FALSE.
                1680         ENDIF
                1681        ENDIF
                1682 
                1683        IF ( iterate4pH ) THEN
                1684 C         Now determine the next iterate zh
                1685         niter_atgen = niter_atgen + 1
                1686 
                1687         IF(ABS(zeqn) .GE. 0.5 _d 0*zeqn_absmin) THEN
                1688 C          if the function evaluation at the current point is
                1689 C          not decreasing faster than with a bisection step (at least linearly)
                1690 C          in absolute value take one bisection step on [ph_min, ph_max]
                1691 C          ph_new = (ph_min + ph_max)/2 _d 0
                1692 C          In terms of [H]_new:
                1693 C            [H]_new = 10**(-ph_new)
                1694 C                   = 10**(-(ph_min + ph_max)/2 _d 0)
                1695 C                   = SQRT(10**(-(ph_min + phmax)))
                1696 C                   = SQRT(zh_max * zh_min)
                1697            zh          = SQRT(zh_max * zh_min)
                1698 C Required to test convergence below
                1699            zh_lnfactor = (zh - zh_prev)/zh_prev
                1700         ELSE
                1701 C          dzeqn/dpH = dzeqn/d[H] * d[H]/dpH
                1702 C                   = -zdeqndh * LOG(10) * [H]
                1703 C          \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10))
                1704 
                1705 C          pH_new = pH_old + \deltapH
                1706 
                1707 C          [H]_new = 10**(-pH_new)
                1708 C                 = 10**(-pH_old - \Delta pH)
                1709 C                 = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10)))
                1710 C                 = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10)))
                1711 C                 = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old))
                1712            zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
                1713 
                1714            IF(ABS(zh_lnfactor) .GT. pz_exp_threshold) THEN
                1715              zh          = zh_prev*EXP(zh_lnfactor)
                1716            ELSE
                1717              zh_delta    = zh_lnfactor*zh_prev
                1718              zh          = zh_prev + zh_delta
                1719            ENDIF
                1720 
                1721 C#if defined(DEBUG_PHSOLVERS)
                1722 C           PRINT*, '[SOLVE_AT_GENERAL] testing zh :', zh, zeqn, zh_lnfactor
                1723 C#endif
                1724 
                1725            IF( zh .LT. zh_min ) THEN
                1726 C              if [H]_new < [H]_min
                1727 C              i.e., if ph_new > ph_max then
                1728 C              take one bisection step on [ph_prev, ph_max]
                1729 C              ph_new = (ph_prev + ph_max)/2 _d 0
                1730 C              In terms of [H]_new:
                1731 C              [H]_new = 10**(-ph_new)
                1732 C                     = 10**(-(ph_prev + ph_max)/2 _d 0)
                1733 C                     = SQRT(10**(-(ph_prev + phmax)))
                1734 C                     = SQRT([H]_old*10**(-ph_max))
                1735 C                     = SQRT([H]_old * zh_min)
                1736              zh               = SQRT(zh_prev * zh_min)
                1737 C Required to test convergence below
                1738              zh_lnfactor      = (zh - zh_prev)/zh_prev
                1739            ENDIF
                1740 
                1741            IF( zh .GT. zh_max ) THEN
                1742 C              if [H]_new > [H]_max
                1743 C              i.e., if ph_new < ph_min, then
                1744 C              take one bisection step on [ph_min, ph_prev]
                1745 C              ph_new = (ph_prev + ph_min)/2 _d 0
                1746 C              In terms of [H]_new:
                1747 C              [H]_new = 10**(-ph_new)
                1748 C                     = 10**(-(ph_prev + ph_min)/2 _d 0)
                1749 C                     = SQRT(10**(-(ph_prev + ph_min)))
                1750 C                     = SQRT([H]_old*10**(-ph_min))
                1751 C                     = SQRT([H]_old * zhmax)
                1752              zh               = SQRT(zh_prev * zh_max)
                1753 C Required to test convergence below
                1754              zh_lnfactor      = (zh - zh_prev)/zh_prev
                1755 
                1756            ENDIF
                1757         ENDIF
                1758 
                1759         zeqn_absmin = MIN( ABS(zeqn), zeqn_absmin)
                1760 C        Stop iterations once |\delta{[H]}/[H]| < rdel
                1761 C        <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
                1762 C        |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
                1763 C
                1764 C        Alternatively:
                1765 C        |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
                1766 C                    ~ 1/LOG(10) * |\Delta [H]|/[H]
                1767 C                    < 1/LOG(10) * rdel
                1768 C
                1769 C        Hence |zeqn/(zdeqndh*zh)| < rdel
                1770 C
                1771 C        rdel <-- pp_rdel_ah_target
                1772 C        l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
                1773 C        IF(l_exitnow) EXIT
                1774          iterate4pH = ( ABS(zh_lnfactor) .GE. pp_rdel_ah_target )
                1775 
                1776        ENDIF
                1777       ENDDO
                1778 
                1779 C      SOLVE_AT_GENERAL = zh
                1780       p_hnew = zh
                1781 
                1782 C      IF(PRESENT(p_val)) THEN
                1783 C           p_val = EQUATION_AT(p_alktot, zh,      p_dictot, p_bortot,
                1784 C     &                       p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1785 C     &                       p_so4tot, p_flutot)
                1786 C        ELSE
                1787 C           p_val = HUGE(1. _d 0)
                1788 C       ENDIF
                1789 C      ENDIF
                1790 #ifdef ALLOW_DEBUG
                1791       IF (debugPrt) CALL DEBUG_LEAVE('SOLVE_AT_GENERAL',myThid)
                1792 #endif
                1793 #endif /* CARBONCHEM_SOLVESAPHE */
                1794       RETURN
                1795       END
                1796 C      END FUNCTION SOLVE_AT_GENERAL
                1797 C      =========================================================================
                1798 
                1799 C      =========================================================================
                1800       SUBROUTINE SOLVE_AT_GENERAL_SEC(p_alktot, p_dictot,
                1801      &                           p_bortot, p_po4tot,
                1802      &                           p_siltot, p_nh4tot,
                1803      &                           p_h2stot, p_so4tot,
                1804      &                           p_flutot, p_hini,
                1805      &                           p_val,    p_hnew,
                1806      &                           i, j, k, bi, bj,
                1807      &                           debugPrt, myIter, myThid )
                1808 
                1809 C      Universal pH solver that converges from any given initial value,
                1810 C      determines upper an lower bounds for the solution if required
                1811       IMPLICIT NONE
                1812 
                1813 C     MITgcm GLobal variables
                1814 #include "SIZE.h"
                1815 #include "EEPARAMS.h"
                1816 #include "DIC_VARS.h"
                1817 
                1818 C     --------------------
                1819 C     Argument variables
                1820 C     --------------------
                1821 
                1822       _RL p_alktot
                1823       _RL p_dictot
                1824       _RL p_bortot
                1825       _RL p_po4tot
                1826       _RL p_siltot
                1827       _RL p_nh4tot
                1828       _RL p_h2stot
                1829       _RL p_so4tot
                1830       _RL p_flutot
                1831       _RL p_hini
                1832       _RL p_hnew
                1833       _RL p_val
                1834       INTEGER i, j, k, bi, bj
                1835       LOGICAL debugPrt
                1836       INTEGER myIter, myThid
                1837 
                1838 #ifdef CARBONCHEM_SOLVESAPHE
                1839 
                1840 C     -----------------
                1841 C     Local variables
                1842 C     -----------------
                1843 
                1844       _RL  zh_ini, zh, zh_1, zh_2, zh_delta
                1845       _RL  p_alknw_inf, p_alknw_sup
                1846       _RL  zh_min, zh_max
                1847       _RL  zeqn, zeqn_1, zeqn_2, zeqn_absmin
                1848       _RL  zdelta, zdeqndh
                1849       _RL pp_rdel_ah_target
                1850       INTEGER niter_atsec
                1851 C      LOGICAL       ::  l_exitnow
                1852       LOGICAL iterate4pH
                1853 
                1854       pp_rdel_ah_target = 1. _d -8
                1855 
                1856 C      IF(PRESENT(p_hini)) THEN
                1857         zh_ini = p_hini
                1858 C      ELSE
                1859 C#if defined(DEBUG_PHSOLVERS)
                1860 C        PRINT*, '[SOLVE_AT_GENERAL_SEC] Calling AHINI_FOR_AT for h_ini'
                1861 C#endif
                1862 C       CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, zh_ini)
                1863 C#if defined(DEBUG_PHSOLVERS)
                1864 C        PRINT*, '[SOLVE_AT_GENERAL_SEC] h_ini :', zh_ini
                1865 C#endif
                1866 C      ENDIF
                1867 
                1868       CALL ANW_INFSUP(p_dictot, p_bortot,
                1869      &               p_po4tot, p_siltot,  p_nh4tot, p_h2stot,
                1870      &               p_so4tot, p_flutot,
                1871      &               p_alknw_inf, p_alknw_sup,
                1872      &               i, j, k, bi, bj, myIter, myThid )
                1873 
                1874       zdelta = (p_alktot-p_alknw_inf)**2 + 4. _d 0*akw(i,j,bi,bj)
                1875      &         /aphscale(i,j,bi,bj)
                1876 
                1877       IF(p_alktot .GE. p_alknw_inf) THEN
                1878         zh_min = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_inf
                1879      &           + SQRT(zdelta) )
                1880       ELSE
                1881         zh_min = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_inf)
                1882      &           + SQRT(zdelta) ) / 2. _d 0
                1883       ENDIF
                1884 
                1885       zdelta = (p_alktot-p_alknw_sup)**2 + 4. _d 0*akw(i,j,bi,bj)
                1886      &         /aphscale(i,j,bi,bj)
                1887 
                1888       IF(p_alktot .LE. p_alknw_sup) THEN
                1889         zh_max = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_sup)
                1890      &           + SQRT(zdelta) ) / 2. _d 0
                1891       ELSE
                1892         zh_max = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_sup
                1893      &           + SQRT(zdelta) )
                1894       ENDIF
                1895 
                1896 C#if defined(DEBUG_PHSOLVERS)
                1897 C           PRINT*, '[SOLVE_AT_GENERAL_SEC] h_min :', zh_min
                1898 C           PRINT*, '[SOLVE_AT_GENERAL_SEC] h_max :', zh_max
                1899 C#endif
                1900 
                1901 C     Uncomment this line for the "safe" initialisation test
                1902       zh = MAX(MIN(zh_max, zh_ini), zh_min)
                1903 C      zh = SQRT(zh_max*zh_min)
                1904 
                1905 C      Reset counters of iterations
                1906       niter_atsec       = 0
                1907 
                1908 C      Prepare the secant iterations: two initial (zh, zeqn) pairs are required
                1909 C      We have the starting value, that needs to be completed by the evaluation
                1910 C      of the equation value it produces.
                1911 C      Complete the initial value with its equation evaluation
                1912 C      (will take the role of the $n-2$ iterate at the first secant evaluation)
                1913 
                1914 C     zh_2 is the initial value
                1915       zh_2   = zh
                1916 C      zeqn_2 = EQUATION_AT(p_alktot, zh_2,     p_dictot, p_bortot,
                1917 C     &                   p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1918 C     &                   p_so4tot, p_flutot)
                1919       CALL EQUATION_AT( p_alktot, zh_2,    p_dictot, p_bortot,
                1920      &                  p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1921      &                  p_so4tot, p_flutot,
                1922      &                  zeqn_2  , zdeqndh,
                1923      &                  i, j, k, bi, bj, myIter, myThid )
                1924 
                1925       zeqn_absmin       = ABS(zeqn_2)
                1926 
                1927 C      Adapt bracketing interval and heuristically set zh_1
                1928 
                1929       IF(zeqn_2 .LT. 0. _d 0) THEN
                1930 C                 If zeqn_2 < 0, then we adjust zh_max:
                1931 C                 we can be sure that zh_min < zh_2 < zh_max.
                1932         zh_max = zh_2
                1933 C                 for zh_1, try 25% (0.1 pH units) below the current zh_max,
                1934 C                 but stay above SQRT(zh_min*zh_max), which would be equivalent
                1935 C                 to a bisection step on [pH@zh_min, pH@zh_max]
                1936         zh_1   = MAX(zh_max/1.25 _d 0, SQRT(zh_min*zh_max))
                1937 
                1938       ELSEIF(zeqn_2 .GT. 0. _d 0) THEN
                1939 C                 If zeqn_2 < 0, then we adjust zh_min:
                1940 C                 we can be sure that zh_min < zh_2 < zh_max.
                1941         zh_min = zh_2
                1942 C                 for zh_1, try 25% (0.1 pH units) above the current zh_min,
                1943 C                 but stay below SQRT(zh_min*zh_max) which would be equivalent
                1944 C                 to a bisection step on [pH@zh_min, pH@zh_max]
                1945         zh_1   = MIN(zh_min*1.25 _d 0, SQRT(zh_min*zh_max))
                1946 
                1947       ELSE
                1948 C       we have got the root; unlikely, but one never knows
                1949 C        SOLVE_AT_GENERAL_SEC = zh_2
                1950          p_hnew = zh_2
                1951 C        IF(PRESENT(p_val)) p_val = zeqn_2
                1952          p_val = zeqn_2
                1953         RETURN
                1954       ENDIF
                1955 
                1956 C      We now have the first pair completed (zh_2, zeqn_2).
                1957 C      Define the second one (zh_1, zeqn_1), which is also the first iterate.
                1958 C      zh_1 has already been set above
                1959 
                1960 C     Update counter of iterations
                1961       niter_atsec = 1
                1962 
                1963 C      zeqn_1 = EQUATION_AT(p_alktot, zh_1,      p_dictot, p_bortot,
                1964 C     &                   p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1965 C     &                   p_so4tot, p_flutot)
                1966       CALL EQUATION_AT( p_alktot, zh_1,     p_dictot, p_bortot,
                1967      &                  p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1968      &                  p_so4tot, p_flutot,
                1969      &                  zeqn_1  , zdeqndh,
                1970      &                  i, j, k, bi, bj, myIter, myThid )
                1971 
                1972 C      Adapt bracketing interval: we know that zh_1 <= zh <= zh_max (if zeqn_1>0)
                1973 C      or zh_min <= zh <= zh_1 (if zeqn_1 < 0), so this can always be done
                1974 
                1975       IF(zeqn_1 .GT. 0. _d 0) THEN
                1976         zh_min = zh_1
                1977       ELSEIF(zeqn_1 .LT. 0. _d 0) THEN
                1978         zh_max = zh_1
                1979       ELSE
                1980 C      zh_1 is the root
                1981 C        SOLVE_AT_GENERAL_SEC = zh_1
                1982         p_hnew = zh_1
                1983 C        IF(PRESENT(p_val)) p_val = zeqn_1
                1984         p_val = zeqn_1
                1985       ENDIF
                1986 
                1987       IF(ABS(zeqn_1) .GT. zeqn_absmin) THEN
                1988 C     Swap zh_2 and zh_1 if ABS(zeqn_2) < ABS(zeqn_1)
                1989 C                 so that zh_2 and zh_1 lead to decreasing equation
                1990 C                 values (in absolute value)
                1991         zh     = zh_1
                1992         zeqn   = zeqn_1
                1993         zh_1   = zh_2
                1994         zeqn_1 = zeqn_2
                1995         zh_2   = zh
                1996         zeqn_2 = zeqn
                1997       ELSE
                1998         zeqn_absmin = ABS(zeqn_1)
                1999       ENDIF
                2000 
                2001 C     Pre-calculate the first secant iterate (this is the second iterate)
                2002       niter_atsec = 2
                2003 
                2004       zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1))
                2005       zh = zh_1 + zh_delta
                2006 
                2007 C      Make sure that zh_min < zh < zh_max (if not,
                2008 C      bisect around zh_1 which is the best estimate)
                2009 
                2010       IF (zh .GT. zh_max) THEN
                2011 C       This can only happen if zh_2 < zh_1
                2012 C                 and zeqn_2 > zeqn_1 > 0
                2013         zh = SQRT(zh_1*zh_max)
                2014       ENDIF
                2015 
                2016       IF (zh .LT. zh_min) THEN
                2017 C       This can only happen if zh_2 > zh_1
                2018 C                and zeqn_2 < zeqn_1 < 0
                2019         zh = SQRT(zh_1*zh_min)
                2020       ENDIF
                2021 
                2022       iterate4pH = .TRUE.
                2023       DO WHILE ( iterate4pH )
                2024 
                2025        IF ( niter_atsec .GE. at_maxniter ) THEN
                2026            zh = -1. _d 0
                2027            iterate4pH = .FALSE.
                2028        ELSE
                2029 
                2030 C        zeqn = EQUATION_AT(p_alktot, zh, p_dictot, p_bortot,
                2031 C      &                  p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                2032 C      &                  p_so4tot, p_flutot)
                2033         CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                2034      &                    p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                2035      &                    p_so4tot, p_flutot,
                2036      &                    zeqn    , zdeqndh,
                2037      &                    i, j, k, bi, bj, myIter, myThid )
                2038 
                2039 C       Adapt bracketing interval: since initially, zh_min <= zh <= zh_max
                2040 C       we are sure that zh will improve either bracket, depending on the sign
                2041 C       of zeqn
                2042         IF(zeqn .GT. 0. _d 0) THEN
                2043           zh_min = zh
                2044         ELSEIF(zeqn .LT. 0. _d 0) THEN
                2045           zh_max = zh
                2046         ELSE
                2047 C         zh is the root
                2048           iterate4pH = .FALSE.
                2049         ENDIF
                2050        ENDIF
                2051 
                2052        IF ( iterate4pH ) THEN
                2053 C      Start calculation of next iterate
                2054         niter_atsec = niter_atsec + 1
                2055 
                2056         zh_2   = zh_1
                2057         zeqn_2 = zeqn_1
                2058         zh_1   = zh
                2059         zeqn_1 = zeqn
                2060 
                2061         IF(ABS(zeqn) .GE. 0.5 _d 0*zeqn_absmin) THEN
                2062 C            if the function evaluation at the current point
                2063 C            is not decreasing faster in absolute value than
                2064 C            we may expect for a bisection step, then take
                2065 C            one bisection step on [ph_min, ph_max]
                2066 C            ph_new = (ph_min + ph_max)/2 _d 0
                2067 C            In terms of [H]_new:
                2068 C            [H]_new = 10**(-ph_new)
                2069 C                   = 10**(-(ph_min + ph_max)/2 _d 0)
                2070 C                   = SQRT(10**(-(ph_min + phmax)))
                2071 C                   = SQRT(zh_max * zh_min)
                2072            zh              = SQRT(zh_max * zh_min)
                2073            zh_delta          = zh - zh_1
                2074         ELSE
                2075 C            \Delta H = -zeqn_1*(h_2 - h_1)/(zeqn_2 - zeqn_1)
                2076 C            H_new = H_1 + \Delta H
                2077            zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1))
                2078            zh      = zh_1 + zh_delta
                2079 
                2080 C#if defined(DEBUG_PHSOLVERS)
                2081 C           PRINT*, '[SOLVE_AT_GENERAL_SEC] testing zh :', zh, zeqn, zh_delta
                2082 C#endif
                2083            IF( zh .LT. zh_min ) THEN
                2084 C              if [H]_new < [H]_min
                2085 C              i.e., if ph_new > ph_max then
                2086 C              take one bisection step on [ph_prev, ph_max]
                2087 C              ph_new = (ph_prev + ph_max)/2 _d 0
                2088 C              In terms of [H]_new:
                2089 C              [H]_new = 10**(-ph_new)
                2090 C                     = 10**(-(ph_prev + ph_max)/2 _d 0)
                2091 C                     = SQRT(10**(-(ph_prev + phmax)))
                2092 C                     = SQRT([H]_old*10**(-ph_max))
                2093 C                     = SQRT([H]_old * zh_min)
                2094              zh              = SQRT(zh_1 * zh_min)
                2095              zh_delta         = zh - zh_1
                2096            ENDIF
                2097 
                2098            IF( zh .GT. zh_max ) THEN
                2099 C              if [H]_new > [H]_max
                2100 C              i.e., if ph_new < ph_min, then
                2101 C              take one bisection step on [ph_min, ph_prev]
                2102 C              ph_new = (ph_prev + ph_min)/2 _d 0
                2103 C              In terms of [H]_new:
                2104 C              [H]_new = 10**(-ph_new)
                2105 C                     = 10**(-(ph_prev + ph_min)/2 _d 0)
                2106 C                     = SQRT(10**(-(ph_prev + ph_min)))
                2107 C                     = SQRT([H]_old*10**(-ph_min))
                2108 C                     = SQRT([H]_old * zhmax)
                2109              zh               = SQRT(zh_1 * zh_max)
                2110              zh_delta          = zh - zh_1
                2111            ENDIF
                2112         ENDIF
                2113 
                2114         zeqn_absmin = MIN(ABS(zeqn), zeqn_absmin)
                2115 
                2116 C       Stop iterations once |([H]-[H_1])/[H_1]| < rdel
                2117 C        l_exitnow = (ABS(zh_delta) < pp_rdel_ah_target*zh_1)
                2118 C        IF(l_exitnow) EXIT
                2119         iterate4pH = ( ABS(zh_delta) .GE. pp_rdel_ah_target*zh_1 )
                2120 
                2121        ENDIF
                2122       ENDDO
                2123 
                2124 C      SOLVE_AT_GENERAL_SEC = zh
                2125       p_hnew = zh
                2126 
                2127 C      IF(PRESENT(p_val)) THEN
                2128        IF(zh .GT. 0. _d 0) THEN
                2129 C           p_val = EQUATION_AT(p_alktot, zh,      p_dictot, p_bortot,
                2130 C     &                       p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                2131 C     &                       p_so4tot, p_flutot)
                2132           CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                2133      &                      p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                2134      &                      p_so4tot, p_flutot,
                2135      &                      p_val   , zdeqndh,
                2136      &                      i, j, k, bi, bj, myIter, myThid )
                2137        ELSE
                2138 c         p_val = HUGE(1. _d 0)
                2139           p_val = 1. _d +65
                2140        ENDIF
                2141 C      ENDIF
                2142 
                2143 #endif /* CARBONCHEM_SOLVESAPHE */
                2144       RETURN
                2145       END
                2146 C      END FUNCTION SOLVE_AT_GENERAL_SEC
                2147 C      =========================================================================
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!