pkg/dic/carbon_chem.F

6d54cf9ca1 Ed H*0001 #include "DIC_OPTIONS.h"
daab022f42 Step*0002
e098791e51 Jean*0003 C--  File carbon_chem.F:
f52bc56573 Jean*0004 C--   Contents
                0005 C--   o CALC_PCO2
                0006 C--   o CALC_PCO2_APPROX
                0007 C--   o CARBON_COEFFS
                0008 C--   o CARBON_COEFFS_PRESSURE_DEP
                0009
                0010 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0011
08536d17ba Step*0012 CBOP
                0013 C !ROUTINE: CALC_PCO2
                0014
                0015 C !INTERFACE: ==========================================================
daab022f42 Step*0016        SUBROUTINE CALC_PCO2(
                0017      I                       donewt,inewtonmax,ibrackmax,
                0018      I                       t,s,diclocal,pt,sit,ta,
                0019      I                       k1local,k2local,
                0020      I                       k1plocal,k2plocal,k3plocal,
                0021      I                       kslocal,kblocal,kwlocal,
                0022      I                       ksilocal,kflocal,
d0092a57ac Step*0023      I                       k0local, fugflocal,
daab022f42 Step*0024      I                       fflocal,btlocal,stlocal,ftlocal,
5625485478 Jean*0025      U                       pHlocal,pCO2surfloc,
e098791e51 Jean*0026      I                       i,j,k,bi,bj,myIter,myThid )
daab022f42 Step*0027
08536d17ba Step*0028 C !DESCRIPTION:
f52bc56573 Jean*0029 C  surface ocean inorganic carbon chemistry to OCMIP2
                0030 C  regulations modified from OCMIP2 code;
                0031 C  Mick Follows, MIT, Oct 1999.
08536d17ba Step*0032
d0092a57ac Step*0033 C Apr 2011: fix vapour bug (following Bennington)
                0034
08536d17ba Step*0035 C !USES: ===============================================================
                0036       IMPLICIT NONE
daab022f42 Step*0037 #include "SIZE.h"
                0038 #include "EEPARAMS.h"
                0039 #include "PARAMS.h"
2ef8966791 Davi*0040 #include "DIC_VARS.h"
daab022f42 Step*0041
                0042 C     == Routine arguments ==
f52bc56573 Jean*0043 C       diclocal = total inorganic carbon (mol/m^3)
daab022f42 Step*0044 C             where 1 T = 1 metric ton = 1000 kg
f52bc56573 Jean*0045 C       ta  = total alkalinity (eq/m^3)
                0046 C       pt  = inorganic phosphate (mol/^3)
                0047 C       sit = inorganic silicate (mol/^3)
daab022f42 Step*0048 C       t   = temperature (degrees C)
ba0b047096 Mart*0049 C       s   = salinity (g/kg)
daab022f42 Step*0050         INTEGER donewt
                0051         INTEGER inewtonmax
                0052         INTEGER ibrackmax
                0053         _RL  t, s, pt, sit, ta
                0054         _RL  pCO2surfloc, diclocal, pHlocal
                0055         _RL  fflocal, btlocal, stlocal, ftlocal
                0056         _RL  k1local, k2local
                0057         _RL  k1plocal, k2plocal, k3plocal
                0058         _RL  kslocal, kblocal, kwlocal, ksilocal, kflocal
d0092a57ac Step*0059         _RL  k0local, fugflocal
e098791e51 Jean*0060         INTEGER i,j,k,bi,bj,myIter
5625485478 Jean*0061         INTEGER myThid
f52bc56573 Jean*0062 CEOP
daab022f42 Step*0063
                0064 C     == Local variables ==
                0065 C INPUT
                0066 C       phlo= lower limit of pH range
                0067 C       phhi= upper limit of pH range
                0068 C       atmpres = atmospheric pressure in atmospheres (1 atm==1013.25mbar)
                0069 C OUTPUT
                0070 C       co2star  = CO2*water (mol/m^3)
                0071 C       pco2surf = oceanic pCO2 (ppmv)
                0072 C ---------------------------------------------------------------------
                0073 C OCMIP NOTE: Some words about units - (JCO, 4/4/1999)
                0074 C     - Models carry tracers in mol/m^3 (on a per volume basis)
f52bc56573 Jean*0075 C     - Conversely, this routine, which was written by
                0076 C       observationalists (C. Sabine and R. Key), passes input
daab022f42 Step*0077 C       arguments in umol/kg (i.e., on a per mass basis)
4e9f2133df Step*0078 C     - I have changed things slightly so that input arguments are in
daab022f42 Step*0079 C       mol/m^3,
f52bc56573 Jean*0080 C     - Thus, all input concentrations (diclocal, ta, pt, and st) should be
                0081 C       given in mol/m^3; output arguments "co2star" and "dco2star"
daab022f42 Step*0082 C       are likewise be in mol/m^3.
                0083 C ---------------------------------------------------------------------
f52bc56573 Jean*0084
daab022f42 Step*0085         _RL  phhi
                0086         _RL  phlo
                0087         _RL  c
                0088         _RL  a
                0089         _RL  a2
                0090         _RL  da
                0091         _RL  b
                0092         _RL  b2
                0093         _RL  db
                0094         _RL  fn
                0095         _RL  df
                0096         _RL  deltax
                0097         _RL  x
                0098         _RL  x2
                0099         _RL  x3
                0100         _RL  xmid
                0101         _RL  ftest
                0102         _RL  htotal
                0103         _RL  htotal2
f52bc56573 Jean*0104         _RL  co2star
daab022f42 Step*0105         _RL  phguess
d0092a57ac Step*0106         _RL  fco2
daab022f42 Step*0107         INTEGER inewton
                0108         INTEGER ibrack
                0109         INTEGER hstep
                0110         _RL  fni(3)
                0111         _RL  xlo
                0112         _RL  xhi
                0113         _RL  xguess
                0114         _RL  k123p
                0115         _RL  k12p
                0116         _RL  k12
                0117 c ---------------------------------------------------------------------
                0118 c import donewt flag
                0119 c set donewt = 1 for newton-raphson iteration
                0120 c set donewt = 0 for bracket and bisection
                0121 c ---------------------------------------------------------------------
                0122 C Change units from the input of mol/m^3 -> mol/kg:
                0123 c (1 mol/m^3)  x (1 m^3/1024.5 kg)
3daafce20b Jean*0124 c where the ocean mean surface density is 1024.5 kg/m^3
f52bc56573 Jean*0125 c Note: mol/kg are actually what the body of this routine uses
                0126 c for calculations.  Units are reconverted back to mol/m^3 at the
daab022f42 Step*0127 c end of this routine.
                0128 c ---------------------------------------------------------------------
f52bc56573 Jean*0129 c To convert input in mol/m^3 -> mol/kg
daab022f42 Step*0130         pt=pt*permil
                0131         sit=sit*permil
                0132         ta=ta*permil
                0133         diclocal=diclocal*permil
                0134 c ---------------------------------------------------------------------
                0135 c set first guess and brackets for [H+] solvers
                0136 c first guess (for newton-raphson)
                0137         phguess = phlocal
                0138
                0139 c bracketing values (for bracket/bisection)
                0140         phhi = 10.0
                0141         phlo = 5.0
                0142 c convert to [H+]...
29ad036528 Step*0143         xguess = 10.0**(-phguess)
                0144         xlo = 10.0**(-phhi)
                0145         xhi = 10.0**(-phlo)
daab022f42 Step*0146         xmid = (xlo + xhi)*0.5
                0147
                0148 c----------------------------------------------------------------
                0149 c iteratively solve for [H+]
f52bc56573 Jean*0150 c (i) Newton-Raphson method with fixed number of iterations,
daab022f42 Step*0151 c     use previous [H+] as first guess
                0152
                0153 c select newton-raphson, inewt=1
                0154 c else select bracket and bisection
                0155
                0156 cQQQQQ
                0157         if( donewt .eq. 1)then
                0158 c.........................................................
                0159 c NEWTON-RAPHSON METHOD
                0160 c.........................................................
                0161           x = xguess
                0162 cdiags
                0163 c         WRITE(0,*)'xguess ',xguess
                0164 cdiags
720e9330bd Step*0165           do inewton = 1, inewtonmax
daab022f42 Step*0166 c set some common combinations of parameters used in
                0167 c the iterative [H+] solvers
                0168             x2=x*x
                0169             x3=x2*x
                0170             k12 = k1local*k2local
                0171             k12p = k1plocal*k2plocal
                0172             k123p = k12p*k3plocal
                0173             c = 1.0 + stlocal/kslocal
                0174             a = x3 + k1plocal*x2 + k12p*x + k123p
                0175             a2=a*a
                0176             da = 3.0*x2 + 2.0*k1plocal*x + k12p
                0177             b = x2 + k1local*x + k12
                0178             b2=b*b
                0179             db = 2.0*x + k1local
                0180
3daafce20b Jean*0181 c Evaluate f([H+]) and f_prime([H+])
daab022f42 Step*0182 c fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree
                0183 c      +hso4+hf+h3po4-ta
                0184             fn = k1local*x*diclocal/b +
                0185      &        2.0*diclocal*k12/b +
                0186      &        btlocal/(1.0 + x/kblocal) +
                0187      &        kwlocal/x +
                0188      &        pt*k12p*x/a +
                0189      &        2.0*pt*k123p/a +
                0190      &        sit/(1.0 + x/ksilocal) -
                0191      &        x/c -
                0192      &        stlocal/(1.0 + kslocal/x/c) -
                0193      &        ftlocal/(1.0 + kflocal/x) -
                0194      &        pt*x3/a -
                0195      &        ta
                0196
                0197 c df = dfn/dx
                0198 cdiags
                0199 c      WRITE(0,*)'values',b2,kblocal,x2,a2,c,x
                0200 cdiags
                0201             df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -
                0202      &        2.0*diclocal*k12*db/b2 -
                0203      &        btlocal/kblocal/(1.0+x/kblocal)**2. -
                0204      &        kwlocal/x2 +
                0205      &        (pt*k12p*(a - x*da))/a2 -
                0206      &        2.0*pt*k123p*da/a2 -
                0207      &        sit/ksilocal/(1.0+x/ksilocal)**2. +
                0208      &        1.0/c +
