function carbon_coeffs(t,s)
global k0 k1 k2 kw kb ks kf k1p k2p k3p ksi ff htotal
global bt st ft

% This routine computes various constants as functions of T and S 
% required by the carbon chemistry solver. These are passed to the 
% chemistry solver by declaring them as global variables. (In the 
% OCMIP fortran code, these were in a common block.)

% Samar Khatiwala (spk@ldeo.columbia.edu)

% SPK 11/8/05: - comments added
%              - deleted some variables (sit pt dic ta) from global declaration

%-----------------------------------------------------------------
% 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
% Note: mol/kg are actually what the body of this routine uses 
% for calculations.  
%---------------------------------------------------------------------

%*************************************************************************
% Calculate all constants needed to convert between various measured
% carbon species. References for each equation are noted in the code. 
% Once calculated, the constants are
% stored and passed in the common block "const". The original version of this
% code was based on the code by Dickson in 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 = 273.15 + t;
tk100 = tk/100.0;
tk1002=tk100*tk100;
invtk=1.0/tk;
dlogtk=log(tk);
is=19.924*s/(1000.-1.005*s);
is2=is*is;
sqrtis=sqrt(is);
s2=s*s;
sqrts=sqrt(s);
s15=s^1.5;
scl=s/1.80655;
%
% 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));
%
% K0 from Weiss 1974
%
k0 = exp(93.4517/tk100 - 60.2409 + 23.3585 * log(tk100) + ...
		s * (.023517 - 0.023656 * tk100 + 0.0047036 * tk1002));

%
% k1 = [H][HCO3]/[H2CO3]
% k2 = [H][CO3]/[HCO3]
%
% Millero p.664 (1995) using Mehrbach et al. data on seawater scale 
%
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));
%
% kb = [H][BO2]/[HBO2]
%
% Millero p.669 (1995) using data from Dickson (1990)
%
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);
%
% k1p = [H][H2PO4]/[H3PO4]
%
% DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
%
k1p = exp(-4576.752*invtk + 115.525 - 18.453 * dlogtk + ...
		(-106.736*invtk + 0.69171) * sqrts + ...
		(-0.65643*invtk - 0.01844) * s);
%
% k2p = [H][HPO4]/[H2PO4]
%
% DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
%
k2p = exp(-8814.715*invtk + 172.0883 - 27.927 * dlogtk + ...
		(-160.340*invtk + 1.3566) * sqrts + ...
		(0.37335*invtk - 0.05778) * s);
%
%------------------------------------------------------------------------
% k3p = [H][PO4]/[HPO4]
%
% DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
%
k3p = exp(-3070.75*invtk - 18.141 + ...
		(17.27039*invtk + 2.81197) * ...
		sqrts + (-44.99486*invtk - 0.09984) * s);
%
%------------------------------------------------------------------------
% ksi = [H][SiO(OH)3]/[Si(OH)4]
%
% Millero p.671 (1995) using data from Yao and Millero (1995)
%
ksi = exp(-8904.2*invtk + 117.385 - 19.334 * dlogtk + ...
		(-458.79*invtk + 3.5913) * sqrtis + ...
		(188.74*invtk - 1.5998) * is + ...
		(-12.1652*invtk + 0.07871) * is2 + ...
		log(1.0-0.001005*s));
%
%------------------------------------------------------------------------
% kw = [H][OH]
%
% Millero p.670 (1995) using composite data
%
kw = exp(-13847.26*invtk + 148.9652 - 23.6521 * dlogtk + ...
		(118.67*invtk - 5.977 + 1.0495 * dlogtk) * ...
		sqrts - 0.01615 * s);
%
%------------------------------------------------------------------------
% ks = [H][SO4]/[HSO4]
%
% Dickson (1990, J. chem. Thermodynamics 22, 113)
%
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 - 0.001005*s));
%
%------------------------------------------------------------------------
% kf = [H][F]/[HF]
%
% Dickson and Riley (1979) -- change pH scale to total
%
kf=exp(1590.2*invtk - 12.641 + 1.525*sqrtis + ...
		log(1.0 - 0.001005*s) + ...
		log(1.0 + (0.1400/96.062)*(scl)/ks));
%
%------------------------------------------------------------------------
% 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;
