C $Header: /u/gcmpack/MITgcm/pkg/bling/bling_carbon_chem.F,v 1.4 2018/01/10 19:36:03 mmazloff Exp $
C $Name:  $

#include "BLING_OPTIONS.h"

C--  File bling_carbon_chem.F:
C--   Contents
C--   o CALC_PCO2
C--   o CALC_PCO2_APPROX
C--   o CARBON_COEFFS
C--   o CARBON_COEFFS_PRESSURE_DEP

CBOP
       subroutine CALC_PCO2(
     I                       donewt,inewtonmax,ibrackmax,
     I                       t,s,diclocal,pt,sit,ta,
     I                       k1local,k2local,
     I                       k1plocal,k2plocal,k3plocal,
     I                       kslocal,kblocal,kwlocal,
     I                       ksilocal,kflocal,
     I                       k0local, fugflocal,
     I                       fflocal,btlocal,stlocal,ftlocal,
     U                       pHlocal,pCO2surfloc,
     I                       i,j,k,bi,bj,myIter,myThid )

C     =================================================================
C     | subroutine bling_carbon_chem
C     | o Calculates ocean inorganic carbon chemistry
C     |   Adapted from pkg/dic/carbon_chem.F
C     |   from OCMIP2 code
C     =================================================================

      implicit none
      
C     === Global variables ===
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "FFIELDS.h"
#include "BLING_VARS.h"

C     == Routine arguments ==
C       diclocal = total inorganic carbon (mol/m^3)
C             where 1 T = 1 metric ton = 1000 kg
C       ta  = total alkalinity (eq/m^3)
C       pt  = inorganic phosphate (mol/^3)
C       sit = inorganic silicate (mol/^3)
C       t   = temperature (degrees C)
C       s   = salinity (PSU)
        INTEGER donewt
        INTEGER inewtonmax
        INTEGER ibrackmax
        _RL  t, s, pt, sit, ta
        _RL  pCO2surfloc, diclocal, pHlocal
        _RL  fflocal, btlocal, stlocal, ftlocal
        _RL  k1local, k2local
        _RL  k1plocal, k2plocal, k3plocal
        _RL  kslocal, kblocal, kwlocal, ksilocal, kflocal
        _RL  k0local, fugflocal
        INTEGER i,j,k,bi,bj,myIter
        INTEGER myThid
CEOP

C     == Local variables ==
C INPUT
C       phlo= lower limit of pH range
C       phhi= upper limit of pH range
C       atmpres = atmospheric pressure in atmospheres (1 atm==1013.25mbar)
C OUTPUT
C       co2star  = CO2*water (mol/m^3)
C       pco2surf = oceanic pCO2 (ppmv)
C ---------------------------------------------------------------------
C OCMIP NOTE: Some words about units - (JCO, 4/4/1999)
C     - Models carry tracers in mol/m^3 (on a per volume basis)
C     - Conversely, this routine, which was written by
C       observationalists (C. Sabine and R. Key), passes input
C       arguments in umol/kg (i.e., on a per mass basis)
C     - I have changed things slightly so that input arguments are in
C       mol/m^3,
C     - Thus, all input concentrations (diclocal, ta, pt, and st) should be
C       given in mol/m^3; output arguments "co2star" and "dco2star"
C       are likewise be in mol/m^3.
C ---------------------------------------------------------------------

        _RL  phhi
        _RL  phlo
        _RL  c
        _RL  a
        _RL  a2
        _RL  da
        _RL  b
        _RL  b2
        _RL  db
        _RL  fn
        _RL  df
        _RL  deltax
        _RL  x
        _RL  x2
        _RL  x3
        _RL  xmid
        _RL  ftest
        _RL  htotal
        _RL  htotal2
        _RL  co2star
        _RL  phguess
        _RL  fco2
        INTEGER inewton
        INTEGER ibrack
        INTEGER hstep
        _RL  fni(3)
        _RL  xlo
        _RL  xhi
        _RL  xguess
        _RL  k123p
        _RL  k12p
        _RL  k12
c ---------------------------------------------------------------------
c import donewt flag
c set donewt = 1 for newton-raphson iteration
c set donewt = 0 for bracket and bisection
c ---------------------------------------------------------------------
C Change units from the input of mol/m^3 -> mol/kg:
c (1 mol/m^3)  x (1 m^3/1024.5 kg)
c where the ocean mean surface density is 1024.5 kg/m^3
c Note: mol/kg are actually what the body of this routine uses
c for calculations.  Units are reconverted back to mol/m^3 at the
c end of this routine.
c ---------------------------------------------------------------------
c To convert input in mol/m^3 -> mol/kg
        pt=pt*permil
        sit=sit*permil
        ta=ta*permil
        diclocal=diclocal*permil
c ---------------------------------------------------------------------
c set first guess and brackets for [H+] solvers
c first guess (for newton-raphson)
        phguess = phlocal


c bracketing values (for bracket/bisection)
        phhi = 10.0
        phlo = 5.0
c convert to [H+]...
        xguess = 10.0**(-phguess)
        xlo = 10.0**(-phhi)
        xhi = 10.0**(-phlo)
        xmid = (xlo + xhi)*0.5


