pkg/darwin/darwin_carbon_chem.F

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