function [dic,co2star,pH,dco2ddic]=diccalc_mex(t,s,pco2,alk,pt,sit)

% [dic,co2star,pH,dco2ddic]=diccalc_mex(t,s,pco2,alk,pt,sit)
% This routine computes DIC, [CO2]aq, ph, and d[CO2]aq/d[DIC] as a function of 
% temperature, salinity, surface ocean pCO2, total alkalinity, inorganic phosphate, 
% and inorganic silicate. 
% 
% USAGE: [dic,co2star,pH,dco2ddic]=diccalc_mex(t,s,pco2,alk,pt,sit)
% INPUTS: t = temperature in deg C
%         s = salinity in psu
%         pCO2 = surface ocean pCO2 in ppm (uatm (same as mole fraction in ppm))
%         alk = total alkalinity in eq/m^3
%	      pt  = inorganic phosphate in mol/m^3
%	      sit = inorganic silicate in mol/m^3 
% 
% OUTPUTS: dic     = total inorganic carbon (mol/m^3)
%          co2star = [CO2]aq (mol/m^3)
%          ph
%          dco2ddic= approximate d([CO2]aq)/d([DIC)

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

% NOTES:
%  - Use this routine if you want to compute [DIC], [CO2]aq, etc in equilibrium 
%    with a known OCEAN pCO2 value.
%  - Use this routine if you want to compute [DIC], [CO2]aq etc in equlilibrium 
%    with a given ATMOSPHERIC pCO2. ([CO2]aq will then be the saturation  
%    [CO2]aq concentration.) The OCEAN pCO2 to be passed to diccalc should be 
%    = pCO2_atm*P/Po, where pCO2_atm is the atmospheric mole fraction of CO2 in 
%    dry air (ppmv), P is the (local) total air pressure at sea level (atm), and 
%    Po is 1 atm. (If you don't know P, set pCO2 = pCO2_atm.)
%  - This function uses co2calc_mex to iterate toward the solution. 
%    Newton-Raphson is used to solve: F(DIC) = pco2(DIC)-pco2_given = 0
%    To make life easier (since what is known is the (approximate) 
%    derivative of [CO2]aq w.r.t. DIC - returned by pco2calc), we solve
%    F(DIC) = co2star(DIC) - co2star_given = 0, where co2star_given
%    is computed from pco2_given via the solubility constant.

% SPK 04/21/06: - created

tol=1e-10;
phlo=6;
phhi=9;
% first guess
dic=repmat(2000*1e-6*1024.5,size(t)); % dic in mol/m^3 (units expected by co2calc)
% Target value of [CO2]aq.
co2star0=pco2.*co2solubility_mex(t,s); % convert pco2 [uatm] to [CO2]aq in mol/m^3
% Compute residual using first guess [DIC]
% Note: The pco2 argument of co2calc_mex is arbitrary. It is only used to 
% compute the output argument co2starair and doesn't influence 
% the chemistry calculation. We don't use co2starair here.
[co2star,co2starair,pCO2surf,dco2ddic,pH]=co2calc_mex(t,s,dic,alk,pt,sit,pco2,phlo,phhi);
d=co2star-co2star0; % residual
while max(abs(d(:)))>tol % iterate until converged
  [co2star,co2starair,pCO2surf,dco2ddic,pH]=co2calc_mex(t,s,dic,alk,pt,sit,pco2,phlo,phhi);
  d=co2star-co2star0;
  dic = dic - d./dco2ddic;
end
[co2star,co2starair,pCO2surf,dco2ddic,pH]=co2calc_mex(t,s,dic,alk,pt,sit,pco2,phlo,phhi);