c----------------------------------------------------------------
c iteratively solve for [H+]
c (i) Newton-Raphson method with fixed number of iterations,
c     use previous [H+] as first guess

c select newton-raphson, inewt=1
c else select bracket and bisection

cQQQQQ
        if( donewt .eq. 1)then
c.........................................................
c NEWTON-RAPHSON METHOD
c.........................................................
          x = xguess
cdiags
c         WRITE(0,*)'xguess ',xguess
cdiags
          do inewton = 1, inewtonmax
c set some common combinations of parameters used in
c the iterative [H+] solvers
            x2=x*x
            x3=x2*x
            k12 = k1local*k2local
            k12p = k1plocal*k2plocal
            k123p = k12p*k3plocal
            c = 1.0 + stlocal/kslocal
            a = x3 + k1plocal*x2 + k12p*x + k123p
            a2=a*a
            da = 3.0*x2 + 2.0*k1plocal*x + k12p
            b = x2 + k1local*x + k12
            b2=b*b
            db = 2.0*x + k1local

c Evaluate f([H+]) and f_prime([H+])
c fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree
c      +hso4+hf+h3po4-ta
            fn = k1local*x*diclocal/b +
     &        2.0*diclocal*k12/b +
     &        btlocal/(1.0 + x/kblocal) +
     &        kwlocal/x +
     &        pt*k12p*x/a +
     &        2.0*pt*k123p/a +
     &        sit/(1.0 + x/ksilocal) -
     &        x/c -
     &        stlocal/(1.0 + kslocal/x/c) -
     &        ftlocal/(1.0 + kflocal/x) -
     &        pt*x3/a -
     &        ta

c df = dfn/dx
cdiags
c      WRITE(0,*)'values',b2,kblocal,x2,a2,c,x
cdiags
            df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -
     &        2.0*diclocal*k12*db/b2 -
     &        btlocal/kblocal/(1.0+x/kblocal)**2. -
     &        kwlocal/x2 +
     &        (pt*k12p*(a - x*da))/a2 -
     &        2.0*pt*k123p*da/a2 -
     &        sit/ksilocal/(1.0+x/ksilocal)**2. +
     &        1.0/c +
     &        stlocal*(1.0 + kslocal/x/c)**(-2.0)*(kslocal/c/x2) +
     &        ftlocal*(1.0 + kflocal/x)**(-2.)*kflocal/x2 -
     &        pt*x2*(3.0*a-x*da)/a2
c evaluate increment in [H+]
            deltax = - fn/df
c update estimate of [H+]
            x = x + deltax
cdiags
c write value of x to check convergence....
c           write(0,*)'inewton, x, deltax ',inewton, x, deltax
c           write(6,*)
cdiags

          end


do c end of newton-raphson method c.................................................... else c.................................................... C BRACKET AND BISECTION METHOD c.................................................... c (ii) If first step use Bracket and Bisection method c with fixed, large number of iterations do ibrack = 1, ibrackmax do hstep = 1,3 if(hstep .eq. 1)x = xhi if(hstep .eq. 2)x = xlo if(hstep .eq. 3)x = xmid c set some common combinations of parameters used in c the iterative [H+] solvers x2=x*x x3=x2*x k12 = k1local*k2local k12p = k1plocal*k2plocal k123p = k12p*k3plocal c = 1.0 + stlocal/kslocal a = x3 + k1plocal*x2 + k12p*x + k123p a2=a*a da = 3.0*x2 + 2.0*k1plocal*x + k12p b = x2 + k1local*x + k12 b2=b*b db = 2.0*x + k1local c evaluate f([H+]) for bracketing and mid-value cases fn = k1local*x*diclocal/b + & 2.0*diclocal*k12/b + & btlocal/(1.0 + x/kblocal) + & kwlocal/x + & pt*k12p*x/a + & 2.0*pt*k123p/a + & sit/(1.0 + x/ksilocal) - & x/c - & stlocal/(1.0 + kslocal/x/c) - & ftlocal/(1.0 + kflocal/x) - & pt*x3/a - & ta fni(hstep) = fn end


do c now bracket solution within two of three ftest = fni(1)/fni(3) if(ftest .gt. 0.0)then xhi = xmid else xlo = xmid end


if xmid = (xlo + xhi)*0.5 cdiags c write value of x to check convergence.... c WRITE(0,*)'bracket-bisection iteration ',ibrack, xmid cdiags end


do c last iteration gives value x = xmid c end of bracket and bisection method c.................................... end