29ad036528 Step*0209      &        stlocal*(1.0 + kslocal/x/c)**(-2.0)*(kslocal/c/x2) +
                0210      &        ftlocal*(1.0 + kflocal/x)**(-2.)*kflocal/x2 -
daab022f42 Step*0211      &        pt*x2*(3.0*a-x*da)/a2
                0212 c evaluate increment in [H+]
                0213             deltax = - fn/df
                0214 c update estimate of [H+]
                0215             x = x + deltax
                0216 cdiags
                0217 c write value of x to check convergence....
                0218 c           write(0,*)'inewton, x, deltax ',inewton, x, deltax
                0219 c           write(6,*)
                0220 cdiags
                0221
720e9330bd Step*0222           end do
daab022f42 Step*0223 c end of newton-raphson method
                0224 c....................................................
f52bc56573 Jean*0225         else
daab022f42 Step*0226 c....................................................
                0227 C BRACKET AND BISECTION METHOD
                0228 c....................................................
                0229 c (ii) If first step use Bracket and Bisection method
                0230 c      with fixed, large number of iterations
720e9330bd Step*0231           do ibrack = 1, ibrackmax
daab022f42 Step*0232             do hstep = 1,3
                0233               if(hstep .eq. 1)x = xhi
f52bc56573 Jean*0234               if(hstep .eq. 2)x = xlo
daab022f42 Step*0235               if(hstep .eq. 3)x = xmid
                0236 c set some common combinations of parameters used in
                0237 c the iterative [H+] solvers
                0238
                0239               x2=x*x
                0240               x3=x2*x
                0241               k12 = k1local*k2local
                0242               k12p = k1plocal*k2plocal
                0243               k123p = k12p*k3plocal
                0244               c = 1.0 + stlocal/kslocal
                0245               a = x3 + k1plocal*x2 + k12p*x + k123p
                0246               a2=a*a
                0247               da = 3.0*x2 + 2.0*k1plocal*x + k12p
                0248               b = x2 + k1local*x + k12
                0249               b2=b*b
                0250               db = 2.0*x + k1local
                0251 c evaluate f([H+]) for bracketing and mid-value cases
                0252               fn = k1local*x*diclocal/b +
                0253      &          2.0*diclocal*k12/b +
                0254      &          btlocal/(1.0 + x/kblocal) +
                0255      &          kwlocal/x +
                0256      &          pt*k12p*x/a +
                0257      &          2.0*pt*k123p/a +
                0258      &          sit/(1.0 + x/ksilocal) -
                0259      &          x/c -
                0260      &          stlocal/(1.0 + kslocal/x/c) -
                0261      &          ftlocal/(1.0 + kflocal/x) -
                0262      &          pt*x3/a -
                0263      &          ta
                0264               fni(hstep) = fn
                0265             end do
                0266 c now bracket solution within two of three
                0267             ftest = fni(1)/fni(3)
                0268             if(ftest .gt. 0.0)then
                0269               xhi = xmid
                0270             else
                0271               xlo = xmid
                0272             end if
                0273             xmid = (xlo + xhi)*0.5
                0274
                0275 cdiags
                0276 c write value of x to check convergence....
                0277 c           WRITE(0,*)'bracket-bisection iteration ',ibrack, xmid
                0278 cdiags
720e9330bd Step*0279           end do
daab022f42 Step*0280 c last iteration gives value
                0281           x = xmid
                0282 c end of bracket and bisection method
                0283 c....................................
                0284         end if
                0285 c iterative [H+] solver finished
                0286 c----------------------------------------------------------------
                0287
                0288 c now determine pCO2 etc...
                0289 c htotal = [H+], hydrogen ion conc
                0290         htotal = x
f52bc56573 Jean*0291 C Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2,
daab022f42 Step*0292 C ORNL/CDIAC-74, dickson and Goyet, eds. (Ch 2 p 10, Eq A.49)
                0293         htotal2=htotal*htotal
4e9f2133df Step*0294         co2star=diclocal*htotal2/(htotal2 + k1local*htotal
daab022f42 Step*0295      &            + k1local*k2local)
                0296         phlocal=-log10(htotal)
                0297
                0298 c ---------------------------------------------------------------
                0299 c Add two output arguments for storing pCO2surf
                0300 c Should we be using K0 or ff for the solubility here?
                0301 c ---------------------------------------------------------------
d0092a57ac Step*0302 #ifdef WATERVAP_BUG
daab022f42 Step*0303         pCO2surfloc = co2star / fflocal
d0092a57ac Step*0304 #else
                0305 c Corrected by Val Bennington (Nov 2010)
                0306         fco2 = co2star / k0local
                0307         pCO2surfloc = fco2/fugflocal
                0308 #endif
daab022f42 Step*0309
                0310 C ----------------------------------------------------------------
                0311 C Reconvert units back to original values for input arguments
                0312 C no longer necessary????
                0313 C ----------------------------------------------------------------
                0314 c       Reconvert from mol/kg -> mol/m^3
                0315         pt=pt/permil
                0316         sit=sit/permil
                0317         ta=ta/permil
                0318         diclocal=diclocal/permil
                0319
f52bc56573 Jean*0320         RETURN
                0321         END
daab022f42 Step*0322
f52bc56573 Jean*0323 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0324
                0325 CBOP
                0326 C !ROUTINE: CALC_PCO2_APPROX
                0327
                0328 C !INTERFACE: ==========================================================
daab022f42 Step*0329        SUBROUTINE CALC_PCO2_APPROX(
                0330      I                       t,s,diclocal,pt,sit,ta,
                0331      I                       k1local,k2local,
                0332      I                       k1plocal,k2plocal,k3plocal,
                0333      I                       kslocal,kblocal,kwlocal,
                0334      I                       ksilocal,kflocal,
d0092a57ac Step*0335      I                       k0local, fugflocal,
daab022f42 Step*0336      I                       fflocal,btlocal,stlocal,ftlocal,
e6a52874f7 Davi*0337      U                       pHlocal,pCO2surfloc,co3local,
e098791e51 Jean*0338      I                       i,j,k,bi,bj,myIter,myThid )
f52bc56573 Jean*0339
                0340 C !DESCRIPTION:
                0341 C     *==========================================================*
daab022f42 Step*0342 C     | SUBROUTINE CALC_PCO2_APPROX                              |
f52bc56573 Jean*0343 C     *==========================================================*
e6a52874f7 Davi*0344 C     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0345 C     C  New efficient pCO2 solver, Mick Follows         CC
                0346 C     C                             Taka Ito             CC
                0347 C     C                             Stephanie Dutkiewicz CC
                0348 C     C  20 April 2003                                   CC
                0349 C     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0350 C      Apr 2011: fix vapour bug (following Bennington)
                0351 C      Oct 2011: add CO3 extimation and pass out
f52bc56573 Jean*0352
                0353 C !USES: ===============================================================
daab022f42 Step*0354       IMPLICIT NONE
                0355
                0356 C     == GLobal variables ==
                0357 #include "SIZE.h"
                0358 #include "EEPARAMS.h"
                0359 #include "PARAMS.h"
2ef8966791 Davi*0360 #include "DIC_VARS.h"
daab022f42 Step*0361
                0362 C     == Routine arguments ==
                0363 C       diclocal = total inorganic carbon (mol/m^3)
                0364 C             where 1 T = 1 metric ton = 1000 kg
                0365 C       ta  = total alkalinity (eq/m^3)
                0366 C       pt  = inorganic phosphate (mol/^3)
                0367 C       sit = inorganic silicate (mol/^3)
                0368 C       t   = temperature (degrees C)
ba0b047096 Mart*0369 C       s   = salinity (g/kg)
daab022f42 Step*0370         _RL  t, s, pt, sit, ta
                0371         _RL  pCO2surfloc, diclocal, pHlocal
                0372         _RL  fflocal, btlocal, stlocal, ftlocal
                0373         _RL  k1local, k2local
                0374         _RL  k1plocal, k2plocal, k3plocal
                0375         _RL  kslocal, kblocal, kwlocal, ksilocal, kflocal
d0092a57ac Step*0376         _RL  k0local, fugflocal
e6a52874f7 Davi*0377         _RL  co3local
e098791e51 Jean*0378         INTEGER i,j,k,bi,bj,myIter
5625485478 Jean*0379         INTEGER myThid
f52bc56573 Jean*0380 CEOP
daab022f42 Step*0381
                0382 C     == Local variables ==
                0383         _RL  phguess
                0384         _RL  cag
                0385         _RL  bohg
                0386         _RL  hguess
                0387         _RL  stuff
                0388         _RL  gamm
                0389         _RL  hnew
                0390         _RL  co2s
f52bc56573 Jean*0391         _RL  h3po4g, h2po4g, hpo4g, po4g
daab022f42 Step*0392         _RL  siooh3g
d0092a57ac Step*0393         _RL  fco2
daab022f42 Step*0394
                0395 c ---------------------------------------------------------------------
                0396 C Change units from the input of mol/m^3 -> mol/kg:
                0397 c (1 mol/m^3)  x (1 m^3/1024.5 kg)
