function [dic,co2star,ph,dco2ddic]=diccalc(t,s,pco2,alk,pt,sit,flag)

% [dic,co2star,ph,dco2ddic]=diccalc(t,s,pco2,alk,pt,sit,flag)
% 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. 
% This is the 'fast' version of diccalc. If flag=1 or missing, it first calls 
% carbon_coeffs(t,s) to precompute various constants. If flag=0, it is assumed 
% that carbon_coeffs(t,s) has already been called.
% 
% USAGE: [dic,co2star,ph,dco2ddic]=diccalc(t,s,pco2,alk,pt,sit,[flag])
% INPUTS: t = temperature in deg C
%         s = salinity in psu
%         pCO2 = surface ocean pCO2 in ppm (uatm (same as mole fraction in ppm))
%         ta = total alkalinity in eq/m^3
%	      pt  = inorganic phosphate in mol/m^3
%	      sit = inorganic silicate in mol/m^3 
%         flag = 0 - don't calculate coeffs; 1 - calculate coeffs internally
% 
% 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 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.
%    This is obviously a hack!

% SPK 11/8/05: - added comments

global ff % solubility in mol/(kg*atm) <- computed in carbon_coeffs

if nargin<7
  flag=1;
end

tol=1e-10;
phlo=6;
phhi=9;
% first guess
dic=2000*1e-6*1024.5; % dic in mol/m^3 (units expected by co2calc)
% precompute
if flag==1
  carbon_coeffs(t,s)
end
% Target value of [CO2]aq.
co2star0=(pco2*1e-6)*(ff*1024.5); % convert pco2 [uatm] to [CO2]aq in mol/m^3
% Compute residual using first guess [DIC]
% Note: The last 2 arguments of co2calc_fast (pCO2 and atmpres) are arbitrary. They 
% are only used to compute the output arguments dco2star and dpco2 and don't influence 
% the chemistry calculation. We don't use dco2star and dpco2 here so their (incorrect) 
% values don't matter.
% ?? Why not pass the given pCO2 value ??
[co2star,dco2star,pCO2,dpco2,ph,dco2ddic]=co2calc_fast(dic,alk,pt,sit,phlo,phhi,280,1);
d=co2star-co2star0; % residual
while abs(d)>tol % iterate until converged
  [co2star,dco2star,pCO2,dpco2,ph,dco2ddic]=co2calc_fast(dic,alk,pt,sit,phlo,phhi,280,1);
  d=co2star-co2star0;
  dic = dic -d/dco2ddic;
end
[co2star,dco2star,pCO2,dpco2,ph,dco2ddic]=co2calc_fast(dic,alk,pt,sit,phlo,phhi,280,1);
