function [k] = ph_tdconst(t,s,p)
% PH_eqconst $Revision: 1.0 $ $Date: 1998/12/07 $
%
% USAGE: [k] = ph_tdconst(t,s,p)
%
% DESCRIPTION:
%   ph_tdconst calculates the values of thermodynamic constants
%   k(i), in seawater for given temperature, salinity, pressure. The
%   calculated constants are on the sws ph scale.
%   these constants assume the concentrations are in molality.
%   the equilibrium constants to be calculated are:
%   k(1) = dissociation constant of h2o
%   k(2) = 1st dissociation constant of co2(aq)
%   k(3) = 2nd dissociation constant of co2(aq)
%   k(4) = dissociation constant of b(oh)3
%   k(5) = 1st dissociation constant of h3po4
%   k(6) = 2nd dissociation constant of h3po4
%   k(7) = 3nd dissociation constant of h3po4
%   k(8) = dissociation constant of hno3
%   k(9) = dissociation constant of nh4
%   k(10) = 1st dissociation constant of h2s
%   k(11) = solubility product of caco3 (aragonite)
%   k(12) = solubility product of caco3 (calcite)
%  correction to in situ pressure uses the basic correction formula:
%        ln(kp/k0) = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp
%  where dv is the change in partial molal volume and dk is the
%  compressibility. pressure pp is water pressure over atmospheric
%  in bars.
%
% INPUT:
%  t     : temperature of seawater at the specific location [deg C]
%  s     : salinity of seawater at the specific location  [PSU]
%  p     : pressure at the specific location  [atm, (with 1 atm
%          at sea surface)]
%
% OUTPUT:
%  k     : equilibrium constants 
%
% ORIGINAL AUTHOR: 
%   Bernard Boudreau 
%   Dept. Oceanography 
%   Dalhousie University
%   Halifax, NS B3H 4J1,Canada 
%   bernie@boudreau.ocean.dal.ca
% MATLAB VERSION FROM
%   Roger Luff 98-12-07
%   GEOMAR Research Centre for Marine Geosciences
%   Wischhofstr. 1-3, D-24148 Kiel, Germany
%   visumod@technologist.com
%
% DISCLAIMER:
%   This software is provided "as it is" without warranty of any kind.
%
% REFERENCES:
%  Boudreau B. P. (1996) "A method-of-lines code for carbon and 
%   nutrient diagenesis in aquatic sediments". Computers & Geosciences
%   22(5), 479-496.
%  Clegg S. L. and Whitefield M. (1995) "A chemical model of seawater 
%   including dissolved ammonia and the stoichiometric dissociation 
%   constant of ammonia in estuarine water and seawater from -2 to 
%   40C." Geochim. Cosmochim. Acta 59(12), 2403-2421.
%  Fofonoff P. and Millard R. C. (1983) "Algorithms for computation 
%   of fundamental properties of seawater." UNESCO Technical Paper 
%   in Marine Science 44, 53.
%  Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN
%   and MATLAB libraries to calculate pH distributions in
%   marine systems", Computers & Geosciences.
%  Millero F. J. (1983) "Influence of pressure on chemical processes 
%   in the sea." In Chemical Oceanography 2nd ed., Vol. 8 (ed. J. P. 
%   Riley and R. Chester), pp. 1-88. Academic Press.
%  Millero F. J. (1995) "Thermodynamics of the carbon dioxide system 
%   in the oceans." Geochim. Cosmochim. Acta 59(4), 661-677.
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% CHECK INPUT ARGUMENTS
%-------------------------------------------------------------------
if nargin ~=3
   error(' ph_eqconst.m: Must pass 3 parameters')
end

tk = t + 273.15d+00;
lntk = log(tk);
s2 = s*s;
sqs = sqrt(s);
s23 = sqrt(s*s*s);
pp = (p-1.0)*1.013d+00;
r = 83.147d+00;
%-------------------------------------------------------------------
%       1) water:
%-------------------------------------------------------------------
lnk1 = 148.9802 - 13847.26/tk - 23.6521*lntk;
lnk1 = lnk1 + (-5.977 + 118.67/tk + 1.0495*lntk)*sqs - 0.01615*s;
dvw = - 20.02 + 0.1119*t - 1.4096d-03*t*t;
dkw = (-5.13 + 0.0794*t)/1000.0;
lnk1 = - dvw/(r*tk)*pp + 0.5*dkw/(r*tk)*pp*pp + lnk1;
k(1) = exp(lnk1);
%-------------------------------------------------------------------
%     2) carbonic acid   (using the combined, i.e. com, equations)
%-------------------------------------------------------------------
lnkc1 = 2.18867 - 2275.036/tk - 1.468591*lntk;
lnkc1 = lnkc1 + (-0.138681 - 9.33291/tk)*sqs;
lnkc1 = lnkc1 + 0.0726483*s - 0.00574938*s23;
a0 = 25.50;
a1 = 0.151;
a2 = - 0.1271;
dv = -(a0 + a1*(s-34.8) + a2*t);
b0 = 3.08;
b1 = 0.578;
b2 = - 0.0877;
dk = - (b0 + b1*(s-34.8) + b2*t)/1000.0;
lnkc1 = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnkc1;
k(2) = exp(lnkc1);

