Back to home page

darwin3
 
 

    

File indexing completed on 2024-12-17 18:34:04 UTC

view on githubraw file Latest commit 45ab1c33 on 2023-04-19 13:45:57 UTC
f0a72e2151 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002 
                0003 C--  File dic_solvesaphe.F:
                0004 C--   Contents:
                0005 C--   o AHINI_FOR_AT
                0006 C--   o ANW_INFSUP
                0007 C--   o EQUATION_AT
                0008 C--   o CALC_PCO2_SOLVESAPHE
                0009 C--   o DARWIN_COEFFS_SURF
                0010 C--   o DARWIN_COEFFS_DEEP
                0011 C--   o SOLVE_AT_FAST
                0012 C--   o SOLVE_AT_GENERAL
                0013 C--   o SOLVE_AT_GENERAL_SEC
                0014 
                0015 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0016 C    SolveSAPHE is free software: you can redistribute it and/or modify
                0017 C    it under the terms of the GNU Lesser General Public License as published by
                0018 C    the Free Software Foundation, either version 3 of the License, or
                0019 C    (at your option) any later version.
                0020 C
                0021 C    SolveSAPHE is distributed in the hope that it will be useful,
                0022 C    but WITHOUT ANY WARRANTY; without even the implied warranty of
                0023 C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
                0024 C    GNU Lesser General Public License for more details.
                0025 C
                0026 C    You should have received a copy of the GNU Lesser General Public License
                0027 C    along with SolveSAPHE.  If not, see <http://www.gnu.org/licenses/>.
                0028 C
                0029 C    Copyright 2013 Guy Munhoven
                0030 C
                0031 C    From the paper: Munhoven, G (2013),
                0032 
** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 33.


C        Mathematics of the total alkalinity–pH equation – pathway to robust
                0033 C        and universal solution algorithms: the SolveSAPHE package v1.0.1.
                0034 
** Warning **
Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 35.


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

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

The LXR team
Valid CSS 2.1! Valid HTML 4.01!