function[pco2,pH,co2,hco3,co3] = co3eq(temp,s,z,alk,dic)
% Function to calculate pCO2 (microatm), bicarbonate ion (micromol/kg), and
% carbonate ion (micromol/kg) concentrations from Alkalinity (micromol/kg) and
% DIC (micromol/kg).  This function requires temperature (degrees C),
% salinity (ppt), and depth (m), and uses the equations in Zeebe and
% Wolf-Gladrow (2000).

t = temp + 273.15;
Pr = z/10;
alk = alk*.000001;
dic = dic*.000001;
R = 83.131;

% Calculate total borate from chlorinity
tbor = .000416 * s / 35.0;

% Calculate Henry's Law coeffieicent (Weiss, 1974)
U1 = -60.2409 + 93.4517 * (100/t) + 23.3585*log(t/100); % note that *log* in matlab is the natural log, ln
U2 = s * (.023517 - .023656 * (t/100) + .0047036 * (t/100)^2);
KH = exp(U1+U2);

% Calculate KB from temp and salinity (Dickson, 1990)
KB = exp((-8966.9 - 2890.53 * s^0.5 - 77.942 * s + 1.728 * s^1.5 - 0.0996 * s^2)/t ...
    + 148.0248 + 137.1942 * s^0.5 + 1.62142 * s - (24.4344 + 25.085 * s^0.5 + ...
    0.2474 * s) * log(t) + 0.053105 * s^0.5 * t);

% Calculate K1 and K2 (Luecker et al., 2000)
K1 = 10^(-(3633.86/t - 61.2172 + 9.67770 * log(t) - 0.011555 * s + 0.0001152 * s^2));
K2 = 10^(-(471.78/t + 25.92990 - 3.16967 * log(t) - 0.01781*s + 0.0001122 * s^2));

% Pressure variation of K1, K2, and KB (Millero, 1995)
dvB = -29.48 + 0.1622 * temp - .002608 * (temp)^2;
dv1 = -25.50 + 0.1271 * temp;
dv2 = -15.82 - 0.0219 * temp;
dkB = -.00284;
dk1 = -.00308 + 0.0000877 * temp;
dk2 = .00113 - .0001475 * temp;
KB = (exp(-(dvB / (R * t)) * Pr + (0.5 * dkB / (R * t)) * Pr^2)) * KB;
K1 = (exp(-(dv1 / (R * t)) * Pr + (0.5 * dk1 / (R * t)) * Pr^2)) * K1;
K2 = (exp(-(dv2 / (R * t)) * Pr + (0.5 * dk2 / (R * t)) * Pr^2)) * K2;

% temperature dependence of KW (DoE, 1994)
KW1 = 148.96502 - 13847.26 / t - 23.65218 * log(t);
KW2 = (118.67 / t - 5.977 + 1.0495 * log(t)) * s^.5 - 0.01615 * s;
KW = exp(KW1 + KW2);

% solve for H ion (Zeebe and Wolf-Gladrow, 2000)
a1 = 1;
a2 = (alk + KB + K1);
a3 = (alk * KB - KB * tbor - KW + alk * K1 + K1 * KB + K1 * K2 - dic * K1);
a4 = (-KW  * KB + alk * KB * K1 - KB * tbor * K1 - KW * K1 + alk * K1 * K2 ...
    + KB * K1 * K2 - dic * KB * K1 - 2 * dic * K1 * K2);
a5 = (-KW * KB * K1 + alk * KB * K1 * K2 - KW * K1 * K2 - KB * tbor * K1 * K2 ...
    - 2 * dic * KB * K1 * K2);
a6 = -KB * KW * K1 * K2;
p = [a1 a2 a3 a4 a5 a6];
r = roots(p);
h = max(real(r));

% calculate bicarbonate, carbonate, and aqueous CO2 using DIC, Alk, and H+
format short g;
hco3 = dic/ (1 + h/K1 + K2/h) * 1e6;
co3 = dic / (1 + h/K2 + h * h / (K1 * K2)) * 1e6;
co2 = dic / (1 + K1/h + K1 * K2 / (h * h)) * 1e6;
pco2 = co2 / KH ;
pH = -log10(h);

% calculate B(OH)4 and OH
BOH4 = KB * tbor / (h + KB);
OH = KW / h;

% recalculate DIC and Alk to check calculations
Ct = (hco3 + co3 + co2) * 1e6;
At = (hco3 + 2*co3 + BOH4 + OH - h) * 1e6;