% _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
%
% 	file: csys.m
%
%
%    	Richard E. Zeebe & Dieter A. Wolf-Gladrow
%
%               			     
%	Alfred Wegener Institute for       
%	Polar and Marine Research		  
%	P.O. Box 12 01 61		       
%	D-27515 Bremerhaven		
%	Germany	      
%	e-mail: rzeebe@awi-bremerhaven.de   wolf@awi-bremerhaven.de  	
%	www     : http://www.awi-bremerhaven.de/
%
% _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
%
% purpose: given two components of CO2 system, calculate all other.
%
%		Compare book by
%
%		Zeebe & Wolf-Gladrow, CO2 in Seawater: 
%		Equilibrium, Kinetics, Isotopes.
% updates:
%
%
%          03.01.01 k -> K
%          01.10.98 include pressure
%          08.05.95 set phflag Total or free
%
% MATLAB 5.3.1
%
% remarks: specify S, T, and P, flag-number for the two variables
%	   pH flag, choose constants (Roy or Mehrbach),
%          and values for two variables.
%
%	   the variable s is equivalent to [CO2]
%
% use:     save file in directory. open matlab. type 'csys'
%
% ---------------------------------------------------------------

%============= INPUT SECTION ===============================%

% flag = 1      h and CO2 given
% flag = 101    h and pco2 given
% flag = 2      CO2 and HCO3 given
% flag = 3      CO2 and CO3 given  
% flag = 4      CO2 and ALK given
% flag = 5      CO2 and DIC given
% flag = 6      h and HCO3 given
% flag = 7      h and CO3 given
% flag = 8      h and ALK given
% flag = 9      h and DIC given
% flag = 10     HCO3 and CO3 given
% flag = 11     HCO3 and ALK given
% flag = 12     HCO3 and DIC given
% flag = 13     CO3 and ALK given
% flag = 14     CO3 and DIC given
% flag = 15     ALK and DIC given

%============= Set flag ======================================%
flag = 15;


%--------------   pH Total or pH free  ------------%
%
%       Specify on which pH scale we are working on.
%       DIC, ALK, pCO2, CO2, CO3, and HCO3 are 
%	physical quantities which do not depend on the pH
%	scale! However, the pH does depend on the scale 
%	used. A certain set of equilibrium
%       constants belongs to a certain pH scale.
%       The scale conversion of the constants is
%       done in equic.m.
%
%       Note that from  [H+]_total > [H+]_free
%       follows           pH_total <   pH_free
%
%	The total scale is recommended.
%

phflag = 0;     % 0: Total scale
                % 1: Free scale

%----- choose K1 and K2: Roy or Mehrbach.

k1k2flag = 0;	% 0: Roy et al. (1993)
		% 1: Mehrbach et al (1973) as 
		%    refit by Lueker et al. (2000) on total scale.


%----- set temperature, salinity and pressure

TC = 25.;
		% deg C
S =  35.;
		
P =  0.;                % bar

%----  load constants
equic;                      
bor = 1.*(416.*(S/35.))* 1.e-6;   % (mol/kg), DOE94

%---- set values for carbonate system
%---- specify two variables according to flag (see above)

ph1   = 8.1;		%  pH
h1    = 10^(-ph1);	%  [H+]
s1    = 10.e-6;		%  [CO2]	(mol/kg)
hco31 = 1750.e-6;	%  [HCO3-]	(mol/kg)
co31  = 250.e-6;	%  [CO3--]	(mol/kg)
alk1  = 2350.e-6;	%  ALK		(mol/kg)
dic1  = 2000.e-6;	%  [DIC]	(mol/kg)
pco21 = 365.e-6;	%  pCO2		(atm)


