! source file: /gxfs_work/geomar/smomw113/UVic_ESCM/2.9_dk_0814/updates/11_valley/source/common/co2calc_TAfCO2.F subroutine co2calc_TAfCO2 (t, s, ta_in, xco2_in, pt_in &, sit_in, atmpres, depth, phlo, phhi &, debug,dicsat) !------------------------------------------------------------------------- cwk PURPOSE: Calculate saturated DIC from TA and atmospheric pCO2, and local t, s, etc. cwk experimental code, wk, 29.07.20 ! Modified from co2calc.f (RCS version 1.8, OCMIP-2) ! Constants are given on "seawater" H scale (hSWS) except for the "free" ! H scale dissociation constants for S and F (necessary to work with hSWS). ! Modified to compute DIC from TA and fCO2 with the purpose to allow for the ! calculation of saturated DIC. ! INPUT: ! t = temperature (degrees C) ! s = salinity (PSU) ! ta_in = total alkalinity (eq/m^3) ! xco2_in = atmospheric CO2 concentration (ppmv) ! pt_in = inorganic phosphate (mol/m^3) ! sit_in = inorganic silicate (mol/m^3) ! atmpres = atmospheric pressure in atmospheres (1 atm==1013.25mbar) ! depth = ocean depth (m) ! phlo = lower limit of pH range ! phhi = upper limit of pH range ! OUTPUT: ! dicsat = total saturated inorganic carbon (mol/m^3) cwk Note; dic is actually dicsat, when returned from this routine cwk the following are also not needed to be returned ! pH = ocean pH (-log10(hSWS)) ! co2star = ocean CO2* (mol/m^3) ! dco2star = delta (atm-ocn) CO2* (mol/m^3) ! hsws = H+ on seawater H scale cwk the following three shall not be returned from _TAfCO2() ! CO3 = carbonate ion concentration ! Omega_c = omega calcite ! Omega_a = omega aragonite cwk ? do I want to keep memory of this one, no I do not think so ! pco2 = oceanic pCO2 (uatm) cwk I have adopted the first lines of this subroutine to not return the latter variables !-------------------------------------------------------------------------- implicit none real permil, pt, pt_in, sit, sit_in, ta, ta_in, dic, di real permeg, xco2, xco2_in, tk, t, tk100, tk1002, dlogtk, s real sqrtis, s2, sqrts, s15, scl, bt, st, ft, ff, x1, phhi real x2, phlo, xacc, hSWS, drtsafe2, hSWS2, co2star, co2starair real atmpres, dco2star, ph, pco2, invtk, is, is2, k0, k1 real k12, k2, kw, kb, ks, kf, k1p, k2p, k3p, ksi, C2K real Kspc, Kspa, ca, co3, omega_c, omega_a, depth, pres real t2, dvc, dva, dk, pitkR, p2itkR, Rgas, fugf cwk real fco2, dicsat0, dicsat real balk, oh, phostop, phosbot, palk, sialk &, freetotot, hfree, hso4, hf, calk logical debug !-------------------------------------------------------------------------- ! Models carry tracers in mol/m^3 but this routine uses umol/kg ! Change units from the input of mol/m^3 -> mol/kg: ! (1 mol/m^3) x (1 m^3/1024.5 kg) ! where the ocean's mean surface density is 1024.5 kg/m^3 permil = 1./1024.5 ! to convert mol/m^3 to mol/kg permeg = 1.e-6 ! to convert uatm to atm. Rgas = 83.1451 ! ideal gas constant R in cm3 bar k-1 C2K = 273.15 ! to cinvert C to K ! Convert input ta = ta_in*permil ! change mol/m^3 to mol/kg xco2 = xco2_in*permeg ! change units from ppmv to parts by volume. pt = pt_in*permil ! change mol/m^3 to mol/kg sit = sit_in*permil ! change mol/m^3 to mol/kg pres = depth*0.1 ! change from meters to decibars ! later I need to do the respective back-conversion! ! dic = dic_in*permil ! change mol/m^3 to mol/kg cwk / unchanged code from here until ..... !-------------------------------------------------------------------------- ! Calculate all constants needed to convert between various measured ! carbon species. References for each equation are noted in the code. ! Version 2 of "Handbook of Methods for the Analysis of the Various ! Parameters of the Carbon Dioxide System in Seawater", DOE, 1994 ! (SOP No. 3, p25-26). ! Derive simple terms used more than once tk = C2K + t tk100 = tk/100. tk1002 = tk100*tk100 invtk = 1./tk dlogtk = log(tk) is = 19.924*s/(1000. - 1.005*s) is2 = is*is sqrtis = sqrt(is) s2 = s*s t2 = t*t sqrts = sqrt(s) s15 = s**1.5 scl = s/1.80655 pitkR = pres/tk/Rgas p2itkR = pres*pitkR !------------------------------------------------------------------------ ! Calculate concentrations for borate, sulfate, and fluoride ! Uppstrom (1974) bt = 0.000232*scl/10.811 ! Morris & Riley (1966) st = 0.14*scl/96.062 ! Riley (1965) ft = 0.000067*scl/18.9984 !------------------------------------------------------------------------ ! f = k0(1-pH2O)*correction term for non-ideality ! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values) ff = exp(-162.8301 + 218.2968/tk100 + 90.9241*log(tk100) & - 1.47696*tk1002 + s*(.025695 - .025225*tk100 & + 0.0049867*tk1002)) !------------------------------------------------------------------------ ! k1 = [H][HCO3]/[H2CO3] ! k2 = [H][CO3]/[HCO3] on hSWS ! Millero p.664 (1995) using Mehrbach et al. data on SEAWATER scale ! (Original reference: Dickson and Millero, DSR, 1987) k1 = 10**(-1.*(3670.7*invtk - 62.008 + 9.7944*dlogtk & - 0.0118*s + 0.000116*s2)) k2 = 10**(-1*(1394.7*invtk + 4.777 - 0.0184*s + 0.000118*s2)) !------------------------------------------------------------------------ ! k1p = [H][H2PO4]/[H3PO4] on hSWS ! Millero p.670 (1995) k1p = exp(-4576.752*invtk + 115.540 - 18.453*dlogtk & + (-106.736*invtk + 0.69171)*sqrts & + (-0.65643*invtk - 0.01844)*s) !------------------------------------------------------------------------ ! k2p = [H][HPO4]/[H2PO4] on hSWS ! Millero p.670 (1995) k2p = exp(-8814.715*invtk + 172.1033 - 27.927*dlogtk & + (-160.340*invtk + 1.3566)*sqrts & + (0.37335*invtk - 0.05778)*s) !------------------------------------------------------------------------ ! k3p = [H][PO4]/[HPO4] on hSWS ! Millero p.670 (1995) k3p = exp(-3070.75*invtk - 18.126 & + (17.27039*invtk + 2.81197)*sqrts & + (-44.99486*invtk - 0.09984)*s) !------------------------------------------------------------------------ ! ksi = [H][SiO(OH)3]/[Si(OH)4] on hSWS ! Millero p.671 (1995) using data from Yao and Millero (1995) ! change to (mol/ kg soln) ! depth dependancy assumed to be the same as boric acid ! typo in Millero 1994 corrected in sign of 0.1622 ksi = exp(-8904.2*invtk + 117.400 - 19.334*dlogtk & + (-458.79*invtk + 3.5913)*sqrtis & + (188.74*invtk - 1.5998)*is & + (-12.1652*invtk + 0.07871)*is2 & + log(1. - 0.001005*s)) !------------------------------------------------------------------------ ! kw = [H][OH] on hSWS ! Millero p.670 (1995) using composite data ! pressure dependancy in Millero 1994 corrected for sea water from ! Millero 1983 kw = exp(-13847.26*invtk + 148.9802 - 23.6521*dlogtk & + (118.67*invtk - 5.977 + 1.0495*dlogtk)*sqrts - 0.01615*s) !------------------------------------------------------------------------ ! ks = [H][SO4]/[HSO4] on free H scale ! Dickson (1990, J. chem. Thermodynamics 22, 113) ! change to (mol/kg soln) ks = 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.001005*s)) !------------------------------------------------------------------------ ! kf = [H][F]/[HF] on free H scale ! Dickson and Riley (1979) ! change to (mol/ kg soln) kf = exp(1590.2*invtk - 12.641 + 1.525*sqrtis & + log(1. - 0.001005*s)) !------------------------------------------------------------------------ ! kb = [H][BO2]/[HBO2] on hSWS ! Dickson p.673 (1990) ! change from htotal to hSWS ! typo in Millero 1994 corrected in sign of 0.1622 kb = 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 & + log((1+(st/ks)+(ft/kf))/(1+(st/ks)))) cwk I moved this up from further below in co2calc_SWS.F ! K0 and fugf from Weiss 1974 Marine Chemistry k0 = exp(93.4517/tk100 - 60.2409 + 23.3585*log(tk100) & + s*(.023517 - 0.023656*tk100 + 0.0047036*tk1002)) ! Fugacity factor (assumes pressure of 1 atm or 1.01325 bars) fugf = exp((-1636.75 + 12.0408*tk - 0.0327957*tk**2 & + 3.16528e-5*tk**3 + 2.*(57.7 - 0.118*tk))*1.01325/(Rgas*tk)) !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2, ! ORNL/CDIAC-74, Dickson and Goyet, eds. (Ch 2 p 10, Eq A.49) ! hSWS2 = hSWS*hSWS ! co2star = (dic*hSWS2/(hSWS2 + k1*hSWS + k1*k2)) co2starair = xco2*ff*atmpres ! dco2star = co2starair - co2star ! Calculate the partial pressures ! pCO2 = co2starair/(k0*fugf)/permeg ! convert atm to uatm cwk compute fco2; fco2 = pco2*fugfac ! fco2 = co2starair /k0/permeg cwk is equal to ! fco2 = xco2*ff*atmpres/k0 fco2 = xco2*ff*atmpres/k0 if (debug) then write(20,*) 'xco2: ', xco2 write(20,*) 'ff: ', ff write(20,*) 'fugf: ', fugf write(20,*) 'atmpres',atmpres write(20,*) 'k0: ',k0 write(20,*) 'permeg ',permeg endif !------------------------------------------------------------------------- ! Calculate [H+] SWS when fCO2 and TA are known at T, S and 1 atm. x1 = 10.**(-phhi) x2 = 10.**(-phlo) xacc = 1.e-10 if (debug) then write (20,*) & 'Input to drtsafe2:' write (20,*) & 'ta (mol/kg): ', ta write (20,*) 'fco2 (atm): ', fco2 write (20,*) 'implicit input:' write (20,*) 'SST (degC): ', t write (20,*) 'SSS (): ', s write (20,*) 'PO4 (mol/kg) ',pt endif hSWS = drtsafe2 (k0,k1, k2, kw, kb, ks, kf, k1p, k2p, k3p &, ksi, bt, st &, ft, sit, pt, fco2, ta, x1, x2, xacc, debug) if (debug) then ph = -log10(hSWS) write(20,*) 'initialy returned pH:', ph endif ! The solution converges to err of xacc. The solution must be within ! the range x1 to x2. ! If DIC and TA are known then either a root finding or iterative method ! must be used to calculate hSWS. In this case we use the Newton-Raphson ! "safe" method taken from "Numerical Recipes" (function "rtsafe.f" with ! error trapping removed). ! As currently set, this procedure iterates about 12 times. The x1 and x2 ! values set below will accomodate ANY oceanographic values. If an initial ! guess of the pH is known, then the number of iterations can be reduced to ! about 5 by narrowing the gap between x1 and x2. It is recommended that ! the first few time steps be run with x1 and x2 set as below. After that, ! set x1 and x2 to the previous value of the pH +/- ~0.5. The current ! setting of xacc will result in co2star accurate to 3 significant figures ! (xx.y). Making xacc bigger will result in faster convergence also, but this ! is not recommended (xacc of 10**-9 drops precision to 2 significant figures). x1 = 10.**(-phhi) x2 = 10.**(-phlo) xacc = 1.e-10 cwk cwk changed call, adding k0 and replace dic by fco2 cwk but why a second time? hSWS = drtsafe2 (k0, k1, k2, kw, kb, ks, kf, k1p, k2p &, k3p, ksi, bt, st &, ft, sit, pt, fco2, ta, x1, x2, xacc, debug) if (debug) then ph = -log10(hSWS) write(20,*) 'pH returned after second call:', ph endif cwk now I need to calculate DICsat from ph cwk adopted from function CalculateTCfromTApH in co2sys ! BAlk = TBF.*KBF./(KBF + H); balk = bt*kb/(kb+hSWS) ! OH = KWF./H; oh = kw/hSWS ! PhosTop = KP1F.*KP2F.*H + 2.*KP1F.*KP2F.*KP3F - H.*H.*H; phostop = k1p*k2p*hSWS + 2*k1p*k2p*k3p - hSWS*hSWS*hSWS ! PhosBot = H.*H.*H + KP1F.*H.*H + KP1F.*KP2F.*H + KP1F.*KP2F.*KP3F; phosbot = hSWS*hSWS*hSWS + k1p*hSWS*hSWS + k1p*k2p*hSWS & + k1p*k2p*k3p ! PAlk = TPF.*PhosTop./PhosBot; palk = pt*phostop/phosbot ! SiAlk = TSiF.*KSiF./(KSiF + H); sialk = sit *ksi/(ksi + hSWS) ! FREEtoTOT = (1 + TSF./KSF); % pH scale conversion factor freetotot = (1. + st/ks) ! Hfree = H./FREEtoTOT; %' for H on the total scale hfree = hSWS/freetotot ! HSO4 = TSF./(1 + KSF./Hfree); %' since KS is on the free scale hso4 = st/(1 + ks/hfree) ! HF = TFF./(1 + KFF./Hfree); %' since KF is on the free scale hf = ft/(1 + kf/hfree) ! CAlk = TAx - BAlk - OH - PAlk - SiAlk + Hfree + HSO4 + HF; calk = ta - balk - oh - palk - sialk + hfree + hso4 + hf ! TCxtemp = CAlk.*(H.*H + K1F.*H + K1F.*K2F)./(K1F.*(H + 2.*K2F)); ! dicsat should be now in mol/kg dicsat0 = calk * (hSWS*hSWS + k1*hSWS + k1*k2) / & (k1 * (hSWS + 2*k2)) ! convert back from mol/kg to mol/m3 dicsat = dicsat0/permil if (debug) then write (20,*) & 'co2calc_TAfCO2.f has found the following ' write (20,*) & ' value of dicsat0 (mol/kg): ',dicsat0 write (20,*) & ' value of dicsat (mol/m3) : ',dicsat write (20,*) & '-------------------------------------' endif cwk the rest of this routine is not needed, I think !------------------------------------------------------------------------- ! Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2, ! ORNL/CDIAC-74, Dickson and Goyet, eds. (Ch 2 p 10, Eq A.49) ! hSWS2 = hSWS*hSWS ! co2star = (dic*hSWS2/(hSWS2 + k1*hSWS + k1*k2)) ! co2starair = xco2*ff*atmpres ! dco2star = co2starair - co2star ! Calculate the partial pressures ! pCO2 = co2star/(k0*fugf)/permeg ! convert atm to uatm cwk compute fco2; fco2 = pco2*fugfac ! fco2 = co2star/k0/permeg !------------------------------------------------------------------------- ! co2star = co2star/permil ! convert from mol/kg -> mol/m^3 ! dco2star = dco2star/permil ! convert from mol/kg -> mol/m^3 return end real function drtsafe2 (k0,k1, k2, kw, kb, ks, kf, k1p, k2p &, k3p, ksi, bt, st, ft, sit, pt, fco2, ta, x1, x2, xacc &, debug_in) !------------------------------------------------------------------------- ! File taken from Numerical Recipes. Modified R.M.Key 4/94 cwk I changed drtsave to drtsave2 cwk I changed all calls to ta_iter_SWS_2 (add k0, change dic into fco2) implicit none integer maxit, j real k1, k2, kw, kb, ks, kf, k1p, k2p, k3p, ksi, bt, st, ft, sit real pt, dic, ta real x1, fl, df, x2, fh, xl, xh, swap, dxold, dx, f, temp, xacc cwk added k0, fco2 real k0, fco2 logical debug_in, debug debug = debug_in ! or debug=.FALSE. maxit = 100 call ta_iter_SWS_2 (k0, k1, k2, kw, kb, ks, kf, k1p, k2p &, k3p, ksi, bt, st, ft, sit, pt, fco2, ta, x1, fl, df) ! ... returns fl for x1 (the lower pH bound if (debug) write (20,*) '------ inside drtsafe2 ######' call ta_iter_SWS_2 (k0, k1, k2, kw, kb, ks, kf, k1p, k2p &, k3p, ksi, bt, st, ft, sit, pt, fco2, ta, x2, fh, df) ! ... returns fh for x2 (the upper pH bound) ! fl and fh are if (debug) write (20,*) 'fl, fh: ', & -log10(fl), -log10(fh) if (fl .lt. 0.) then xl = x1 xh = x2 else xh = x1 xl = x2 swap = fl fl = fh fh = swap endif drtsafe2 = 0.5*(x1 + x2) if (debug) write (20,*) 'first pH-estimate: ', -log10(drtsafe2) dxold = abs(x2 - x1) dx = dxold call ta_iter_SWS_2 (k0, k1, k2, kw, kb, ks, kf, k1p, k2p &, k3p, ksi, bt, st, ft, sit, pt, fco2, ta, drtsafe2, f, df) do j=1,maxit if (((drtsafe2 - xh)*df - f)*((drtsafe2 - xl)*df - f) .ge. 0. & .or. abs(2.*f) .gt. abs(dxold*df)) then dxold = dx dx = 0.5*(xh - xl) drtsafe2 = xl + dx if (xl .eq. drtsafe2) return else dxold = dx dx = f/df temp = drtsafe2 drtsafe2 = drtsafe2 - dx if (temp .eq. drtsafe2) return endif if (abs(dx) .lt. xacc) return call ta_iter_SWS_2 (k0, k1, k2, kw, kb, ks, kf, k1p, k2p &, k3p, ksi,bt, st, ft, sit, pt, fco2, ta, drtsafe2, f, df) if (f .lt. 0.) then xl = drtsafe2 fl = f else xh = drtsafe2 fh = f endif if (debug) write (20,*) & 'j, f, df, ph: ',j, f, df,-log10(drtsafe2) enddo if (debug) write (20,*) 'ending drtsafe2 loop ', & 'and returning to co2calc_TAfCO2' return end cwk I added k0 as first parameter and replaced dic by fco2 in the second line subroutine ta_iter_SWS_2 (k0, k1, k2, kw, kb, ks, kf, k1p, k2p &, k3p, ksi, bt, st, ft, sit, pt, fco2, ta, x, fn, df) !------------------------------------------------------------------------- ! This routine expresses TA as a function of DIC, hSWS and constants. ! It also calculates the derivative of this function with respect to ! hSWS. It is used in the iterative solution for hSWS. In the call ! "x" is the input value for hSWS, "fn" is the calculated value for TA ! and "df" is the value for dTA/dhSWS ! Modified from ta_iter_1.f (RCS version 1.2, OCMIP-2) implicit none cwk next block still unchanged real k1, k2, kw, kb, ks, kf, k1p, k2p, k3p, ksi, bt, st, ft, sit real pt, dic, ta, x, fn, df real rx, x2, x3, k12, k12p, k123p, c, rc, a, ra, da, b, rb, db real rkb, rksi, kwrx, ptra, xrc, crx, crx2, t1, t2, t3, t4, t5, t6 cwk real fco2,k0 real hco3,co3,balk,oh logical debug ! the following block is from the original ta_iter_SWS; so far I leave it unchanged rx = 1./x x2 = x*x x3 = x2*x k12 = k1*k2 k12p = k1p*k2p k123p = k12p*k3p c = 1. + st/ks + ft/kf rc = 1./c a = x3 + k1p*x2 + k12p*x + k123p ra = 1./a da = 3.*x2 + 2.*k1p*x + k12p b = x2 + k1*x + k12 rb = 1./b db = 2.*x + k1 rkb = 1./kb rksi = 1./ksi kwrx = kw*rx ptra = pt*ra xrc = x*rc crx = c*rx crx2 = crx*rx cwk ... except for these two ... ! t1 = k1*dic*rb ! t2 = 2.*dic*k12*rb t3 = 1. + x*rkb t4 = 1. + x*rksi t5 = 1. + ks*crx t6 = 1. + kf*crx ! I adopt from CO2SYS: ! HCO3 = K0F.*K1F.*fCO2i./H; hco3 = k0*k1*fco2/x ! CO3 = K0F.*K1F.*K2F.*fCO2i./(H.*H); co3 = k0*k1*k2*fco2/(x*x) ! BAlk = TBF.*KBF./(KBF + H); balk = bt*kb/(kb+x) ! OH = KWF./H; oh = kw/x !------------------------------------------------------------------------- ! fn = hco3+co3+borate+oh+hpo4+2*po4+silicate-hfree-hso4-hf-h3po4-ta cwk Translating fn into co2sys terminology: cwk RESIDUAL = TA - Calk - Balk - OH - PAlk - SiAlk + Hfree + HSO4 + HF cwk fn and RESIDUAL are related by *(-1), hope that that does not cause problems .. cwk the following is th RESIDUAL equation from the original ta_iter_SWS. ! Comparing the co2sys equations CalculatepHfromTATC and CalculatepHfromTAfCO2 ! I find that the residual equations are overall identical ... ! fn = hco3 + 2*co3 + balk + oh + hpo4 + 2.* po4 !! fn = t1*x + t2 + bt/t3 + kwrx + ptra*k12p*x + 2.*ptra*k123p ! + SiAlk - Hfree - HSO4 - HF - ?? - ta !! & + sit/t4 - xrc - st/t5 - ft/t6 - ptra*x3 - ta ! ... but the individual terms are partly calculated differently (since dic is not known) ! i.e hco3, 2*co3 are differnt, while balk and oh only have different names here fn = hco3 + 2*co3 + balk + oh + ptra*k12p*x + 2.*ptra*k123p & + sit/t4 - xrc - st/t5 - ft/t6 - ptra*x3 - ta ! df = dfn/dx cwk the following is the SLOPE dTA/dpH still from the routine ta_iter_SWS ! df = t1 - t1*x*db*rb - t2*db*rb - bt*rkb/t3**2 - kwrx*rx ! & + (ptra*k12p*(a - x*da))*ra - 2.*ptra*k123p*da*ra ! & - sit*rksi/t4**2 - rc - st*t5**(-2)*(ks*crx2) ! & - ft*t6**(-2)*(kf*crx2) - ptra*x2*(3.*a - x*da)*ra cwk I have to replace it by the slope taken from the co2sys function CalculatepHfromTAfCO2(TAi,fCO2i) ! from co2sys: find Slope dTA/dpH ! (this is not exact, but keeps all important terms): ! SLOPE = log(10.)*(HCO3 + 4.*CO3 + BALK.*H./(KBF + H) + OH + H) ! note: 'log' in fortran and matlab is the natural logarithm df = log(10.)*(hco3 + 4*co3 + balk*x/(kb+x) + oh + x) return end