3daafce20b Jean*0398 c where the ocean mean surface density is 1024.5 kg/m^3
daab022f42 Step*0399 c Note: mol/kg are actually what the body of this routine uses
                0400 c for calculations.  Units are reconverted back to mol/m^3 at the
                0401 c end of this routine.
                0402 c To convert input in mol/m^3 -> mol/kg
                0403         pt=pt*permil
                0404         sit=sit*permil
                0405         ta=ta*permil
                0406         diclocal=diclocal*permil
                0407 c ---------------------------------------------------------------------
                0408 c set first guess and brackets for [H+] solvers
                0409 c first guess (for newton-raphson)
                0410         phguess = phlocal
                0411 cmick - new approx method
                0412 cmick - make estimate of htotal (hydrogen ion conc) using
                0413 cmick   appromate estimate of CA, carbonate alkalinity
29ad036528 Step*0414         hguess = 10.0**(-phguess)
daab022f42 Step*0415 cmick - first estimate borate contribution using guess for [H+]
                0416         bohg = btlocal*kblocal/(hguess+kblocal)
                0417
                0418 cmick - first estimate of contribution from phosphate
                0419 cmick based on Dickson and Goyet
                0420         stuff = hguess*hguess*hguess
                0421      &           + (k1plocal*hguess*hguess)
                0422      &           + (k1plocal*k2plocal*hguess)
                0423      &           + (k1plocal*k2plocal*k3plocal)
                0424         h3po4g = (pt*hguess*hguess*hguess) / stuff
                0425         h2po4g = (pt*k1plocal*hguess*hguess) / stuff
                0426         hpo4g  = (pt*k1plocal*k2plocal*hguess) / stuff
                0427         po4g   = (pt*k1plocal*k2plocal*k3plocal) / stuff
                0428
                0429 cmick - estimate contribution from silicate
                0430 cmick based on Dickson and Goyet
                0431         siooh3g = sit*ksilocal / (ksilocal + hguess)
                0432
                0433 cmick - now estimate carbonate alkalinity
                0434         cag = ta - bohg - (kwlocal/hguess) + hguess
75e97f1e14 Davi*0435      &           - hpo4g - 2.0 _d 0*po4g + h3po4g
daab022f42 Step*0436      &           - siooh3g
                0437
                0438 cmick - now evaluate better guess of hydrogen ion conc
                0439 cmick   htotal = [H+], hydrogen ion conc
                0440         gamm  = diclocal/cag
75e97f1e14 Davi*0441         stuff = (1.0 _d 0-gamm)*(1.0 _d 0-gamm)*k1local*k1local
                0442      &          - 4.0 _d 0*k1local*k2local*(1.0 _d 0-2.0 _d 0*gamm)
                0443         hnew  = 0.5 _d 0*( (gamm-1.0 _d 0)*k1local + sqrt(stuff) )
daab022f42 Step*0444 cmick - now determine [CO2*]
                0445         co2s  = diclocal/
75e97f1e14 Davi*0446      &   (1.0 _d 0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew)))
daab022f42 Step*0447 cmick - return update pH to main routine
                0448         phlocal = -log10(hnew)
                0449
4e9f2133df Step*0450 c NOW EVALUATE CO32-, carbonate ion concentration
                0451 c used in determination of calcite compensation depth
                0452 c Karsten Friis & Mick - Sep 2004
e6a52874f7 Davi*0453         co3local = k1local*k2local*diclocal /
                0454      &         (hnew*hnew + k1local*hnew + k1local*k2local)
4e9f2133df Step*0455
daab022f42 Step*0456 c ---------------------------------------------------------------
                0457 c surface pCO2 (following Dickson and Goyet, DOE...)
d0092a57ac Step*0458 #ifdef WATERVAP_BUG
daab022f42 Step*0459         pCO2surfloc = co2s/fflocal
d0092a57ac Step*0460 #else
                0461 c bug fix by Bennington
                0462         fco2 = co2s/k0local
                0463         pco2surfloc = fco2/fugflocal
                0464 #endif
daab022f42 Step*0465
                0466 C ----------------------------------------------------------------
                0467 c Reconvert from mol/kg -> mol/m^3
                0468         pt=pt/permil
                0469         sit=sit/permil
                0470         ta=ta/permil
                0471         diclocal=diclocal/permil
                0472
f52bc56573 Jean*0473         RETURN
                0474         END
                0475
                0476 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0477 CBOP
                0478 C !ROUTINE: CARBON_COEFFS
                0479
                0480 C !INTERFACE: ==========================================================
daab022f42 Step*0481       SUBROUTINE CARBON_COEFFS(
                0482      I                   ttemp,stemp,
5625485478 Jean*0483      I                   bi,bj,iMin,iMax,jMin,jMax,myThid)
f52bc56573 Jean*0484
                0485 C !DESCRIPTION:
                0486 C     *==========================================================*
daab022f42 Step*0487 C     | SUBROUTINE CARBON_COEFFS                                 |
                0488 C     | determine coefficients for surface carbon chemistry      |
                0489 C     | adapted from OCMIP2:  SUBROUTINE CO2CALC                 |
                0490 C     | mick follows, oct 1999                                   |
f52bc56573 Jean*0491 C     | minor changes to tidy, swd aug 2002                      |
                0492 C     *==========================================================*
daab022f42 Step*0493 C INPUT
f52bc56573 Jean*0494 C       diclocal = total inorganic carbon (mol/m^3)
daab022f42 Step*0495 C             where 1 T = 1 metric ton = 1000 kg
f52bc56573 Jean*0496 C       ta  = total alkalinity (eq/m^3)
                0497 C       pt  = inorganic phosphate (mol/^3)
                0498 C       sit = inorganic silicate (mol/^3)
daab022f42 Step*0499 C       t   = temperature (degrees C)
ba0b047096 Mart*0500 C       s   = salinity (g/kg)
daab022f42 Step*0501 C OUTPUT
                0502 C IMPORTANT: Some words about units - (JCO, 4/4/1999)
f52bc56573 Jean*0503 C     - Models carry tracers in mol/m^3 (on a per volume basis)
                0504 C     - Conversely, this routine, which was written by observationalists
                0505 C       (C. Sabine and R. Key), passes input arguments in umol/kg
                0506 C       (i.e., on a per mass basis)
                0507 C     - I have changed things slightly so that input arguments are in mol/m^3,
                0508 C     - Thus, all input concentrations (diclocal, ta, pt, and st) should be
                0509 C       given in mol/m^3; output arguments "co2star" and "dco2star"
                0510 C       are likewise be in mol/m^3.
d0092a57ac Step*0511 C
                0512 C Apr 2011: fix vapour bug (following Bennington)
daab022f42 Step*0513 C--------------------------------------------------------------------------
f52bc56573 Jean*0514
                0515 C !USES: ===============================================================
daab022f42 Step*0516         IMPLICIT NONE
                0517 C     == GLobal variables ==
                0518 #include "SIZE.h"
                0519 #include "EEPARAMS.h"
                0520 #include "PARAMS.h"
                0521 #include "GRID.h"
2ef8966791 Davi*0522 #include "DIC_VARS.h"
f0dfb321df Jean*0523
daab022f42 Step*0524 C     == Routine arguments ==
                0525 C ttemp and stemp are local theta and salt arrays
                0526 C dont really need to pass T and S in, could use theta, salt in
                0527 C common block in DYNVARS.h, but this way keeps subroutine more
                0528 C general
8bf2c0e0ad Step*0529         _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0530         _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
daab022f42 Step*0531         INTEGER bi,bj,iMin,iMax,jMin,jMax
5625485478 Jean*0532         INTEGER myThid
f52bc56573 Jean*0533 CEOP
4e9f2133df Step*0534
                0535 C LOCAL VARIABLES
daab022f42 Step*0536         _RL  t
                0537         _RL  s
                0538         _RL  tk
                0539         _RL  tk100
                0540         _RL  tk1002
                0541         _RL  dlogtk
                0542         _RL  sqrtis
                0543         _RL  sqrts
                0544         _RL  s15
                0545         _RL  scl
                0546         _RL  s2
                0547         _RL  invtk
                0548         _RL  is
                0549         _RL  is2
d0092a57ac Step*0550 c add Bennington
                0551         _RL  P1atm
                0552         _RL  Rgas
                0553         _RL  RT
                0554         _RL  delta
                0555         _RL  B1
                0556         _RL  B
425561518f Jona*0557 #ifdef CARBONCHEM_TOTALPHSCALE
                0558 c pH scale converstions
                0559         _RL total2free
8d230ec051 Jona*0560         _RL free2total
92031a2908 Jean*0561         _RL free2sw
8d230ec051 Jona*0562         _RL sw2free
425561518f Jona*0563         _RL total2sw
8d230ec051 Jona*0564         _RL sw2total
425561518f Jona*0565 #endif
daab022f42 Step*0566         INTEGER i
                0567         INTEGER j
                0568
                0569 C.....................................................................
                0570 C OCMIP note:
                0571 C Calculate all constants needed to convert between various measured
4e9f2133df Step*0572 C carbon species. References for each equation are noted in the code.
daab022f42 Step*0573 C Once calculated, the constants are
f52bc56573 Jean*0574 C stored and passed in the common block "const". The original version
                0575 C of this code was based on the code by dickson in Version 2 of
                0576 C  Handbook of Methods C for the Analysis of the Various Parameters of
                0577 C the Carbon Dioxide System in Seawater , DOE, 1994 (SOP No. 3, p25-26).
