#include "DIC_OPTIONS.h" C-- File carbon_chem.F: C-- Contents C-- o CALC_PCO2 C-- o CALC_PCO2_APPROX C-- o CARBON_COEFFS C-- o CARBON_COEFFS_PRESSURE_DEP C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: CALC_PCO2 C !INTERFACE: ========================================================== 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 !DESCRIPTION: C surface ocean inorganic carbon chemistry to OCMIP2 C regulations modified from OCMIP2 code; C Mick Follows, MIT, Oct 1999. C Apr 2011: fix vapour bug (following Bennington) C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "DYNVARS.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "DIC_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 --------------------------------------------------------------- #ifdef WATERVAP_BUG pCO2surfloc = co2star / fflocal #else c Corrected by Val Bennington (Nov 2010) fco2 = co2star / k0local pCO2surfloc = fco2/fugflocal #endif 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 C !ROUTINE: CALC_PCO2_APPROX C !INTERFACE: ========================================================== 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 C !USES: =============================================================== IMPLICIT NONE C == GLobal variables == #include "SIZE.h" #include "DYNVARS.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "DIC_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) 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...) #ifdef WATERVAP_BUG pCO2surfloc = co2s/fflocal #else c bug fix by Bennington fco2 = co2s/k0local pco2surfloc = fco2/fugflocal #endif 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 C !ROUTINE: CARBON_COEFFS C !INTERFACE: ========================================================== 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-------------------------------------------------------------------------- C !USES: =============================================================== IMPLICIT NONE C == GLobal variables == #include "SIZE.h" #include "DYNVARS.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "DIC_VARS.h" 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 #ifdef CARBONCHEM_TOTALPHSCALE c pH scale converstions _RL total2free _RL free2sw _RL total2sw #endif 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 for: C temperature tk = 273.15 _d 0 + t tk100 = tk/100. _d 0 tk1002=tk100*tk100 invtk=1.0 _d 0/tk dlogtk=log(tk) C ionic strength is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s) is2=is*is sqrtis=sqrt(is) C salinity 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) fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) * P1atm / RT) 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 C Independent of pH scale 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 C Both on seawater pH 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) C On total pH scale 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) C The constant 115.525 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 115.540 constant to stay on the seawater C pH scale (everything else is the same). 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)) C The constant 172.0833 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 172.1033 constant to stay on the seawater C pH scale (everything else is the same). 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) C The constant 18.141 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 18.126 constant to stay on the seawater C pH scale (everything else is the same). 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) C The constant 117.385 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 117.400 constant to stay on the seawater C pH scale (everything else is the same). C Note: converted here from mol/kg-H2O to mol/kg-SW 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 C The constant 148.9652 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 148.9802 constant to stay on the seawater C pH scale (everything else is the same). 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) C On the free pH scale C Note: converted here from mol/kg-H2O to mol/kg-SW 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 (following Dickson & Goyet, 1994) C Note: converted here from mol/kg-H2O to mol/kg-SW 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------------------------------------------------------------------------ #ifdef CARBONCHEM_TOTALPHSCALE C Convert between pH scales, we want everything to end up on the total scale C ak1 and ak2 are on seawater scale and aks is on the free scale total2free = 1. _d 0/ & (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj)) free2sw = 1. _d 0 & + st(i,j,bi,bj)/ aks(i,j,bi,bj) & + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free) total2sw = total2free * free2sw ak1(i,j,bi,bj) = ak1(i,j,bi,bj)/total2sw ak2(i,j,bi,bj) = ak2(i,j,bi,bj)/total2sw aks(i,j,bi,bj) = aks(i,j,bi,bj)/total2free #endif 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 C !ROUTINE: CARBON_COEFFS_PRESSURE_DEP C !INTERFACE: ========================================================== 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-------------------------------------------------------------------------- C !USES: =============================================================== IMPLICIT NONE C == GLobal variables == #include "SIZE.h" #include "DYNVARS.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "FFIELDS.h" #include "DIC_VARS.h" 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 INTEGER i INTEGER j INTEGER k _RL bdepth _RL cdepth _RL pressc _RL Ksp_T_Calc _RL xvalue _RL zdum _RL tmpa1 _RL tmpa2 _RL tmpa3 _RL logKspc _RL dv _RL dk _RL pfactor _RL bigR #ifdef CARBONCHEM_TOTALPHSCALE c pH scale converstions _RL total2free_surf _RL free2sw_surf _RL total2sw_surf _RL total2free _RL free2sw _RL total2sw #endif 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. _d 0 cdepth = 0. _d 0 pressc = 1. _d 0 do k = 1,Klevel cdepth = bdepth + 0.5 _d 0*drF(k) bdepth = bdepth + drF(k) pressc = 1. _d 0 + 0.1 _d 0*cdepth end do cmick................................................... c write(6,*)'depth,pressc ',cdepth,pressc cmick.................................................... do i=imin,imax do j=jmin,jmax if (hFacC(i,j,Klevel,bi,bj).gt.0. _d 0) then t = ttemp(i,j) s = stemp(i,j) C terms used more than once for: C temperature tk = 273.15 _d 0 + t tk100 = tk/100. _d 0 tk1002=tk100*tk100 invtk=1.0 _d 0/tk dlogtk=log(tk) C ionic strength is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s) is2=is*is sqrtis=sqrt(is) C salinity s2=s*s sqrts=sqrt(s) s15=s**1.5 _d 0 scl=s/1.80655 _d 0 bigR = 83.14472 _d 0 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 C Independent of pH scale 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 C Both on seawater pH 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*(1394.7 _d 0*invtk + 4.777 _d 0 - & 0.0184 _d 0*s + 0.000118 _d 0*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 _d 0-0.085 _d 0*t)*(pressc-1.0 _d 0)/(83.143 _d 0*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 _d 0-0.040 _d 0*t)*(pressc-1.0 _d 0)/(83.143 _d 0*tk) ) C------------------------------------------------------------------------ C kb = [H][BO2]/[HBO2] C Millero p.669 (1995) using data from dickson (1990) C On total pH scale 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) #ifdef CARBONCHEM_TOTALPHSCALE C JML Not yet, need to account for pH scale conversions C (pressure correct on the seawater scale, not the total scale) 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) C bigR = 83.145 C dv = -29.48 + 0.1622*t + 2.608d-3*t*t C dk = -2.84d-3 C pfactor = - (dv/(bigR*tk))*pressc C & + (0.5*dk/(bigR*tk))*pressc*pressc C akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor) #else 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 _d 0 dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t dk = -2.84 _d -3 pfactor = - (dv/(bigR*tk))*pressc & + (0.5 _d 0*dk/(bigR*tk))*pressc*pressc akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor) #endif C------------------------------------------------------------------------ C k1p = [H][H2PO4]/[H3PO4] C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974) C The constant 115.525 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 115.540 constant to stay on the seawater C pH scale (everything else is the same). 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) C The constant 172.0833 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 172.1033 constant to stay on the seawater C pH scale (everything else is the same). ak2p(i,j,bi,bj) = exp(-8814.715 _d 0*invtk + 172.0883 _d 0 - & 27.927 _d 0*dlogtk + & (-160.34 _d 00*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) C The constant 18.141 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 18.126 constant to stay on the seawater C pH scale (everything else is the same). 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) C The constant 117.385 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 117.400 constant to stay on the seawater C pH scale (everything else is the same). C Note: converted here from mol/kg-H2O to mol/kg-SW 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 C The constant 148.9652 is an approximation to convert to the total pH scale C (Dickson et al., 2007). Use Millero's 148.9802 constant to stay on the seawater C pH scale (everything else is the same). 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) C On the free pH scale C Note: converted here from mol/kg-H2O to mol/kg-SW 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 (following Dickson & Goyet, 1994) C Note: converted here from mol/kg-H2O to mol/kg-SW 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------------------------------------------------------------------------ #ifdef CARBONCHEM_TOTALPHSCALE C JML Convert between pH scales, we want everything to end up on the total scale C ak1 and ak2 are on seawater scale and are pressure corrected, so they just need converting. C aks is on the free scale, it needs pressure correcting then converting to the total scale. C akf needs pressure correcting too, but is on the right (total) pH scale. C akb, ak1p, ak2p, ak3p, aksi and akw are on the total scale. Convert to the seawater C scale before pressure correction, and then convert back to the total scale. C All the terms for pressure correction are from Table 9, Millero (1995), p675 C Info about pH conversions from Munhoven (2013) and Orr and Epitalon (2015) total2free_surf = 1. _d 0/ & (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj)) free2sw_surf = 1. _d 0 & + st(i,j,bi,bj)/ aks(i,j,bi,bj) & + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free_surf) total2sw_surf = total2free_surf * free2sw_surf C------------------------------------------------------------------------ C Pressure correct aks on native free pH scale dv=-18.03 _d 0 + 0.0466 _d 0*t + 0.316 _d -3 *t*t dk=-4.53 _d -3 + 0.09 _d -3 *t pfactor = -(dv/(bigR*tk))*pressc & +(0.5*dk/(bigR*tk))*pressc*pressc aks(i,j,bi,bj) = aks(i,j,bi,bj)*exp(pfactor) total2free = 1. _d 0/ & (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj)) C free2sw has an additional component from fluoride free2sw = 1. _d 0 & + st(i,j,bi,bj)/ aks(i,j,bi,bj) aks(i,j,bi,bj) = aks(i,j,bi,bj)/total2free C------------------------------------------------------------------------ C Pressure correct akf on the free scale before converting back to total akf(i,j,bi,bj) = akf(i,j,bi,bj) * total2free_surf dv = -9.78 _d 0 -0.0090 _d 0 *t -0.942 _d -3 *t*t dk = -3.91 _d -3 + 0.054 _d -3 *t pfactor = -(dv/(bigR*tk))*pressc & +(0.5*dk/(bigR*tk))*pressc*pressc akf(i,j,bi,bj) = akf(i,j,bi,bj) * exp(pfactor) akf(i,j,bi,bj) = akf(i,j,bi,bj)/total2free free2sw = free2sw & + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free) total2sw = total2free * free2sw C------------------------------------------------------------------------ C Convert pressure corrected ak1 and ak2 from the seawater pH scale to the total scale ak1(i,j,bi,bj) = ak1(i,j,bi,bj)/total2sw ak2(i,j,bi,bj) = ak1(i,j,bi,bj)/total2sw C------------------------------------------------------------------------ C Pressure correct akb on the seawater scale before converting back to total scale dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t dk = -2.84 _d -3 pfactor = - (dv/(bigR*tk))*pressc & + (0.5*dk/(bigR*tk))*pressc*pressc akb(i,j,bi,bj) = total2sw_surf*akb(i,j,bi,bj)*exp(pfactor) akb(i,j,bi,bj) = akb(i,j,bi,bj)/total2sw C------------------------------------------------------------------------ C Pressure correct ak1p on the seawater scale before converting back to total scale dv = -14.51 _d 0 + 0.1211 _d 0*t -0.321 _d -3*t*t dk = -2.67 _d -3 + 0.0427 _d -3*t pfactor = - (dv/(bigR*tk))*pressc & + (0.5*dk/(bigR*tk))*pressc*pressc ak1p(i,j,bi,bj) = total2sw_surf*ak1p(i,j,bi,bj)*exp(pfactor) ak1p(i,j,bi,bj) = ak1p(i,j,bi,bj)/total2sw C------------------------------------------------------------------------ C Pressure correct ak2p on the seawater scale before converting back to total scale dv = -23.12 _d 0 + 0.1758 _d 0*t -2.647 _d -3*t*t dk = -5.15 _d -3 + 0.09 _d -3*t pfactor = - (dv/(bigR*tk))*pressc & + (0.5*dk/(bigR*tk))*pressc*pressc ak2p(i,j,bi,bj) = total2sw_surf*ak2p(i,j,bi,bj)*exp(pfactor) ak2p(i,j,bi,bj) = ak2p(i,j,bi,bj)/total2sw C------------------------------------------------------------------------ C Pressure correct ak3p on the seawater scale before converting back to total scale dv = -26.57 _d 0 + 0.2020 _d 0*t -3.042 _d -3*t*t dk = -4.08 _d -3 + 0.0714 _d -3*t pfactor = - (dv/(bigR*tk))*pressc & + (0.5*dk/(bigR*tk))*pressc*pressc ak3p(i,j,bi,bj) = total2sw_surf*ak3p(i,j,bi,bj)*exp(pfactor) ak3p(i,j,bi,bj) = ak3p(i,j,bi,bj)/total2sw C------------------------------------------------------------------------ C Pressure correct akw on the seawater scale before converting back to total scale dv = -20.02 _d 0 + 0.1119 _d 0*t -1.409 _d -3*t*t dk = -5.13 _d -3 + 0.0794 _d -3 *t pfactor = - (dv/(bigR*tk))*pressc & + (0.5*dk/(bigR*tk))*pressc*pressc akw(i,j,bi,bj) = total2sw_surf*akw(i,j,bi,bj)*exp(pfactor) akw(i,j,bi,bj) = akw(i,j,bi,bj)/total2sw C------------------------------------------------------------------------ C Pressure correct aksi on the seawater scale before converting back to total scale C Estimated from the effect on akb dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t dk = -2.84 _d -3 pfactor = - (dv/(bigR*tk))*pressc & + (0.5*dk/(bigR*tk))*pressc*pressc aksi(i,j,bi,bj) = total2sw_surf*aksi(i,j,bi,bj)*exp(pfactor) aksi(i,j,bi,bj) = aksi(i,j,bi,bj)/total2sw #endif C------------------------------------------------------------------------ C Solubility product for calcite 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 _d 0 - (0.077993 _d 0*tk) & + (2839.319 _d 0/tk) + (71.595 _d 0*log10(tk)) tmpa2 = +(-0.77712 _d 0 + (0.0028426 _d 0*tk) & + (178.34 _d 0/tk) )*sqrts tmpa3 = -(0.07711 _d 0*s) + (0.0041249 _d 0*s15) logKspc = tmpa1 + tmpa2 + tmpa3 Ksp_T_Calc = 10.0 _d 0**logKspc c write(6,*)i,j,k,tmpa1,tmpa2,tmpa3,logkspc,Ksp_T_Calc 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. _d 0 - 10. _d 0)/10. _d 0 xvalue = ( (48.8 _d 0 - 0.53 _d 0*t)*zdum & + (-0.00588 _d 0 + 0.0001845 _d 0*t)*zdum*zdum) & / (188.93 _d 0*(t + 273.15 _d 0)) Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*10**(xvalue) C------------------------------------------------------------------------ else 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 Ksp_TP_Calc(i,j,bi,bj) = 0. _d 0 endif end do end do return end