% 
function [co2star,dco2star,pCO2surf,dpCO2,ph,dco2ddic]=co2calc(t,s,dic_in,ta_in,pt_in,...
                    sit_in,phlo,phhi,xco2_in,atmpres)

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

%
%-------------------------------------------------------------------------
% SUBROUTINE CO2CALC
%
% PURPOSE
%	Calculate delta co2* from total alkalinity and total CO2 at
% temperature (t), salinity (s) and "atmpres" atmosphere total pressure. 
%
% USAGE
%       call co2calc(t,s,dic_in,ta_in,pt_in,sit_in
%    &                  ,phlo,phhi,ph,xco2_in,atmpres
%    &                  ,co2star,dco2star,pCO2surf,dpco2)
%
% INPUT
%	dic_in = total inorganic carbon (mol/m^3) 
%                where 1 T = 1 metric ton = 1000 kg
%	ta_in  = total alkalinity (eq/m^3) 
%	pt_in  = inorganic phosphate (mol/m^3) 
%	sit_in = inorganic silicate (mol/m^3) 
%	t      = temperature (degrees C)
%	s      = salinity (PSU)
%	phlo   = lower limit of pH range
%	phhi   = upper limit of pH range
%	xco2_in=atmospheric mole fraction CO2 in dry air (ppmv) 
%	atmpres= atmospheric pressure in atmospheres (1 atm==1013.25mbar)
%
%       Note: arguments dic_in, ta_in, pt_in, sit_in, and xco2_in are 
%             used to initialize variables dic, ta, pt, sit, and xco2.
%             * Variables dic, ta, pt, and sit are in the common block 
%               "species".
%             * Variable xco2 is a local variable.
%             * Variables with "_in" suffix have different units 
%               than those without.

% OUTPUT
%	co2star  = CO2*water (mol/m^3)
%	dco2star = delta CO2 (mol/m^3)
%       pco2surf = oceanic pCO2 (ppmv)
%       dpco2    = Delta pCO2, i.e, pCO2ocn - pCO2atm (ppmv)
%       ph
%       dco2ddic = approximate d([CO2]aq)/d([DIC)
%
% IMPORTANT: Some words about units - (JCO, 4/4/1999)
%     - Models carry tracers in mol/m^3 (on a per volume basis)
%     - Conversely, this routine, which was written by observationalists 
%       (C. Sabine and R. Key), passes input arguments in umol/kg  
%       (i.e., on a per mass basis)
%     - I have changed things slightly so that input arguments are in mol/m^3,
%     - Thus, all input concentrations (dic_in, ta_in, pt_in, and st_in) 
%       should be given in mol/m^3; output arguments "co2star" and "dco2star"  
%       are likewise in mol/m^3.

% FILES and PROGRAMS NEEDED
%	drtsafe
%	ta_iter_1
%
%--------------------------------------------------------------------------
%
global k0 k1 k2 kw kb ks kf k1p k2p k3p ksi ff htotal
global bt st ft sit pt dic ta
%

%---------------------------------------------------------------------
% 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.  
%---------------------------------------------------------------------
permil = 1.0 / 1024.5;
% To convert input in mol/m^3 -> mol/kg 
pt=pt_in*permil;
sit=sit_in*permil;
ta=ta_in*permil;
dic=dic_in*permil;

%       ---------------------------------------------------------------------
%       Change units from uatm to atm. That is, atm is what the body of 
%       this routine uses for calculations.
%       ---------------------------------------------------------------------
permeg=1.e-6;
% To convert input in uatm -> atm
xco2=xco2_in*permeg;
%       ---------------------------------------------------------------------
%
%*************************************************************************
% 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;
%*************************************************************************
%
% Calculate [H+] total when DIC and TA are known at T, S and 1 atm.
% 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 htotal. 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).
%
% Parentheses added around negative exponents (Keith Lindsay)
%
x1 = 10.0^(-phhi);
x2 = 10.0^(-phlo);
%xacc = 1.e-10;
xacc = 1.e-12; % SPK setting
htotal = drtsafe(@ta_iter_1,x1,x2,xacc);
%
% 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)
%
htotal2=htotal*htotal;
co2star=dic*htotal2/(htotal2 + k1*htotal + k1*k2);
co2starair=xco2*ff*atmpres;
dco2star=co2starair-co2star;
ph=-log10(htotal);