daab022f42 Step*0578 C....................................................................
                0579
51c3bf0077 Step*0580         do i=imin,imax
                0581          do j=jmin,jmax
f0dfb321df Jean*0582           IF ( maskC(i,j,1,bi,bj).EQ.oneRS ) THEN
8bf2c0e0ad Step*0583            t = ttemp(i,j)
                0584            s = stemp(i,j)
425561518f Jona*0585 C terms used more than once for:
                0586 C temperature
75e97f1e14 Davi*0587            tk = 273.15 _d 0 + t
                0588            tk100 = tk/100. _d 0
daab022f42 Step*0589            tk1002=tk100*tk100
75e97f1e14 Davi*0590            invtk=1.0 _d 0/tk
daab022f42 Step*0591            dlogtk=log(tk)
425561518f Jona*0592 C ionic strength
75e97f1e14 Davi*0593            is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)
daab022f42 Step*0594            is2=is*is
                0595            sqrtis=sqrt(is)
425561518f Jona*0596 C salinity
daab022f42 Step*0597            s2=s*s
                0598            sqrts=sqrt(s)
75e97f1e14 Davi*0599            s15=s**1.5 _d 0
                0600            scl=s/1.80655 _d 0
d0092a57ac Step*0601 C -----------------------------------------------------------------------
                0602 C added by Val Bennington Nov 2010
                0603 C Fugacity Factor needed for non-ideality in ocean
                0604 C ff used for atmospheric correction for water vapor and pressure
                0605 C Weiss (1974) Marine Chemistry
                0606            P1atm = 1.01325 _d 0 ! bars
                0607            Rgas = 83.1451 _d 0 ! bar*cm3/(mol*K)
                0608            RT = Rgas*tk
                0609            delta = (57.7 _d 0 - 0.118 _d 0*tk)
                0610            B1 = -1636.75 _d 0 + 12.0408 _d 0*tk - 0.0327957 _d 0*tk*tk
                0611            B = B1 + 3.16528 _d 0*tk*tk*tk*(0.00001 _d 0)
7122734d29 Davi*0612            fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * P1atm / RT)
daab022f42 Step*0613 C------------------------------------------------------------------------
                0614 C f = k0(1-pH2O)*correction term for non-ideality
                0615 C Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
75e97f1e14 Davi*0616            ff(i,j,bi,bj) = exp(-162.8301 _d 0 + 218.2968 _d 0/tk100  +
                0617      &          90.9241 _d 0*log(tk100) - 1.47696 _d 0*tk1002 +
                0618      &          s * (.025695 _d 0 - .025225 _d 0*tk100 +
                0619      &          0.0049867 _d 0*tk1002))
daab022f42 Step*0620 C------------------------------------------------------------------------
                0621 C K0 from Weiss 1974
425561518f Jona*0622 C Independent of pH scale
75e97f1e14 Davi*0623            ak0(i,j,bi,bj) = exp(93.4517 _d 0/tk100 - 60.2409 _d 0 +
                0624      &        23.3585 _d 0 * log(tk100) +
                0625      &        s * (0.023517 _d 0 - 0.023656 _d 0*tk100 +
                0626      &        0.0047036 _d 0*tk1002))
daab022f42 Step*0627 C------------------------------------------------------------------------
                0628 C k1 = [H][HCO3]/[H2CO3]
                0629 C k2 = [H][CO3]/[HCO3]
92031a2908 Jean*0630 C Millero p.664 (1995) using Mehrbach et al. data
425561518f Jona*0631 C Both on seawater pH scale
75e97f1e14 Davi*0632            ak1(i,j,bi,bj)=10.**(-1. _d 0*(3670.7 _d 0*invtk -
                0633      &          62.008 _d 0 + 9.7944 _d 0*dlogtk -
                0634      &          0.0118 _d 0 * s + 0.000116 _d 0*s2))
                0635            ak2(i,j,bi,bj)=10.**(-1. _d 0*(1394.7 _d 0*invtk+ 4.777 _d 0-
                0636      &          0.0184 _d 0*s + 0.000118 _d 0*s2))
daab022f42 Step*0637 C------------------------------------------------------------------------
                0638 C kb = [H][BO2]/[HBO2]
                0639 C Millero p.669 (1995) using data from dickson (1990)
425561518f Jona*0640 C On total pH scale
75e97f1e14 Davi*0641            akb(i,j,bi,bj)=exp((-8966.90 _d 0- 2890.53 _d 0*sqrts -
                0642      &          77.942 _d 0*s + 1.728 _d 0*s15 - 0.0996 _d 0*s2)*invtk +
                0643      &          (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) +
                0644      &          (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) *
                0645      &          dlogtk + 0.053105 _d 0*sqrts*tk)
daab022f42 Step*0646 C------------------------------------------------------------------------
                0647 C k1p = [H][H2PO4]/[H3PO4]
                0648 C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
92031a2908 Jean*0649 C The constant 115.525 is an approximation to convert to the total pH scale
425561518f Jona*0650 C  (Dickson et al., 2007). Use Millero's 115.540 constant to stay on the seawater
                0651 C  pH scale (everything else is the same).
75e97f1e14 Davi*0652            ak1p(i,j,bi,bj) = exp(-4576.752 _d 0*invtk + 115.525 _d 0 -
                0653      &          18.453 _d 0*dlogtk +
                0654      &          (-106.736 _d 0*invtk + 0.69171 _d 0)*sqrts +
                0655      &          (-0.65643 _d 0*invtk - 0.01844 _d 0)*s)
daab022f42 Step*0656 C------------------------------------------------------------------------
                0657 C k2p = [H][HPO4]/[H2PO4]
                0658 C DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
92031a2908 Jean*0659 C The constant 172.0833 is an approximation to convert to the total pH scale
425561518f Jona*0660 C  (Dickson et al., 2007). Use Millero's 172.1033 constant to stay on the seawater
                0661 C  pH scale (everything else is the same).
75e97f1e14 Davi*0662            ak2p(i,j,bi,bj) = exp(-8814.715 _d 0*invtk + 172.0883 _d 0 -
                0663      &          27.927 _d 0*dlogtk +
                0664      &          (-160.340 _d 0*invtk + 1.3566 _d 0) * sqrts +
                0665      &          (0.37335 _d 0*invtk - 0.05778 _d 0) * s)
daab022f42 Step*0666 C------------------------------------------------------------------------
                0667 C k3p = [H][PO4]/[HPO4]
                0668 C DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
92031a2908 Jean*0669 C The constant 18.141 is an approximation to convert to the total pH scale
425561518f Jona*0670 C  (Dickson et al., 2007). Use Millero's 18.126 constant to stay on the seawater
                0671 C  pH scale (everything else is the same).
75e97f1e14 Davi*0672            ak3p(i,j,bi,bj) = exp(-3070.75 _d 0*invtk - 18.141 _d 0 +
                0673      &          (17.27039 _d 0*invtk + 2.81197 _d 0) *
                0674      &          sqrts + (-44.99486 _d 0*invtk - 0.09984 _d 0) * s)
daab022f42 Step*0675 C------------------------------------------------------------------------
                0676 C ksi = [H][SiO(OH)3]/[Si(OH)4]
                0677 C Millero p.671 (1995) using data from Yao and Millero (1995)
92031a2908 Jean*0678 C The constant 117.385 is an approximation to convert to the total pH scale
425561518f Jona*0679 C  (Dickson et al., 2007). Use Millero's 117.400 constant to stay on the seawater
                0680 C  pH scale (everything else is the same).
80d7ca01bb Jona*0681 C  Note: converted here from mol/kg-H2O to mol/kg-SW
75e97f1e14 Davi*0682            aksi(i,j,bi,bj) = exp(-8904.2 _d 0*invtk + 117.385 _d 0 -
                0683      &          19.334 _d 0*dlogtk +
                0684      &          (-458.79 _d 0*invtk + 3.5913 _d 0) * sqrtis +
                0685      &          (188.74 _d 0*invtk - 1.5998 _d 0) * is +
                0686      &          (-12.1652 _d 0*invtk + 0.07871 _d 0) * is2 +
                0687      &          log(1.0 _d 0-0.001005 _d 0*s))
daab022f42 Step*0688 C------------------------------------------------------------------------
                0689 C kw = [H][OH]
                0690 C Millero p.670 (1995) using composite data
92031a2908 Jean*0691 C The constant 148.9652 is an approximation to convert to the total pH scale
425561518f Jona*0692 C  (Dickson et al., 2007). Use Millero's 148.9802 constant to stay on the seawater
                0693 C  pH scale (everything else is the same).
75e97f1e14 Davi*0694            akw(i,j,bi,bj) = exp(-13847.26 _d 0*invtk + 148.9652 _d 0 -
                0695      &          23.6521 _d 0*dlogtk +
                0696      &          (118.67 _d 0*invtk - 5.977 _d 0 + 1.0495 _d 0 * dlogtk)
                0697      &          * sqrts - 0.01615 _d 0 * s)
daab022f42 Step*0698 C------------------------------------------------------------------------
                0699 C ks = [H][SO4]/[HSO4]
                0700 C dickson (1990, J. chem. Thermodynamics 22, 113)
425561518f Jona*0701 C On the free pH scale
80d7ca01bb Jona*0702 C Note: converted here from mol/kg-H2O to mol/kg-SW
75e97f1e14 Davi*0703            aks(i,j,bi,bj)=exp(-4276.1 _d 0*invtk + 141.328 _d 0 -
                0704      &          23.093 _d 0*dlogtk +
                0705      &   (-13856. _d 0*invtk + 324.57 _d 0 - 47.986 _d 0*dlogtk)*sqrtis+
                0706      &   (35474. _d 0*invtk - 771.54 _d 0 + 114.723 _d 0*dlogtk)*is -
e18333c42b Davi*0707      &          2698. _d 0*invtk*is**1.5 _d 0 + 1776. _d 0*invtk*is2 +
75e97f1e14 Davi*0708      &          log(1.0 _d 0 - 0.001005 _d 0*s))
daab022f42 Step*0709 C------------------------------------------------------------------------
                0710 C kf = [H][F]/[HF]