if c iterative [H+] solver finished c---------------------------------------------------------------- c now determine pCO2 etc... c htotal = [H+], hydrogen ion conc htotal = x C Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2, C ORNL/CDIAC-74, dickson and Goyet, eds. (Ch 2 p 10, Eq A.49) htotal2=htotal*htotal co2star=diclocal*htotal2/(htotal2 + k1local*htotal & + k1local*k2local) phlocal=-log10(htotal) c --------------------------------------------------------------- c Add two output arguments for storing pCO2surf c Should we be using K0 or ff for the solubility here? c --------------------------------------------------------------- fco2 = co2star / k0local pCO2surfloc = fco2/fugflocal C ---------------------------------------------------------------- C Reconvert units back to original values for input arguments C no longer necessary???? C ---------------------------------------------------------------- c Reconvert from mol/kg -> mol/m^3 pt=pt/permil sit=sit/permil ta=ta/permil diclocal=diclocal/permil RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP subroutine CALC_PCO2_APPROX( I t,s,diclocal,pt,sit,ta, I k1local,k2local, I k1plocal,k2plocal,k3plocal, I kslocal,kblocal,kwlocal, I ksilocal,kflocal, I k0local, fugflocal, I fflocal,btlocal,stlocal,ftlocal, U pHlocal,pCO2surfloc,co3local, I i,j,k,bi,bj,myIter,myThid ) C !DESCRIPTION: C *==========================================================* C | SUBROUTINE CALC_PCO2_APPROX | C *==========================================================* C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C New efficient pCO2 solver, Mick Follows CC C C Taka Ito CC C C Stephanie Dutkiewicz CC C C 20 April 2003 CC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Apr 2011: fix vapour bug (following Bennington) C Oct 2011: add CO3 extimation and pass out implicit none C === Global variables === #include "SIZE.h" #include "DYNVARS.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "BLING_VARS.h" C == Routine arguments == C diclocal = total inorganic carbon (mol/m^3) C where 1 T = 1 metric ton = 1000 kg C ta = total alkalinity (eq/m^3) C pt = inorganic phosphate (mol/^3) C sit = inorganic silicate (mol/^3) C t = temperature (degrees C) C s = salinity (PSU) _RL t, s, pt, sit, ta _RL pCO2surfloc, diclocal, pHlocal _RL fflocal, btlocal, stlocal, ftlocal _RL k1local, k2local _RL k1plocal, k2plocal, k3plocal _RL kslocal, kblocal, kwlocal, ksilocal, kflocal _RL k0local, fugflocal _RL co3local INTEGER i,j,k,bi,bj,myIter INTEGER myThid CEOP C == Local variables == _RL phguess _RL cag _RL bohg _RL hguess _RL stuff _RL gamm _RL hnew _RL co2s _RL h3po4g, h2po4g, hpo4g, po4g _RL siooh3g _RL fco2 c --------------------------------------------------------------------- C Change units from the input of mol/m^3 -> mol/kg: c (1 mol/m^3) x (1 m^3/1024.5 kg) c where the ocean mean surface density is 1024.5 kg/m^3 c Note: mol/kg are actually what the body of this routine uses c for calculations. Units are reconverted back to mol/m^3 at the c end of this routine. c To convert input in mol/m^3 -> mol/kg pt=pt*permil sit=sit*permil ta=ta*permil diclocal=diclocal*permil c --------------------------------------------------------------------- c set first guess and brackets for [H+] solvers c first guess (for newton-raphson) phguess = phlocal cmick - new approx method cmick - make estimate of htotal (hydrogen ion conc) using cmick appromate estimate of CA, carbonate alkalinity hguess = 10.0**(-phguess) cmick - first estimate borate contribution using guess for [H+] bohg = btlocal*kblocal/(hguess+kblocal) cmick - first estimate of contribution from phosphate cmick based on Dickson and Goyet stuff = hguess*hguess*hguess & + (k1plocal*hguess*hguess) & + (k1plocal*k2plocal*hguess) & + (k1plocal*k2plocal*k3plocal) h3po4g = (pt*hguess*hguess*hguess) / stuff h2po4g = (pt*k1plocal*hguess*hguess) / stuff hpo4g = (pt*k1plocal*k2plocal*hguess) / stuff po4g = (pt*k1plocal*k2plocal*k3plocal) / stuff cmick - estimate contribution from silicate cmick based on Dickson and Goyet siooh3g = sit*ksilocal / (ksilocal + hguess) cmm - rare, but ice melt can cause prescribed silicate alkalinity to be cmm larger than total alkalinity. Crop here to account for this IF (siooh3g.GT.0.2 _d 0*ta) siooh3g = 0.2 _d 0*ta cmick - now estimate carbonate alkalinity cag = ta - bohg - (kwlocal/hguess) + hguess & - hpo4g - 2.0 _d 0*po4g + h3po4g & - siooh3g cmick - now evaluate better guess of hydrogen ion conc cmick htotal = [H+], hydrogen ion conc gamm = diclocal/cag stuff = (1.0 _d 0-gamm)*(1.0 _d 0-gamm)*k1local*k1local & - 4.0 _d 0*k1local*k2local*(1.0 _d 0-2.0 _d 0*gamm) hnew = 0.5 _d 0*( (gamm-1.0 _d 0)*k1local + sqrt(stuff) ) cmick - now determine [CO2*] co2s = diclocal/ & (1.0 _d 0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew))) cmick - return update pH to main routine phlocal = -log10(hnew) c NOW EVALUATE CO32-, carbonate ion concentration c used in determination of calcite compensation depth c Karsten Friis & Mick - Sep 2004 co3local = k1local*k2local*diclocal / & (hnew*hnew + k1local*hnew + k1local*k2local) c --------------------------------------------------------------- c surface pCO2 (following Dickson and Goyet, DOE...) fco2 = co2s/k0local pco2surfloc = fco2/fugflocal C ---------------------------------------------------------------- c Reconvert from mol/kg -> mol/m^3 pt=pt/permil sit=sit/permil ta=ta/permil diclocal=diclocal/permil RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP subroutine CARBON_COEFFS( I ttemp,stemp, I bi,bj,iMin,iMax,jMin,jMax,myThid) C !DESCRIPTION: C *==========================================================* C | SUBROUTINE CARBON_COEFFS | C | determine coefficients for surface carbon chemistry | C | adapted from OCMIP2: SUBROUTINE CO2CALC | C | mick follows, oct 1999 | C | minor changes to tidy, swd aug 2002 | C *==========================================================* C INPUT C diclocal = total inorganic carbon (mol/m^3) C where 1 T = 1 metric ton = 1000 kg C ta = total alkalinity (eq/m^3) C pt = inorganic phosphate (mol/^3) C sit = inorganic silicate (mol/^3) C t = temperature (degrees C) C s = salinity (PSU) C OUTPUT C IMPORTANT: Some words about units - (JCO, 4/4/1999) C - Models carry tracers in mol/m^3 (on a per volume basis) C - Conversely, this routine, which was written by observationalists C (C. Sabine and R. Key), passes input arguments in umol/kg C (i.e., on a per mass basis) C - I have changed things slightly so that input arguments are in mol/m^3, C - Thus, all input concentrations (diclocal, ta, pt, and st) should be C given in mol/m^3; output arguments "co2star" and "dco2star" C are likewise be in mol/m^3. C C Apr 2011: fix vapour bug (following Bennington) C-------------------------------------------------------------------------- implicit none C === Global variables === #include "SIZE.h" #include "DYNVARS.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "BLING_VARS.h" CMM( #ifdef ALLOW_EXF # include "EXF_FIELDS.h" #endif CMM) C == Routine arguments == C ttemp and stemp are local theta and salt arrays C dont really need to pass T and S in, could use theta, salt in C common block in DYNVARS.h, but this way keeps subroutine more C general _RL ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER bi,bj,iMin,iMax,jMin,jMax INTEGER myThid CEOP C LOCAL VARIABLES _RL t _RL s _RL tk _RL tk100 _RL tk1002 _RL dlogtk _RL sqrtis _RL sqrts _RL s15 _RL scl _RL s2 _RL invtk _RL is _RL is2 c add Bennington _RL P1atm _RL Rgas _RL RT _RL delta _RL B1 _RL B INTEGER i INTEGER j C..................................................................... C OCMIP note: C Calculate all constants needed to convert between various measured C carbon species. References for each equation are noted in the code. C Once calculated, the constants are C stored and passed in the common block "const". The original version C of this code was based on the code by dickson in Version 2 of C Handbook of Methods C for the Analysis of the Various Parameters of C the Carbon Dioxide System in Seawater , DOE, 1994 (SOP No. 3, p25-26). C.................................................................... do i=imin,imax do j=jmin,jmax if (hFacC(i,j,1,bi,bj).gt.0. _d 0) then t = ttemp(i,j) s = stemp(i,j) C terms used more than once tk = 273.15 _d 0 + t tk100 = tk/100. _d 0 tk1002=tk100*tk100 invtk=1.0 _d 0/tk dlogtk=log(tk) is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s) is2=is*is sqrtis=sqrt(is) s2=s*s sqrts=sqrt(s) s15=s**1.5 _d 0 scl=s/1.80655 _d 0 C ----------------------------------------------------------------------- C added by Val Bennington Nov 2010 C Fugacity Factor needed for non-ideality in ocean C ff used for atmospheric correction for water vapor and pressure C Weiss (1974) Marine Chemistry P1atm = 1.01325 _d 0 ! bars Rgas = 83.1451 _d 0 ! bar*cm3/(mol*K) RT = Rgas*tk delta = (57.7 _d 0 - 0.118 _d 0*tk) B1 = -1636.75 _d 0 + 12.0408 _d 0*tk - 0.0327957 _d 0*tk*tk B = B1 + 3.16528 _d 0*tk*tk*tk*(0.00001 _d 0) #ifdef USE_EXF_ATMPRES C Atm pressure in N/m2, convert to bars fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * & (apressure (i,j,bi,bj) * 0.00001) & / RT) #else fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * P1atm / RT) #endif C------------------------------------------------------------------------ C f = k0(1-pH2O)*correction term for non-ideality C Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values) ff(i,j,bi,bj) = exp(-162.8301 _d 0 + 218.2968 _d 0/tk100 + & 90.9241 _d 0*log(tk100) - 1.47696 _d 0*tk1002 + & s * (.025695 _d 0 - .025225 _d 0*tk100 + & 0.0049867 _d 0*tk1002)) C------------------------------------------------------------------------ C K0 from Weiss 1974 ak0(i,j,bi,bj) = exp(93.4517 _d 0/tk100 - 60.2409 _d 0 + & 23.3585 _d 0 * log(tk100) + & s * (0.023517 _d 0 - 0.023656 _d 0*tk100 + & 0.0047036 _d 0*tk1002)) C------------------------------------------------------------------------ C k1 = [H][HCO3]/[H2CO3] C k2 = [H][CO3]/[HCO3] C Millero p.664 (1995) using Mehrbach et al. data on seawater scale ak1(i,j,bi,bj)=10.**(-1. _d 0*(3670.7 _d 0*invtk - & 62.008 _d 0 + 9.7944 _d 0*dlogtk - & 0.0118 _d 0 * s + 0.000116 _d 0*s2)) ak2(i,j,bi,bj)=10.**(-1. _d 0*(1394.7 _d 0*invtk+ 4.777 _d 0- & 0.0184 _d 0*s + 0.000118 _d 0*s2)) C------------------------------------------------------------------------ C kb = [H][BO2]/[HBO2] C Millero p.669 (1995) using data from dickson (1990) akb(i,j,bi,bj)=exp((-8966.90 _d 0- 2890.53 _d 0*sqrts - & 77.942 _d 0*s + 1.728 _d 0*s15 - 0.0996 _d 0*s2)*invtk + & (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) + & (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) * & dlogtk + 0.053105 _d 0*sqrts*tk) C------------------------------------------------------------------------ C k1p = [H][H2PO4]/[H3PO4] C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974) ak1p(i,j,bi,bj) = exp(-4576.752 _d 0*invtk + 115.525 _d 0 - & 18.453 _d 0*dlogtk + & (-106.736 _d 0*invtk + 0.69171 _d 0)*sqrts + & (-0.65643 _d 0*invtk - 0.01844 _d 0)*s) C------------------------------------------------------------------------ C k2p = [H][HPO4]/[H2PO4] C DOE(1994) eq 7.2.23 with footnote using data from Millero (1974)) ak2p(i,j,bi,bj) = exp(-8814.715 _d 0*invtk + 172.0883 _d 0 - & 27.927 _d 0*dlogtk + & (-160.340 _d 0*invtk + 1.3566 _d 0) * sqrts + & (0.37335 _d 0*invtk - 0.05778 _d 0) * s) C------------------------------------------------------------------------ C k3p = [H][PO4]/[HPO4] C DOE(1994) eq 7.2.26 with footnote using data from Millero (1974) ak3p(i,j,bi,bj) = exp(-3070.75 _d 0*invtk - 18.141 _d 0 + & (17.27039 _d 0*invtk + 2.81197 _d 0) * & sqrts + (-44.99486 _d 0*invtk - 0.09984 _d 0) * s) C------------------------------------------------------------------------ C ksi = [H][SiO(OH)3]/[Si(OH)4] C Millero p.671 (1995) using data from Yao and Millero (1995) aksi(i,j,bi,bj) = exp(-8904.2 _d 0*invtk + 117.385 _d 0 - & 19.334 _d 0*dlogtk + & (-458.79 _d 0*invtk + 3.5913 _d 0) * sqrtis + & (188.74 _d 0*invtk - 1.5998 _d 0) * is + & (-12.1652 _d 0*invtk + 0.07871 _d 0) * is2 + & log(1.0 _d 0-0.001005 _d 0*s)) C------------------------------------------------------------------------ C kw = [H][OH] C Millero p.670 (1995) using composite data akw(i,j,bi,bj) = exp(-13847.26 _d 0*invtk + 148.9652 _d 0 - & 23.6521 _d 0*dlogtk + & (118.67 _d 0*invtk - 5.977 _d 0 + 1.0495 _d 0 * dlogtk) & * sqrts - 0.01615 _d 0 * s) C------------------------------------------------------------------------ C ks = [H][SO4]/[HSO4] C dickson (1990, J. chem. Thermodynamics 22, 113) aks(i,j,bi,bj)=exp(-4276.1 _d 0*invtk + 141.328 _d 0 - & 23.093 _d 0*dlogtk + & (-13856. _d 0*invtk + 324.57 _d 0 - 47.986 _d 0*dlogtk)*sqrtis+ & (35474. _d 0*invtk - 771.54 _d 0 + 114.723 _d 0*dlogtk)*is - & 2698. _d 0*invtk*is**1.5 _d 0 + 1776. _d 0*invtk*is2 + & log(1.0 _d 0 - 0.001005 _d 0*s)) C------------------------------------------------------------------------ C kf = [H][F]/[HF] C dickson and Riley (1979) -- change pH scale to total akf(i,j,bi,bj)=exp(1590.2 _d 0*invtk - 12.641 _d 0 + & 1.525 _d 0*sqrtis + log(1.0 _d 0 - 0.001005 _d 0*s) + & log(1.0 _d 0 + (0.1400 _d 0/96.062 _d 0)*(scl)/aks(i,j,bi,bj))) C------------------------------------------------------------------------ C Calculate concentrations for borate, sulfate, and fluoride C Uppstrom (1974) bt(i,j,bi,bj) = 0.000232 _d 0 * scl/10.811 _d 0 C Morris & Riley (1966) st(i,j,bi,bj) = 0.14 _d 0 * scl/96.062 _d 0 C Riley (1965) ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0 C------------------------------------------------------------------------ else c add Bennington fugf(i,j,bi,bj)=0. _d 0 ff(i,j,bi,bj)=0. _d 0 ak0(i,j,bi,bj)= 0. _d 0 ak1(i,j,bi,bj)= 0. _d 0 ak2(i,j,bi,bj)= 0. _d 0 akb(i,j,bi,bj)= 0. _d 0 ak1p(i,j,bi,bj) = 0. _d 0 ak2p(i,j,bi,bj) = 0. _d 0 ak3p(i,j,bi,bj) = 0. _d 0 aksi(i,j,bi,bj) = 0. _d 0 akw(i,j,bi,bj) = 0. _d 0 aks(i,j,bi,bj)= 0. _d 0 akf(i,j,bi,bj)= 0. _d 0 bt(i,j,bi,bj) = 0. _d 0 st(i,j,bi,bj) = 0. _d 0 ft(i,j,bi,bj) = 0. _d 0 endif end