lnkc2 = -0.84226 - 3741.1288/tk - 1.437139*lntk;
lnkc2 = lnkc2 + (-0.128417 - 24.41239/tk)*sqs;
lnkc2 = lnkc2 + 0.1195308*s - 0.0091284*s23;
a0 = 15.82;
a1 = - 0.321;
a2 = 0.0219;
dv = -(a0 + a1*(s-34.8) + a2*t);
b0 = - 1.13;
b1 = 0.314;
b2 = 0.1475;
dk = - (b0 + b1*(s-34.8) + b2*t)/1000.0;
lnkc2 = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnkc2;
k(3) = exp(lnkc2);
%-------------------------------------------------------------------
%     3) boric acid
%-------------------------------------------------------------------
lnkb = (-8966.90 -2890.51*sqs -77.942*s +1.726*s23 -0.0993*s2)/tk;
lnkb = lnkb + (148.0248 + 137.194*sqs + 1.62247*s);
lnkb = lnkb + (-24.4344 - 25.085*sqs - 0.2474*s)*lntk;
lnkb = lnkb + 0.053105*sqs*tk;
a0 = 29.48;
a1 = - 0.295;
a2 = - 0.1622;
a3 = 2.608d-03;
dv = -(a0 + a1*(s-34.8) + a2*t + a3*t*t);
b0 = 2.84;
b1 = - 0.354;
dk = - (b0 + b1*(s-34.8))/1000.0;
lnkb = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnkb;
k(4) = exp(lnkb);
%-------------------------------------------------------------------
%     4) phosphoric acid
%-------------------------------------------------------------------
lnkp1 = 115.54 - 4576.752/tk - 18.453*lntk;
lnkp1 = lnkp1 + (0.69171 - 106.736/tk)*sqs;
lnkp1 = lnkp1 + (-0.01844 - 0.65643/tk)*s;
dv = - 14.51 + 0.1211*t - 0.321d-03*t*t;
dk = (-2.67 + 0.0427*t)/1000.0;
lnkp1 = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnkp1;
k(5) = exp(lnkp1);

lnkp2 = 172.1033 - 8814.715/tk - 27.927*lntk;
lnkp2 = lnkp2 + (1.3566 - 160.34/tk)*sqs;
lnkp2 = lnkp2 + (-0.05778 + 0.37335/tk)*s;
dv = - 23.12 + 0.1758*t - 2.647d-03*t*t;
dk = (-5.15 + 0.0900*t)/1000.0;
lnkp2 = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnkp2;
k(6) = exp(lnkp2);

lnkp3 = -18.126 - 3070.75/tk;
lnkp3 = lnkp3 + (2.81197 + 17.27039/tk)*sqs;
lnkp3 = lnkp3 + (-0.09984 - 44.99486/tk)*s;
dv = - 26.57 + 0.2020*t - 3.042d-03*t*t;
dk = (-4.08 + 0.0714*t)/1000.0;
lnkp3 = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnkp3;
k(7) = exp(lnkp3);
%-------------------------------------------------------------------
%     5) nitric acid
%-------------------------------------------------------------------
k(8) = 23.44;
%-------------------------------------------------------------------
%     6) ammonia
%-------------------------------------------------------------------
lnk9 =-6285.33/tk+0.0001635*tk-0.25444+(0.46532-123.7184/tk)*sqs+(-0.01992+3.17556/tk)*s;
%-------------------------------------------------------------------
%     note: millero (1983) gives the dv and dk for the reaction
%                   nh3 + h2o = nh4 + oh
%     while the value calculated by this program is for the reaction
%                      nh4 = nh3 + h
%     to get the desired value, the dv for water is subtracted from
%     the millero dv and the the negative of this resulting number is
%     used.
%-------------------------------------------------------------------
dv = - 26.43 + 0.0889*t - 0.905-03*t*t;
dk = (-5.03 + 0.0814*t)/1000.0;
dv = - (dv - dvw);
dk = - (dk - dkw);
lnk9 = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnk9;
k(9) = exp(lnk9);
%-------------------------------------------------------------------
%     7) hydrogen sulfide
%-------------------------------------------------------------------
lnks = 225.838 - 13275/tk - 34.6435*lntk + 0.3449*sqs - 0.0274*s;
dv = - 11.07 + 0.009*t - 0.942-03*t*t;
dk = (-2.89 + 0.054*t)/1000.0;
lnks = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp + lnks;
k(10) = exp(lnks);
%-------------------------------------------------------------------
%     8) calcite & aragonite
%-------------------------------------------------------------------
logk12 = - 171.9065 - 0.077993*tk + 2839.319/tk + 71.595*log10(tk);
logk12 = logk12+(-0.77712+2.8426d-03*tk+178.34/tk)*sqs-0.07711*s+4.1249d-03*s23;
k(12) = 10^(logk12);
lnk12 = log(k(12));
dvc = - (48.76 - 0.5304*t);
dkc = - (11.76 - 0.3692*t)/1000.0;
lnk12 = - dvc/(r*tk)*pp + 0.5*dkc/(r*tk)*pp*pp + lnk12;
k(12) = exp(lnk12);
logk11 = - 171.945 - 0.077993*tk + 2903.293/tk + 71.595*log10(tk);
logk11 = logk11 + (-0.068393 + 1.7276d-03*tk + 88.135/tk)*sqs-0.10018*s + 5.9415d-03*s23;
k(11) = 10^(logk11);
dva = dvc + 2.8;
dka = dkc;
lnk11 = log(k(11));
lnk11 = - dvc/(r*tk)*pp + 0.5*dkc/(r*tk)*pp*pp + lnk11;
k(11) = exp(lnk11);

k=k*1.0d+03;
k(1) = k(1)*1.0d+03;
k(11) = k(11)*1.0d+03;
k(12) = k(12)*1.0d+03;