425561518f Jona*0711 C dickson and Riley (1979) -- change pH scale to total (following Dickson & Goyet, 1994)
80d7ca01bb Jona*0712 C Note: converted here from mol/kg-H2O to mol/kg-SW
75e97f1e14 Davi*0713            akf(i,j,bi,bj)=exp(1590.2 _d 0*invtk - 12.641 _d 0 +
f52bc56573 Jean*0714      &   1.525 _d 0*sqrtis + log(1.0 _d 0 - 0.001005 _d 0*s) +
75e97f1e14 Davi*0715      &   log(1.0 _d 0 + (0.1400 _d 0/96.062 _d 0)*(scl)/aks(i,j,bi,bj)))
daab022f42 Step*0716 C------------------------------------------------------------------------
                0717 C Calculate concentrations for borate, sulfate, and fluoride
                0718 C Uppstrom (1974)
75e97f1e14 Davi*0719            bt(i,j,bi,bj) = 0.000232 _d 0 * scl/10.811 _d 0
daab022f42 Step*0720 C Morris & Riley (1966)
75e97f1e14 Davi*0721            st(i,j,bi,bj) = 0.14 _d 0 * scl/96.062 _d 0
daab022f42 Step*0722 C Riley (1965)
75e97f1e14 Davi*0723            ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0
daab022f42 Step*0724 C------------------------------------------------------------------------
425561518f Jona*0725 #ifdef CARBONCHEM_TOTALPHSCALE
                0726 C Convert between pH scales, we want everything to end up on the total scale
                0727 C ak1 and ak2 are on seawater scale and aks is on the free scale
                0728            total2free = 1. _d 0/
8d230ec051 Jona*0729      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
425561518f Jona*0730
92031a2908 Jean*0731            free2total = (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                0732
                0733            free2sw = 1. _d 0
                0734      &                 + st(i,j,bi,bj)/ aks(i,j,bi,bj)
                0735      &                 + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free)
                0736
8d230ec051 Jona*0737            sw2free = 1. _d 0 / free2sw
92031a2908 Jean*0738
                0739            total2sw = total2free * free2sw
                0740
8d230ec051 Jona*0741            sw2total = 1. _d 0 / total2sw
425561518f Jona*0742
8d230ec051 Jona*0743            ak1(i,j,bi,bj)  = ak1(i,j,bi,bj)*sw2total
                0744            ak2(i,j,bi,bj)  = ak2(i,j,bi,bj)*sw2total
92031a2908 Jean*0745            aks(i,j,bi,bj)  = aks(i,j,bi,bj)*free2total
425561518f Jona*0746 #endif
f0dfb321df Jean*0747           ELSE
d0092a57ac Step*0748 c add Bennington
c5830e1a8b Jona*0749            fugf(i,j,bi,bj)=0. _d 0
                0750            ff(i,j,bi,bj)=0. _d 0
                0751            ak0(i,j,bi,bj)= 0. _d 0
                0752            ak1(i,j,bi,bj)= 0. _d 0
                0753            ak2(i,j,bi,bj)= 0. _d 0
                0754            akb(i,j,bi,bj)= 0. _d 0
                0755            ak1p(i,j,bi,bj) = 0. _d 0
                0756            ak2p(i,j,bi,bj) = 0. _d 0
                0757            ak3p(i,j,bi,bj) = 0. _d 0
                0758            aksi(i,j,bi,bj) = 0. _d 0
                0759            akw(i,j,bi,bj) = 0. _d 0
                0760            aks(i,j,bi,bj)= 0. _d 0
                0761            akf(i,j,bi,bj)= 0. _d 0
                0762            bt(i,j,bi,bj) = 0. _d 0
                0763            st(i,j,bi,bj) = 0. _d 0
                0764            ft(i,j,bi,bj) = 0. _d 0
f0dfb321df Jean*0765           ENDIF
                0766          enddo
                0767         enddo
daab022f42 Step*0768
f52bc56573 Jean*0769         RETURN
                0770         END
                0771
                0772 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
daab022f42 Step*0773
f52bc56573 Jean*0774 CBOP
                0775 C !ROUTINE: CARBON_COEFFS_PRESSURE_DEP
                0776
                0777 C !INTERFACE: ==========================================================
4e9f2133df Step*0778       SUBROUTINE CARBON_COEFFS_PRESSURE_DEP(
                0779      I                   ttemp,stemp,
                0780      I                   bi,bj,iMin,iMax,jMin,jMax,
5625485478 Jean*0781      I                   Klevel,myThid)
f52bc56573 Jean*0782
                0783 C !DESCRIPTION:
                0784 C     *==========================================================*
4e9f2133df Step*0785 C     | SUBROUTINE CARBON_COEFFS                                 |
                0786 C     | determine coefficients for surface carbon chemistry      |
                0787 C     | adapted from OCMIP2:  SUBROUTINE CO2CALC                 |
                0788 C     | mick follows, oct 1999                                   |
f52bc56573 Jean*0789 C     | minor changes to tidy, swd aug 2002                      |
                0790 C     | MODIFIED FOR PRESSURE DEPENDENCE                         |
                0791 C     | Karsten Friis and Mick Follows 2004                      |
                0792 C     *==========================================================*
4e9f2133df Step*0793 C INPUT
                0794 C       diclocal = total inorganic carbon (mol/m^3)
                0795 C             where 1 T = 1 metric ton = 1000 kg
                0796 C       ta  = total alkalinity (eq/m^3)
                0797 C       pt  = inorganic phosphate (mol/^3)
                0798 C       sit = inorganic silicate (mol/^3)
                0799 C       t   = temperature (degrees C)
ba0b047096 Mart*0800 C       s   = salinity (g/kg)
4e9f2133df Step*0801 C OUTPUT
                0802 C IMPORTANT: Some words about units - (JCO, 4/4/1999)
f52bc56573 Jean*0803 C     - Models carry tracers in mol/m^3 (on a per volume basis)
                0804 C     - Conversely, this routine, which was written by observationalists
                0805 C       (C. Sabine and R. Key), passes input arguments in umol/kg
                0806 C       (i.e., on a per mass basis)
                0807 C     - I have changed things slightly so that input arguments are in mol/m^3,
                0808 C     - Thus, all input concentrations (diclocal, ta, pt, and st) should be
                0809 C       given in mol/m^3; output arguments "co2star" and "dco2star"
                0810 C       are likewise be in mol/m^3.
                0811 C
                0812 C NOW INCLUDES:
                0813 C PRESSURE DEPENDENCE of K1, K2, solubility product of calcite
                0814 C based on Takahashi, GEOSECS Atlantic Report, Vol. 1 (1981)
                0815 C
4e9f2133df Step*0816 C--------------------------------------------------------------------------
f52bc56573 Jean*0817
                0818 C !USES: ===============================================================
4e9f2133df Step*0819         IMPLICIT NONE
                0820 C     == GLobal variables ==
                0821 #include "SIZE.h"
                0822 #include "EEPARAMS.h"
                0823 #include "PARAMS.h"
                0824 #include "GRID.h"
2ef8966791 Davi*0825 #include "DIC_VARS.h"
f0dfb321df Jean*0826
4e9f2133df Step*0827 C     == Routine arguments ==
                0828 C ttemp and stemp are local theta and salt arrays
                0829 C dont really need to pass T and S in, could use theta, salt in
                0830 C common block in DYNVARS.h, but this way keeps subroutine more
                0831 C general
8bf2c0e0ad Step*0832         _RL  ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0833         _RL  stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4e9f2133df Step*0834         INTEGER bi,bj,iMin,iMax,jMin,jMax
f52bc56573 Jean*0835 C K is depth index
4e9f2133df Step*0836         INTEGER Klevel
5625485478 Jean*0837         INTEGER myThid
f52bc56573 Jean*0838 CEOP
4e9f2133df Step*0839
                0840 C LOCAL VARIABLES
                0841         _RL  t
                0842         _RL  s
                0843         _RL  tk
                0844         _RL  tk100
                0845         _RL  tk1002
                0846         _RL  dlogtk
                0847         _RL  sqrtis
                0848         _RL  sqrts
                0849         _RL  s15
                0850         _RL  scl
                0851         _RL  s2
                0852         _RL  invtk
                0853         _RL  is
                0854         _RL  is2
                0855         INTEGER i
                0856         INTEGER j
                0857         INTEGER k
                0858         _RL  bdepth
                0859         _RL  cdepth
                0860         _RL  pressc
                0861         _RL  Ksp_T_Calc
                0862         _RL  xvalue
                0863         _RL  zdum
                0864         _RL  tmpa1
                0865         _RL  tmpa2
                0866         _RL  tmpa3
                0867         _RL  logKspc
                0868         _RL  dv
                0869         _RL  dk
                0870         _RL  pfactor
                0871         _RL  bigR
425561518f Jona*0872 #ifdef CARBONCHEM_TOTALPHSCALE
                0873 c pH scale converstions
                0874         _RL total2free_surf
92031a2908 Jean*0875         _RL free2sw_surf
425561518f Jona*0876         _RL total2sw_surf
                0877         _RL total2free