%
%---------------------------------------------------------------
% Add two output arguments for storing pCO2surf
% Should we be using K0 or ff for the solubility here?
%---------------------------------------------------------------
pCO2surf = co2star / ff;
%sprintf('%10.7f\n',ff)
dpCO2    = pCO2surf - xco2*atmpres;

dco2ddic = (htotal/k2)/(dic/co2star); % approximate d(co2star)/d(dic)
%  Convert units of output arguments
%      Note: co2star and dco2star are calculated in mol/kg within this routine 
%      Thus Convert now from mol/kg -> mol/m^3
co2star  = co2star / permil;
dco2star = dco2star / permil;
% Note: pCO2surf and dpCO2 are calculated in atm above. 
% Thus convert now to uatm
pCO2surf = pCO2surf / permeg;
dpCO2    = dpCO2 / permeg;

%-------------------------

function drtsafe_out=drtsafe(funcd,x1,x2,xacc)
%
%	File taken from Numerical Recipes. Modified  R.M.Key 4/94
%
maxit=100;
[fl,df]=funcd(x1);
[fh,df]=funcd(x2);
if (fl < 0.0) 
  xl=x1;
  xh=x2;
else
  xh=x1;
  xl=x2;
  swap=fl;
  fl=fh;
  fh=swap;
end
drtsafe_out=.5*(x1+x2);
dxold=abs(x2-x1);
dx=dxold;
[f,df]=funcd(drtsafe_out);
for j=1:maxit
  if (((drtsafe_out-xh)*df-f)*((drtsafe_out-xl)*df-f) >= 0.0 | ...
      abs(2.0*f) > abs(dxold*df)) 
    dxold=dx;
    dx=0.5*(xh-xl);
    drtsafe_out=xl+dx;
    if (xl == drtsafe_out) 
      return
    end
  else
    dxold=dx;
    dx=f/df;
    temp=drtsafe_out;
    drtsafe_out=drtsafe_out-dx;
    if (temp == drtsafe_out)
      return
    end
  end
  if(abs(dx) < xacc)
    return
   end
   [f,df]=funcd(drtsafe_out);
   if (f < 0.0)
     xl=drtsafe_out;
     fl=f;
   else
     xh=drtsafe_out;
     fh=f;
   end
end

function [fn,df]=ta_iter_1(x)

%
% This routine expresses TA as a function of DIC, htotal and constants.
% It also calculates the derivative of this function with respect to 
% htotal. It is used in the iterative solution for htotal. In the call
% "x" is the input value for htotal, "fn" is the calculated value for TA
% and "df" is the value for dTA/dhtotal
%
global k0 k1 k2 kw kb ks kf k1p k2p k3p ksi ff htotal
global bt st ft sit pt dic ta

x2=x*x;
x3=x2*x;
k12 = k1*k2;
k12p = k1p*k2p;
k123p = k12p*k3p;
c = 1.0 + st/ks;
a = x3 + k1p*x2 + k12p*x + k123p;
a2=a*a;
da = 3.0*x2 + 2.0*k1p*x + k12p;
b = x2 + k1*x + k12;
b2=b*b;
db = 2.0*x + k1;
%
% fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree+hso4+hf+h3po4-ta
%
fn = k1*x*dic/b + ...
 	     2.0*dic*k12/b + ...
 	     bt/(1.0 + x/kb) + ...
 	     kw/x + ...
 	     pt*k12p*x/a + ...
 	     2.0*pt*k123p/a + ...
 	     sit/(1.0 + x/ksi) - ...
 	     x/c - ...
 	     st/(1.0 + ks/x/c) - ...
 	     ft/(1.0 + kf/x) - ...
 	     pt*x3/a - ...
 	     ta;
%
%	df = dfn/dx
%
df = ((k1*dic*b) - k1*x*dic*db)/b2 	- ...
	     2.0*dic*k12*db/b2 - ...
 	     bt/kb/(1.0+x/kb)^2. - ...
 	     kw/x2 + ...
 	     (pt*k12p*(a - x*da))/a2 - ...
 	     2.0*pt*k123p*da/a2 - ...
 	     sit/ksi/(1.0+x/ksi)^2. - ...
 	     1.0/c + ...
 	     st*(1.0 + ks/x/c)^(-2.0)*(ks/c/x2)  + ...
 	     ft*(1.0 + kf/x)^(-2.)*kf/x2 - ...
 	     pt*x2*(3.0*a-x*da)/a2;
