pkg/bling/bling_solvesaphe.F

e0f9a7ba0b Matt*0001 #include "BLING_OPTIONS.h"
                0002
                0003 C--  File bling_solvesaphe.F:
                0004 C--   Contents:
                0005 C--   o AHINI_FOR_AT
                0006 C--   o ANW_INFSUP
                0007 C--   o CALC_PCO2_SOLVESAPHE
                0008 C--   o DIC_COEFFS_SURF
                0009 C--   o DIC_COEFFS_DEEP
15ec8028f7 aver*0010 C--   o EQUATION_AT
e0f9a7ba0b Matt*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 "BLING_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 "BLING_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
                0229 C INTERFACE: ==========================================================
                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 "BLING_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)
e0f9a7ba0b Matt*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 "BLING_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
                0412       _RL inv_t_k
                0413       _RL ion_st, is_2, sqrtis
                0414       _RL s_2, sqrts, s_15, scl, s35
                0415       _RL log_fw2sw
                0416       _RL B, B1
                0417       _RL delta
                0418       _RL P1atm
                0419       _RL RT
                0420       _RL Rgas
                0421 C conversions for different pH scales
                0422       _RL total2free
                0423       _RL free2sw
                0424       _RL total2sw
                0425
                0426         do i=imin,imax
                0427          do j=jmin,jmax
                0428           if (hFacC(i,j,1,bi,bj).gt.0. _d 0) then
                0429            t = ttemp(i,j)
                0430            s = stemp(i,j)
                0431 C terms used more than once for:
                0432 C temperature
                0433            t_k         = 273.15 _d 0 + t
                0434            t_k_o_100   = t_k/100. _d 0
                0435            t_k_o_100_2 = t_k_o_100*t_k_o_100
                0436            inv_t_k=1.0 _d 0/t_k
                0437            dlog_t_k=LOG(t_k)
                0438 C ionic strength (converted to kgsw)
                0439            ion_st=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)
                0440            is_2=ion_st*ion_st
                0441            sqrtis=sqrt(ion_st)
                0442 C salinity
                0443            s_2  = s*s
                0444            sqrts=sqrt(s)
                0445            s_15 = s*sqrts
                0446            scl  = s/1.80655 _d 0
                0447            s35  = s/35. _d 0
                0448
                0449            log_fw2sw = LOG( 1. _d 0 - 0.001005 _d 0*s )
                0450
                0451       IF ( selectBTconst.EQ.1 ) THEN
                0452 C -----------------------------------------------------------------------
                0453 C       Calculate the total borate concentration in mol/kg-SW
                0454 C       given the salinity of a sample
                0455 C       Ref: Uppström (1974), cited by  Dickson et al. (2007, chapter 5, p 10)
                0456 C            Millero (1982) cited in Millero (1995)
                0457 C       pH scale  : N/A
                0458            bt(i,j,bi,bj) = 0.000232 _d 0* scl/10.811 _d 0
                0459       ELSEIF ( selectBTconst.EQ.2 ) THEN
                0460 C -----------------------------------------------------------------------
                0461 C       Calculate the total borate concentration in mol/kg-SW
                0462 C       given the salinity of a sample
                0463 C       References: New formulation from Lee et al (2010)
                0464 C       pH scale  : N/A
                0465            bt(i,j,bi,bj) = 0.0002414 _d 0* scl/10.811 _d 0
                0466       ENDIF
                0467
                0468       IF ( selectFTconst.EQ.1 ) THEN
                0469 C -----------------------------------------------------------------------
                0470 C       Calculate the total fluoride concentration in mol/kg-SW
                0471 C       given the salinity of a sample
                0472 C       References: Riley (1965)
                0473 C       pH scale  : N/A
                0474            ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0
                0475       ELSEIF ( selectFTconst.EQ.2 ) THEN
                0476 C -----------------------------------------------------------------------
                0477 C       Calculate the total fluoride concentration in mol/kg-SW
                0478 C       given the salinity of a sample
                0479 C       References: Culkin (1965) (???)
                0480 C       pH scale  : N/A
                0481            ft(i,j,bi,bj) = 0.000068 _d 0*(s35)
                0482       ENDIF
                0483
                0484 C -----------------------------------------------------------------------
                0485 C       Calculate the total sulfate concentration in mol/kg-SW
                0486 C       given the salinity of a sample
                0487 C       References: Morris & Riley (1966) quoted in Handbook (2007)
                0488 C       pH scale  : N/A
                0489            st(i,j,bi,bj) = 0.14 _d 0 * scl/96.062 _d 0
                0490
                0491 C -----------------------------------------------------------------------
                0492 C       Calculate the total calcium concentration in mol/kg-SW
                0493 C       given the salinity of a sample
                0494 C       References: Culkin (1965) (???)
                0495 C       pH scale  : N/A
                0496            cat(i,j,bi,bj) = 0.010282 _d 0*(s35)
                0497
                0498 C -----------------------------------------------------------------------
                0499 C       Calculate K0 in (mol/kg-SW)/atmosphere
                0500 C       References: Weiss (1979) [(mol/kg-SW)/atm]
                0501 C       pH scale  : N/A
                0502 C       Note      : currently no pressure correction
                0503            ak0(i,j,bi,bj)  = EXP( 93.4517 _d 0/t_k_o_100 - 60.2409 _d 0
                0504      &                 + 23.3585 _d 0*LOG(t_k_o_100)
                0505      &                 + s * (0.023517 _d 0 - 0.023656 _d 0*t_k_o_100
                0506      &                 + 0.0047036 _d 0*t_k_o_100_2))
                0507
                0508 C------------------------------------------------------------------------
                0509 C       Calculate f = k0(1-pH2O)*correction term for non-ideality
                0510 C       References: Weiss & Price (1980, Mar. Chem., 8, 347-359
                0511 C                   Eq 13 with table 6 values)
                0512 C       pH scale  : N/A
                0513            ff(i,j,bi,bj) = exp(-162.8301 _d 0 + 218.2968 _d 0/t_k_o_100
                0514      &          + 90.9241 _d 0*log(t_k_o_100) - 1.47696 _d 0*t_k_o_100_2
                0515      &          + s * (.025695 _d 0 - .025225 _d 0*t_k_o_100
                0516      &          + 0.0049867 _d 0*t_k_o_100_2))
                0517
                0518 C------------------------------------------------------------------------
                0519 C       Calculate Fugacity Factor needed for non-ideality in ocean
                0520 C       References: Weiss (1974) Marine Chemistry
                0521 C       pH scale  : N/A
                0522            P1atm = 1.01325 _d 0 ! bars
                0523            Rgas = 83.1451 _d 0  ! bar*cm3/(mol*K)
                0524            RT = Rgas*t_k
                0525            delta = (57.7 _d 0 - 0.118 _d 0*t_k)
                0526            B1 = -1636.75 _d 0 + 12.0408 _d 0*t_k - 0.0327957 _d 0*t_k*t_k
                0527            B  = B1 + 3.16528 _d 0*t_k*t_k*t_k*(0.00001 _d 0)
                0528
                0529 C   "x2" term often neglected (assumed=1) in applications of Weiss's (1974) eq.9
                0530 C    x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite)
                0531            fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * P1atm / RT)
                0532       IF ( selectK1K2const.EQ.1 ) THEN
                0533 C -----------------------------------------------------------------------
                0534 C       Calculate first dissociation constant of carbonic acid
                0535 C       in mol/kg-SW on the Seawater pH-scale.
                0536 C       References: Millero (1995, eq 35 -- pK1(MEHR));
                0537 C                   Mehrbach et al. (1973) data
                0538 C       pH scale:   Seawater
                0539            ak1(i,j,bi,bj)=10.**(-1. _d 0*(3670.7 _d 0*inv_t_k
                0540      &          - 62.008 _d 0 + 9.7944 _d 0*dlog_t_k
                0541      &          - 0.0118 _d 0 * s + 0.000116 _d 0*s_2))
                0542 C       Calculate second dissociation constant of carbonic acid
                0543 C       in mol/kg-SW on the Seawater pH-scale.
                0544 C       References: Millero (1995, eq 36 -- pK2(MEHR))
                0545 C                   Mehrbach et al. (1973) data
                0546 C       pH scale:   Seawater
                0547            ak2(i,j,bi,bj)=10.**(-1. _d 0*(1394.7 _d 0*inv_t_k
                0548      &           + 4.777 _d 0 - 0.0184 _d 0*s + 0.000118 _d 0*s_2))
                0549       ELSEIF ( selectK1K2const.EQ.2 ) THEN
                0550 C -----------------------------------------------------------------------
                0551 C       Calculate first dissociation constant of carbonic acid
                0552 C       in mol/kg-SW on the Total pH-scale.
                0553 C       References: Roy et al. (1993) -- also Handbook (1994)
                0554 C       pH scale  : Total
                0555 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0556            ak1(i,j,bi,bj)  =  EXP(-2307.1255 _d 0*inv_t_k + 2.83655 _d 0
                0557      &              - 1.5529413 _d 0*dlog_t_k
                0558      &              + (-4.0484 _d 0*inv_t_k - 0.20760841)*sqrts
                0559      &              + 0.08468345*s
                0560      &              - 0.00654208*s_15
                0561      &              + log_fw2sw )
                0562 C       Calculate second dissociation constant of carbonic acid
                0563 C       in mol/kg-SW on the Total pH-scale.
                0564 C       References: Roy et al. (1993) -- also Handbook (1994)
                0565 C       pH scale  : Total
                0566 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0567            ak2(i,j,bi,bj) = EXP(-3351.6106 _d 0*inv_t_k - 9.226508 _d 0
                0568      &              - 0.2005743 _d 0*dlog_t_k
                0569      &              + ( -23.9722 _d 0*inv_t_k - 0.106901773 _d 0)*sqrts
                0570      &              + 0.1130822*s - 0.00846934 _d 0*s_15
                0571      &              + log_fw2sw )
                0572       ELSEIF ( selectK1K2const.EQ.3 ) THEN
                0573 C -----------------------------------------------------------------------
                0574 C       Calculate first dissociation constant of carbonic acid
                0575 C       in mol/kg-SW on the Seawater pH-scale.
                0576 C       References: Millero (1995, eq 50 -- ln K1(COM))
                0577 C       pH scale:   Seawater
                0578            ak1(i,j,bi,bj)  = EXP(2.18867 _d 0 - 2275.0360 _d 0*inv_t_k
                0579      &              - 1.468591 _d 0*dlog_t_k
                0580      &              + ( -0.138681 _d 0 - 9.33291 _d 0*inv_t_k)*sqrts
                0581      &              + 0.0726483 _d 0*s - 0.00574938 _d 0*s_15)
                0582 C       Calculate second dissociation constant of carbonic acid
                0583 C       in mol/kg-SW on the Seawater pH-scale.
                0584 C       References: Millero (1995, eq 51 -- ln K2(COM))
                0585 C       pH scale:   Seawater
                0586            ak2(i,j,bi,bj)  = EXP(-0.84226 _d 0 - 3741.1288 _d 0*inv_t_k
                0587      &              -  1.437139 _d 0*dlog_t_k
                0588      &              + (-0.128417 _d 0 - 24.41239 _d 0*inv_t_k)*sqrts
                0589      &              +  0.1195308 _d 0*s - 0.00912840 _d 0*s_15)
                0590       ELSEIF ( selectK1K2const.EQ.4 ) THEN
                0591 C -----------------------------------------------------------------------
                0592 C       Calculate first dissociation constant of carbonic acid
                0593 C       in mol/kg-SW on the Total pH-scale.
                0594 C       Suitable when 2 < T < 35 and 19 < S < 43
                0595 C       References: Luecker et al. (2000) -- also Handbook (2007)
                0596 C       pH scale:   Total
                0597            ak1(i,j,bi,bj)  = 10. _d 0**( 61.2172 _d 0
                0598      &                 - 3633.86 _d 0*inv_t_k - 9.67770 _d 0*dlog_t_k
                0599      &                 + s*(0.011555 - s*0.0001152 _d 0))
                0600 C       Calculate second dissociation constant of carbonic acid
                0601 C       in mol/kg-SW on the Total pH-scale.
                0602 C       Suitable when 2 < T < 35 and 19 < S < 43
                0603 C       References: Luecker et al. (2000) -- also Handbook (2007)
                0604 C       pH scale:   Total
                0605            ak2(i,j,bi,bj)  = 10. _d 0**(-25.9290 _d 0
                0606      &                 - 471.78 _d 0*inv_t_k + 3.16967 _d 0*dlog_t_k
                0607      &                 + s*(0.01781 _d 0 - s*0.0001122 _d 0))
                0608       ELSEIF ( selectK1K2const.EQ.5 ) THEN
                0609 C -----------------------------------------------------------------------
                0610 C       Calculate first dissociation constant of carbonic acid
                0611 C       in mol/kg-SW on the Seawater pH-scale.
                0612 C       Suitable when 0 < T < 50 and 1 < S < 50
                0613 C       References: Millero (2010, Mar. Fresh Wat. Res.)
                0614 C                   Millero (1979) pressure correction
                0615 C       pH scale:   Seawater
                0616            ak1(i,j,bi,bj) = 10.0 _d 0**(-1*(6320.813 _d 0*inv_t_k
                0617      &      + 19.568224 _d 0*dlog_t_k -126.34048 _d 0
                0618      &      + 13.4038 _d 0*sqrts + 0.03206 _d 0*s - (5.242 _d -5)*s_2
                0619      &      + (-530.659 _d 0*sqrts - 5.8210 _d 0*s)*inv_t_k
                0620      &      -2.0664 _d 0*sqrts*dlog_t_k))
                0621 C       Calculate second dissociation constant of carbonic acid
                0622 C       in mol/kg-SW on the Seawater pH-scale.
                0623 C       Suitable when 0 < T < 50 and 1 < S < 50
                0624 C       References: Millero (2010, Mar. Fresh Wat. Res.)
                0625 C                   Millero (1979) pressure correction
                0626 C       pH scale:   Seawater
                0627            ak2(i,j,bi,bj) = 10.0 _d 0**(-1*(5143.692 _d 0*inv_t_k
                0628      &     + 14.613358 _d 0*dlog_t_k -90.18333 _d 0
                0629      &     + 21.3728 _d 0*sqrts + 0.1218 _d 0*s - (3.688 _d -4)*s_2
                0630      &     + (-788.289 _d 0*sqrts - 19.189 _d 0*s)*inv_t_k
                0631      &     -3.374 _d 0*sqrts*dlog_t_k))
                0632       ELSEIF ( selectK1K2const.EQ.6 ) THEN
                0633 C -----------------------------------------------------------------------
                0634 C       Calculate first dissociation constant of carbonic acid
                0635 C       in mol/kg-SW on the Seawater pH-scale.
                0636 C       Suitable when 0 < T < 50 and 1 < S < 50
                0637 C       References: Waters, Millero, Woosley (Mar. Chem., 165, 66-67, 2014)
                0638 C                   Millero (1979) pressure correction
                0639 C       pH scale:   Seawater
                0640            ak1(i,j,bi,bj) = 10.0 _d 0**(-1.*(6320.813 _d 0*inv_t_k
                0641      &     + 19.568224 _d 0*dlog_t_k -126.34048 _d 0
                0642      &     + 13.409160 _d 0*sqrts + 0.031646 _d 0*s - (5.1895 _d -5)*s_2
                0643      &     + (-531.3642 _d 0*sqrts - 5.713 _d 0*s)*inv_t_k
                0644      &     -2.0669166 _d 0*sqrts*dlog_t_k))
                0645 C       Calculate second dissociation constant of carbonic acid
                0646 C       in mol/kg-SW on the Seawater pH-scale.
                0647 C       Suitable when 0 < T < 50 and 1 < S < 50
                0648 C       References: Waters, Millero, Woosley (Mar. Chem., 165, 66-67, 2014)
                0649 C                   Millero (1979) pressure correction
                0650 C       pH scale:   Seawater
                0651            ak2(i,j,bi,bj) = 10.0 _d 0**(-1.*
                0652      &     ( 5143.692 _d 0*inv_t_k + 14.613358 _d 0*dlog_t_k
                0653      &      - 90.18333 _d 0 + 21.225890 _d 0*sqrts + 0.12450870 _d 0*s
                0654      &      - (3.7243 _d -4)*s_2
                0655      &      + (-779.3444 _d 0*sqrts - 19.91739 _d 0*s)*inv_t_k
                0656      &      - 3.3534679 _d 0*sqrts*dlog_t_k ) )
                0657       ENDIF /* selectK1K2const */
                0658 C -----------------------------------------------------------------------
                0659 C       Calculate boric acid dissociation constant KB
                0660 C       in mol/kg-SW on the total pH-scale.
                0661 C       References: Dickson (1990, eq. 23) -- also Handbook (2007, eq. 37)
                0662 C       pH scale  : total
                0663            akb(i,j,bi,bj)  = EXP(( -8966.90 _d 0 - 2890.53 _d 0*sqrts
                0664      &      -77.942 _d 0*s + 1.728 _d 0*s_15 - 0.0996 _d 0*s_2 )*inv_t_k
                0665      &      + (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s)
                0666      &      + (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s)*
                0667      &      dlog_t_k + 0.053105 _d 0*sqrts*t_k )
                0668 C -----------------------------------------------------------------------
                0669 C       Calculate the first dissociation constant
                0670 C       of phosphoric acid (H3PO4) in seawater
                0671 C       References: Yao and Millero (1995)
                0672 C       pH scale  : Seawater
                0673            ak1p(i,j,bi,bj) = EXP(115.54 _d 0
                0674      &              - 4576.752 _d 0*inv_t_k - 18.453 _d 0*dlog_t_k
                0675      &              + ( 0.69171 _d 0 -  106.736 _d 0*inv_t_k)*sqrts
                0676      &              + (-0.01844 _d 0 -  0.65643 _d 0*inv_t_k)*s)
                0677 C -----------------------------------------------------------------------
                0678 C       Calculate the second dissociation constant
                0679 C       of phosphoric acid (H3PO4) in seawater
                0680 C       References: Yao and Millero (1995)
                0681 C       pH scale  : Seawater
                0682            ak2p(i,j,bi,bj) = EXP( 172.1033 _d 0
                0683      &              - 8814.715 _d 0*inv_t_k
                0684      &              -   27.927 _d 0*dlog_t_k
                0685      &              + (  1.3566 _d 0 -  160.340 _d 0*inv_t_k)*sqrts
                0686      &              + (-0.05778 _d 0 +  0.37335 _d 0*inv_t_k)*s)
                0687 C -----------------------------------------------------------------------
                0688 C       Calculate the third dissociation constant
                0689 C       of phosphoric acid (H3PO4) in seawater
                0690 C       References: Yao and Millero (1995)
                0691 C       pH scale  : Seawater
                0692            ak3p(i,j,bi,bj) = EXP(-18.126 _d 0  -  3070.75 _d 0*inv_t_k
                0693      &                + ( 2.81197 _d 0 + 17.27039 _d 0*inv_t_k)*sqrts
                0694      &                + (-0.09984 _d 0 - 44.99486 _d 0*inv_t_k)*s)
                0695 C -----------------------------------------------------------------------
                0696 C       Calculate the first dissociation constant
                0697 C       of silicic acid (H4SiO4) in seawater
                0698 C       References: Yao and Millero (1995) cited by Millero (1995)
                0699 C       pH scale  : Seawater (according to Dickson et al, 2007)
                0700 C       Note      : converted here from mol/kg-H2O to mol/kg-sw
                0701            aksi(i,j,bi,bj) = EXP(
                0702      &           117.40 _d 0 - 8904.2 _d 0*inv_t_k
                0703      &         - 19.334 _d 0 * dlog_t_k
                0704      &         + ( 3.5913 _d 0 -  458.79 _d 0*inv_t_k) * sqrtis
                0705      &         + (-1.5998 _d 0 +  188.74 _d 0*inv_t_k) * ion_st
                0706      &         + (0.07871 _d 0 - 12.1652 _d 0*inv_t_k) * ion_st*ion_st
                0707      &         + log_fw2sw )
                0708 C -----------------------------------------------------------------------
                0709 C       Calculate the dissociation constant
                0710 C       of ammonium in sea-water [mol/kg-SW]
                0711 C       References: Yao and Millero (1995)
                0712 C       pH scale  : Seawater
                0713            akn(i,j,bi,bj)  = EXP(-0.25444 _d 0 -  6285.33 _d 0*inv_t_k
                0714      &                + 0.0001635 _d 0 * t_k
                0715      &                + ( 0.46532 _d 0 - 123.7184 _d 0*inv_t_k)*sqrts
                0716      &                + (-0.01992 _d 0 +  3.17556 _d 0*inv_t_k)*s)
                0717 C -----------------------------------------------------------------------
                0718 C       Calculate the dissociation constant of hydrogen sulfide in sea-water
                0719 C       References: Millero et al. (1988) (cited by Millero (1995)
                0720 C       pH scale  : - Seawater (according to Yao and Millero, 1995,
                0721 C                               p. 82: "refitted if necessary")
                0722 C                   - Total (according to Lewis and Wallace, 1998)
                0723 C       Note      : we stick to Seawater here for the time being
                0724 C       Note      : the fits from Millero (1995) and Yao and Millero (1995)
                0725 C                   derive from Millero et al. (1998), with all the coefficients
                0726 C                   multiplied by -ln(10)
                0727            akhs(i,j,bi,bj) = EXP( 225.838 _d 0 - 13275.3 _d 0*inv_t_k
                0728      &               - 34.6435 _d 0 * dlog_t_k
                0729      &               +  0.3449 _d 0*sqrts -  0.0274 _d 0*s)
                0730 C -----------------------------------------------------------------------
                0731 C       Calculate the dissociation constant of hydrogen sulfate (bisulfate)
                0732 C       References: Dickson (1990) -- also Handbook (2007)
                0733 C       pH scale  : free
                0734 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0735            aks(i,j,bi,bj) = EXP(141.328 _d 0
                0736      &                -   4276.1 _d 0*inv_t_k -  23.093 _d 0*dlog_t_k
                0737      &                + ( 324.57 _d 0 - 13856. _d 0*inv_t_k
                0738      &                -   47.986 _d 0*dlog_t_k) * sqrtis
                0739      &                + (-771.54 _d 0 + 35474. _d 0*inv_t_k
                0740      &                +  114.723 _d 0*dlog_t_k) * ion_st
                0741      &                - 2698. _d 0*inv_t_k*ion_st**1.5 _d 0
                0742      &                + 1776. _d 0*inv_t_k*ion_st*ion_st
                0743      &                + log_fw2sw )
                0744       IF ( selectHFconst.EQ.1 ) THEN
                0745 C -----------------------------------------------------------------------
                0746 C       Calculate the dissociation constant \beta_{HF} [(mol/kg-SW)^{-1}]
                0747 C       in (mol/kg-SW)^{-1}, where
                0748 C         \beta_{HF} = \frac{ [HF] }{ [H^{+}] [F^{-}] }
                0749 C       References: Dickson and Riley (1979)
                0750 C       pH scale  : free
                0751 C       Note      : converted here from mol/kg-H2O to mol/kg-SW
                0752            akf(i,j,bi,bj) = EXP(1590.2 _d 0*inv_t_k - 12.641 _d 0
                0753      &                 + 1.525 _d 0*sqrtis
                0754      &                 + log_fw2sw )
                0755       ELSEIF ( selectHFconst.EQ.2 ) THEN
                0756 C -----------------------------------------------------------------------
                0757 C       Calculate the dissociation constant for hydrogen fluoride
                0758 C       in mol/kg-SW
                0759 C       References: Perez and Fraga (1987)
                0760 C       pH scale  : Total (according to Handbook, 2007)
                0761            akf(i,j,bi,bj) = EXP( 874. _d 0*inv_t_k - 9.68 _d 0
                0762      &                           + 0.111 _d 0*sqrts )
                0763       ENDIF
                0764 C -----------------------------------------------------------------------
                0765 C       Calculate the water dissociation constant Kw in (mol/kg-SW)^2
                0766 C       References: Millero (1995) for value at pressc = 0
                0767 C       pH scale  : Seawater
                0768            akw(i,j,bi,bj) =  EXP(148.9802 _d 0 - 13847.26 _d 0*inv_t_k
                0769      &              -   23.6521 _d 0*dlog_t_k
                0770      &              + (-5.977 _d 0 + 118.67 _d 0*inv_t_k
                0771      &              + 1.0495 _d 0*dlog_t_k)*sqrts
                0772      &              -   0.01615 _d 0*s )
                0773 C -----------------------------------------------------------------------
                0774 C pH scale conversion factors and conversions. Everything on total pH scale
                0775            total2free = 1. _d 0/
                0776      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                0777            free2sw = 1. _d 0
                0778      &                 + (st(i,j,bi,bj)/ aks(i,j,bi,bj))
                0779      &                 + (ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free))
                0780            total2sw = total2free * free2sw
                0781            aphscale(i,j,bi,bj) = 1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj)
                0782       IF ( selectK1K2const.NE.2 .AND. selectK1K2const.NE.4 ) THEN
                0783 C Convert to the total pH scale
                0784            ak1(i,j,bi,bj)  = ak1(i,j,bi,bj)/total2sw
                0785            ak2(i,j,bi,bj)  = ak2(i,j,bi,bj)/total2sw
                0786       ENDIF
                0787            ak1p(i,j,bi,bj) = ak1p(i,j,bi,bj)/total2sw
                0788            ak2p(i,j,bi,bj) = ak2p(i,j,bi,bj)/total2sw
                0789            ak3p(i,j,bi,bj) = ak3p(i,j,bi,bj)/total2sw
                0790            aksi(i,j,bi,bj) = aksi(i,j,bi,bj)/total2sw
                0791            akn (i,j,bi,bj) = akn (i,j,bi,bj)/total2sw
                0792            akhs(i,j,bi,bj) = akhs(i,j,bi,bj)/total2sw
                0793            aks (i,j,bi,bj) = aks (i,j,bi,bj)/total2free
                0794       IF ( selectHFconst.EQ.1 ) THEN
                0795            akf (i,j,bi,bj) = akf (i,j,bi,bj)/total2free
                0796       ENDIF
                0797            akw (i,j,bi,bj) = akw (i,j,bi,bj)/total2free
                0798 C -----------------------------------------------------------------------
                0799 C       Calculate the stoichiometric solubility product
                0800 C       of calcite in seawater
                0801 C       References: Mucci (1983)
                0802 C       pH scale  : N/A
                0803 C       Units     : (mol/kg-SW)^2
                0804            Ksp_TP_Calc(i,j,bi,bj) = 10. _d 0**(-171.9065 _d 0
                0805      &             - 0.077993 _d 0*t_k
15ec8028f7 aver*0806      &             + 2839.319 _d 0*inv_t_k + 71.595 _d 0*LOG10(t_k)
e0f9a7ba0b Matt*0807      &             + ( -0.77712 _d 0 + 0.0028426 _d 0*t_k
                0808      &             + 178.34 _d 0*inv_t_k)*sqrts
                0809      &             - 0.07711 _d 0*s + 0.0041249 _d 0*s_15)
                0810 C -----------------------------------------------------------------------
                0811 C       Calculate the stoichiometric solubility product
                0812 C       of aragonite in seawater
                0813 C       References: Mucci (1983)
                0814 C       pH scale  : N/A
                0815 C       Units     : (mol/kg-SW)^2
                0816            Ksp_TP_Arag(i,j,bi,bj) = 10. _d 0**(-171.945 _d 0
                0817      &             - 0.077993 _d 0*t_k
15ec8028f7 aver*0818      &             + 2903.293 _d 0*inv_t_k + 71.595 _d 0*LOG10(t_k)
e0f9a7ba0b Matt*0819      &             + ( -0.068393 _d 0 + 0.0017276 _d 0*t_k
                0820      &             + 88.135 _d 0*inv_t_k)*sqrts
                0821      &             - 0.10018 _d 0*s + 0.0059415 _d 0*s_15)
                0822          else
                0823            bt(i,j,bi,bj)  = 0. _d 0
                0824            st(i,j,bi,bj)  = 0. _d 0
                0825            ft(i,j,bi,bj)  = 0. _d 0
                0826            cat(i,j,bi,bj) = 0. _d 0
                0827            fugf(i,j,bi,bj)= 0. _d 0
                0828            ff(i,j,bi,bj)  = 0. _d 0
                0829            ak0(i,j,bi,bj) = 0. _d 0
                0830            ak1(i,j,bi,bj) = 0. _d 0
                0831            ak2(i,j,bi,bj) = 0. _d 0
                0832            akb(i,j,bi,bj) = 0. _d 0
                0833            ak1p(i,j,bi,bj)= 0. _d 0
                0834            ak2p(i,j,bi,bj)= 0. _d 0
                0835            ak3p(i,j,bi,bj)= 0. _d 0
                0836            aksi(i,j,bi,bj)= 0. _d 0
                0837            akw(i,j,bi,bj) = 0. _d 0
                0838            aks(i,j,bi,bj) = 0. _d 0
                0839            akf(i,j,bi,bj) = 0. _d 0
                0840            akn(i,j,bi,bj) = 0. _d 0
                0841            akhs(i,j,bi,bj)= 0. _d 0
                0842            Ksp_TP_Calc(i,j,bi,bj) = 0. _d 0
                0843            Ksp_TP_Arag(i,j,bi,bj) = 0. _d 0
                0844            aphscale(i,j,bi,bj)    = 0. _d 0
                0845          endif
                0846          end do
                0847         end do
                0848 #endif /* CARBONCHEM_SOLVESAPHE */
                0849       RETURN
                0850       END
                0851 C     END SUBROUTINE DIC_COEFFS_SURF
                0852 C -----------------------------------------------------------------------
                0853 C -----------------------------------------------------------------------
                0854       SUBROUTINE DIC_COEFFS_DEEP(
                0855      &                   ttemp,stemp,
                0856      &                   bi,bj,iMin,iMax,jMin,jMax,
                0857      &                   Klevel,myThid)
                0858 C     Add depth dependence to carbon chemistry coefficients loaded into
                0859 C     common block. Corrections occur on the Seawater pH scale and are
                0860 C     converted back to the total scale
                0861       IMPLICIT NONE
                0862 C     MITgcm GLobal variables
                0863 #include "SIZE.h"
                0864 #include "EEPARAMS.h"
                0865 #include "PARAMS.h"
                0866 #include "GRID.h"
                0867 #include "BLING_VARS.h"
                0868       _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0869       _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0870       INTEGER bi,bj,iMin,iMax,jMin,jMax,Klevel,myThid
                0871 #ifdef CARBONCHEM_SOLVESAPHE
                0872 C LOCAL VARIABLES
                0873       INTEGER i, j, k
                0874       _RL  bdepth
                0875       _RL  cdepth
                0876       _RL  pressc
                0877       _RL t
                0878       _RL s
                0879       _RL zds
                0880       _RL t_k
                0881       _RL t_k_o_100
                0882       _RL t_k_o_100_2
                0883       _RL dlog_t_k
                0884       _RL sqrtis
                0885       _RL sqrts
                0886       _RL s_2
                0887       _RL s_15
                0888       _RL inv_t_k
                0889       _RL ion_st
                0890       _RL is_2
                0891       _RL scl
                0892       _RL zrt
                0893       _RL B1
                0894       _RL B
                0895       _RL delta
                0896       _RL Pzatm
                0897       _RL zdvi
                0898       _RL zdki
                0899       _RL pfac
                0900 C pH scale converstions
                0901       _RL total2free_surf
                0902       _RL free2sw_surf
                0903       _RL total2sw_surf
                0904       _RL total2free
                0905       _RL free2sw
                0906       _RL total2sw
                0907 c determine pressure (bar) from depth
                0908 c 1 BAR at z=0m (atmos pressure)
                0909 c use UPPER surface of cell so top layer pressure = 0 bar
                0910 c for surface exchange coeffs
                0911 cmick..............................
                0912 c        write(6,*)'Klevel ',klevel
                0913         bdepth = 0. _d 0
                0914         cdepth = 0. _d 0
                0915         pressc = 1. _d 0
                0916         do k = 1,Klevel
                0917             cdepth = bdepth + 0.5 _d 0*drF(k)
                0918             bdepth = bdepth + drF(k)
                0919             pressc = 1. _d 0 + 0.1 _d 0*cdepth
                0920         end do
                0921 cmick...................................................
                0922 c        write(6,*)'depth,pressc ',cdepth,pressc
                0923 cmick....................................................
                0924         do i=imin,imax
                0925          do j=jmin,jmax
                0926           if (hFacC(i,j,Klevel,bi,bj).gt.0. _d 0) then
                0927            t = ttemp(i,j)
                0928            s = stemp(i,j)
                0929 C terms used more than once for:
                0930 C temperature
                0931            t_k = 273.15 _d 0 + t
                0932            zrt= 83.14472 _d 0 * t_k
                0933            t_k_o_100 = t_k/100. _d 0
                0934            t_k_o_100_2=t_k_o_100*t_k_o_100
                0935            inv_t_k=1.0 _d 0/t_k
                0936            dlog_t_k=log(t_k)
                0937 C ionic strength
                0938            ion_st=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)
                0939            is_2=ion_st*ion_st
                0940            sqrtis=sqrt(ion_st)
                0941 C salinity
                0942            s_2=s*s
                0943            sqrts=sqrt(s)
                0944            s_15=s**1.5 _d 0
                0945            scl=s/1.80655 _d 0
                0946            zds = s - 34.8 _d 0
                0947            total2free_surf = 1. _d 0/
                0948      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                0949            free2sw_surf = 1. _d 0
                0950      &                 + st(i,j,bi,bj)/ aks(i,j,bi,bj)
                0951      &                 + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free_surf)
                0952            total2sw_surf = total2free_surf * free2sw_surf
                0953 C------------------------------------------------------------------------
                0954 C       Recalculate Fugacity Factor needed for non-ideality in ocean
                0955 C       with pressure dependence.
                0956 C       Reference : Weiss (1974) Marine Chemistry
                0957 C       pH scale  : N/A
                0958            Pzatm = 1.01325 _d 0+pressc ! bars
                0959            delta = (57.7 _d 0 - 0.118 _d 0*t_k)
                0960            B1 = -1636.75 _d 0 + 12.0408 _d 0*t_k -0.0327957 _d 0*t_k*t_k
                0961            B  = B1 + 3.16528 _d 0*t_k*t_k*t_k*(0.00001 _d 0)
                0962 C   "x2" term often neglected (assumed=1) in applications of Weiss's (1974) eq.9
                0963 C    x2 = 1 - x1 = 1 - xCO2 (it is very close to 1, but not quite)
                0964            fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * Pzatm / zrt)
                0965
                0966 C -----------------------------------------------------------------------
                0967 C       Apply pressure dependence to the dissociation constant of hydrogen
                0968 C       sulfate (bisulfate).  Ref: Millero (1995) for pressure correction
                0969            zdvi         =  -18.03 _d 0 + t*(0.0466 _d 0 + t*0.316 _d -3)
                0970            zdki         = ( -4.53 _d 0 + t*0.0900 _d 0)*1. _d -3
                0971            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                0972            aks(i,j,bi,bj) = total2free_surf*aks(i,j,bi,bj) * exp(pfac)
                0973
                0974            total2free = 1. _d 0/
                0975      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                0976
                0977 C free2sw has an additional component from fluoride
                0978            free2sw    = 1. _d 0
                0979      &                 + st(i,j,bi,bj)/ aks(i,j,bi,bj)
                0980
                0981            aks(i,j,bi,bj) = aks(i,j,bi,bj)/total2free
                0982
                0983 C -----------------------------------------------------------------------
                0984 C       Apply pressure dependence to dissociation constant for hydrogen fluoride
                0985 C       References: Millero (1995) for pressure correction
                0986            zdvi   =   -9.78 _d 0 + t*(-0.0090 _d 0 - t*0.942 _d -3)
                0987            zdki   = ( -3.91 _d 0 + t*0.054 _d 0)*1. _d -3
                0988            pfac   = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                0989            akf(i,j,bi,bj) = total2free_surf*akf(i,j,bi,bj) * exp(pfac)
                0990
                0991            free2sw = free2sw
                0992      &               + ft(i,j,bi,bj)/akf(i,j,bi,bj)
                0993            total2sw = total2free * free2sw
                0994
                0995            akf(i,j,bi,bj) = akf(i,j,bi,bj)/total2free
                0996
                0997 C -----------------------------------------------------------------------
                0998 C       Apply pressure dependence to 1rst dissociation constant of carbonic acid
                0999 C       References: Millero (1982) pressure correction
                1000            zdvi =  -25.50 _d 0 -0.151 _d 0*zds + 0.1271 _d 0*t
                1001            zdki = ( -3.08 _d 0 -0.578 _d 0*zds + 0.0877 _d 0*t)*1. _d -3
                1002            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1003            ak1(i,j,bi,bj) = (total2free_surf*ak1(i,j,bi,bj)
                1004      &                      * exp(pfac))/total2sw
                1005
                1006 C -----------------------------------------------------------------------
                1007 C       Apply pressure dependence to 2nd dissociation constant of carbonic acid
                1008 C       References: Millero (1979) pressure correction
                1009            zdvi = -15.82 _d 0 + 0.321 _d 0*zds - 0.0219 _d 0*t
                1010            zdki = ( 1.13 _d 0 - 0.314 _d 0*zds - 0.1475 _d 0*t)*1. _d -3
                1011            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1012            ak2(i,j,bi,bj) = (total2free_surf*ak2(i,j,bi,bj)
                1013      &                      * exp(pfac))/total2sw
                1014
                1015 C -----------------------------------------------------------------------
                1016 C       Apply pressure dependence to the boric acid dissociation constant KB
                1017 C       References: Millero (1979) pressure correction
                1018            zdvi      = -29.48 _d 0 + 0.295 _d 0*zds + 0.1622 _d 0*t
                1019      &                 - 0.002608 _d 0*t*t
                1020            zdki      = (-2.84 _d 0 + 0.354 _d 0*zds)*1. _d -3
                1021            pfac =  (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1022            akb(i,j,bi,bj) = (total2free_surf*akb(i,j,bi,bj)
                1023      &                      * exp(pfac))/total2sw
                1024
                1025 C -----------------------------------------------------------------------
                1026 C       Apply pressure dependence to water dissociation constant Kw in
                1027 C       (mol/kg-SW)^2. Ref.: Millero (pers. comm. 1996) for pressure correction
                1028            zdvi    =  -20.02 _d 0 + 0.1119 _d 0*t - 0.1409 _d -2*t*t
                1029            zdki    = ( -5.13 _d 0 + 0.0794 _d 0*t)*1. _d -3
                1030            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1031            akw(i,j,bi,bj) = (total2free_surf*akw(i,j,bi,bj)
                1032      &                      * exp(pfac))/total2sw
                1033
                1034 C -----------------------------------------------------------------------
                1035 C       Apply pressure dependence to the first dissociation constant
                1036 C       of phosphoric acid (H3PO4) in seawater
                1037 C       References: Millero (1995) for pressure correction
                1038            zdvi      =  -14.51 _d 0 + 0.1211 _d 0*t - 0.321 _d -3*t*t
                1039            zdki      = ( -2.67 _d 0 + 0.0427 _d 0*t)*1. _d -3
                1040            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1041            ak1p(i,j,bi,bj) = (total2free_surf*ak1p(i,j,bi,bj)
                1042      &                       * exp(pfac))/total2sw
                1043
                1044 C -----------------------------------------------------------------------
                1045 C       Apply pressure dependence to the second dissociation constant
                1046 C       of phosphoric acid (H3PO4) in seawater
                1047 C       References: Millero (1995) for pressure correction
                1048            zdvi       =  -23.12 _d 0 + 0.1758 _d 0*t -2.647 _d -3*t*t
                1049            zdki       = ( -5.15 _d 0 +   0.09 _d 0*t)*1. _d -3
                1050            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1051            ak2p(i,j,bi,bj) = (total2free_surf*ak2p(i,j,bi,bj)
                1052      &                       * exp(pfac))/total2sw
                1053
                1054 C -----------------------------------------------------------------------
                1055 C       Apply pressure dependence to the third dissociation constant
                1056 C       of phosphoric acid (H3PO4) in seawater
                1057 C       References: Millero (1995) for pressure correction
                1058            zdvi      =  -26.57 _d 0 + 0.2020 _d 0*t -3.042 _d -3*t*t
                1059            zdki      = ( -4.08 _d 0 + 0.0714 _d 0*t)*1. _d -3
                1060            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1061            ak3p(i,j,bi,bj) = (total2free_surf*ak3p(i,j,bi,bj)
                1062      &                       * exp(pfac))/total2sw
                1063
                1064 C -----------------------------------------------------------------------
                1065 C       Apply pressure dependence to the first dissociation constant
                1066 C       of silicic acid (H4SiO4) in seawater
                1067 C       References: Millero (1979) pressure correction. Note: Pressure
                1068 C        correction estimated to be the same as borate (Millero, 1995)
                1069            zdvi      = -29.48 _d 0 + 0.295 _d 0*zds + 0.1622 _d 0*t
                1070      &                 - 0.002608 _d 0*t*t
                1071            zdki      = (-2.84 _d 0 + 0.354 _d 0*zds)*1. _d -3
                1072            pfac =  (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1073            aksi(i,j,bi,bj) = (total2free_surf*aksi(i,j,bi,bj)
                1074      &                       * exp(pfac))/total2sw
                1075
                1076 C -----------------------------------------------------------------------
                1077 C       Apply pressure dependence to the dissociation constant of hydrogen
                1078 C       sulfide in sea-water
                1079 C       References: Millero (1995) for pressure correction
                1080            zdvi         =  -14.80 _d 0 + t*(0.0020 _d 0 - t*0.400 _d -3)
                1081            zdki         = (  2.89 _d 0 + t*0.054 _d 0)*1. _d -3
                1082            pfac  = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1083            akhs(i,j,bi,bj) = (total2free_surf*akhs(i,j,bi,bj)
                1084      &                       * exp(pfac))/total2sw
                1085
                1086 C -----------------------------------------------------------------------
                1087 C       Apply pressure dependence to the dissociation constant
                1088 C       of ammonium in sea-water [mol/kg-SW]
                1089 C       References: Millero (1995) for pressure correction
                1090            zdvi         =  -26.43 _d 0 + t*(0.0889 _d 0 - t*0.905 _d -3)
                1091            zdki         = ( -5.03 _d 0 + t*0.0814 _d 0)*1. _d -3
                1092            pfac  = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1093            akn(i,j,bi,bj) = (total2free_surf*akn(i,j,bi,bj)
                1094      &                      * exp(pfac))/total2sw
                1095
                1096 C -----------------------------------------------------------------------
                1097 C       Apply pressure dependence to the stoichiometric solubility product
                1098 C       of calcite in seawater
                1099 C       References: Millero (1995) for pressure correction
                1100            zdvi      =  -48.76 _d 0 + 0.5304 _d 0*t
                1101            zdki      = (-11.76 _d 0 + 0.3692 _d 0*t)*1. _d -3
                1102            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1103            Ksp_TP_Calc(i,j,bi,bj) = Ksp_TP_Calc(i,j,bi,bj) * exp(pfac)
                1104
                1105 C -----------------------------------------------------------------------
                1106 C       Apply pressure dependence to the stoichiometric solubility product
                1107 C       of aragonite in seawater
                1108 C       References: Millero (1979) for pressure correction
                1109            zdvi      =  -48.76 _d 0 + 0.5304 _d 0*t  + 2.8 _d 0
                1110            zdki      = (-11.76 _d 0 + 0.3692 _d 0*t)*1. _d -3
                1111            pfac = (-zdvi + zdki*pressc/2. _d 0)*pressc/zrt
                1112            Ksp_TP_Arag(i,j,bi,bj) = Ksp_TP_Arag(i,j,bi,bj) * exp(pfac)
                1113          else
                1114            bt(i,j,bi,bj)  = 0. _d 0
                1115            st(i,j,bi,bj)  = 0. _d 0
                1116            ft(i,j,bi,bj)  = 0. _d 0
                1117            cat(i,j,bi,bj) = 0. _d 0
                1118            fugf(i,j,bi,bj)= 0. _d 0
                1119            ff(i,j,bi,bj)  = 0. _d 0
                1120            ak0(i,j,bi,bj) = 0. _d 0
                1121            ak1(i,j,bi,bj) = 0. _d 0
                1122            ak2(i,j,bi,bj) = 0. _d 0
                1123            akb(i,j,bi,bj) = 0. _d 0
                1124            ak1p(i,j,bi,bj)= 0. _d 0
                1125            ak2p(i,j,bi,bj)= 0. _d 0
                1126            ak3p(i,j,bi,bj)= 0. _d 0
                1127            aksi(i,j,bi,bj)= 0. _d 0
                1128            akw(i,j,bi,bj) = 0. _d 0
                1129            aks(i,j,bi,bj) = 0. _d 0
                1130            akf(i,j,bi,bj) = 0. _d 0
                1131            akn(i,j,bi,bj) = 0. _d 0
                1132            akhs(i,j,bi,bj)= 0. _d 0
                1133            Ksp_TP_Calc(i,j,bi,bj) = 0. _d 0
                1134            Ksp_TP_Arag(i,j,bi,bj) = 0. _d 0
                1135            aphscale(i,j,bi,bj)    = 0. _d 0
                1136          endif
                1137        enddo
                1138       enddo
                1139 #endif /* CARBONCHEM_SOLVESAPHE */
                1140
                1141       RETURN
                1142       END
                1143 C     END SUBROUTINE DIC_COEFFS_DEEP
                1144 C      =========================================================================
                1145
                1146 C      =========================================================================
                1147       SUBROUTINE EQUATION_AT(p_alktot, p_h,      p_dictot,
                1148      &                       p_bortot, p_po4tot, p_siltot,
                1149      &                       p_nh4tot, p_h2stot, p_so4tot,
                1150      &                       p_flutot, p_eqn   , p_deriveqn,
                1151      &                       i, j, k, bi, bj, myIter,
                1152      &                       myThid )
                1153
                1154       IMPLICIT NONE
                1155
                1156 C     MITgcm GLobal variables
                1157 #include "SIZE.h"
                1158 #include "EEPARAMS.h"
                1159 #include "BLING_VARS.h"
                1160
                1161 C      --------------------
                1162 C      Argument variables
                1163 C      --------------------
                1164
                1165       _RL p_alktot
                1166       _RL p_h
                1167       _RL p_dictot
                1168       _RL p_bortot
                1169       _RL p_po4tot
                1170       _RL p_siltot
                1171       _RL p_nh4tot
                1172       _RL p_h2stot
                1173       _RL p_so4tot
                1174       _RL p_flutot
                1175       _RL p_deriveqn
                1176       _RL p_eqn
                1177       INTEGER i,j,k,bi,bj,myIter,myThid
                1178
                1179 #ifdef CARBONCHEM_SOLVESAPHE
                1180
                1181 C      -----------------
                1182 C      Local variables
                1183 C      -----------------
                1184
                1185       _RL znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
                1186       _RL znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
                1187       _RL znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
                1188       _RL znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
                1189       _RL znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4
                1190       _RL znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s
                1191       _RL znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
                1192       _RL znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
                1193       _RL                                      zalk_wat
                1194
                1195 C      H2CO3 - HCO3 - CO3 : n=2, m=0
                1196       znumer_dic = 2. _d 0*ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1197      &             + p_h*  ak1(i,j,bi,bj)
                1198       zdenom_dic =      ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1199      &             + p_h*(ak1(i,j,bi,bj) + p_h)
                1200       zalk_dic   = p_dictot * (znumer_dic/zdenom_dic)
                1201
                1202 C      B(OH)3 - B(OH)4 : n=1, m=0
                1203       znumer_bor =      akb(i,j,bi,bj)
                1204       zdenom_bor =      akb(i,j,bi,bj) + p_h
                1205       zalk_bor   = p_bortot * (znumer_bor/zdenom_bor)
                1206
                1207 C      H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
                1208       znumer_po4 =
                1209      &        3. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)*ak3p(i,j,bi,bj)
                1210      &      + p_h*(2. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1211      &               + p_h*ak1p(i,j,bi,bj))
                1212       zdenom_po4 =    ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)*ak3p(i,j,bi,bj)
                1213      &           + p_h*(ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1214      &           + p_h*(ak1p(i,j,bi,bj) + p_h))
                1215       zalk_po4   = p_po4tot * (znumer_po4/zdenom_po4 - 1. _d 0)
                1216 C       Zero level of H3PO4 = 1
                1217
                1218 C      H4SiO4 - H3SiO4 : n=1, m=0
                1219       znumer_sil =      aksi(i,j,bi,bj)
                1220       zdenom_sil =      aksi(i,j,bi,bj) + p_h
                1221       zalk_sil   = p_siltot * (znumer_sil/zdenom_sil)
                1222
                1223 C      NH4 - NH3 : n=1, m=0
                1224       znumer_nh4 =      akn(i,j,bi,bj)
                1225       zdenom_nh4 =      akn(i,j,bi,bj) + p_h
                1226       zalk_nh4   = p_nh4tot * (znumer_nh4/zdenom_nh4)
                1227
                1228 C      H2S - HS : n=1, m=0
                1229       znumer_h2s =      akhs(i,j,bi,bj)
                1230       zdenom_h2s =      akhs(i,j,bi,bj) + p_h
                1231       zalk_h2s   = p_h2stot * (znumer_h2s/zdenom_h2s)
                1232
                1233 C      HSO4 - SO4 : n=1, m=1
                1234       znumer_so4 =      aks(i,j,bi,bj)
                1235       zdenom_so4 =      aks(i,j,bi,bj) + p_h
                1236       zalk_so4   = p_so4tot * (znumer_so4/zdenom_so4 - 1. _d 0)
                1237
                1238 C      HF - F : n=1, m=1
                1239       znumer_flu =      akf(i,j,bi,bj)
                1240       zdenom_flu =      akf(i,j,bi,bj) + p_h
                1241       zalk_flu   = p_flutot * (znumer_flu/zdenom_flu - 1. _d 0)
                1242
                1243 C      H2O - OH
                1244       zalk_wat   = akw(i,j,bi,bj)/p_h - p_h/aphscale(i,j,bi,bj)
                1245
                1246 C      EQUATION_AT =    zalk_dic + zalk_bor + zalk_po4 + zalk_sil &
                1247 C              + zalk_nh4 + zalk_h2s + zalk_so4 + zalk_flu &
                1248 C                  + zalk_wat - p_alktot
                1249       p_eqn =   zalk_dic + zalk_bor + zalk_po4 + zalk_sil
                1250      &         + zalk_nh4 + zalk_h2s + zalk_so4 + zalk_flu
                1251      &         + zalk_wat - p_alktot
                1252
                1253 C      IF(PRESENT(p_deriveqn)) THEN
                1254
                1255 C      H2CO3 - HCO3 - CO3 : n=2
                1256         zdnumer_dic = ak1(i,j,bi,bj)*ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1257      &                + p_h*(4. _d 0*ak1(i,j,bi,bj)*ak2(i,j,bi,bj)
                1258      &                               + p_h*      ak1(i,j,bi,bj))
                1259         zdalk_dic   = -p_dictot*(zdnumer_dic/zdenom_dic**2)
                1260
                1261 C      B(OH)3 - B(OH)4 : n=1
                1262         zdnumer_bor = akb(i,j,bi,bj)
                1263         zdalk_bor   = -p_bortot*(zdnumer_bor/zdenom_bor**2)
                1264
                1265 C      H3PO4 - H2PO4 - HPO4 - PO4 : n=3
                1266         zdnumer_po4 = ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1267      &               *ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1268      &               *ak3p(i,j,bi,bj)
                1269      & + p_h*(4. _d 0*ak1p(i,j,bi,bj)*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1270      &               *ak3p(i,j,bi,bj)
                1271      & + p_h*(9. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)*ak3p(i,j,bi,bj)
                1272      &              + ak1p(i,j,bi,bj)*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1273      & + p_h*(4. _d 0*ak1p(i,j,bi,bj)*ak2p(i,j,bi,bj)
                1274      & + p_h*         ak1p(i,j,bi,bj))))
                1275         zdalk_po4   = -p_po4tot * (zdnumer_po4/zdenom_po4**2)
                1276
                1277 C      H4SiO4 - H3SiO4 : n=1
                1278         zdnumer_sil = aksi(i,j,bi,bj)
                1279         zdalk_sil   = -p_siltot * (zdnumer_sil/zdenom_sil**2)
                1280
                1281 C      NH4 - NH3 : n=1
                1282         zdnumer_nh4 = akn(i,j,bi,bj)
                1283         zdalk_nh4   = -p_nh4tot * (zdnumer_nh4/zdenom_nh4**2)
                1284
                1285 C      H2S - HS : n=1
                1286         zdnumer_h2s = akhs(i,j,bi,bj)
                1287         zdalk_h2s   = -p_h2stot * (zdnumer_h2s/zdenom_h2s**2)
                1288
                1289 C      HSO4 - SO4 : n=1
                1290         zdnumer_so4 = aks(i,j,bi,bj)
                1291         zdalk_so4   = -p_so4tot * (zdnumer_so4/zdenom_so4**2)
                1292
                1293 C      HF - F : n=1
                1294         zdnumer_flu = akf(i,j,bi,bj)
                1295         zdalk_flu   = -p_flutot * (zdnumer_flu/zdenom_flu**2)
                1296
                1297         p_deriveqn = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil
                1298      &             + zdalk_nh4 + zdalk_h2s + zdalk_so4 + zdalk_flu
                1299      &             - akw(i,j,bi,bj)/p_h**2 - 1. _d 0/aphscale(i,j,bi,bj)
                1300 C      ENDIF
                1301
                1302 #endif /* CARBONCHEM_SOLVESAPHE */
                1303       RETURN
                1304       END
                1305 C     END SUBROUTINE EQUATION_AT
                1306 C      =========================================================================
                1307
                1308 C      =========================================================================
                1309       SUBROUTINE SOLVE_AT_FAST(p_alktot, p_dictot, p_bortot,
                1310      &                         p_po4tot, p_siltot, p_nh4tot,
                1311      &                         p_h2stot, p_so4tot, p_flutot,
                1312      &                         p_hini,   p_val,    p_hnew,
                1313      &                         i, j, k, bi, bj,
                1314      &                         debugPrt, myIter, myThid )
                1315 C      =========================================================================
                1316 C      Fast version of SOLVE_AT_GENERAL, without any bounds checking.
                1317
                1318       IMPLICIT NONE
                1319
                1320 C     MITgcm GLobal variables
                1321 #include "SIZE.h"
                1322 #include "EEPARAMS.h"
                1323 #include "BLING_VARS.h"
                1324
                1325 C     _RL SOLVE_AT_FAST
                1326
                1327 C     --------------------
                1328 C     Argument variables
                1329 C     --------------------
                1330
                1331       _RL p_alktot
                1332       _RL p_dictot
                1333       _RL p_bortot
                1334       _RL p_po4tot
                1335       _RL p_siltot
                1336       _RL p_nh4tot
                1337       _RL p_h2stot
                1338       _RL p_so4tot
                1339       _RL p_flutot
                1340       _RL p_hini
                1341       _RL p_val
                1342       _RL p_hnew
                1343       INTEGER i, j, k, bi, bj
                1344       LOGICAL debugPrt
                1345       INTEGER myIter, myThid
                1346
                1347 #ifdef CARBONCHEM_SOLVESAPHE
                1348
                1349 C     --------------------
                1350 C     Local variables
                1351 C     --------------------
                1352
                1353       _RL zh_ini, zh, zh_prev, zh_lnfactor
                1354       _RL zhdelta
                1355       _RL zeqn, zdeqndh
                1356
                1357 C     LOGICAL l_exitnow
                1358       LOGICAL iterate4pH
                1359       _RL pz_exp_threshold
                1360       _RL pp_rdel_ah_target
                1361       INTEGER niter_atfast
                1362
                1363       pp_rdel_ah_target = 1. _d -8
                1364       pz_exp_threshold = 1.0 _d 0
                1365 C      =========================================================================
                1366
                1367 C     IF(PRESENT(p_hini)) THEN
                1368       zh_ini = p_hini
                1369 C     ELSE
                1370 C     #if defined(DEBUG_PHSOLVERS)
                1371 C        PRINT*, '[SOLVE_AT_FAST] Calling AHINI_FOR_AT for h_ini'
                1372 C     #endif
                1373 C
                1374 C        CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, zh_ini)
                1375 C
                1376 C     #if defined(DEBUG_PHSOLVERS)
                1377 C        PRINT*, '[SOLVE_AT_FAST] h_ini :', zh_ini
                1378 C     #endif
                1379 C     ENDIF
                1380
                1381       zh = zh_ini
                1382 C Reset counters of iterations
                1383       niter_atfast    = 0
                1384
                1385       iterate4pH = .TRUE.
                1386       DO WHILE ( iterate4pH )
                1387
                1388        niter_atfast = niter_atfast + 1
                1389        IF ( niter_atfast .GT. at_maxniter ) THEN
                1390          zh = -1. _d 0
                1391          iterate4pH = .FALSE.
                1392        ELSE
                1393
                1394          zh_prev = zh
                1395
                1396 C        zeqn = EQUATION_AT(p_alktot, zh,       p_dictot, p_bortot,
                1397 C     &                      p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1398 C     &                      p_so4tot, p_flutot, P_DERIVEQN = zdeqndh)
                1399
                1400         CALL EQUATION_AT( p_alktot, zh      , p_dictot, p_bortot,
                1401      &                    p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1402      &                    p_so4tot, p_flutot,
                1403      &                    zeqn    , zdeqndh,
                1404      &                    i, j, k, bi, bj, myIter, myThid )
                1405 C zh is the root
                1406         iterate4pH = ( zeqn .NE. 0. _d 0 )
                1407        ENDIF
                1408
                1409        IF ( iterate4pH ) THEN
                1410         zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
                1411         IF(ABS(zh_lnfactor) .GT. pz_exp_threshold) THEN
                1412            zh          = zh_prev*EXP(zh_lnfactor)
                1413         ELSE
                1414            zhdelta     = zh_lnfactor*zh_prev
                1415            zh          = zh_prev + zhdelta
                1416         ENDIF
                1417
                1418 C     #if defined(DEBUG_PHSOLVERS)
                1419 C        PRINT*, '[SOLVE_AT_FAST] testing zh :', zh, zeqn, zh_lnfactor
                1420 C     #endif
                1421 C
                1422 C      ! Stop iterations once |\delta{[H]}/[H]| < rdel
                1423 C      ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
                1424 C      ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
                1425 C
                1426 C      ! Alternatively:
                1427 C      ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
                1428 C      !             ~ 1/LOG(10) * |\Delta [H]|/[H]
                1429 C      !             < 1/LOG(10) * rdel
                1430 C
                1431 C      ! Hence |zeqn/(zdeqndh*zh)| < rdel
                1432 C
                1433 C      ! rdel <- pp_rdel_ah_target
                1434
                1435 C        l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
                1436 C        IF(l_exitnow) EXIT
                1437         iterate4pH = ( ABS(zh_lnfactor) .GE. pp_rdel_ah_target )
                1438
                1439        ENDIF
                1440       ENDDO
                1441
                1442 C     SOLVE_AT_FAST = zh
                1443       p_hnew = zh
                1444
                1445 C     IF(PRESENT(p_val)) THEN
                1446         IF(zh .GT. 0. _d 0) THEN
                1447 C           p_val = EQUATION_AT(p_alktot, zh,       p_dictot, p_bortot,
                1448 C     &                          p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1449 C     &                          p_so4tot, p_flutot)
                1450            CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                1451      &                       p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1452      &                       p_so4tot, p_flutot,
                1453      &                       p_val   , zdeqndh,
                1454      &                       i, j, k, bi, bj, myIter, myThid )
                1455         ELSE
                1456 c          p_val = HUGE(1. _d 0)
                1457            p_val = 1. _d +65
                1458         ENDIF
                1459 C     ENDIF
                1460
                1461 #endif /* CARBONCHEM_SOLVESAPHE */
                1462       RETURN
                1463       END
                1464 C END FUNCTION SOLVE_AT_FAST
                1465 C      =========================================================================
                1466
                1467 C      =========================================================================
                1468       SUBROUTINE SOLVE_AT_GENERAL(p_alktot, p_dictot, p_bortot,
                1469      &                         p_po4tot, p_siltot, p_nh4tot,
                1470      &                         p_h2stot, p_so4tot, p_flutot,
                1471      &                         p_hini,   p_val,    p_hnew,
                1472      &                         i, j, k, bi, bj,
                1473      &                         debugPrt, myIter, myThid )
                1474 C      Universal pH solver that converges from any given initial value,
                1475 C      determines upper an lower bounds for the solution if required
                1476
                1477       IMPLICIT NONE
                1478
                1479 C     MITgcm GLobal variables
                1480 #include "SIZE.h"
                1481 #include "EEPARAMS.h"
                1482 #include "BLING_VARS.h"
                1483
                1484 C      _RL SOLVE_AT_GENERAL
                1485
                1486 C     --------------------
                1487 C     Argument variables
                1488 C     --------------------
                1489
                1490        _RL p_alktot
                1491        _RL p_dictot
                1492        _RL p_bortot
                1493        _RL p_po4tot
                1494        _RL p_siltot
                1495        _RL p_nh4tot
                1496        _RL p_h2stot
                1497        _RL p_so4tot
                1498        _RL p_flutot
                1499        _RL p_hini
                1500        _RL p_hnew
                1501        _RL p_val
                1502       INTEGER i, j, k, bi, bj
                1503       LOGICAL debugPrt
                1504       INTEGER myIter, myThid
                1505
                1506 #ifdef CARBONCHEM_SOLVESAPHE
                1507
                1508 C     -----------------
                1509 C     Local variables
                1510 C     -----------------
                1511
                1512       _RL zh_ini, zh, zh_prev, zh_lnfactor
                1513       _RL p_alknw_inf, p_alknw_sup
                1514       _RL zh_min, zh_max
                1515       _RL zdelta, zh_delta
                1516       _RL zeqn, zdeqndh, zeqn_absmin
                1517       INTEGER niter_atgen
                1518
                1519 C      LOGICAL       :: l_exitnow
                1520       LOGICAL iterate4pH
                1521       _RL pz_exp_threshold
                1522       _RL pp_rdel_ah_target
                1523
                1524       pp_rdel_ah_target = 1. _d -8
                1525       pz_exp_threshold = 1. _d 0
                1526
                1527 C      IF(PRESENT(p_hini)) THEN
                1528         zh_ini = p_hini
                1529 C      ELSE
                1530 C#if defined(DEBUG_PHSOLVERS)
                1531 C        PRINT*, '[SOLVE_AT_GENERAL] Calling AHINI_FOR_AT for h_ini'
                1532 C#endif
                1533 C        CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, zh_ini)
                1534 C#if defined(DEBUG_PHSOLVERS)
                1535 C        PRINT*, '[SOLVE_AT_GENERAL] h_ini :', zh_ini
                1536 C#endif
                1537 C      ENDIF
                1538 #ifdef ALLOW_DEBUG
                1539       IF (debugPrt) CALL DEBUG_CALL('ANW_INFSUP',myThid)
                1540 #endif
                1541       CALL ANW_INFSUP(p_dictot, p_bortot,
                1542      &                p_po4tot, p_siltot,  p_nh4tot, p_h2stot,
                1543      &                p_so4tot, p_flutot,
                1544      &                p_alknw_inf, p_alknw_sup,
                1545      &                i, j, k, bi, bj, myIter, myThid )
                1546
                1547       zdelta = (p_alktot-p_alknw_inf)**2 + 4. _d 0*akw(i,j,bi,bj)
                1548      &         /aphscale(i,j,bi,bj)
                1549
                1550       IF(p_alktot .GE. p_alknw_inf) THEN
                1551         zh_min = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_inf
                1552      &           + SQRT(zdelta) )
                1553       ELSE
                1554         zh_min = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_inf)
                1555      &           + SQRT(zdelta) ) / 2. _d 0
                1556       ENDIF
                1557
                1558       zdelta = (p_alktot-p_alknw_sup)**2 + 4. _d 0*akw(i,j,bi,bj)
                1559      &         /aphscale(i,j,bi,bj)
                1560
                1561       IF(p_alktot .LE. p_alknw_sup) THEN
                1562         zh_max = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_sup)
                1563      &           + SQRT(zdelta) ) / 2. _d 0
                1564       ELSE
                1565         zh_max = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_sup
                1566      &           + SQRT(zdelta) )
                1567       ENDIF
                1568
                1569 C#if defined(DEBUG_PHSOLVERS)
                1570 C           PRINT*, '[SOLVE_AT_GENERAL] h_min :', zh_min
                1571 C           PRINT*, '[SOLVE_AT_GENERAL] h_max :', zh_max
                1572 C#endif
                1573
                1574       zh = MAX(MIN(zh_max, zh_ini), zh_min)
                1575 C     Uncomment this line for the "safe" initialisation test
                1576 C      zh = SQRT(zh_max*zh_min)
                1577
                1578 C     Reset counters of iterations
                1579       niter_atgen       = 0
                1580 c     zeqn_absmin       = HUGE(1. _d 0)
                1581       zeqn_absmin       = 1. _d +65
                1582
                1583       iterate4pH = .TRUE.
                1584       DO WHILE ( iterate4pH )
                1585
                1586        IF ( niter_atgen .GE. at_maxniter ) THEN
                1587            zh = -1. _d 0
                1588            iterate4pH = .FALSE.
                1589        ELSE
                1590
                1591         zh_prev = zh
                1592 #ifdef ALLOW_DEBUG
                1593         IF (debugPrt) CALL DEBUG_CALL('EQUATION_AT',myThid)
                1594 #endif
                1595 C        zeqn = EQUATION_AT(p_alktot, zh,      p_dictot, p_bortot,
                1596 C     &                   p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1597 C     &                   p_so4tot, p_flutot, P_DERIVEQN = zdeqndh)
                1598         CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                1599      &                    p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1600      &                    p_so4tot, p_flutot,
                1601      &                    zeqn   , zdeqndh,
                1602      &                    i, j, k, bi, bj, myIter, myThid )
                1603
                1604 C         Adapt bracketing interval
                1605         IF(zeqn .GT. 0. _d 0) THEN
                1606            zh_min = zh_prev
                1607         ELSEIF(zeqn .LT. 0. _d 0) THEN
                1608            zh_max = zh_prev
                1609         ELSE
                1610 C          zh is the root; unlikely but, one never knows
                1611            iterate4pH = .FALSE.
                1612         ENDIF
                1613        ENDIF
                1614
                1615        IF ( iterate4pH ) THEN
                1616 C         Now determine the next iterate zh
                1617         niter_atgen = niter_atgen + 1
                1618
                1619         IF(ABS(zeqn) .GE. 0.5 _d 0*zeqn_absmin) THEN
                1620 C          if the function evaluation at the current point is
                1621 C          not decreasing faster than with a bisection step (at least linearly)
                1622 C          in absolute value take one bisection step on [ph_min, ph_max]
                1623 C          ph_new = (ph_min + ph_max)/2 _d 0
                1624 C          In terms of [H]_new:
                1625 C            [H]_new = 10**(-ph_new)
                1626 C                   = 10**(-(ph_min + ph_max)/2 _d 0)
                1627 C                   = SQRT(10**(-(ph_min + phmax)))
                1628 C                   = SQRT(zh_max * zh_min)
                1629            zh          = SQRT(zh_max * zh_min)
                1630 C Required to test convergence below
                1631            zh_lnfactor = (zh - zh_prev)/zh_prev
                1632         ELSE
                1633 C          dzeqn/dpH = dzeqn/d[H] * d[H]/dpH
                1634 C                   = -zdeqndh * LOG(10) * [H]
                1635 C          \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10))
                1636
                1637 C          pH_new = pH_old + \deltapH
                1638
                1639 C          [H]_new = 10**(-pH_new)
                1640 C                 = 10**(-pH_old - \Delta pH)
                1641 C                 = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10)))
                1642 C                 = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10)))
                1643 C                 = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old))
                1644            zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
                1645
                1646            IF(ABS(zh_lnfactor) .GT. pz_exp_threshold) THEN
                1647              zh          = zh_prev*EXP(zh_lnfactor)
                1648            ELSE
                1649              zh_delta    = zh_lnfactor*zh_prev
                1650              zh          = zh_prev + zh_delta
                1651            ENDIF
                1652
                1653 C#if defined(DEBUG_PHSOLVERS)
                1654 C           PRINT*, '[SOLVE_AT_GENERAL] testing zh :', zh, zeqn, zh_lnfactor
                1655 C#endif
                1656
                1657            IF( zh .LT. zh_min ) THEN
                1658 C              if [H]_new < [H]_min
                1659 C              i.e., if ph_new > ph_max then
                1660 C              take one bisection step on [ph_prev, ph_max]
                1661 C              ph_new = (ph_prev + ph_max)/2 _d 0
                1662 C              In terms of [H]_new:
                1663 C              [H]_new = 10**(-ph_new)
                1664 C                     = 10**(-(ph_prev + ph_max)/2 _d 0)
                1665 C                     = SQRT(10**(-(ph_prev + phmax)))
                1666 C                     = SQRT([H]_old*10**(-ph_max))
                1667 C                     = SQRT([H]_old * zh_min)
                1668              zh               = SQRT(zh_prev * zh_min)
                1669 C Required to test convergence below
                1670              zh_lnfactor      = (zh - zh_prev)/zh_prev
                1671            ENDIF
                1672
                1673            IF( zh .GT. zh_max ) THEN
                1674 C              if [H]_new > [H]_max
                1675 C              i.e., if ph_new < ph_min, then
                1676 C              take one bisection step on [ph_min, ph_prev]
                1677 C              ph_new = (ph_prev + ph_min)/2 _d 0
                1678 C              In terms of [H]_new:
                1679 C              [H]_new = 10**(-ph_new)
                1680 C                     = 10**(-(ph_prev + ph_min)/2 _d 0)
                1681 C                     = SQRT(10**(-(ph_prev + ph_min)))
                1682 C                     = SQRT([H]_old*10**(-ph_min))
                1683 C                     = SQRT([H]_old * zhmax)
                1684              zh               = SQRT(zh_prev * zh_max)
                1685 C Required to test convergence below
                1686              zh_lnfactor      = (zh - zh_prev)/zh_prev
                1687
                1688            ENDIF
                1689         ENDIF
                1690
                1691         zeqn_absmin = MIN( ABS(zeqn), zeqn_absmin)
                1692 C        Stop iterations once |\delta{[H]}/[H]| < rdel
                1693 C        <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
                1694 C        |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
                1695 C
                1696 C        Alternatively:
                1697 C        |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
                1698 C                    ~ 1/LOG(10) * |\Delta [H]|/[H]
                1699 C                    < 1/LOG(10) * rdel
                1700 C
                1701 C        Hence |zeqn/(zdeqndh*zh)| < rdel
                1702 C
                1703 C        rdel <-- pp_rdel_ah_target
                1704 C        l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
                1705 C        IF(l_exitnow) EXIT
                1706          iterate4pH = ( ABS(zh_lnfactor) .GE. pp_rdel_ah_target )
                1707
                1708        ENDIF
                1709       ENDDO
                1710
                1711 C      SOLVE_AT_GENERAL = zh
                1712       p_hnew = zh
                1713
                1714 C      IF(PRESENT(p_val)) THEN
                1715 C           p_val = EQUATION_AT(p_alktot, zh,      p_dictot, p_bortot,
                1716 C     &                       p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1717 C     &                       p_so4tot, p_flutot)
                1718 C        ELSE
                1719 C           p_val = HUGE(1. _d 0)
                1720 C       ENDIF
                1721 C      ENDIF
                1722 #ifdef ALLOW_DEBUG
                1723       IF (debugPrt) CALL DEBUG_LEAVE('SOLVE_AT_GENERAL',myThid)
                1724 #endif
                1725 #endif /* CARBONCHEM_SOLVESAPHE */
                1726       RETURN
                1727       END
                1728 C      END FUNCTION SOLVE_AT_GENERAL
                1729 C      =========================================================================
                1730
                1731 C      =========================================================================
                1732       SUBROUTINE SOLVE_AT_GENERAL_SEC(p_alktot, p_dictot,
                1733      &                           p_bortot, p_po4tot,
                1734      &                           p_siltot, p_nh4tot,
                1735      &                           p_h2stot, p_so4tot,
                1736      &                           p_flutot, p_hini,
                1737      &                           p_val,    p_hnew,
                1738      &                           i, j, k, bi, bj,
                1739      &                           debugPrt, myIter, myThid )
                1740
                1741 C      Universal pH solver that converges from any given initial value,
                1742 C      determines upper an lower bounds for the solution if required
                1743       IMPLICIT NONE
                1744
                1745 C     MITgcm GLobal variables
                1746 #include "SIZE.h"
                1747 #include "EEPARAMS.h"
                1748 #include "BLING_VARS.h"
                1749
                1750 C     --------------------
                1751 C     Argument variables
                1752 C     --------------------
                1753
                1754       _RL p_alktot
                1755       _RL p_dictot
                1756       _RL p_bortot
                1757       _RL p_po4tot
                1758       _RL p_siltot
                1759       _RL p_nh4tot
                1760       _RL p_h2stot
                1761       _RL p_so4tot
                1762       _RL p_flutot
                1763       _RL p_hini
                1764       _RL p_hnew
                1765       _RL p_val
                1766       INTEGER i, j, k, bi, bj
                1767       LOGICAL debugPrt
                1768       INTEGER myIter, myThid
                1769
                1770 #ifdef CARBONCHEM_SOLVESAPHE
                1771
                1772 C     -----------------
                1773 C     Local variables
                1774 C     -----------------
                1775
                1776       _RL  zh_ini, zh, zh_1, zh_2, zh_delta
                1777       _RL  p_alknw_inf, p_alknw_sup
                1778       _RL  zh_min, zh_max
                1779       _RL  zeqn, zeqn_1, zeqn_2, zeqn_absmin
                1780       _RL  zdelta, zdeqndh
                1781       _RL pp_rdel_ah_target
                1782       INTEGER niter_atsec
                1783 C      LOGICAL       ::  l_exitnow
                1784       LOGICAL iterate4pH
                1785
                1786       pp_rdel_ah_target = 1. _d -8
                1787
                1788 C      IF(PRESENT(p_hini)) THEN
                1789         zh_ini = p_hini
                1790 C      ELSE
                1791 C#if defined(DEBUG_PHSOLVERS)
                1792 C        PRINT*, '[SOLVE_AT_GENERAL_SEC] Calling AHINI_FOR_AT for h_ini'
                1793 C#endif
                1794 C       CALL AHINI_FOR_AT(p_alktot, p_dictot, p_bortot, zh_ini)
                1795 C#if defined(DEBUG_PHSOLVERS)
                1796 C        PRINT*, '[SOLVE_AT_GENERAL_SEC] h_ini :', zh_ini
                1797 C#endif
                1798 C      ENDIF
                1799
                1800       CALL ANW_INFSUP(p_dictot, p_bortot,
                1801      &               p_po4tot, p_siltot,  p_nh4tot, p_h2stot,
                1802      &               p_so4tot, p_flutot,
                1803      &               p_alknw_inf, p_alknw_sup,
                1804      &               i, j, k, bi, bj, myIter, myThid )
                1805
                1806       zdelta = (p_alktot-p_alknw_inf)**2 + 4. _d 0*akw(i,j,bi,bj)
                1807      &         /aphscale(i,j,bi,bj)
                1808
                1809       IF(p_alktot .GE. p_alknw_inf) THEN
                1810         zh_min = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_inf
                1811      &           + SQRT(zdelta) )
                1812       ELSE
                1813         zh_min = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_inf)
                1814      &           + SQRT(zdelta) ) / 2. _d 0
                1815       ENDIF
                1816
                1817       zdelta = (p_alktot-p_alknw_sup)**2 + 4. _d 0*akw(i,j,bi,bj)
                1818      &         /aphscale(i,j,bi,bj)
                1819
                1820       IF(p_alktot .LE. p_alknw_sup) THEN
                1821         zh_max = aphscale(i,j,bi,bj)*(-(p_alktot-p_alknw_sup)
                1822      &           + SQRT(zdelta) ) / 2. _d 0
                1823       ELSE
                1824         zh_max = 2. _d 0*akw(i,j,bi,bj) /( p_alktot-p_alknw_sup
                1825      &           + SQRT(zdelta) )
                1826       ENDIF
                1827
                1828 C#if defined(DEBUG_PHSOLVERS)
                1829 C           PRINT*, '[SOLVE_AT_GENERAL_SEC] h_min :', zh_min
                1830 C           PRINT*, '[SOLVE_AT_GENERAL_SEC] h_max :', zh_max
                1831 C#endif
                1832
                1833 C     Uncomment this line for the "safe" initialisation test
                1834       zh = MAX(MIN(zh_max, zh_ini), zh_min)
                1835 C      zh = SQRT(zh_max*zh_min)
                1836
                1837 C      Reset counters of iterations
                1838       niter_atsec       = 0
                1839
                1840 C      Prepare the secant iterations: two initial (zh, zeqn) pairs are required
                1841 C      We have the starting value, that needs to be completed by the evaluation
                1842 C      of the equation value it produces.
                1843 C      Complete the initial value with its equation evaluation
                1844 C      (will take the role of the $n-2$ iterate at the first secant evaluation)
                1845
                1846 C     zh_2 is the initial value
                1847       zh_2   = zh
                1848 C      zeqn_2 = EQUATION_AT(p_alktot, zh_2,     p_dictot, p_bortot,
                1849 C     &                   p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1850 C     &                   p_so4tot, p_flutot)
                1851       CALL EQUATION_AT( p_alktot, zh_2,    p_dictot, p_bortot,
                1852      &                  p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1853      &                  p_so4tot, p_flutot,
                1854      &                  zeqn_2  , zdeqndh,
                1855      &                  i, j, k, bi, bj, myIter, myThid )
                1856
                1857       zeqn_absmin       = ABS(zeqn_2)
                1858
                1859 C      Adapt bracketing interval and heuristically set zh_1
                1860
                1861       IF(zeqn_2 .LT. 0. _d 0) THEN
                1862 C                 If zeqn_2 < 0, then we adjust zh_max:
                1863 C                 we can be sure that zh_min < zh_2 < zh_max.
                1864         zh_max = zh_2
                1865 C                 for zh_1, try 25% (0.1 pH units) below the current zh_max,
                1866 C                 but stay above SQRT(zh_min*zh_max), which would be equivalent
                1867 C                 to a bisection step on [pH@zh_min, pH@zh_max]
                1868         zh_1   = MAX(zh_max/1.25 _d 0, SQRT(zh_min*zh_max))
                1869
                1870       ELSEIF(zeqn_2 .GT. 0. _d 0) THEN
                1871 C                 If zeqn_2 < 0, then we adjust zh_min:
                1872 C                 we can be sure that zh_min < zh_2 < zh_max.
                1873         zh_min = zh_2
                1874 C                 for zh_1, try 25% (0.1 pH units) above the current zh_min,
                1875 C                 but stay below SQRT(zh_min*zh_max) which would be equivalent
                1876 C                 to a bisection step on [pH@zh_min, pH@zh_max]
                1877         zh_1   = MIN(zh_min*1.25 _d 0, SQRT(zh_min*zh_max))
                1878
                1879       ELSE
                1880 C       we have got the root; unlikely, but one never knows
                1881 C        SOLVE_AT_GENERAL_SEC = zh_2
                1882          p_hnew = zh_2
                1883 C        IF(PRESENT(p_val)) p_val = zeqn_2
                1884          p_val = zeqn_2
                1885         RETURN
                1886       ENDIF
                1887
                1888 C      We now have the first pair completed (zh_2, zeqn_2).
                1889 C      Define the second one (zh_1, zeqn_1), which is also the first iterate.
                1890 C      zh_1 has already been set above
                1891
                1892 C     Update counter of iterations
                1893       niter_atsec = 1
                1894
                1895 C      zeqn_1 = EQUATION_AT(p_alktot, zh_1,      p_dictot, p_bortot,
                1896 C     &                   p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1897 C     &                   p_so4tot, p_flutot)
                1898       CALL EQUATION_AT( p_alktot, zh_1,     p_dictot, p_bortot,
                1899      &                  p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1900      &                  p_so4tot, p_flutot,
                1901      &                  zeqn_1  , zdeqndh,
                1902      &                  i, j, k, bi, bj, myIter, myThid )
                1903
                1904 C      Adapt bracketing interval: we know that zh_1 <= zh <= zh_max (if zeqn_1>0)
                1905 C      or zh_min <= zh <= zh_1 (if zeqn_1 < 0), so this can always be done
                1906
                1907       IF(zeqn_1 .GT. 0. _d 0) THEN
                1908         zh_min = zh_1
                1909       ELSEIF(zeqn_1 .LT. 0. _d 0) THEN
                1910         zh_max = zh_1
                1911       ELSE
                1912 C      zh_1 is the root
                1913 C        SOLVE_AT_GENERAL_SEC = zh_1
                1914         p_hnew = zh_1
                1915 C        IF(PRESENT(p_val)) p_val = zeqn_1
                1916         p_val = zeqn_1
                1917       ENDIF
                1918
                1919       IF(ABS(zeqn_1) .GT. zeqn_absmin) THEN
                1920 C     Swap zh_2 and zh_1 if ABS(zeqn_2) < ABS(zeqn_1)
                1921 C                 so that zh_2 and zh_1 lead to decreasing equation
                1922 C                 values (in absolute value)
                1923         zh     = zh_1
                1924         zeqn   = zeqn_1
                1925         zh_1   = zh_2
                1926         zeqn_1 = zeqn_2
                1927         zh_2   = zh
                1928         zeqn_2 = zeqn
                1929       ELSE
                1930         zeqn_absmin = ABS(zeqn_1)
                1931       ENDIF
                1932
                1933 C     Pre-calculate the first secant iterate (this is the second iterate)
                1934       niter_atsec = 2
                1935
                1936       zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1))
                1937       zh = zh_1 + zh_delta
                1938
                1939 C      Make sure that zh_min < zh < zh_max (if not,
                1940 C      bisect around zh_1 which is the best estimate)
                1941
                1942       IF (zh .GT. zh_max) THEN
                1943 C       This can only happen if zh_2 < zh_1
                1944 C                 and zeqn_2 > zeqn_1 > 0
                1945         zh = SQRT(zh_1*zh_max)
                1946       ENDIF
                1947
                1948       IF (zh .LT. zh_min) THEN
                1949 C       This can only happen if zh_2 > zh_1
                1950 C                and zeqn_2 < zeqn_1 < 0
                1951         zh = SQRT(zh_1*zh_min)
                1952       ENDIF
                1953
                1954       iterate4pH = .TRUE.
                1955       DO WHILE ( iterate4pH )
                1956
                1957        IF ( niter_atsec .GE. at_maxniter ) THEN
                1958            zh = -1. _d 0
                1959            iterate4pH = .FALSE.
                1960        ELSE
                1961
                1962 C        zeqn = EQUATION_AT(p_alktot, zh, p_dictot, p_bortot,
                1963 C      &                  p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1964 C      &                  p_so4tot, p_flutot)
                1965         CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                1966      &                    p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                1967      &                    p_so4tot, p_flutot,
                1968      &                    zeqn    , zdeqndh,
                1969      &                    i, j, k, bi, bj, myIter, myThid )
                1970
                1971 C       Adapt bracketing interval: since initially, zh_min <= zh <= zh_max
                1972 C       we are sure that zh will improve either bracket, depending on the sign
                1973 C       of zeqn
                1974         IF(zeqn .GT. 0. _d 0) THEN
                1975           zh_min = zh
                1976         ELSEIF(zeqn .LT. 0. _d 0) THEN
                1977           zh_max = zh
                1978         ELSE
                1979 C         zh is the root
                1980           iterate4pH = .FALSE.
                1981         ENDIF
                1982        ENDIF
                1983
                1984        IF ( iterate4pH ) THEN
                1985 C      Start calculation of next iterate
                1986         niter_atsec = niter_atsec + 1
                1987
                1988         zh_2   = zh_1
                1989         zeqn_2 = zeqn_1
                1990         zh_1   = zh
                1991         zeqn_1 = zeqn
                1992
                1993         IF(ABS(zeqn) .GE. 0.5 _d 0*zeqn_absmin) THEN
                1994 C            if the function evaluation at the current point
                1995 C            is not decreasing faster in absolute value than
                1996 C            we may expect for a bisection step, then take
                1997 C            one bisection step on [ph_min, ph_max]
                1998 C            ph_new = (ph_min + ph_max)/2 _d 0
                1999 C            In terms of [H]_new:
                2000 C            [H]_new = 10**(-ph_new)
                2001 C                   = 10**(-(ph_min + ph_max)/2 _d 0)
                2002 C                   = SQRT(10**(-(ph_min + phmax)))
                2003 C                   = SQRT(zh_max * zh_min)
                2004            zh              = SQRT(zh_max * zh_min)
                2005            zh_delta          = zh - zh_1
                2006         ELSE
                2007 C            \Delta H = -zeqn_1*(h_2 - h_1)/(zeqn_2 - zeqn_1)
                2008 C            H_new = H_1 + \Delta H
                2009            zh_delta = -zeqn_1/((zeqn_2-zeqn_1)/(zh_2 - zh_1))
                2010            zh      = zh_1 + zh_delta
                2011
                2012 C#if defined(DEBUG_PHSOLVERS)
                2013 C           PRINT*, '[SOLVE_AT_GENERAL_SEC] testing zh :', zh, zeqn, zh_delta
                2014 C#endif
                2015            IF( zh .LT. zh_min ) THEN
                2016 C              if [H]_new < [H]_min
                2017 C              i.e., if ph_new > ph_max then
                2018 C              take one bisection step on [ph_prev, ph_max]
                2019 C              ph_new = (ph_prev + ph_max)/2 _d 0
                2020 C              In terms of [H]_new:
                2021 C              [H]_new = 10**(-ph_new)
                2022 C                     = 10**(-(ph_prev + ph_max)/2 _d 0)
                2023 C                     = SQRT(10**(-(ph_prev + phmax)))
                2024 C                     = SQRT([H]_old*10**(-ph_max))
                2025 C                     = SQRT([H]_old * zh_min)
                2026              zh              = SQRT(zh_1 * zh_min)
                2027              zh_delta         = zh - zh_1
                2028            ENDIF
                2029
                2030            IF( zh .GT. zh_max ) THEN
                2031 C              if [H]_new > [H]_max
                2032 C              i.e., if ph_new < ph_min, then
                2033 C              take one bisection step on [ph_min, ph_prev]
                2034 C              ph_new = (ph_prev + ph_min)/2 _d 0
                2035 C              In terms of [H]_new:
                2036 C              [H]_new = 10**(-ph_new)
                2037 C                     = 10**(-(ph_prev + ph_min)/2 _d 0)
                2038 C                     = SQRT(10**(-(ph_prev + ph_min)))
                2039 C                     = SQRT([H]_old*10**(-ph_min))
                2040 C                     = SQRT([H]_old * zhmax)
                2041              zh               = SQRT(zh_1 * zh_max)
                2042              zh_delta          = zh - zh_1
                2043            ENDIF
                2044         ENDIF
                2045
                2046         zeqn_absmin = MIN(ABS(zeqn), zeqn_absmin)
                2047
                2048 C       Stop iterations once |([H]-[H_1])/[H_1]| < rdel
                2049 C        l_exitnow = (ABS(zh_delta) < pp_rdel_ah_target*zh_1)
                2050 C        IF(l_exitnow) EXIT
                2051         iterate4pH = ( ABS(zh_delta) .GE. pp_rdel_ah_target*zh_1 )
                2052
                2053        ENDIF
                2054       ENDDO
                2055
                2056 C      SOLVE_AT_GENERAL_SEC = zh
                2057       p_hnew = zh
                2058
                2059 C      IF(PRESENT(p_val)) THEN
                2060        IF(zh .GT. 0. _d 0) THEN
                2061 C           p_val = EQUATION_AT(p_alktot, zh,      p_dictot, p_bortot,
                2062 C     &                       p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                2063 C     &                       p_so4tot, p_flutot)
                2064           CALL EQUATION_AT( p_alktot, zh,       p_dictot, p_bortot,
                2065      &                      p_po4tot, p_siltot, p_nh4tot, p_h2stot,
                2066      &                      p_so4tot, p_flutot,
                2067      &                      p_val   , zdeqndh,
                2068      &                      i, j, k, bi, bj, myIter, myThid )
                2069        ELSE
                2070 c         p_val = HUGE(1. _d 0)
                2071           p_val = 1. _d +65
                2072        ENDIF
                2073 C      ENDIF
                2074
                2075 #endif /* CARBONCHEM_SOLVESAPHE */
                2076       RETURN
                2077       END
                2078 C      END FUNCTION SOLVE_AT_GENERAL_SEC
ba0b047096 Mart*2079 C      =========================================================================