do end


do RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP subroutine CARBON_COEFFS_PRESSURE_DEP( I ttemp,stemp, I bi,bj,iMin,iMax,jMin,jMax, I Klevel,myThid) C !DESCRIPTION: C *==========================================================* C | SUBROUTINE CARBON_COEFFS | C | determine coefficients for surface carbon chemistry | C | adapted from OCMIP2: SUBROUTINE CO2CALC | C | mick follows, oct 1999 | C | minor changes to tidy, swd aug 2002 | C | MODIFIED FOR PRESSURE DEPENDENCE | C | Karsten Friis and Mick Follows 2004 | C *==========================================================* C INPUT C diclocal = total inorganic carbon (mol/m^3) C where 1 T = 1 metric ton = 1000 kg C ta = total alkalinity (eq/m^3) C pt = inorganic phosphate (mol/^3) C sit = inorganic silicate (mol/^3) C t = temperature (degrees C) C s = salinity (PSU) C OUTPUT C IMPORTANT: Some words about units - (JCO, 4/4/1999) C - Models carry tracers in mol/m^3 (on a per volume basis) C - Conversely, this routine, which was written by observationalists C (C. Sabine and R. Key), passes input arguments in umol/kg C (i.e., on a per mass basis) C - I have changed things slightly so that input arguments are in mol/m^3, C - Thus, all input concentrations (diclocal, ta, pt, and st) should be C given in mol/m^3; output arguments "co2star" and "dco2star" C are likewise be in mol/m^3. C C NOW INCLUDES: C PRESSURE DEPENDENCE of K1, K2, solubility product of calcite C based on Takahashi, GEOSECS Atlantic Report, Vol. 1 (1981) C C-------------------------------------------------------------------------- implicit none C === Global variables === #include "SIZE.h" #include "DYNVARS.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "BLING_VARS.h" CMM( #ifdef ALLOW_EXF # include "EXF_FIELDS.h" #endif CMM) C == Routine arguments == C ttemp and stemp are local theta and salt arrays C dont really need to pass T and S in, could use theta, salt in C common block in DYNVARS.h, but this way keeps subroutine more C general _RL ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER bi,bj,iMin,iMax,jMin,jMax C K is depth index INTEGER Klevel INTEGER myThid CEOP C LOCAL VARIABLES _RL t _RL s _RL tk _RL tk100 _RL tk1002 _RL dlogtk _RL sqrtis _RL sqrts _RL s15 _RL scl _RL s2 _RL invtk _RL is _RL is2 c add Bennington _RL P1atm _RL Rgas _RL RT _RL delta _RL B1 _RL B INTEGER i INTEGER j INTEGER k _RL bdepth _RL cdepth _RL pressc _RL Ksp_T_Calc _RL Ksp_T_Arag _RL xvalue _RL zdum _RL tmpa1 _RL tmpa2 _RL tmpa3 _RL logKspc _RL logKspa _RL dv _RL dk _RL pfactor _RL bigR C..................................................................... C OCMIP note: C Calculate all constants needed to convert between various measured C carbon species. References for each equation are noted in the code. C Once calculated, the constants are C stored and passed in the common block "const". The original version C of this code was based on the code by dickson in Version 2 of C Handbook of Methods C for the Analysis of the Various Parameters of C the Carbon Dioxide System in Seawater , DOE, 1994 (SOP No. 3, p25-26). C.................................................................... c determine pressure (bar) from depth c 1 BAR at z=0m (atmos pressure) c use UPPER surface of cell so top layer pressure = 0 bar c for surface exchange coeffs cmick.............................. c write(6,*)'Klevel ',klevel bdepth = 0.0d0 cdepth = 0.0d0 pressc = 1.0d0 do k = 1,Klevel cdepth = bdepth + 0.5d0*drF(k) bdepth = bdepth + drF(k) pressc = 1.0d0 + 0.1d0*cdepth end