8d230ec051 Jona*0878         _RL free2total
92031a2908 Jean*0879         _RL free2sw
8d230ec051 Jona*0880         _RL sw2free
425561518f Jona*0881         _RL total2sw
8d230ec051 Jona*0882         _RL sw2total
425561518f Jona*0883 #endif
4e9f2133df Step*0884 C.....................................................................
                0885 C OCMIP note:
                0886 C Calculate all constants needed to convert between various measured
                0887 C carbon species. References for each equation are noted in the code.
                0888 C Once calculated, the constants are
                0889 C stored and passed in the common block "const". The original version
                0890 C of this code was based on the code by dickson in Version 2 of
3daafce20b Jean*0891 C  Handbook of Methods C for the Analysis of the Various Parameters of
                0892 C the Carbon Dioxide System in Seawater , DOE, 1994 (SOP No. 3, p25-26).
4e9f2133df Step*0893 C....................................................................
                0894
                0895 c determine pressure (bar) from depth
                0896 c 1 BAR at z=0m (atmos pressure)
                0897 c use UPPER surface of cell so top layer pressure = 0 bar
                0898 c for surface exchange coeffs
                0899
                0900 cmick..............................
                0901 c        write(6,*)'Klevel ',klevel
                0902
c5830e1a8b Jona*0903         bdepth = 0. _d 0
                0904         cdepth = 0. _d 0
                0905         pressc = 1. _d 0
4e9f2133df Step*0906         do k = 1,Klevel
c5830e1a8b Jona*0907             cdepth = bdepth + 0.5 _d 0*drF(k)
4e9f2133df Step*0908             bdepth = bdepth + drF(k)
c5830e1a8b Jona*0909             pressc = 1. _d 0 + 0.1 _d 0*cdepth
f0dfb321df Jean*0910         enddo
4e9f2133df Step*0911 cmick...................................................
                0912 c        write(6,*)'depth,pressc ',cdepth,pressc
                0913 cmick....................................................
                0914
51c3bf0077 Step*0915         do i=imin,imax
                0916          do j=jmin,jmax
f0dfb321df Jean*0917           IF ( maskC(i,j,Klevel,bi,bj).EQ.oneRS ) THEN
8bf2c0e0ad Step*0918            t = ttemp(i,j)
                0919            s = stemp(i,j)
c5830e1a8b Jona*0920 C terms used more than once for:
                0921 C temperature
                0922            tk = 273.15 _d 0 + t
                0923            tk100 = tk/100. _d 0
4e9f2133df Step*0924            tk1002=tk100*tk100
c5830e1a8b Jona*0925            invtk=1.0 _d 0/tk
4e9f2133df Step*0926            dlogtk=log(tk)
c5830e1a8b Jona*0927 C ionic strength
                0928            is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)
4e9f2133df Step*0929            is2=is*is
                0930            sqrtis=sqrt(is)
c5830e1a8b Jona*0931 C salinity
4e9f2133df Step*0932            s2=s*s
                0933            sqrts=sqrt(s)
c5830e1a8b Jona*0934            s15=s**1.5 _d 0
                0935            scl=s/1.80655 _d 0
4e9f2133df Step*0936
c5830e1a8b Jona*0937            bigR = 83.14472 _d 0
4e9f2133df Step*0938 C------------------------------------------------------------------------
                0939 C f = k0(1-pH2O)*correction term for non-ideality
c5830e1a8b Jona*0940 C Weiss & Price (1980, Mar. Chem., 8, 347-359;  Eq 13 with table 6 values)
                0941            ff(i,j,bi,bj) = exp(-162.8301 _d 0 + 218.2968 _d 0/tk100  +
                0942      &          90.9241 _d 0*log(tk100) - 1.47696 _d 0*tk1002 +
                0943      &          s * (.025695 _d 0 - .025225 _d 0*tk100 +
                0944      &          0.0049867 _d 0*tk1002))
4e9f2133df Step*0945 C------------------------------------------------------------------------
                0946 C K0 from Weiss 1974
425561518f Jona*0947 C Independent of pH scale
c5830e1a8b Jona*0948            ak0(i,j,bi,bj) = exp(93.4517 _d 0/tk100 - 60.2409 _d 0 +
                0949      &        23.3585 _d 0 * log(tk100) +
                0950      &        s * (0.023517 _d 0 - 0.023656 _d 0*tk100 +
                0951      &        0.0047036 _d 0*tk1002))
4e9f2133df Step*0952 C------------------------------------------------------------------------
                0953 C k1 = [H][HCO3]/[H2CO3]
                0954 C k2 = [H][CO3]/[HCO3]
92031a2908 Jean*0955 C Millero p.664 (1995) using Mehrbach et al. data
425561518f Jona*0956 C Both on seawater pH scale
c5830e1a8b Jona*0957            ak1(i,j,bi,bj)=10**(-1 _d 0*(3670.7 _d 0*invtk -
                0958      &          62.008 _d 0 + 9.7944 _d 0*dlogtk -
                0959      &          0.0118 _d 0 * s + 0.000116 _d 0*s2))
                0960            ak2(i,j,bi,bj)=10**(-1*(1394.7 _d 0*invtk + 4.777 _d 0 -
                0961      &          0.0184 _d 0*s + 0.000118 _d 0*s2))
4e9f2133df Step*0962 C NOW PRESSURE DEPENDENCE:
f52bc56573 Jean*0963 c Following Takahashi (1981) GEOSECS report - quoting Culberson and
4e9f2133df Step*0964 c Pytkowicz (1968)
                0965 c pressc = pressure in bars
                0966            ak1(i,j,bi,bj) = ak1(i,j,bi,bj)*
b5953f02df Jona*0967      &         exp( (24.2 _d 0-0.085 _d 0*t)
                0968      &             *(pressc-1.0 _d 0)/(83.143 _d 0*tk) )
4e9f2133df Step*0969 c FIRST GO FOR K2: According to GEOSECS (1982) report
                0970 c          ak2(i,j,bi,bj) = ak2(i,j,bi,bj)*
                0971 c    &             exp( (26.4-0.040*t)*(pressc-1.0)/(83.143*tk) )
                0972 c SECOND GO FOR K2: corrected coeff according to CO2sys documentation
                0973 c                   E. Lewis and D. Wallace (1998) ORNL/CDIAC-105
                0974            ak2(i,j,bi,bj) = ak2(i,j,bi,bj)*
b5953f02df Jona*0975      &         exp( (16.4 _d 0-0.040 _d 0*t)
                0976      &             *(pressc-1.0 _d 0)/(83.143 _d 0*tk) )
4e9f2133df Step*0977 C------------------------------------------------------------------------
                0978 C kb = [H][BO2]/[HBO2]
                0979 C Millero p.669 (1995) using data from dickson (1990)
425561518f Jona*0980 C On total pH scale
c5830e1a8b Jona*0981            akb(i,j,bi,bj)=exp((-8966.90 _d 0- 2890.53 _d 0*sqrts -
                0982      &          77.942 _d 0*s + 1.728 _d 0*s15 - 0.0996 _d 0*s2)*invtk +
                0983      &          (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) +
                0984      &          (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) *
                0985      &          dlogtk + 0.053105 _d 0*sqrts*tk)
425561518f Jona*0986 #ifdef CARBONCHEM_TOTALPHSCALE
                0987 C JML Not yet, need to account for pH scale conversions
                0988 C    (pressure correct on the seawater scale, not the total scale)
                0989 C Mick and Karsten - Dec 04
                0990 C ADDING pressure dependence based on Millero (1995), p675
                0991 C with additional info from CO2sys documentation (E. Lewis and
                0992 C D. Wallace, 1998 - see endnotes for commentary on Millero, 95)
                0993 C           bigR = 83.145
                0994 C           dv = -29.48 + 0.1622*t + 2.608d-3*t*t
                0995 C           dk = -2.84d-3
                0996 C           pfactor = - (dv/(bigR*tk))*pressc
                0997 C     &               + (0.5*dk/(bigR*tk))*pressc*pressc
                0998 C           akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor)
                0999 #else
4e9f2133df Step*1000 C Mick and Karsten - Dec 04
                1001 C ADDING pressure dependence based on Millero (1995), p675
f52bc56573 Jean*1002 C with additional info from CO2sys documentation (E. Lewis and
4e9f2133df Step*1003 C D. Wallace, 1998 - see endnotes for commentary on Millero, 95)
c5830e1a8b Jona*1004            bigR = 83.145 _d 0
                1005            dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t
                1006            dk = -2.84 _d -3
f52bc56573 Jean*1007            pfactor = - (dv/(bigR*tk))*pressc
c5830e1a8b Jona*1008      &               + (0.5 _d 0*dk/(bigR*tk))*pressc*pressc
4e9f2133df Step*1009            akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor)
425561518f Jona*1010 #endif
4e9f2133df Step*1011 C------------------------------------------------------------------------
                1012 C k1p = [H][H2PO4]/[H3PO4]
                1013 C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
92031a2908 Jean*1014 C The constant 115.525 is an approximation to convert to the total pH scale
425561518f Jona*1015 C  (Dickson et al., 2007). Use Millero's 115.540 constant to stay on the seawater
                1016 C  pH scale (everything else is the same).
c5830e1a8b Jona*1017            ak1p(i,j,bi,bj) = exp(-4576.752 _d 0*invtk + 115.525 _d 0 -
                1018      &          18.453 _d 0*dlogtk +
                1019      &          (-106.736 _d 0*invtk + 0.69171 _d 0)*sqrts +
                1020      &          (-0.65643 _d 0*invtk - 0.01844 _d 0)*s)
4e9f2133df Step*1021 C------------------------------------------------------------------------
                1022 C k2p = [H][HPO4]/[H2PO4]
