pkg/bling/bling_carbon_chem.F

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