do do i=imin,imax do j=jmin,jmax if (hFacC(i,j,Klevel,bi,bj).gt.0.d0) then t = ttemp(i,j) s = stemp(i,j) C terms used more than once tk = 273.15 + t tk100 = tk/100.0 tk1002=tk100*tk100 invtk=1.0/tk dlogtk=log(tk) is=19.924*s/(1000.-1.005*s) is2=is*is sqrtis=sqrt(is) s2=s*s sqrts=sqrt(s) s15=s**1.5 scl=s/1.80655 C----------------------------------------------------------------------- C Ariane adding this - copied from non-pressure-dep case C added by Val Bennington Nov 2010 C Fugacity Factor needed for non-ideality in ocean C ff used for atmospheric correction for water vapor and pressure C Weiss (1974) Marine Chemistry P1atm = 1.01325 _d 0 ! bars Rgas = 83.1451 _d 0 ! bar*cm3/(mol*K) RT = Rgas*tk delta = (57.7 _d 0 - 0.118 _d 0*tk) B1 = -1636.75 _d 0 + 12.0408 _d 0*tk - 0.0327957 _d 0*tk*tk B = B1 + 3.16528 _d 0*tk*tk*tk*(0.00001 _d 0) CMM( #ifdef USE_EXF_ATMPRES CMM apresure in N/m2. move to bars fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * & (apressure (i,j,bi,bj) * 0.00001) & / RT) #else fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * P1atm / RT) #endif CMM) C------------------------------------------------------------------------ C f = k0(1-pH2O)*correction term for non-ideality C Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values) ff(i,j,bi,bj) = exp(-162.8301 + 218.2968/tk100 + & 90.9241*log(tk100) - 1.47696*tk1002 + & s * (.025695 - .025225*tk100 + & 0.0049867*tk1002)) C------------------------------------------------------------------------ C K0 from Weiss 1974 ak0(i,j,bi,bj) = exp(93.4517/tk100 - 60.2409 + & 23.3585 * log(tk100) + & s * (0.023517 - 0.023656*tk100 + & 0.0047036*tk1002)) C------------------------------------------------------------------------ C k1 = [H][HCO3]/[H2CO3] C k2 = [H][CO3]/[HCO3] C Millero p.664 (1995) using Mehrbach et al. data on seawater scale ak1(i,j,bi,bj)=10**(-1*(3670.7*invtk - & 62.008 + 9.7944*dlogtk - & 0.0118 * s + 0.000116*s2)) ak2(i,j,bi,bj)=10**(-1*(1394.7*invtk + 4.777 - & 0.0184*s + 0.000118*s2)) C NOW PRESSURE DEPENDENCE: c Following Takahashi (1981) GEOSECS report - quoting Culberson and c Pytkowicz (1968) c pressc = pressure in bars ak1(i,j,bi,bj) = ak1(i,j,bi,bj)* & exp( (24.2-0.085*t)*(pressc-1.0)/(83.143*tk) ) c FIRST GO FOR K2: According to GEOSECS (1982) report c ak2(i,j,bi,bj) = ak2(i,j,bi,bj)* c & exp( (26.4-0.040*t)*(pressc-1.0)/(83.143*tk) ) c SECOND GO FOR K2: corrected coeff according to CO2sys documentation c E. Lewis and D. Wallace (1998) ORNL/CDIAC-105 ak2(i,j,bi,bj) = ak2(i,j,bi,bj)* & exp( (16.4-0.040*t)*(pressc-1.0)/(83.143*tk) ) C------------------------------------------------------------------------ C kb = [H][BO2]/[HBO2] C Millero p.669 (1995) using data from dickson (1990) akb(i,j,bi,bj)=exp((-8966.90 - 2890.53*sqrts - 77.942*s + & 1.728*s15 - 0.0996*s2)*invtk + & (148.0248 + 137.1942*sqrts + 1.62142*s) + & (-24.4344 - 25.085*sqrts - 0.2474*s) * & dlogtk + 0.053105*sqrts*tk) C Mick and Karsten - Dec 04 C ADDING pressure dependence based on Millero (1995), p675 C with additional info from CO2sys documentation (E. Lewis and C D. Wallace, 1998 - see endnotes for commentary on Millero, 95) bigR = 83.145 dv = -29.48 + 0.1622*t + 2.608d-3*t*t dk = -2.84d-3 pfactor = - (dv/(bigR*tk))*pressc & + (0.5*dk/(bigR*tk))*pressc*pressc akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor) C------------------------------------------------------------------------ C k1p = [H][H2PO4]/[H3PO4] C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974) ak1p(i,j,bi,bj) = exp(-4576.752*invtk + 115.525 - & 18.453*dlogtk + & (-106.736*invtk + 0.69171)*sqrts + & (-0.65643*invtk - 0.01844)*s) C------------------------------------------------------------------------ C k2p = [H][HPO4]/[H2PO4] C DOE(1994) eq 7.2.23 with footnote using data from Millero (1974)) ak2p(i,j,bi,bj) = exp(-8814.715*invtk + 172.0883 - & 27.927*dlogtk + & (-160.340*invtk + 1.3566) * sqrts + & (0.37335*invtk - 0.05778) * s) C------------------------------------------------------------------------ C k3p = [H][PO4]/[HPO4] C DOE(1994) eq 7.2.26 with footnote using data from Millero (1974) ak3p(i,j,bi,bj) = exp(-3070.75*invtk - 18.141 + & (17.27039*invtk + 2.81197) * & sqrts + (-44.99486*invtk - 0.09984) * s) C------------------------------------------------------------------------ C ksi = [H][SiO(OH)3]/[Si(OH)4] C Millero p.671 (1995) using data from Yao and Millero (1995) aksi(i,j,bi,bj) = exp(-8904.2*invtk + 117.385 - & 19.334*dlogtk + & (-458.79*invtk + 3.5913) * sqrtis + & (188.74*invtk - 1.5998) * is + & (-12.1652*invtk + 0.07871) * is2 + & log(1.0-0.001005*s)) C------------------------------------------------------------------------ C kw = [H][OH] C Millero p.670 (1995) using composite data akw(i,j,bi,bj) = exp(-13847.26*invtk + 148.9652 - & 23.6521*dlogtk + & (118.67*invtk - 5.977 + 1.0495 * dlogtk) * & sqrts - 0.01615 * s) C------------------------------------------------------------------------ C ks = [H][SO4]/[HSO4] C dickson (1990, J. chem. Thermodynamics 22, 113) aks(i,j,bi,bj)=exp(-4276.1*invtk + 141.328 - & 23.093*dlogtk + & (-13856*invtk + 324.57 - 47.986*dlogtk)*sqrtis + & (35474*invtk - 771.54 + 114.723*dlogtk)*is - & 2698*invtk*is**1.5 + 1776*invtk*is2 + & log(1.0 - 0.001005*s)) C------------------------------------------------------------------------ C kf = [H][F]/[HF] C dickson and Riley (1979) -- change pH scale to total akf(i,j,bi,bj)=exp(1590.2*invtk - 12.641 + 1.525*sqrtis + & log(1.0 - 0.001005*s) + & log(1.0 + (0.1400/96.062)*(scl)/aks(i,j,bi,bj))) C------------------------------------------------------------------------ C Calculate concentrations for borate, sulfate, and fluoride C Uppstrom (1974) bt(i,j,bi,bj) = 0.000232 * scl/10.811 C Morris & Riley (1966) st(i,j,bi,bj) = 0.14 * scl/96.062 C Riley (1965) ft(i,j,bi,bj) = 0.000067 * scl/18.9984 C------------------------------------------------------------------------ C solubility product for calcite caxx and aragonite C c Following Takahashi (1982) GEOSECS handbook C NOT SURE THIS IS WORKING??? C Ingle et al. (1973) c Ksp_T_Calc = ( -34.452 - 39.866*(s**0.333333) c & + 110.21*log(s) - 7.5752d-6 * (tk**2.0) c & ) * 1.0d-7 c with pressure dependence Culberson and Pytkowicz (1968) c xvalue = (36-0.20*t)*(pressc-1.0)/(83.143*tk) c Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*exp(xvalue) c c C Following Mucci (1983) - from Zeebe/Wolf-Gladrow equic.m tmpa1 = - 171.9065 - (0.077993*tk) + (2839.319/tk) & + (71.595*log10(tk)) tmpa2 = +(-0.77712 + (0.0028426*tk) + (178.34/tk) )*sqrts tmpa3 = -(0.07711*s) + (0.0041249*s15) logKspc = tmpa1 + tmpa2 + tmpa3 Ksp_T_Calc = 10.0**logKspc C tmpa1 = - 171.945 - (0.077993*tk) + (2903.293/tk) & + (71.595*log10(tk)) tmpa2 = +(-0.068393 + (0.0017276*tk) + (88.135/tk) )*sqrts tmpa3 = -(0.10018*s) + (0.0059415*s15) logKspa = tmpa1 + tmpa2 + tmpa3 Ksp_T_Arag = 10.0**logKspa c with pressure dependence Culberson and Pytkowicz (1968) c xvalue = (36.0-0.20*t)*(pressc-1.0)/(83.143*tk) c Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*exp(xvalue) c alternative pressure depdendence c following Millero (1995) but using info from Appendix A11 of c Zeebe and Wolf-Gladrow (2001) book c dv = -48.6 - 0.5304*t c dk = -11.76d-3 - 0.3692*t c xvalue = - (dv/(bigR*tk))*pressc c & + (0.5*dk/(bigR*tk))*pressc*pressc c Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*exp(xvalue) c alternative pressure dependence from Ingle (1975) zdum = (pressc*10.0d0 - 10.0d0)/10.0d0 xvalue = ( (48.8d0 - 0.53d0*t)*zdum & + (-0.00588d0 + 0.0001845d0*t)*zdum*zdum) & / (188.93d0*(t + 273.15d0)) Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*10**(xvalue) Ksp_TP_Arag(i,j,bi,bj) = Ksp_T_Arag*10**(xvalue) C------------------------------------------------------------------------ else fugf(i,j,bi,bj)=0.d0 ff(i,j,bi,bj)=0.d0 ak0(i,j,bi,bj)= 0.d0 ak1(i,j,bi,bj)= 0.d0 ak2(i,j,bi,bj)= 0.d0 akb(i,j,bi,bj)= 0.d0 ak1p(i,j,bi,bj) = 0.d0 ak2p(i,j,bi,bj) = 0.d0 ak3p(i,j,bi,bj) = 0.d0 aksi(i,j,bi,bj) = 0.d0 akw(i,j,bi,bj) = 0.d0 aks(i,j,bi,bj)= 0.d0 akf(i,j,bi,bj)= 0.d0 bt(i,j,bi,bj) = 0.d0 st(i,j,bi,bj) = 0.d0 ft(i,j,bi,bj) = 0.d0 Ksp_TP_Calc(i,j,bi,bj) = 0.d0 Ksp_TP_Arag(i,j,bi,bj) = 0.d0 endif end


do end


do return end