425561518f Jona*1023 C DOE(1994) eq 7.2.23 with footnote using data from Millero (1974)
92031a2908 Jean*1024 C The constant 172.0833 is an approximation to convert to the total pH scale
425561518f Jona*1025 C  (Dickson et al., 2007). Use Millero's 172.1033 constant to stay on the seawater
                1026 C  pH scale (everything else is the same).
c5830e1a8b Jona*1027            ak2p(i,j,bi,bj) = exp(-8814.715 _d 0*invtk + 172.0883 _d 0 -
                1028      &          27.927 _d 0*dlogtk +
                1029      &          (-160.34 _d 00*invtk + 1.3566 _d 0) * sqrts +
                1030      &          (0.37335 _d 0*invtk - 0.05778 _d 0) * s)
4e9f2133df Step*1031 C------------------------------------------------------------------------
                1032 C k3p = [H][PO4]/[HPO4]
                1033 C DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
92031a2908 Jean*1034 C The constant 18.141 is an approximation to convert to the total pH scale
425561518f Jona*1035 C  (Dickson et al., 2007). Use Millero's 18.126 constant to stay on the seawater
                1036 C  pH scale (everything else is the same).
c5830e1a8b Jona*1037            ak3p(i,j,bi,bj) = exp(-3070.75 _d 0*invtk - 18.141 _d 0 +
                1038      &          (17.27039 _d 0*invtk + 2.81197 _d 0) *
                1039      &          sqrts + (-44.99486 _d 0*invtk - 0.09984 _d 0) * s)
4e9f2133df Step*1040 C------------------------------------------------------------------------
                1041 C ksi = [H][SiO(OH)3]/[Si(OH)4]
                1042 C Millero p.671 (1995) using data from Yao and Millero (1995)
92031a2908 Jean*1043 C The constant 117.385 is an approximation to convert to the total pH scale
425561518f Jona*1044 C  (Dickson et al., 2007). Use Millero's 117.400 constant to stay on the seawater
                1045 C  pH scale (everything else is the same).
80d7ca01bb Jona*1046 C Note: converted here from mol/kg-H2O to mol/kg-SW
c5830e1a8b Jona*1047            aksi(i,j,bi,bj) = exp(-8904.2 _d 0*invtk + 117.385 _d 0 -
                1048      &          19.334 _d 0*dlogtk +
                1049      &          (-458.79 _d 0*invtk + 3.5913 _d 0) * sqrtis +
                1050      &          (188.74 _d 0*invtk - 1.5998 _d 0) * is +
                1051      &          (-12.1652 _d 0*invtk + 0.07871 _d 0) * is2 +
                1052      &          log(1.0 _d 0-0.001005 _d 0*s))
4e9f2133df Step*1053 C------------------------------------------------------------------------
                1054 C kw = [H][OH]
                1055 C Millero p.670 (1995) using composite data
92031a2908 Jean*1056 C The constant 148.9652 is an approximation to convert to the total pH scale
425561518f Jona*1057 C  (Dickson et al., 2007). Use Millero's 148.9802 constant to stay on the seawater
                1058 C  pH scale (everything else is the same).
c5830e1a8b Jona*1059            akw(i,j,bi,bj) = exp(-13847.26 _d 0*invtk + 148.9652 _d 0 -
b5953f02df Jona*1060      &          23.6521 _d 0*dlogtk + (118.67 _d 0*invtk
                1061      &          - 5.977 _d 0 + 1.0495 _d 0 * dlogtk) *
c5830e1a8b Jona*1062      &          sqrts - 0.01615 _d 0 * s)
4e9f2133df Step*1063 C------------------------------------------------------------------------
                1064 C ks = [H][SO4]/[HSO4]
                1065 C dickson (1990, J. chem. Thermodynamics 22, 113)
425561518f Jona*1066 C On the free pH scale
80d7ca01bb Jona*1067 C Note: converted here from mol/kg-H2O to mol/kg-SW
c5830e1a8b Jona*1068            aks(i,j,bi,bj)=exp(-4276.1 _d 0*invtk + 141.328 _d 0 -
                1069      &          23.093 _d 0*dlogtk +
b5953f02df Jona*1070      &          (-13856 _d 0*invtk + 324.57 _d 0 - 47.986 _d 0*dlogtk)
                1071      &          *sqrtis +
                1072      &          (35474 _d 0*invtk - 771.54 _d 0 + 114.723 _d 0*dlogtk)
                1073      &          *is -
c5830e1a8b Jona*1074      &          2698 _d 0*invtk*is**1.5 _d 0 + 1776 _d 0*invtk*is2 +
                1075      &          log(1.0 _d 0 - 0.001005 _d 0*s))
4e9f2133df Step*1076 C------------------------------------------------------------------------
                1077 C kf = [H][F]/[HF]
425561518f Jona*1078 C dickson and Riley (1979) -- change pH scale to total (following Dickson & Goyet, 1994)
80d7ca01bb Jona*1079 C Note: converted here from mol/kg-H2O to mol/kg-SW
92031a2908 Jean*1080            akf(i,j,bi,bj)=exp(1590.2 _d 0*invtk - 12.641 _d 0
                1081      &          + 1.525 _d 0*sqrtis
c5830e1a8b Jona*1082      &          + log(1.0 _d 0 - 0.001005 _d 0*s)
b5953f02df Jona*1083      &          + log(1.0 _d 0 + (0.1400 _d 0/96.062 _d 0)
                1084      &          *(scl)/aks(i,j,bi,bj)))
4e9f2133df Step*1085 C------------------------------------------------------------------------
                1086 C Calculate concentrations for borate, sulfate, and fluoride
                1087 C Uppstrom (1974)
c5830e1a8b Jona*1088            bt(i,j,bi,bj) = 0.000232 _d 0 * scl/10.811 _d 0
4e9f2133df Step*1089 C Morris & Riley (1966)
c5830e1a8b Jona*1090            st(i,j,bi,bj) = 0.14 _d 0 * scl/96.062 _d 0
4e9f2133df Step*1091 C Riley (1965)
c5830e1a8b Jona*1092            ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0
4e9f2133df Step*1093 C------------------------------------------------------------------------
425561518f Jona*1094
                1095 #ifdef CARBONCHEM_TOTALPHSCALE
                1096 C JML Convert between pH scales, we want everything to end up on the total scale
                1097 C ak1 and ak2 are on seawater scale and are pressure corrected, so they just need converting.
                1098 C aks is on the free scale, it needs pressure correcting then converting to the total scale.
                1099 C akf needs pressure correcting too, but is on the right (total) pH scale.
                1100 C akb, ak1p, ak2p, ak3p, aksi and akw are on the total scale. Convert to the seawater
                1101 C  scale before pressure correction, and then convert back to the total scale.
                1102
                1103 C All the terms for pressure correction are from Table 9, Millero (1995), p675
                1104 C Info about pH conversions from Munhoven (2013) and Orr and Epitalon (2015)
c5830e1a8b Jona*1105            total2free_surf = 1. _d 0/
92031a2908 Jean*1106      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                1107
                1108            free2sw_surf = 1. _d 0
                1109      &                 + st(i,j,bi,bj)/ aks(i,j,bi,bj)
                1110      &                 + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free_surf)
425561518f Jona*1111
92031a2908 Jean*1112            total2sw_surf = total2free_surf * free2sw_surf
80d7ca01bb Jona*1113
92031a2908 Jean*1114 C------------------------------------------------------------------------
425561518f Jona*1115 C Pressure correct aks on native free pH scale
c5830e1a8b Jona*1116            dv=-18.03 _d 0 + 0.0466 _d 0*t + 0.316 _d -3 *t*t
                1117            dk=-4.53 _d -3 + 0.09 _d -3 *t
425561518f Jona*1118            pfactor = -(dv/(bigR*tk))*pressc
                1119      &            +(0.5*dk/(bigR*tk))*pressc*pressc