%----------- These values are test values, set test= 0/1 -------%
test = 0;
if test == 1;
  TC = 25.;
  S = 35.;
  P =  0.;
  equic;
  bor = (416.*(S/35.))* 1.e-6;   % (mol/kg), DOE94   
  ph1   = 8.0902;       %       pH (8.0902 Total, 8.1979 Free)
  h1    = 10^(-ph1);    %       [H+]
  s1    = 10.1304e-6;   %       [CO2]
  hco31 = 1735.9e-06;   %       [HCO3-]
  co31  = 253.9924e-6;  %       [CO3--]
  dic1  = 2000.e-6;     %       [DIC]
  alk1  = 2350.e-6;     %       Alk
  pco21 = 356.8058e-6;	%  	pCO2
end;


%============= END INPUT SECTION ===============================%

%
% rest of file contains the numerical routines.

%==================================================================%

% ----------------- case 1.) h and s given
if flag == 1
disp('flag = 1, pH and CO2 given');
h = h1;
s = s1;
dic = s*(1.+K1/h1+K1*K2/h1/h1);
hco3 = dic/(1+h1/K1+K2/h1);
co3 = dic/(1+h1/K2+h1*h1/K1/K2);
alk = s*(K1/h1+2.*K1*K2/h1/h1)+Kb*bor/(Kb+h1)+Kw/h1-h1;                        
pco2 = s/Kh;
% ----------- change units: mumol/kg
CO2     = 1.e6*s1
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
pco2    = pco2*1.e6
%alktest = (2*co31 + hco31 + bor/(1+h1/Kb) + Kw/h1 - h1)*1.e6
end;
% ----------------- case 101.) h and pco2 given
if flag == 101
disp('flag = 101, pH and pCO2 given');
h = h1;
pco2 = pco21;
s = Kh*pco2;
dic = s*(1.+K1/h1+K1*K2/h1/h1);
hco3 = dic/(1+h1/K1+K2/h1);
co3 = dic/(1+h1/K2+h1*h1/K1/K2);
alk = s*(K1/h1+2.*K1*K2/h1/h1)+Kb*bor/(Kb+h1)+Kw/h1-h1;
% ----------- change units: mumol/kg
CO2     = 1.e6*s                       
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
pco2    = pco2*1.e6
%alktest = (2*co31 + hco31 + bor/(1+h1/Kb) + Kw/h1 - h1)*1.e6
end;
% ------------ s and HCO3 given ------------------
if flag == 2
disp('flag = 2, CO2 and HCO3 given');
s = s1;
hco3 = hco31;
p3 = -hco3/K1;
p2 = s - hco3;
p1 = s*K1 - hco3*K2;
p0 = s*K1*K2;
p = [p3 p2 p1 p0];
r = roots(p);
h = max(real(r));
h*1.e12;
dic = s*(1.+K1/h+K1*K2/h/h);
%       hco3 = dic/(1+h/K1+K2/h); 
co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ s and CO3 given ------------------
if flag == 3
disp('flag = 3, CO2 and CO3 given');
s = s1;
co3 = co31;
p4 = -co3/K1/K2;
p3 = -co3/K2;
p2 = s-co3;
p1 = s*K1;
p0 = s*K1*K2;
p = [p4 p3 p2 p1 p0];      
r = roots(p);
h = max(real(r));
h*1.e12;
dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
%    co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ s and ALK given ------------------
if flag == 4
disp('flag = 4, CO2 and ALK given');
s = s1;
alk = alk1;
p4 = 1.;              
p3 = Kb+alk;
p2 = alk*Kb-s*K1-Kb*bor-Kw;
p1 = -s*Kb*K1-s*2.*K1*K2-Kw*Kb;
p0 = -2.*s*Kb*K1*K2;
p = [p4 p3 p2 p1 p0];
r = roots(p);
h = max(real(r));
h*1.e12;
dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
%    alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ s and DIC given ------------------
if flag == 5
disp('flag = 5, CO2 and DIC given');
s = s1;
dic = dic1;
p2 = dic - s;
p1 = -s*K1;
p0 = -s*K1*K2;
p = [p2 p1 p0];
r = roots(p);
h = max(real(r));
h*1.e12;
%    dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6       
ALK     = alk*1.e6
end
% ------------ h and HCO3 given ------------------
if flag == 6
disp('flag = 6, pH and HCO3 given');
h = h1;
hco3 = hco31;
dic = hco3 * (1+h/K1+K2/h);
s = dic / (1.+K1/h+K1*K2/h/h);
h*1.e12;
%    dic = s*(1.+K1/h+K1*K2/h/h);
%    hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6        
end
% ------------ h and CO3 given ------------------
if flag == 7
disp('flag = 7, pH and CO3 given');
h = h1;
co3 = co31;
dic = co3 * (1+h/K2+h*h/K1/K2);
s = dic / (1.+K1/h+K1*K2/h/h);
h*1.e12;
%    dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
%    co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end                    
% ------------ h and ALK given ------------------
if flag == 8
disp('flag = 8, pH and ALK given');
h = h1;
alk = alk1;
s = (alk-Kw/h+h-Kb*bor/(Kb+h)) / (K1/h+2.*K1*K2/h/h);
h*1.e12;
dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
%    alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ h and DIC given ------------------
if flag == 9
disp('flag = 9, pH and DIC given');
h = h1;
dic = dic1;
s = dic / (1.+K1/h+K1*K2/h/h);
h*1.e12;
%    dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ HCO3 and CO3 given ------------------
if flag == 10
disp('flag = 10, HCO3 and CO3 given');
hco3 = hco31;  
co3 = co31;
p3 = -co3/K1/K2;
p2 = -co3/K2 + hco3/K1;
p1 = -co3 + hco3;
p0 = hco3*K2;
p = [p3 p2 p1 p0];
r = roots(p);
h = max(real(r));
h*1.e12;
dic = hco3 * (1+h/K1+K2/h);
s = dic / (1.+K1/h+K1*K2/h/h);
%    hco3 = dic/(1+h/K1+K2/h);
%    co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6     
end
% ------------ HCO3 and ALK given ------------------
%
%       don't split the lines,
%       matlab does not understand it.
%
if flag == 11
disp('flag = 11, HCO3 and ALK given');
hco3 = hco31;
alk = alk1;
p5  = 1.;
p4  = alk - hco3 + K1 + Kb;
p3  = alk*(Kb+K1)-hco3*(K1+Kb+2.*K2)-Kw+K1*Kb+K1*K2-Kb*bor;
tmp = alk*(Kb*K1+K1*K2)-hco3*((Kb+2.*K2)*K1+2.*Kb*K2+K1*K2);
p2  = tmp +(-K1*Kb*bor-Kw*Kb-K1*Kw+K1*K2*Kb);
tmp = alk*Kb*K1*K2-hco3*(2.*Kb*K1*K2+K2*K1*(Kb+2.*K2));
p1  = tmp +(-K1*K2*Kb*bor-K1*Kw*Kb-K1*K2*Kw);
p0  = -hco3*2.*K2*Kb*K1*K2-K1*K2*Kw*Kb;
p   = [p5 p4 p3 p2 p1 p0];
r   = roots(p);
h   = max(real(r));
h*1.e12;               
dic = hco3 * (1+h/K1+K2/h);
s = dic / (1.+K1/h+K1*K2/h/h);
%   dic = s*(1.+K1/h+K1*K2/h/h);
%   hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ HCO3 and DIC given ------------------
if flag == 12
disp('flag = 12, HCO3 and DIC given');
hco3 = hco31;
dic = dic1;
p2 = hco3/K1;
p1 = hco3-dic;    
p0 = hco3*K2;
p = [p2 p1 p0];
r = roots(p);
h = min(real(r));        % min instead of max !!!!!
h*1.e12;
s = dic / (1.+K1/h+K1*K2/h/h);
dic = s*(1.+K1/h+K1*K2/h/h);
%   hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ CO3 and ALK given ------------------
if flag == 13
disp('flag = 13, CO3 and ALK given');    
co3 = co31;
alk = alk1;
p5  = -co3/K2+1.;
p4  = alk - co3*(K1/K2+(Kb+2.*K2)/K2) + Kb + K1;
tmp = alk*(Kb+K1)-co3*(K1+K1*(Kb+2.*K2)/K2+2.*Kb);
p3  = tmp+(-Kb*bor-Kw+K1*Kb+K1*K2);
tmp = alk*(Kb*K1+K1*K2)-co3*(K1*(Kb+2.*K2)+2.*Kb*K1);
p2  = tmp+(-Kw*Kb-K1*Kb*bor-K1*Kw+K1*K2*Kb);
tmp = alk*Kb*K1*K2-co3*2.*Kb*K1*K2-K1*Kw*Kb;
p1  = tmp+(-K1*K2*Kb*bor-K1*K2*Kw);
p0  = -K1*K2*Kw*Kb;
p   = [p5 p4 p3 p2 p1 p0];
r   = roots(p);
h   = max(real(r));
h*1.e12;
%
% NOTE : calculate dic from dic = co3*(1+h/K2+h^2/K1/K2);
%                  not from dic = hco3*(1+h/K1+K2/h);
%
dic = co3 * (1+h/K2+h^2/K1/K2);
s = dic / (1.+K1/h+K1*K2/h/h);
%   dic = s*(1.+K1/h+K1*K2/h/h);    
hco3 = dic/(1+h/K1+K2/h);
%   co3 = dic/(1+h/K2+h*h/K1/K2);
%   alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1);
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ CO3 and DIC given ------------------
if flag == 14
disp('flag = 14, CO3 and DIC given');
co3 = co31;
dic = dic1;
p2 = co3/K1/K2;
p1 = co3/K2;
p0 = co3-dic;
p = [p2 p1 p0];
r = roots(p);        
h = max(real(r));
h*1.e12;
s = dic / (1.+K1/h+K1*K2/h/h);
%   dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
%   co3 = dic/(1+h/K2+h*h/K1/K2);
alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1); %        6.3096e-09
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6
end
% ------------ ALK and DIC given ------------------
if flag == 15
disp('flag = 15, ALK and DIC given');
alk = alk1;
dic = dic1;
p5  = -1.;        
p4  = -alk-Kb-K1;
p3  = dic*K1-alk*(Kb+K1)+Kb*bor+Kw-Kb*K1-K1*K2;
tmp = dic*(Kb*K1+2.*K1*K2)-alk*(Kb*K1+K1*K2)+Kb*bor*K1;
p2  = tmp+(Kw*Kb+Kw*K1-Kb*K1*K2);
tmp = 2.*dic*Kb*K1*K2-alk*Kb*K1*K2+Kb*bor*K1*K2;
p1  = tmp+(+Kw*Kb*K1+Kw*K1*K2);
p0  = Kw*Kb*K1*K2;
p   = [p5 p4 p3 p2 p1 p0];
r   = roots(p);
h   = max(real(r));
%   test = p5*h^5+p4*h^4+p3*h^3+p2*h^2+p1*h+p0
h*1.e12;
s = dic / (1.+K1/h+K1*K2/h/h);
%   dic = s*(1.+K1/h+K1*K2/h/h);
hco3 = dic/(1+h/K1+K2/h);
co3 = dic/(1+h/K2+h*h/K1/K2);
%   alk = s*(K1/h+2.*K1*K2/h/h)+Kb*bor/(Kb+h)+Kw/h-h;
% ----------- change units: mumol/kg
%h1 = 10^(-ph1); %
CO2     = s*1.e6
pCO2    = s*1.e6/Kh
HCO3    = hco3*1.e6     
CO3     = co3*1.e6
DIC     = dic*1.e6
ALK     = alk*1.e6

end
;

% ======================================================
if phflag == 0;
   disp('Total pH scale used.');
   phtotal = -log10(h)
   oh = Kw./h;
   phfree  = -log10(h/total2free)
end;

if phflag == 1;
   disp('Free pH scale used.');
   phfree  = -log10(h)
   oh = Kw./h;   
   phtotal = -log10(h*total2free)
 end;

return;

