[Seacarb-commits] r2 - in pkg: . seacarb seacarb/R seacarb/data seacarb/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 24 17:00:56 CET 2008
Author: phgrosjean
Date: 2008-12-24 17:00:56 +0100 (Wed, 24 Dec 2008)
New Revision: 2
Added:
pkg/seacarb/
pkg/seacarb/ChangeLog
pkg/seacarb/DESCRIPTION
pkg/seacarb/R/
pkg/seacarb/R/K1.R
pkg/seacarb/R/K1p.R
pkg/seacarb/R/K2.R
pkg/seacarb/R/K2p.R
pkg/seacarb/R/K3p.R
pkg/seacarb/R/Kb.R
pkg/seacarb/R/Kf.R
pkg/seacarb/R/Kh.R
pkg/seacarb/R/Khs.R
pkg/seacarb/R/Kn.R
pkg/seacarb/R/Ks.R
pkg/seacarb/R/Ksi.R
pkg/seacarb/R/Kspa.R
pkg/seacarb/R/Kspc.R
pkg/seacarb/R/Kw.R
pkg/seacarb/R/bjerrum.R
pkg/seacarb/R/bor.R
pkg/seacarb/R/buffer.R
pkg/seacarb/R/carb.R
pkg/seacarb/R/kconv.R
pkg/seacarb/R/pCa.R
pkg/seacarb/R/pHconv.R
pkg/seacarb/R/pTA.R
pkg/seacarb/R/pgas.R
pkg/seacarb/R/phinsi.R
pkg/seacarb/R/pmix.R
pkg/seacarb/R/ppH.R
pkg/seacarb/R/rho.R
pkg/seacarb/R/speciation.R
pkg/seacarb/README
pkg/seacarb/data/
pkg/seacarb/data/seacarb_test.txt
pkg/seacarb/man/
pkg/seacarb/man/K1.Rd
pkg/seacarb/man/K1p.Rd
pkg/seacarb/man/K2.Rd
pkg/seacarb/man/K2p.Rd
pkg/seacarb/man/K3p.Rd
pkg/seacarb/man/Kb.Rd
pkg/seacarb/man/Kf.Rd
pkg/seacarb/man/Kh.Rd
pkg/seacarb/man/Khs.Rd
pkg/seacarb/man/Kn.Rd
pkg/seacarb/man/Ks.Rd
pkg/seacarb/man/Ksi.Rd
pkg/seacarb/man/Kspa.Rd
pkg/seacarb/man/Kspc.Rd
pkg/seacarb/man/Kw.Rd
pkg/seacarb/man/bjerrum.Rd
pkg/seacarb/man/bor.Rd
pkg/seacarb/man/buffer.Rd
pkg/seacarb/man/carb.Rd
pkg/seacarb/man/kconv.Rd
pkg/seacarb/man/pCa.Rd
pkg/seacarb/man/pHconv.Rd
pkg/seacarb/man/pTA.Rd
pkg/seacarb/man/pgas.Rd
pkg/seacarb/man/phinsi.Rd
pkg/seacarb/man/pmix.Rd
pkg/seacarb/man/ppH.Rd
pkg/seacarb/man/rho.Rd
pkg/seacarb/man/seacarb_test.Rd
pkg/seacarb/man/speciation.Rd
Log:
seacarb 2.0.5 added
Added: pkg/seacarb/ChangeLog
===================================================================
--- pkg/seacarb/ChangeLog (rev 0)
+++ pkg/seacarb/ChangeLog 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,108 @@
+--------------------------
+version 2.0.5, 2008-12-08
+
+Changes made thanks to comments kindly provided by Andreas Hofman (A.Hofmann at nioo.knaw.nl):
+- K2 modified to fix a bug which provided warnings when several pressures were used.
+- Khs modified to fix a bug (in coefficient a2).
+- The value of the R constant was updated according to Dickson et al. (2007). The new value is 8.314475 J/(K*mol).
+- kconv was modified in order to include pressure correction.
+
+--------------------------
+version 2.0.4, 2008-11-06
+
+- pH, pTA, pmix, pCa, and pgas now accept Sit (concentration of silicate) and Pt (concentration of total phosphate) as input arguments.
+- Various style changes in the help files.
+
+--------------------------
+version 2.0.3, 2008-11-04
+
+- The code of all functions which calculate constants was modified in order to be able to use vectors as input arguments.
+- The man files were changed accordingly by adding a "details" section.
+
+--------------------------
+version 2.0.1, 2008-10-28
+
+- added the data set seacarb_test
+- added and example to use carb with a data frame
+- various cosmetic changes in the help files
+
+--------------------------
+version 2.0, 2008-10-27
+
+- added function "pCa" which calculates the changes in the saturation states of aragonite and calcite resulting from the manipulation of the calcium concentration
+- added function "pgas" which calculates the carbonate chemistry after changes in pCO2 generated by gas bubbling
+- added function "pmix" which calculates the carbonate chemistry after mixing of two water samples with different pCO2
+- added function "ppH" which calculates the carbonate chemistry after pH manipulations through addition of acid or base
+- added function "pTA" which calculates the carbonate chemistry following addition of CO3 or HCO3
+- carb function modified to return only S, T, P, pH, pCO2, fCO2, HCO3, CO3, DIC, ALK, OmegaAragonite and OmegaCalcite.
+- added function "buffer" which calculates the buffer parameters of the seawater carbonate system.
+- carb function modified in order to closely follow the recommendations of the "Guide to Best Practices for Ocean CO2 Measurements" (Dickson et al., 2007). Phosphate ans ilicate concentrations are now taken into account.
+- K1 function uses the equation of Lueker et al. (2000) by default
+- K2 function uses the equation the method of Lueker et al. (2000) by default
+- Kf function uses, by default, the equation of Perez and Fraga (1987). The equation of Dickson and Roy can be used with the argument kf = 'dg'.
+
+--------------------------
+version 1.2.3, 2007-11-27
+
+- carb function modified to allow NA values in the input data. A warning is returned if this happens but all other data rows are processed. Change made by Bernard Gentili.
+
+--------------------------
+version 1.2.2, 2007-09-09
+
+- correction of the definition of PiH and PhiH in the "carb" documentation file
+- BetaD is the Revelle factor (the "carb" documentation has been updated
+- cosmetic changes to some documentation files
+
+--------------------------
+version 1.2.1, 2007-08-21
+
+- cosmetic changes to some documentation files
+
+--------------------------
+version 1.2, 2007-08-21
+
+- added function "speciation" which estimates the concentration of the various ionic forms of a molecule as a function of pH (contributed by Karline Soetaert)
+- added function "kconv" which provides conversion factors to change the pH scale of dissociation constants (contributed by Karline Soetaert)
+- added function "pHconv" which provides conversion factors for changing the pH scale (contributed by Karline Soetaert)
+- added function "Kn" which provides the ammonium dissociation constant (contributed by Karline Soetaert)
+- added function "Khs" which provides the dissociation constant of hydrogen sulfide (contributed by Karline Soetaert)
+- added function "Ksi" which provides the Si(OH)4 dissociation constant (contributed by Karline Soetaert)
+- added function "bjerrum" which makes a bjerrum plot (contributed by Karline Soetaert)
+
+--------------------------
+version 1.1.1, 2007-08-14
+
+- bor(), K1(), K1p() and K2() now return the values
+- return values are now set for all variables
+
+--------------------------
+version 1.1, 2007-07-24
+
+- carb() has now 5 more flags (code contributed by Jim Orr):
+ # flag = 21 pH-pCO2 given
+ # flag = 22 pCO2-HCO3 given
+ # flag = 23 pCO2-CO3 given
+ # flag = 24 pCO2-ALK given
+ # flag = 25 pCO2-DIC given
+
+--------------------------
+version 1.0.5, 2007-07-23
+
+- minor changes to the documentation (tried to have the superscripts right)
+
+--------------------------
+version 1.0, 2007-01-08
+
+- added function phinsi which calculates the pH at in situ temperature from pH values measured in the laboratory and other ancillay data
+
+--------------------------
+version 1.0, 2007-01-08
+
+- a test data file is now provided
+- various cosmetic changes
+
+--------------------------
+version 0.98, 2006-05-07
+
+- data entry via a file is NO LONGER possible
+- results of the carb command are no longer automatically saved in a file
Added: pkg/seacarb/DESCRIPTION
===================================================================
--- pkg/seacarb/DESCRIPTION (rev 0)
+++ pkg/seacarb/DESCRIPTION 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,10 @@
+Package: seacarb
+Title: Calculates parameters of the seawater carbonate system
+Version: 2.0.5
+Date: 2008-12-08
+Author: Heloise Lavigne and Aurelien Proye and Jean-Pierre Gattuso. Portions of code and/or corrections were contributed by Jean-Marie Epitalon, Bernard Gentili, Jim Orr and Karline Soetaert
+Description: Calculates parameters of the seawater carbonate system
+Maintainer: Jean-Pierre Gattuso <gattuso at obs-vlfr.fr>
+URL: http://www.obs-vlfr.fr/~gattuso/seacarb.php
+License: GPL version 2 or newer
+Packaged: Mon Dec 8 15:36:17 2008; gattuso
Added: pkg/seacarb/R/K1.R
===================================================================
--- pkg/seacarb/R/K1.R (rev 0)
+++ pkg/seacarb/R/K1.R 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,129 @@
+# Copyright (C) 2008 Jean-Pierre Gattuso and Héloïse Lavigne and Aurélien Proye
+#
+# This file is part of seacarb.
+#
+# Seacarb is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or any later version.
+#
+# Seacarb is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with seacarb; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+"K1" <-
+function(S=35,T=25,P=0,k1k2='l')
+{
+
+nK <- max(length(S), length(T), length(P), length(k1k2))
+
+##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
+
+if(length(S)!=nK){S <- rep(S[1], nK)}
+if(length(T)!=nK){T <- rep(T[1], nK)}
+if(length(P)!=nK){P <- rep(P[1], nK)}
+if(length(k1k2)!=nK){k1k2 <- rep(k1k2[1], nK)}
+
+
+#-------Constantes----------------
+
+#---- issues de equic----
+tk = 273.15; # [K] (for conversion [deg C] <-> [K])
+TK = T + tk; # TC [C]; T[K]
+Cl = S / 1.80655; # Cl = chlorinity; S = salinity (per mil)
+cl3 = Cl^(1/3);
+ION = 0.00147 + 0.03592 * Cl + 0.000068 * Cl * Cl; # ionic strength
+iom0 = 19.924*S/(1000-1.005*S);
+ST = 0.14/96.062/1.80655*S; # (mol/kg soln) total sulfate
+
+
+bor = (416.*(S/35.))* 1e-6; # (mol/kg), DOE94
+
+# --------------------- K1 ---------------------------------------
+# first acidity constant:
+# [H^+] [HCO_3^-] / [CO2] = K_1
+#
+# Mehrbach et al (1973) refit by Lueker et al. (2000).
+#
+#(Lueker et al., 2000 in Guide to the Best Practices for Ocean CO2 Measurements
+# Dickson, Sabin and Christian , 2007, Chapter 5, p. 13)
+#
+# pH-scale: 'total'. mol/kg-soln
+
+logK1lue <- (-3633.86)/TK + 61.2172 - 9.67770*log(TK) + 0.011555*S - 0.0001152*S*S
+K1lue <- 10^logK1lue
+
+# --------------------- K1 ---------------------------------------
+# first acidity constant:
+# [H^+] [HCO_3^-] / [CO2] = K_1
+#
+# (Roy et al., 1993 in Dickson and Goyet, 1994, Chapter 5, p. 14)
+# pH-scale: 'total'. mol/kg-soln
+
+tmp1 = 2.83655 - 2307.1266 / TK - 1.5529413 * log(TK);
+tmp2 = - (0.20760841 + 4.0484 / TK) * sqrt(S);
+tmp3 = + 0.08468345 * S - 0.00654208 * S * sqrt(S);
+tmp4 = + log(1 - 0.001005 * S);
+
+lnK1roy = tmp1 + tmp2 + tmp3 + tmp4;
+
+ K1roy = exp(lnK1roy);
+
+
+# ---------- Choice between methods (Lueker or Roy) ----------
+
+K1 <- K1lue
+
+for(i in (1:nK)){
+if(k1k2[i]=='l'){K1[i] <- K1lue[i] }
+if(k1k2[i]=='r'){K1[i] <- K1roy[i] }
+}
+
+# ------------------- Pression effect --------------------------------
+for(i in (1:nK)){
+if (P[i] > 0.0)
+{
+
+ RGAS = 8.314510; # J mol-1 deg-1 (perfect Gas)
+ R = 83.14472; # mol bar deg-1
+ # conversion cm3 -> m3 *1.e-6
+ # bar -> Pa = N m-2 *1.e+5
+ # => *1.e-1 or *1/10
+
+
+ # index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8,
+ # K1P 9, K2P 10, K3P 11
+
+ #----- note: there is an error in Table 9 of Millero, 1995.
+ #----- The coefficients -b0 and b1
+ #----- have to be multiplied by 1.e-3!
+
+ #----- there are some more errors!
+ #----- the signs (+,-) of coefficients in Millero 95 do not
+ #----- agree with Millero 79
+
+
+
+ a0 = c(-25.5, -15.82, -29.48, -25.60, -18.03, -9.78, -48.76, -46, -14.51, -23.12, -26.57);
+ a1 = c(0.1271, -0.0219, 0.1622, 0.2324, 0.0466, -0.0090, 0.5304, 0.5304, 0.1211, 0.1758, 0.2020);
+ a2 = c(0.0, 0.0, 2.608*1e-3, -3.6246*1e-3, 0.316*1e-3, -0.942*1e-3, 0.0, 0.0, -0.321*1e-3, -2.647*1e-3, -3.042*1e-3);
+ b0 = c(-3.08*1e-3, 1.13*1e-3, -2.84*1e-3, -5.13*1e-3, -4.53*1e-3, -3.91*1e-3, -11.76*1e-3, -11.76*1e-3, -2.67*1e-3, -5.15*1e-3, -4.08*1e-3);
+ b1 = c(0.0877*1e-3, -0.1475*1e-3, 0.0, 0.0794*1e-3, 0.09*1e-3, 0.054*1e-3, 0.3692*1e-3, 0.3692*1e-3, 0.0427*1e-3, 0.09*1e-3, 0.0714*1e-3);
+ b2 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ deltav = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ deltak = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ lnkpok0 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ for (ipc in 1:length(a0))
+ {
+ deltav[ipc] = a0[ipc] + a1[ipc] *T[i] + a2[ipc] *T[i]*T[i];
+ deltak[ipc] = (b0[ipc] + b1[ipc] *T[i] + b2[ipc] *T[i]*T[i]);
+ lnkpok0[ipc] = -(deltav[ipc] /(R*TK[i]))*P[i] + (0.5*deltak[ipc] /(R*TK[i]))*P[i]*P[i];
+ }
+
+ K1[i] = K1[i]*exp(lnkpok0[1]);
+
+}
+}
+attr(K1,"unit") = "mol/kg-soln"
+attr(K1,"pH scale") = "total hydrogen ion concentration"
+return(K1)
+}
Added: pkg/seacarb/R/K1p.R
===================================================================
--- pkg/seacarb/R/K1p.R (rev 0)
+++ pkg/seacarb/R/K1p.R 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,99 @@
+# Copyright (C) 2003 Jean-Pierre Gattuso and Aurelien Proye
+#
+# This file is part of seacarb.
+#
+# Seacarb is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or any later version.
+#
+# Seacarb is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with seacarb; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+#
+"K1p" <-
+function(S=35,T=25,P=0){
+
+nK <- max(length(S), length(T), length(P))
+
+##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
+
+if(length(S)!=nK){S <- rep(S[1], nK)}
+if(length(T)!=nK){T <- rep(T[1], nK)}
+if(length(P)!=nK){P <- rep(P[1], nK)}
+
+
+#-------Constantes----------------
+
+#---- issues de equic----
+tk = 273.15; # [K] (for conversion [deg C] <-> [K])
+TC = T + tk; # TC [C]; T[K]
+Cl = S / 1.80655; # Cl = chlorinity; S = salinity (per mille)
+cl3 = Cl^(1/3);
+ION = 0.00147 + 0.03592 * Cl + 0.000068 * Cl * Cl; # ionic strength
+iom0 = 19.924*S/(1000-1.005*S);
+ST = 0.14/96.062/1.80655*S; # (mol/kg soln) total sulfate
+
+
+bor = (416.*(S/35.))* 1e-6; # (mol/kg), DOE94
+
+ # --------------------- Phosphoric acid ---------------------
+ #
+ # (DOE, 1994) (Dickson and Goyet): pH_T, mol/(kg-soln)
+ # Ch.5 p. 16
+ #
+
+ lnK1P = -4576.752 / TC + 115.525 - 18.453*log(TC) + (-106.736 / TC + 0.69171) * sqrt(S) + (-0.65643 / TC - 0.01844) * S;
+
+ K1P = exp(lnK1P);
+
+for(i in (1:nK)){
+ if (P[i] > 0.0)
+ {
+
+ RGAS = 8.314510; # J mol-1 deg-1 (perfect Gas)
+ R = 83.14472; # mol bar deg-1
+ # conversion cm3 -> m3 *1.e-6
+ # bar -> Pa = N m-2 *1.e+5
+ # => *1.e-1 or *1/10
+
+
+ # index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8,
+ # K1P 9, K2P 10, K3P 11
+
+ #----- note: there is an error in Table 9 of Millero, 1995.
+ #----- The coefficients -b0 and b1
+ #----- have to be multiplied by 1.e-3!
+
+ #----- there are some more errors!
+ #----- the signs (+,-) of coefficients in Millero 95 do not
+ #----- agree with Millero 79
+
+
+
+ a0 = c(-25.5, -15.82, -29.48, -25.60, -18.03, -9.78, -48.76, -46, -14.51, -23.12, -26.57);
+ a1 = c(0.1271, -0.0219, 0.1622, 0.2324, 0.0466, -0.0090, 0.5304, 0.5304, 0.1211, 0.1758, 0.2020);
+ a2 = c(0.0, 0.0, 2.608*1e-3, -3.6246*1e-3, 0.316*1e-3, -0.942*1e-3, 0.0, 0.0, -0.321*1e-3, -2.647*1e-3, -3.042*1e-3);
+ b0 = c(-3.08*1e-3, 1.13*1e-3, -2.84*1e-3, -5.13*1e-3, -4.53*1e-3, -3.91*1e-3, -11.76*1e-3, -11.76*1e-3, -2.67*1e-3, -5.15*1e-3, -4.08*1e-3);
+ b1 = c(0.0877*1e-3, -0.1475*1e-3, 0.0, 0.0794*1e-3, 0.09*1e-3, 0.054*1e-3, 0.3692*1e-3, 0.3692*1e-3, 0.0427*1e-3, 0.09*1e-3, 0.0714*1e-3);
+ b2 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ deltav = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ deltak = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ lnkpok0 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ for (ipc in 1:length(a0))
+ {
+ deltav[ipc] = a0[ipc] + a1[ipc] *T[i] + a2[ipc] *T[i]*T[i];
+ deltak[ipc] = (b0[ipc] + b1[ipc] *T[i] + b2[ipc] *T[i]*T[i]);
+ lnkpok0[ipc] = -(deltav[ipc] /(R*TC[i]))*P[i] + (0.5*deltak[ipc] /(R*TC[i]))*P[i]*P[i];
+ }
+
+ K1P[i] = K1P[i]*exp(lnkpok0[9]);
+
+
+ }
+}
+
+ attr(K1P,"unit") = "mol/kg-soln"
+ attr(K1P,"pH scale") = "total hydrogen ion concentration"
+ return(K1P)
+}
Added: pkg/seacarb/R/K2.R
===================================================================
--- pkg/seacarb/R/K2.R (rev 0)
+++ pkg/seacarb/R/K2.R 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,137 @@
+# Copyright (C) 2008 Jean-Pierre Gattuso and Héloïse Lavigne and Aurélien Proye
+#
+# This file is part of seacarb.
+#
+# Seacarb is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or any later version.
+#
+# Seacarb is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with seacarb; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+#
+"K2" <-
+function(S=35,T=25,P=0,k1k2='l'){
+
+nK <- max(length(S), length(T), length(P), length(k1k2))
+
+##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
+
+if(length(S)!=nK){S <- rep(S[1], nK)}
+if(length(T)!=nK){T <- rep(T[1], nK)}
+if(length(P)!=nK){P <- rep(P[1], nK)}
+if(length(k1k2)!=nK){k1k2 <- rep(k1k2[1], nK)}
+
+
+#-------Constantes----------------
+
+#---- issues de equic----
+tk = 273.15; # [K] (for conversion [deg C] <-> [K])
+TK = T + tk; # TC [C]; T[K]
+Cl = S / 1.80655; # Cl = chlorinity; S = salinity (per mille)
+cl3 = Cl^(1/3);
+ION = 0.00147 + 0.03592 * Cl + 0.000068 * Cl * Cl; # ionic strength
+iom0 = 19.924*S/(1000-1.005*S);
+ST = 0.14/96.062/1.80655*S; # (mol/kg soln) total sulfate
+
+bor = (416.*(S/35.))* 1e-6; # (mol/kg), DOE94
+
+
+
+# --------------------- K2 lueker et al. 2000 ------------------------------
+#
+# second acidity constant:
+# [H^+] [CO_3^--] / [HCO_3^-] = K_2
+#
+# Mehrbach et al. (1973) refit by Lueker et al. (2000).
+#
+# (Lueker et al., 2000 in Guide to the Best Practices for Ocean CO2 Measurements
+# Dickson, Sabin and Christian , 2007, Chapter 5, p. 14)
+#
+# pH-scale: 'total'. mol/kg-soln
+
+logK2lue <- -471.78/TK - 25.9290 + 3.16967*log(TK) + 0.01781*S - 0.0001122*S*S
+
+K2lue <- 10^(logK2lue)
+
+
+
+# --------------------- K2 Roy et al. 1993----------------------------------------
+#
+# second acidity constant:
+# [H^+] [CO_3^--] / [HCO_3^-] = K_2
+#
+# (Roy et al., 1993 in Dickson and Goyet, 1994, Chapter 5, p. 15)
+# pH-scale: 'total'. mol/kg-soln
+
+tmp1 = -9.226508 - 3351.6106 / TK - 0.2005743 * log(TK);
+tmp2 = (-0.106901773 - 23.9722 / TK) * sqrt(S);
+tmp3 = 0.1130822 * S - 0.00846934 * S^1.5 + log(1 - 0.001005 * S);
+
+lnK2roy = tmp1 + tmp2 + tmp3;
+
+K2roy <- exp(lnK2roy)
+
+# ---------- Choice between methods (Lueker or Roy) ----------
+
+K2 <- K2lue
+
+for(i in (1:nK)){
+if(k1k2[i]=='l'){K2[i] <- K2lue[i] }
+if(k1k2[i]=='r'){K2[i] <- K2roy[i] }
+}
+
+# ------------------- Pression effect --------------------------------
+for(i in (1:nK)){
+if (P[i] > 0.0)
+{
+
+ RGAS = 8.314510; # J mol-1 deg-1 (perfect Gas)
+ R = 83.14472; # mol bar deg-1
+ # conversion cm3 -> m3 *1.e-6
+ # bar -> Pa = N m-2 *1.e+5
+ # => *1.e-1 or *1/10
+
+
+ # index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8,
+ # K1P 9, K2P 10, K3P 11
+
+ #----- note: there is an error in Table 9 of Millero, 1995.
+ #----- The coefficients -b0 and b1
+ #----- have to be multiplied by 1.e-3!
+
+ #----- there are some more errors!
+ #----- the signs (+,-) of coefficients in Millero 95 do not
+ #----- agree with Millero 79
+
+
+
+ a0 = c(-25.5, -15.82, -29.48, -25.60, -18.03, -9.78, -48.76, -46, -14.51, -23.12, -26.57);
+ a1 = c(0.1271, -0.0219, 0.1622, 0.2324, 0.0466, -0.0090, 0.5304, 0.5304, 0.1211, 0.1758, 0.2020);
+ a2 = c(0.0, 0.0, 2.608*1e-3, -3.6246*1e-3, 0.316*1e-3, -0.942*1e-3, 0.0, 0.0, -0.321*1e-3, -2.647*1e-3, -3.042*1e-3);
+ b0 = c(-3.08*1e-3, 1.13*1e-3, -2.84*1e-3, -5.13*1e-3, -4.53*1e-3, -3.91*1e-3, -11.76*1e-3, -11.76*1e-3, -2.67*1e-3, -5.15*1e-3, -4.08*1e-3);
+ b1 = c(0.0877*1e-3, -0.1475*1e-3, 0.0, 0.0794*1e-3, 0.09*1e-3, 0.054*1e-3, 0.3692*1e-3, 0.3692*1e-3, 0.0427*1e-3, 0.09*1e-3, 0.0714*1e-3);
+ b2 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ deltav = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ deltak = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ lnkpok0 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ for (ipc in 1:length(a0))
+ {
+ deltav[ipc] = a0[ipc] + a1[ipc] *T[i] + a2[ipc] *T[i]*T[i];
+ deltak[ipc] = (b0[ipc] + b1[ipc] *T[i] + b2[ipc] *T[i]*T[i]);
+ lnkpok0[ipc] = -(deltav[ipc] /(R*TK[i]))*P[i] + (0.5*deltak[ipc] /(R*TK[i]))*P[i]*P[i];
+ }
+
+ K2[i] <- K2[i]*exp(lnkpok0[2])
+}
+}
+
+attr(K2,"unit") = "mol/kg-soln"
+attr(K2,"pH scale") = "total hydrogen ion concentration"
+return(K2)
+}
+
+
+
+
Added: pkg/seacarb/R/K2p.R
===================================================================
--- pkg/seacarb/R/K2p.R (rev 0)
+++ pkg/seacarb/R/K2p.R 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,103 @@
+# Copyright (C) 2008 Jean-Pierre Gattuso and Héloïse Lavigne and Aurélien Proye
+#
+# This file is part of seacarb.
+#
+# Seacarb is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or any later version.
+#
+# Seacarb is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with seacarb; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+#
+"K2p" <-
+function(S=35,T=25,P=0){
+
+nK <- max(length(S), length(T), length(P))
+
+##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
+
+if(length(S)!=nK){S <- rep(S[1], nK)}
+if(length(T)!=nK){T <- rep(T[1], nK)}
+if(length(P)!=nK){P <- rep(P[1], nK)}
+
+
+#-------Constantes----------------
+
+#---- issues de equic----
+tk = 273.15; # [K] (for conversion [deg C] <-> [K])
+TK = T + tk; # T [C]; TK[K]
+Cl = S / 1.80655; # Cl = chlorinity; S = salinity (per mille)
+cl3 = Cl^(1/3);
+ION = 0.00147 + 0.03592 * Cl + 0.000068 * Cl * Cl; # ionic strength
+iom0 = 19.924*S/(1000-1.005*S);
+ST = 0.14/96.062/1.80655*S; # (mol/kg soln) total sulfate
+
+
+bor = (416.*(S/35.))* 1e-6; # (mol/kg), DOE94
+
+ # --------------------- Phosphoric acid ---------------------
+ #
+ #
+ # Guide to Best Practices in Ocean CO2 Measurements 2007 Chap 5 p 15
+ # (Dickson and Goyet): pH_T, mol/(kg-soln)
+ #
+ #
+
+ lnK2P = -8814.715 / TK + 172.0883 - 27.927 * log(TK) + (-160.34 / TK + 1.3566) * sqrt(S) + (0.37335 / TK - 0.05778) * S;
+
+ K2P = exp(lnK2P);
+
+ for(i in (1:nK)){
+ if (P[i] > 0.0)
+ {
+
+ RGAS = 8.314510; # J mol-1 deg-1 (perfect Gas)
+ R = 83.14472; # mol bar deg-1
+ # conversion cm3 -> m3 *1.e-6
+ # bar -> Pa = N m-2 *1.e+5
+ # => *1.e-1 or *1/10
+
+
+ # index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8,
+ # K1P 9, K2P 10, K3P 11
+
+ #----- note: there is an error in Table 9 of Millero, 1995.
+ #----- The coefficients -b0 and b1
+ #----- have to be multiplied by 1.e-3!
+
+ #----- there are some more errors!
+ #----- the signs (+,-) of coefficients in Millero 95 do not
+ #----- agree with Millero 79
+
+
+
+ a0 = c(-25.5, -15.82, -29.48, -25.60, -18.03, -9.78, -48.76, -46, -14.51, -23.12, -26.57);
+ a1 = c(0.1271, -0.0219, 0.1622, 0.2324, 0.0466, -0.0090, 0.5304, 0.5304, 0.1211, 0.1758, 0.2020);
+ a2 = c(0.0, 0.0, 2.608*1e-3, -3.6246*1e-3, 0.316*1e-3, -0.942*1e-3, 0.0, 0.0, -0.321*1e-3, -2.647*1e-3, -3.042*1e-3);
+ b0 = c(-3.08*1e-3, 1.13*1e-3, -2.84*1e-3, -5.13*1e-3, -4.53*1e-3, -3.91*1e-3, -11.76*1e-3, -11.76*1e-3, -2.67*1e-3, -5.15*1e-3, -4.08*1e-3);
+ b1 = c(0.0877*1e-3, -0.1475*1e-3, 0.0, 0.0794*1e-3, 0.09*1e-3, 0.054*1e-3, 0.3692*1e-3, 0.3692*1e-3, 0.0427*1e-3, 0.09*1e-3, 0.0714*1e-3);
+ b2 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ deltav = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ deltak = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ lnkpok0 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ for (ipc in 1:length(a0))
+ {
+ deltav[ipc] = a0[ipc] + a1[ipc] *T[i] + a2[ipc] *T[i]*T[i];
+ deltak[ipc] = (b0[ipc] + b1[ipc] *T[i] + b2[ipc] *T[i]*T[i]);
+ lnkpok0[ipc] = -(deltav[ipc] /(R*TK[i]))*P[i] + (0.5*deltak[ipc] /(R*TK[i]))*P[i]*P[i];
+ }
+
+
+ K2P[i] = K2P[i]*exp(lnkpok0[10]);
+
+
+ }
+
+}
+
+attr(K2p,"unit") = "mol/kg-soln"
+attr(K2p,"pH scale") = "total hydrogen ion concentration"
+return(K2P)
+}
Added: pkg/seacarb/R/K3p.R
===================================================================
--- pkg/seacarb/R/K3p.R (rev 0)
+++ pkg/seacarb/R/K3p.R 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,98 @@
+# Copyright (C) 2008 Jean-Pierre Gattuso and Héloïse Lavigne and Aurélien Proye
+#
+# This file is part of seacarb.
+#
+# Seacarb is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or any later version.
+#
+# Seacarb is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with seacarb; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+#
+"K3p" <-
+function(S=35,T=25,P=0){
+
+nK <- max(length(S), length(T), length(P))
+
+##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
+
+if(length(S)!=nK){S <- rep(S[1], nK)}
+if(length(T)!=nK){T <- rep(T[1], nK)}
+if(length(P)!=nK){P <- rep(P[1], nK)}
+
+
+
+#-------Constantes----------------
+
+#---- issues de equic----
+tk = 273.15; # [K] (for conversion [deg C] <-> [K])
+TK = T + tk; # T [C]; TK[K]
+Cl = S / 1.80655; # Cl = chlorinity; S = salinity (per mille)
+cl3 = Cl^(1/3);
+ION = 0.00147 + 0.03592 * Cl + 0.000068 * Cl * Cl; # ionic strength
+iom0 = 19.924*S/(1000-1.005*S);
+ST = 0.14/96.062/1.80655*S; # (mol/kg soln) total sulfate
+
+
+bor = (416.*(S/35.))* 1e-6; # (mol/kg), DOE94
+
+ # --------------------- Phosphoric acid ---------------------
+ #
+ # Guide to Best Practices in Ocean CO2 Measurements 2007 Chap 5 p 15
+ #
+ # Ch.5 p. 15
+ #
+
+ lnK3P = -3070.75 / TK - 18.141 + (17.27039 / TK + 2.81197) * sqrt(S) + (-44.99486 / TK - 0.09984) * S;
+
+ K3P = exp(lnK3P);
+
+ for(i in (1:nK)){
+ if (P[i] > 0.0)
+ {
+
+ RGAS = 8.314510; # J mol-1 deg-1 (perfect Gas)
+ R = 83.14472; # mol bar deg-1
+ # conversion cm3 -> m3 *1.e-6
+ # bar -> Pa = N m-2 *1.e+5
+ # => *1.e-1 or *1/10
+
+
+ # index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8,
+ # K1P 9, K2P 10, K3P 11
+
+ #----- note: there is an error in Table 9 of Millero, 1995.
+ #----- The coefficients -b0 and b1
+ #----- have to be multiplied by 1.e-3!
+
+ #----- there are some more errors!
+ #----- the signs (+,-) of coefficients in Millero 95 do not
+ #----- agree with Millero 79
+
+ a0 = c(-25.5, -15.82, -29.48, -25.60, -18.03, -9.78, -48.76, -46, -14.51, -23.12, -26.57);
+ a1 = c(0.1271, -0.0219, 0.1622, 0.2324, 0.0466, -0.0090, 0.5304, 0.5304, 0.1211, 0.1758, 0.2020);
+ a2 = c(0.0, 0.0, 2.608*1e-3, -3.6246*1e-3, 0.316*1e-3, -0.942*1e-3, 0.0, 0.0, -0.321*1e-3, -2.647*1e-3, -3.042*1e-3);
+ b0 = c(-3.08*1e-3, 1.13*1e-3, -2.84*1e-3, -5.13*1e-3, -4.53*1e-3, -3.91*1e-3, -11.76*1e-3, -11.76*1e-3, -2.67*1e-3, -5.15*1e-3, -4.08*1e-3);
+ b1 = c(0.0877*1e-3, -0.1475*1e-3, 0.0, 0.0794*1e-3, 0.09*1e-3, 0.054*1e-3, 0.3692*1e-3, 0.3692*1e-3, 0.0427*1e-3, 0.09*1e-3, 0.0714*1e-3);
+ b2 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ deltav = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ deltak = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ lnkpok0 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ for (ipc in 1:length(a0))
+ {
+ deltav[ipc] = a0[ipc] + a1[ipc] *T[i] + a2[ipc] *T[i]*T[i];
+ deltak[ipc] = (b0[ipc] + b1[ipc] *T[i] + b2[ipc] *T[i]*T[i]);
+ lnkpok0[ipc] = -(deltav[ipc] /(R*TK[i]))*P[i] + (0.5*deltak[ipc] /(R*TK[i]))*P[i]*P[i];
+ }
+
+
+ K3P[i] = K3P[i]*exp(lnkpok0[11]);
+
+ }
+}
+ attr(K3p,"unit") = "mol/kg-soln"
+ attr(K3p,"pH scale") = "total hydrogen ion concentration"
+ return(K3P)
+}
Added: pkg/seacarb/R/Kb.R
===================================================================
--- pkg/seacarb/R/Kb.R (rev 0)
+++ pkg/seacarb/R/Kb.R 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,106 @@
+# Copyright (C) 2008 Jean-Pierre Gattuso and Héloïse Lavigne and Aurélien Proye
+#
+# This file is part of seacarb.
+#
+# Seacarb is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or any later version.
+#
+# Seacarb is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with seacarb; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+#
+
+"Kb" <-
+function(S=35,T=25,P=0){
+
+nK <- max(length(S), length(T), length(P))
+
+##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
+
+if(length(S)!=nK){S <- rep(S[1], nK)}
+if(length(T)!=nK){T <- rep(T[1], nK)}
+if(length(P)!=nK){P <- rep(P[1], nK)}
+
+
+#-------Constantes----------------
+
+#---- issues de equic----
+tk = 273.15; # [K] (for conversion [deg C] <-> [K])
+TK = T + tk; # T [C]; TK [K]
+Cl = S / 1.80655; # Cl = chlorinity; S = salinity (per mille)
+cl3 = Cl^(1/3);
+ION = 0.00147 + 0.03592 * Cl + 0.000068 * Cl * Cl; # ionic strength
+iom0 = 19.924*S/(1000-1.005*S);
+ST = 0.14/96.062/1.80655*S; # (mol/kg soln) total sulfate
+
+
+bor = (416.*(S/35.))* 1e-6; # (mol/kg), DOE94
+
+
+ #---------------------------------------------------------------------
+ # --------------------- Kb --------------------------------------------
+ # Kbor = [H+][B(OH)4-]/[B(OH)3]
+ #
+ # (Dickson, 1990 in Guide to Best Practices in Ocean CO2 Measurements 2007)
+ # pH-scale: 'total'. mol/kg-soln
+
+
+ tmp1 = (-8966.90-2890.53*sqrt(S)-77.942*S+1.728*S^(3/2)-0.0996*S*S);
+ tmp2 = +148.0248+137.1942*sqrt(S)+1.62142*S;
+ tmp3 = +(-24.4344-25.085*sqrt(S)-0.2474*S)*log(TK);
+
+ lnKb = tmp1 / TK + tmp2 + tmp3 + 0.053105*sqrt(S)*TK;
+ Kb <- exp(lnKb)
+
+ # ------------------- Pression effect --------------------------------
+ for(i in (1:nK)){
+
+ if (P[i] > 0.0){
+
+ RGAS = 8.314510; # J mol-1 deg-1 (perfect Gas)
+ R = 83.14472; # mol bar deg-1
+ # conversion cm3 -> m3 *1.e-6
+ # bar -> Pa = N m-2 *1.e+5
+ # => *1.e-1 or *1/10
+
+
+ # index: K1 1, K2 2, Kb 3, Kw 4, Ks 5, Kf 6, Kspc 7, Kspa 8,
+ # K1P 9, K2P 10, K3P 11
+
+ #----- note: there is an error in Table 9 of Millero, 1995.
+ #----- The coefficients -b0 and b1
+ #----- have to be multiplied by 1.e-3!
+
+ #----- there are some more errors!
+ #----- the signs (+,-) of coefficients in Millero 95 do not
+ #----- agree with Millero 79
+
+
+
+ a0 = c(-25.5, -15.82, -29.48, -25.60, -18.03, -9.78, -48.76, -46, -14.51, -23.12, -26.57);
+ a1 = c(0.1271, -0.0219, 0.1622, 0.2324, 0.0466, -0.0090, 0.5304, 0.5304, 0.1211, 0.1758, 0.2020);
+ a2 = c(0.0, 0.0, 2.608*1e-3, -3.6246*1e-3, 0.316*1e-3, -0.942*1e-3, 0.0, 0.0, -0.321*1e-3, -2.647*1e-3, -3.042*1e-3);
+ b0 = c(-3.08*1e-3, 1.13*1e-3, -2.84*1e-3, -5.13*1e-3, -4.53*1e-3, -3.91*1e-3, -11.76*1e-3, -11.76*1e-3, -2.67*1e-3, -5.15*1e-3, -4.08*1e-3);
+ b1 = c(0.0877*1e-3, -0.1475*1e-3, 0.0, 0.0794*1e-3, 0.09*1e-3, 0.054*1e-3, 0.3692*1e-3, 0.3692*1e-3, 0.0427*1e-3, 0.09*1e-3, 0.0714*1e-3);
+ b2 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ deltav = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ deltak = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+ lnkpok0 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ for (ipc in 1:length(a0))
+ {
+ deltav[ipc] = a0[ipc] + a1[ipc] *T[i] + a2[ipc] *T[i]*T[i];
+ deltak[ipc] = (b0[ipc] + b1[ipc] *T[i] + b2[ipc] *T[i]*T[i]);
+ lnkpok0[ipc] = -(deltav[ipc] /(R*TK[i]))*P[i] + (0.5*deltak[ipc] /(R*TK[i]))*P[i]*P[i];
+ }
+
+
+ Kb[i] = Kb[i]*exp(lnkpok0[3])
+
+ }
+}
+ attr(Kb,"unit") = "mol/kg-soln"
+ attr(Kb,"pH scale") = "total hydrogen ion concentration"
+ return(Kb)
+}
Added: pkg/seacarb/R/Kf.R
===================================================================
--- pkg/seacarb/R/Kf.R (rev 0)
+++ pkg/seacarb/R/Kf.R 2008-12-24 16:00:56 UTC (rev 2)
@@ -0,0 +1,131 @@
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/seacarb -r 2
More information about the Seacarb-commits
mailing list