80d7ca01bb Jona*1120            aks(i,j,bi,bj) = aks(i,j,bi,bj)*exp(pfactor)
425561518f Jona*1121
c5830e1a8b Jona*1122            total2free = 1. _d 0/
92031a2908 Jean*1123      &                 (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
                1124
8d230ec051 Jona*1125            free2total = (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
92031a2908 Jean*1126
                1127 C free2sw has an additional component from fluoride
                1128            free2sw    = 1. _d 0
                1129      &                 + st(i,j,bi,bj)/ aks(i,j,bi,bj)
                1130
8d230ec051 Jona*1131            aks(i,j,bi,bj) = aks(i,j,bi,bj)*free2total
92031a2908 Jean*1132
425561518f Jona*1133 C------------------------------------------------------------------------
                1134 C Pressure correct akf on the free scale before converting back to total
                1135            akf(i,j,bi,bj) = akf(i,j,bi,bj) * total2free_surf
c5830e1a8b Jona*1136            dv  =  -9.78 _d 0 -0.0090 _d 0 *t -0.942 _d -3 *t*t
                1137            dk  =  -3.91 _d -3 + 0.054 _d -3 *t
425561518f Jona*1138            pfactor = -(dv/(bigR*tk))*pressc
                1139      &               +(0.5*dk/(bigR*tk))*pressc*pressc
                1140            akf(i,j,bi,bj) = akf(i,j,bi,bj) * exp(pfactor)
92031a2908 Jean*1141            akf(i,j,bi,bj) = akf(i,j,bi,bj)*free2total
                1142
                1143 C free2sw has an additional component from fluoride, add it here
                1144            free2sw = free2sw
                1145      &               + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free)
b5953f02df Jona*1146
8d230ec051 Jona*1147            sw2free  = 1. _d 0 / free2sw
92031a2908 Jean*1148
                1149            total2sw = total2free * free2sw
425561518f Jona*1150
8d230ec051 Jona*1151            sw2total = 1. _d 0 / total2sw
92031a2908 Jean*1152
425561518f Jona*1153 C------------------------------------------------------------------------
                1154 C Convert pressure corrected ak1 and ak2 from the seawater pH scale to the total scale
8d230ec051 Jona*1155            ak1(i,j,bi,bj)  = ak1(i,j,bi,bj)*sw2total
                1156            ak2(i,j,bi,bj)  = ak2(i,j,bi,bj)*sw2total
425561518f Jona*1157
                1158 C------------------------------------------------------------------------
                1159 C Pressure correct akb on the seawater scale before converting back to total scale
c5830e1a8b Jona*1160             dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t
                1161             dk = -2.84 _d -3
425561518f Jona*1162             pfactor = - (dv/(bigR*tk))*pressc
                1163      &                + (0.5*dk/(bigR*tk))*pressc*pressc
                1164             akb(i,j,bi,bj) = total2sw_surf*akb(i,j,bi,bj)*exp(pfactor)
8d230ec051 Jona*1165             akb(i,j,bi,bj) = akb(i,j,bi,bj)*sw2total
425561518f Jona*1166
                1167 C------------------------------------------------------------------------
                1168 C Pressure correct ak1p on the seawater scale before converting back to total scale
c5830e1a8b Jona*1169             dv = -14.51 _d 0 + 0.1211 _d 0*t -0.321 _d -3*t*t
                1170             dk = -2.67 _d -3 + 0.0427 _d -3*t
425561518f Jona*1171             pfactor = - (dv/(bigR*tk))*pressc
                1172      &                + (0.5*dk/(bigR*tk))*pressc*pressc
                1173             ak1p(i,j,bi,bj) = total2sw_surf*ak1p(i,j,bi,bj)*exp(pfactor)
8d230ec051 Jona*1174             ak1p(i,j,bi,bj) = ak1p(i,j,bi,bj)*sw2total
425561518f Jona*1175
                1176 C------------------------------------------------------------------------
                1177 C Pressure correct ak2p on the seawater scale before converting back to total scale
c5830e1a8b Jona*1178             dv = -23.12 _d 0 + 0.1758 _d 0*t -2.647 _d -3*t*t
                1179             dk = -5.15 _d -3 + 0.09 _d -3*t
425561518f Jona*1180             pfactor = - (dv/(bigR*tk))*pressc
                1181      &                + (0.5*dk/(bigR*tk))*pressc*pressc
                1182             ak2p(i,j,bi,bj) = total2sw_surf*ak2p(i,j,bi,bj)*exp(pfactor)
8d230ec051 Jona*1183             ak2p(i,j,bi,bj) = ak2p(i,j,bi,bj)*sw2total
92031a2908 Jean*1184
425561518f Jona*1185 C------------------------------------------------------------------------
                1186 C Pressure correct ak3p on the seawater scale before converting back to total scale
c5830e1a8b Jona*1187             dv = -26.57 _d 0 + 0.2020 _d 0*t -3.042 _d -3*t*t
                1188             dk = -4.08 _d -3 + 0.0714 _d -3*t
425561518f Jona*1189             pfactor = - (dv/(bigR*tk))*pressc
                1190      &                + (0.5*dk/(bigR*tk))*pressc*pressc
                1191             ak3p(i,j,bi,bj) = total2sw_surf*ak3p(i,j,bi,bj)*exp(pfactor)
8d230ec051 Jona*1192             ak3p(i,j,bi,bj) = ak3p(i,j,bi,bj)*sw2total
92031a2908 Jean*1193
425561518f Jona*1194 C------------------------------------------------------------------------
                1195 C Pressure correct akw on the seawater scale before converting back to total scale
c5830e1a8b Jona*1196             dv = -20.02 _d 0 + 0.1119 _d 0*t -1.409 _d -3*t*t
                1197             dk = -5.13 _d -3 + 0.0794 _d -3 *t
425561518f Jona*1198             pfactor = - (dv/(bigR*tk))*pressc
                1199      &                + (0.5*dk/(bigR*tk))*pressc*pressc
                1200             akw(i,j,bi,bj) = total2sw_surf*akw(i,j,bi,bj)*exp(pfactor)
8d230ec051 Jona*1201             akw(i,j,bi,bj) = akw(i,j,bi,bj)*sw2total
425561518f Jona*1202
                1203 C------------------------------------------------------------------------
                1204 C Pressure correct aksi on the seawater scale before converting back to total scale
                1205 C Estimated from the effect on akb
c5830e1a8b Jona*1206             dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t
                1207             dk = -2.84 _d -3
425561518f Jona*1208             pfactor = - (dv/(bigR*tk))*pressc
                1209      &                + (0.5*dk/(bigR*tk))*pressc*pressc
                1210             aksi(i,j,bi,bj) = total2sw_surf*aksi(i,j,bi,bj)*exp(pfactor)
8d230ec051 Jona*1211             aksi(i,j,bi,bj) = aksi(i,j,bi,bj)*sw2total
425561518f Jona*1212 #endif
                1213
                1214 C------------------------------------------------------------------------
                1215 C Solubility product for calcite
4e9f2133df Step*1216 C
                1217 c Following Takahashi (1982) GEOSECS handbook
                1218 C NOT SURE THIS IS WORKING???
                1219 C Ingle et al. (1973)
                1220 c          Ksp_T_Calc = ( -34.452 - 39.866*(s**0.333333)
                1221 c    &                  + 110.21*log(s) - 7.5752d-6 * (tk**2.0)
                1222 c    &                  ) * 1.0d-7
                1223 c with pressure dependence Culberson and Pytkowicz (1968)
                1224 c          xvalue  =  (36-0.20*t)*(pressc-1.0)/(83.143*tk)
                1225 c          Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*exp(xvalue)
                1226 c
                1227 c
                1228 C Following Mucci (1983) - from Zeebe/Wolf-Gladrow equic.m
92031a2908 Jean*1229             tmpa1 = - 171.9065 _d 0 - (0.077993 _d 0*tk)
c5830e1a8b Jona*1230      &            + (2839.319 _d 0/tk) + (71.595 _d 0*log10(tk))
92031a2908 Jean*1231             tmpa2 = +(-0.77712 _d 0 + (0.0028426 _d 0*tk)
c5830e1a8b Jona*1232      &           + (178.34 _d 0/tk) )*sqrts
                1233             tmpa3 = -(0.07711 _d 0*s) + (0.0041249 _d 0*s15)
                1234             logKspc = tmpa1 + tmpa2 + tmpa3
                1235             Ksp_T_Calc = 10.0 _d 0**logKspc
4e9f2133df Step*1236 c        write(6,*)i,j,k,tmpa1,tmpa2,tmpa3,logkspc,Ksp_T_Calc
                1237 c with pressure dependence Culberson and Pytkowicz (1968)
                1238 c        xvalue  =  (36.0-0.20*t)*(pressc-1.0)/(83.143*tk)
                1239 c        Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*exp(xvalue)
                1240
                1241 c alternative pressure depdendence
                1242 c following Millero (1995) but using info from Appendix A11 of
                1243 c Zeebe and Wolf-Gladrow (2001) book
                1244 c          dv = -48.6 - 0.5304*t
                1245 c          dk = -11.76d-3 - 0.3692*t
                1246 c          xvalue = - (dv/(bigR*tk))*pressc
                1247 c    &               + (0.5*dk/(bigR*tk))*pressc*pressc
                1248 c          Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*exp(xvalue)
                1249
                1250 c alternative pressure dependence from Ingle (1975)
                1251
c5830e1a8b Jona*1252             zdum   = (pressc*10. _d 0 - 10. _d 0)/10. _d 0
                1253             xvalue = ( (48.8 _d 0 - 0.53 _d 0*t)*zdum
                1254      &                 + (-0.00588 _d 0 + 0.0001845 _d 0*t)*zdum*zdum)
                1255      &            / (188.93 _d 0*(t + 273.15 _d 0))
4e9f2133df Step*1256
c5830e1a8b Jona*1257             Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*10**(xvalue)
4e9f2133df Step*1258
                1259 C------------------------------------------------------------------------
f0dfb321df Jean*1260           ELSE
c5830e1a8b Jona*1261             ff(i,j,bi,bj)=0. _d 0
                1262             ak0(i,j,bi,bj)= 0. _d 0
                1263             ak1(i,j,bi,bj)= 0. _d 0
                1264             ak2(i,j,bi,bj)= 0. _d 0
                1265             akb(i,j,bi,bj)= 0. _d 0
                1266             ak1p(i,j,bi,bj) = 0. _d 0
                1267             ak2p(i,j,bi,bj) = 0. _d 0
                1268             ak3p(i,j,bi,bj) = 0. _d 0
                1269             aksi(i,j,bi,bj) = 0. _d 0
                1270             akw(i,j,bi,bj) = 0. _d 0
                1271             aks(i,j,bi,bj)= 0. _d 0
                1272             akf(i,j,bi,bj)= 0. _d 0
                1273             bt(i,j,bi,bj) = 0. _d 0
                1274             st(i,j,bi,bj) = 0. _d 0
                1275             ft(i,j,bi,bj) = 0. _d 0
                1276             Ksp_TP_Calc(i,j,bi,bj) = 0. _d 0
f0dfb321df Jean*1277           ENDIF
                1278          enddo
                1279         enddo
4e9f2133df Step*1280
f0dfb321df Jean*1281         RETURN
                1282         END