% _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
%
% 	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:
%
%
%          10.24.04 user interface implemented
%                   data reading from file external implemented
%          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'
%
% ---------------------------------------------------------------

clear all;
runUserInterface = 1;
    %  = 0 -> will NOT run user interface
    % >= 1 -> will run user interface (get file data not functioning)
% ---------------------------------------------------------------
% MULTIPLE INPUT VALUES SECTION
%
% To use multiple input values, change the "runUserInterface" 
% variable (above) to 0 (Zero) and place the appropriate values 
% in the variables below.
%
%
% 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
% ---------------------------------------------------------------
if(runUserInterface == 0)
    numOfInputs = 3;                % Number of values to calculate
    flag = 11;                       % Given Variables (see chart)
    phflag = 0;                     % Total Scale (0) or Free Scale (1)
    k1k2flag = 1;                   % Roy (0) or Mehrbach (1)
    TC = 25;                        % Temperature   (Celcius)
    S =  35;                        % Salinity
    P =  0;                         % Pressure      (bar)
    ph1   = [8.2, 8.3, 8.4];        % pH
    s1    = [10, 11, 12];           % [CO2]         (µmol/kg)
    hco31 = [1700, 1750, 1800];     % [HCO3-]       (µmol/kg)
    co31  = [200, 250, 300];        % [CO3--]       (µmol/kg)
    alk1  = [2350, 2400, 2450];     % ALK           (µmol/kg)
    dic1  = [1950, 2000, 2050];     % [DIC]         (µmol/kg)
    pco21 = [355, 365, 375];        % pCO2          (µatm)
end





% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% ==================== END OF INPUT SECTION =====================
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++




debug = 1;
if(runUserInterface >= 1)
    fprintf('\n');
% ---------------------------------------------------------------
% Variables From File? section
%   - If data is manually inputted, runUserInterface = 1
%   - If data is read from a file, runUserInterface = 2
% ---------------------------------------------------------------       
    valid = 0;
    while(valid == 0)
        from_file = input('Will data be retrieved from a file [Y/N]? ', 's');
        
        % Error checking\
        if(~isempty(from_file))
            if(from_file == 'n' | from_file == 'N')
                runUserInterface = 1;
                valid = 1;
            elseif(from_file == 'y' | from_file == 'Y')
                runUserInterface = 2;
                valid = 1;
            else
                fprintf('\n\n*** INVALID INPUT ***\n\n');            
            end
        else
            fprintf('Default Value Taken: [N]\n\n');
            runUserInterface = 1;
            valid = 1;
        end
    end
    
    
    
% ---------------------------------------------------------------
% Select Values Given section
%   - flag will be set with a numeric value either corresponding 
%     to the users input or the default value
% ---------------------------------------------------------------   
    valid = 0;
    while(valid == 0)
        fprintf('\nSelect Values Given\n');
        fprintf('----------------------\n');
        fprintf('1 - pH and C02 given\t\t[Default]\n');
        fprintf('101 - pH and pCO2 given\n');
        fprintf('2 - CO2 and HCO3 given\n');
        fprintf('3 - CO2 and CO3 given\n');
        fprintf('4 - CO2 and ALK given\n');
        fprintf('5 - CO2 and DIC given\n');
        fprintf('6 - pH and HCO3 given\n');
        fprintf('7 - pH and CO3 given\n');
        fprintf('8 - pH and ALK given\n');
        fprintf('9 - pH and DIC given\n');
        fprintf('10 - HCO3 and CO3 given\n');
        fprintf('11 - HCO3 and ALK given\n');
        fprintf('12 - HCO3 and DIC given\n');
        fprintf('13 - CO3 and ALK given\n');
        fprintf('14 - CO3 and DIC given\n');
        fprintf('15 - ALK and DIC given\n');
        
        % Gets user input.  If the user does not enter a selection, the
        % default (1) will be used.
        str_flag = input('Selection: ', 's');
        if(~isempty(str_flag))
            flag = str2num(str_flag);
            % Error checking
            if(isnumeric(flag))
                if((flag >= 1 & flag <= 15) | flag == 101)
                    valid = 1;
                else
                    fprintf('\n\n*** INVALID INPUT ***\n\n');
                end
            else
                fprintf('\n\n*** INVALID INPUT ***\n\n');              
            end
        else
            fprintf('Default Value Taken: [1]\n\n');
            flag = 1;
            valid = 1;
        end
    end





% ---------------------------------------------------------------
% Scale Section
%   - phflag will be set with a numeric value either 
%     corresponding to the users input or the default value
% ---------------------------------------------------------------
    %--------------   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.
    %--------------------------------------------------%
    valid = 0;
    while(valid == 0)
        fprintf('\nSelect Scale\n');
        fprintf('1 - Total scale\t\t[Default]\n');
        fprintf('2 - Free scale\n');
        
        % Gets user input.  If the user does not enter a selection, the
        % default (0) will be used.
        str_phflag = input('Selection: ', 's');
        if(~isempty(str_phflag))
            phflag = str2num(str_phflag) - 1;
            % Error checking
            if(isnumeric(phflag))
                if(phflag == 0 | phflag == 1)
                    valid = 1;
                else
                    fprintf('\n\n*** INVALID INPUT ***\n\n');
                end
            else
                fprintf('\n\n*** INVALID INPUT ***\n\n');
            end
        else
            fprintf('Default Value Taken [1]\n\n');
            phflag = 0;
            valid = 1;
        end
    end
    
    



% ---------------------------------------------------------------
% Equilibrium Constants selection
%   - k1k2flag will be set with a numeric value either 
%     corresponding to the users input or the default value
% ---------------------------------------------------------------
    valid = 0;
    while(valid == 0)
        fprintf('\nSelect Equilibrium Constants\n');
        fprintf('1 - Roy et al. (1993)\n');
        fprintf('2 - Mehrbach et al (1973) as refit by Lueker et al. (2000)\n');
        fprintf('    on total scale (Recommended)\n');
        fprintf('    [Default]\n');
        
        % Gets user input.  If the user does not enter a selection, the
        % default (1) will be used.
        str_k1k2flag = input('Selection: ', 's');
        if(~isempty(str_k1k2flag))
            k1k2flag = str2num(str_k1k2flag) - 1;
            % Error checking
            if(isnumeric(k1k2flag))
                if(k1k2flag == 0 | k1k2flag == 1)
                    valid = 1;
                else
                    fprintf('\n\n*** INVALID INPUT ***\n\n');
                end
            else
                fprintf('\n\n*** INVALID INPUT ***\n\n');
            end
        else
            fprintf('Default Value Taken [2]\n\n');
            k1k2flag = 1;
            valid = 1;
        end
    end





% ---------------------------------------------------------------
% Temperature, salinity, and pressure section
%   - TC, S, & P will be set with a numeric value either 
%     corresponding to the users input(s) or the default value(s)
% ---------------------------------------------------------------   
    % Gets user input.  If the user does not enter a selection, the
    % default (25) will be used.
    valid = 0;
    while(valid == 0)
        str_TC = input('\nEnter temperature (Celcius): ', 's');
        if(~isempty(str_TC))
            TC = str2num(str_TC);
            if(isnumeric(TC))
                valid = 1;
            else
                fprintf('\n\n*** INVALID INPUT ***\n\n');
            end
        else
            fprintf('Default Value Taken [25]\n\n');            
            TC = 25;
            valid = 1;
        end
    end

    % Gets user input.  If the user does not enter a selection, the
    % default (35) will be used.
    valid = 0;
    while(valid == 0)
        str_S = input('Enter salinity: ', 's');
        if(~isempty(str_S))
            S = str2num(str_S);
            if(isnumeric(S))
                valid = 1;
            else
                fprintf('\n\n*** INVALID INPUT ***\n\n');
            end                
        else
            fprintf('Default Value Taken [35]\n\n');
            S = 35;
            valid = 1;
        end
    end
    
    % Gets user input.  If the user does not enter a selection, the
    % default (0) will be used.
    valid = 0;
    while(valid == 0)
        str_P = input('Enter pressure (bar): ', 's');
        if(~isempty(str_P))
            P = str2num(str_P);
            if(isnumeric(P))
                valid = 1;
            else
                fprintf('\n\n*** INVALID INPUT ***\n\n');
            end                
        else
            fprintf('Default Value Taken [0]\n\n');
            P = 0;
            valid = 1;
        end
    end
    
    
    
    
    
    
    numOfInputs = 1;
    i = 1;
% ---------------------------------------------------------------
% Input File section
%   - filename will be set with a string value from the user's
%     input
%   - delim_type will be set with a character value from the 
%     user's input
% ---------------------------------------------------------------   
    if(runUserInterface == 2)
        valid = 0;
        while(valid == 0)
            filename = input('Enter Input Filename: ', 's');
            
            if(isempty(filename))
                 fprintf('\n\n*** FILE NAME NEEDED ***\n\n');
            else
                valid = 1;
            end
        end
        
        valid = 0;
        while(valid == 0)
            fprintf('\nDelimiter Used\n');
            fprintf('--------------\n');
            fprintf('1 - Comma [Default]\n');
            fprintf('2 - Space\n');
            fprintf('3 - Tab\n');
            
            delim_type = input('Selection: ', 's');
            
            % Error checking
            if(~isempty(delim_type))
                if(delim_type < '1' | delim_type > '3')
                    fprintf('\n\n*** INVALID INPUT ***\n\n');
                else
                    valid = 1;
                end
            else
                fprintf('Default Value Taken [1]\n\n');
                delim_type = '1';
                valid = 1;
            end
        end
    
        % Reads from a file depending on the specified delimiter
        if(delim_type == '1')
            input_data = csvread(filename);
        elseif(delim_type == '2')
            input_data = dlmread(filename, ' ');
        else
            input_data = dlmread(filename, '\t');
        end
            
        numOfInputs = length(input_data);
    end
    
    
    

    
    while(i <= numOfInputs)
        % ---------------------------------------------------------------
        % Variables based on input from Select Variables Given section
        % ("flag" variable)
        % ---------------------------------------------------------------
    
        % Gets user input.  If the user does not enter a selection, the
        % default (8.2) will be used.
        if(flag == 1 | flag == 101 | flag == 6 | flag == 7 | flag == 8 | flag == 9)
            if(runUserInterface == 1)
                valid = 0;
                while(valid == 0)
                    str_ph1 = input('Enter pH: ', 's');
                    if(~isempty(str_ph1))
                        ph1 = str2num(str_ph1);
                        if(isnumeric(ph1))
                            valid = 1;
                        else
                            fprintf('\n\n*** INVALID INPUT ***\n\n');
                        end                
                    else
                        fprintf('Default Value Taken [8.2]\n\n');
                        ph1 = 8.2;
                        valid = 1;
                    end
                end
            else
                ph1(i) = input_data(i, 1);
            end
        end
        
        
        % Gets user input.  If the user does not enter a selection, the
        % default (10) will be used.
        if(flag == 1 | flag == 2 | flag == 3 | flag == 4 | flag == 5)
            if(runUserInterface == 1)
                valid = 0;
                while(valid == 0)
                    str_s1 = input('Enter CO2 (µmol/kg): ', 's');
                    if(~isempty(str_s1))
                        s1 = str2num(str_s1);
                        if(isnumeric(s1))
                            valid = 1;
                        else
                            fprintf('\n\n*** INVALID INPUT ***\n\n');
                        end                
                    else
                        fprintf('Default Value Taken [10]\n\n');
                        s1 = 10;
                        valid = 1;
                    end
                end
            elseif(flag > 1)
                s1(i) = input_data(i, 1);
            else
                s1(i) = input_data(i, 2);
            end               
        end
        
        
        % Gets user input.  If the user does not enter a selection, the
        % default (1750) will be used.
        if(flag == 2 | flag == 6 | flag == 10 | flag == 11 | flag == 12)
            if(runUserInterface == 1)
                valid = 0;
                while(valid == 0)
                    str_hco31 = input('Enter HCO3- (µmol/kg): ', 's');
                    if(~isempty(str_hco31))
                        hco31 = str2num(str_hco31);
                        if(isnumeric(hco31))
                            valid = 1;
                        else
                            fprintf('\n\n*** INVALID INPUT ***\n\n');
                        end                
                    else
                        fprintf('Default Value Taken [1750]\n\n');
                        hco31 = 1750;
                        valid = 1;
                    end
                end
            elseif(flag >= 10)
                hco31(i) = input_data(i, 1);
            else
                hco31(i) = input_data(i, 2);
            end
        end
        
        
        % Gets user input.  If the user does not enter a selection, the
        % default (250) will be used.
        if(flag == 3 | flag == 7 | flag == 10 | flag == 13 | flag == 14)
            if(runUserInterface == 1)
                valid = 0;
                while(valid == 0)
                    str_co31 = input('Enter CO3-- (µmol/kg): ', 's');
                    if(~isempty(str_co31))
                        co31 = str2num(str_co31);
                        if(isnumeric(co31))
                            valid = 1;
                        else
                            fprintf('\n\n*** INVALID INPUT ***\n\n');
                        end                
                    else
                        fprintf('Default Value Taken [250]\n\n');
                        co31 = 250.0;
                        valid = 1;
                    end
                end
            elseif(flag >= 13)
                co31(i) = input_data(i, 1);
            else
                co31(i) = input_data(i, 2);
            end
        end
        
        
        % Gets user input.  If the user does not enter a selection, the
        % default (2400) will be used.
        if(flag == 4 | flag == 8 | flag == 11 | flag == 13 | flag == 15)
            if(runUserInterface == 1)
                valid = 0;
                while(valid == 0)
                    str_alk1 = input('Enter ALK (µmol/kg): ', 's');
                    if(~isempty(str_alk1))
                        alk1 = str2num(str_alk1);
                        if(isnumeric(alk1))
                            valid = 1;
                        else
                            fprintf('\n\n*** INVALID INPUT ***\n\n');
                        end                
                    else
                        fprintf('Default Value Taken [2400]\n\n');
                        alk1 = 2400;
                        valid = 1;
                    end
                end
            elseif(flag == 15)
                alk1(i) = input_data(i, 1);
            else
                alk1(i) = input_data(i, 2);
            end
        end
        
        
        % Gets user input.  If the user does not enter a selection, the
        % default (2000) will be used.
        if(flag == 5 | flag == 9 | flag == 12 | flag == 14 | flag ==  15)
            if(runUserInterface == 1)
                valid = 0;
                while(valid == 0)
                    str_dic1 = input('Enter DIC (µmol/kg): ', 's');
                    if(~isempty(str_dic1))
                        dic1 = str2num(str_dic1);
                        if(isnumeric(dic1))
                            valid = 1;
                        else
                            fprintf('\n\n*** INVALID INPUT ***\n\n');
                        end                
                    else
                        fprintf('Default Value Taken [8.2]\n\n');
                        dic1 = 2000;
                        valid = 1;
                    end
                end
            else
                dic1(i) = input_data(i, 2);
            end
        end
        
        
        % Gets user input.  If the user does not enter a selection, the
        % default (365) will be used.
        if(flag == 101)
            if(runUserInterface == 1)
                valid = 0;
                while(valid == 0)
                    str_pco21 = input('Enter pCO2 (µatm): ', 's');
                    if(~isempty(str_pco21))
                        pco21 = str2num(str_pco21);
                        if(isnumeric(pco21))
                            valid = 1;
                        else
                            fprintf('\n\n*** INVALID INPUT ***\n\n');
                        end                
                    else
                        fprintf('Default Value Taken [365]\n\n');
                        pco21 = 365;
                        valid = 1;
                    end
                end
            else
                pco21(i) = input_data(i, 2);
            end
        end
        
        
        if(runUserInterface == 1)
            fprintf('\n');
        end
        
        i = i + 1;
    end
    
    fclose('all');
end





% ---------------------------------------------------------------
% µmol/kg to mol/kg and µatm to atm conversion and other 
% declarations
% ---------------------------------------------------------------
if(flag == 1 | flag == 2 | flag == 3 | flag == 4 | flag == 5)
    s1 = s1 / 1000000;
end

if(flag == 2 | flag == 6 | flag == 10 | flag == 11 | flag == 12)
    hco31 = hco31 / 1000000;
end

if(flag == 3 | flag == 7 | flag == 10 | flag == 13 | flag == 14)
    co31 = co31 / 1000000;
end

if(flag == 4 | flag == 8 | flag == 11 | flag == 13 | flag == 15)
    alk1 = alk1 / 1000000;
end

if(flag == 5 | flag == 9 | flag == 12 | flag == 14 | flag ==  15)
    dic1 = dic1 / 1000000;
end

if(flag == 101)
    pco21 = pco21 / 1000000;
end

if(flag == 1 | flag == 101 | flag == 6 | flag == 7 | flag == 8 | flag == 9)
    h1 = 10.^(-ph1);	%  [H+]    
end





% ---------------------------------------------------------------
% Constants load section
% ---------------------------------------------------------------
equic;                      
bor = 1.*(416.*(S/35.))* 1.e-6;   % (mol/kg), DOE94
i = 1;  % For loops in Calculation functions section
CO2_r  = zeros(size(numOfInputs));
HCO3_r = zeros(size(numOfInputs));
CO3_r  = zeros(size(numOfInputs));
ALK_r  = zeros(size(numOfInputs));
DIC_r  = zeros(size(numOfInputs));
PC02_r = zeros(size(numOfInputs));




% ---------------------------------------------------------------
% Calculation functions section
% ---------------------------------------------------------------
% Case 1  (h and s given)
while (flag == 1 & i <= numOfInputs)
    h = h1(i);
    s = s1(i);
    
    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;                        
    pco2 = s/Kh;
    
    % ----------- change units: mumol/kg
    CO2_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = pco2*1.e6;
    H_r(i)      = h;
    %alktest = (2*co31 + hco31 + bor/(1+h/Kb) + Kw/h - h)*1.e6
    
    i = i + 1;
end


% Case 101  (h and pco2 given)
while (flag == 101 & i <= numOfInputs)
    h = h1(i);
    pco2 = pco21(i);
    
    s = Kh*pco2;
    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
    CO2_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = pco2*1.e6;
    H_r(i)      = h;
    %alktest = (2*co31 + hco31 + bor/(1+h/Kb) + Kw/h - h)*1.e6
    
    i = i + 1;
end


% Case 2  (s and HCO3 given)
while (flag == 2 & i <= numOfInputs)
    s = s1(i);
    hco3 = hco31(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 3  (s and CO3 given)
while (flag == 3 & i <= numOfInputs)
    s = s1(i);
    co3 = co31(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 4  (s and ALK given)
while (flag == 4 & i <= numOfInputs)
    s = s1(i);
    alk = alk1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 5  (s and DIC given)
while (flag == 5 & i <= numOfInputs)
    s = s1(i);
    dic = dic1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 6  (h and HCO3 given)
while (flag == 6 & i <= numOfInputs)
    h = h1(i);
    hco3 = hco31(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 7  (h and CO3 given)
while (flag == 7 & i <= numOfInputs)
    h = h1(i);
    co3 = co31(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end       



% Case 8  (h and ALK given)
while (flag == 8 & i <= numOfInputs)
    h = h1(i);
    alk = alk1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 9  (h and DIC given)
while (flag == 9 & i <= numOfInputs)
    h = h1(i);
    dic = dic1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 10  (HCO3 and CO3 given)
while (flag == 10 & i <= numOfInputs)
    hco3 = hco31(i);  
    co3 = co31(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 11  (HCO3 and ALK given)
%
%       don't split the lines,
%       matlab does not understand it.
%
while (flag == 11 & i <= numOfInputs)
    hco3 = hco31(i);
    alk = alk1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 12  (HCO3 and DIC given)
while (flag == 12 & i <= numOfInputs)
    hco3 = hco31(i);
    dic = dic1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;

    i = i + 1;
end



% Case 13  (CO3 and ALK given)
while (flag == 13 & i <= numOfInputs)
    co3 = co31(i);
    alk = alk1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 14  (CO3 and DIC given)
while (flag == 14 & i <= numOfInputs)
    co3 = co31(i);
    dic = dic1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end



% Case 15  (ALK and DIC given)
while (flag == 15 & i <= numOfInputs)
    alk = alk1(i);
    dic = dic1(i);
    
    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_r(i)    = 1.e6*s;
    HCO3_r(i)   = hco3*1.e6;
    CO3_r(i)    = co3*1.e6;
    DIC_r(i)    = dic*1.e6;
    ALK_r(i)    = alk*1.e6;
    PC02_r(i)   = s*1.e6/Kh;
    H_r(i)      = h;
    
    i = i + 1;
end





% ---------------------------------------------------------------
% Displays Results section
% ---------------------------------------------------------------

% Prints to the screen
if(runUserInterface ~= 2)
    fprintf('\nResults:');
    
    % C02
    fprintf('\n\tCO2  = ');
    for i=1 : length(CO2_r)
        fprintf('%i\t', CO2_r(i));
    end
    
    % HCO3
    fprintf('\n\tHCO3 = ');
    for i=1 : length(HCO3_r)
        fprintf('%i\t', HCO3_r(i));
    end
    
    % CO3
    fprintf('\n\tCO3  = ');
    for i=1 : length(CO3_r)
        fprintf('%i\t', CO3_r(i));
    end
    
    % DIC
    fprintf('\n\tDIC  = ');
    for i=1 : length(DIC_r)
        fprintf('%i\t', DIC_r(i));
    end
    
    % ALK
    fprintf('\n\tALK  = ');
    for i=1 : length(ALK_r)
        fprintf('%i\t', ALK_r(i));
    end
    
    % PCO2
    fprintf('\n\tPCO2 = ');
    for i=1 : length(PC02_r)
        fprintf('%i\t', PC02_r(i));
    end
    
    % pH
    phtotal = zeros(size(H_r));
    phfree = zeros(size(H_r));
    
    i = 1;
    if(phflag == 0)
        fprintf('\nTotal pH scale used:');
        while (i <= length(H_r))
            phtotal(i) = -log10(H_r(i));
            oh = Kw./H_r(i);
            phfree(i) = -log10(H_r(i)/total2free);
            i = i + 1;
        end
        
        i = 1;
        fprintf('\n\tpH Total = ');
        while (i <= length(phtotal))
            fprintf('%i\t', phtotal(i));
            i = i + 1;
        end
        
        i = 1;
        fprintf('\n\tpH Free  = ');
        while (i <= length(phfree))
            fprintf('%i\t', phfree(i));
            i = i + 1;
        end
    end
    
    
    if(phflag == 1)
        fprintf('\nFree pH scale used:');
        while(i <= length(H_r))
            phfree(i)  = -log10(H_r(i));
            oh = Kw./H_r(i);   
            phtotal(i) = -log10(H_r(i)*total2free);
            i = i + 1;   
        end
        fprintf('\n\tpH Total = ');
        for i = 1 : length(phtotal)
            fprintf('%i\t', phtotal(i));
        end
        
        fprintf('\n\tpH Free  = ');
        for i = 1 : length(phfree)
            fprintf('%i\t', phfree(i));
        end
    end
    

% Prints to a file
else
    valid = 0;
    while(valid == 0)
        fprintf('\n');
        output_file = input('Enter Output Filename: ', 's');
        
        if(isempty(output_file))
             fprintf('\n\n*** FILE NAME NEEDED ***\n\n');
        else
            valid = 1;
            fid = fopen(output_file, 'w');
        end
    end
        
    if(delim_type == '1')
        delim = ',';
        fprintf(fid, 'C02');
        fprintf(fid, delim);
        fprintf(fid, 'HCO3-');
        fprintf(fid, delim);
        fprintf(fid, 'CO3--');
        fprintf(fid, delim);
        fprintf(fid, 'DIC');
        fprintf(fid, delim);
        fprintf(fid, 'ALK');
        fprintf(fid, delim);
        fprintf(fid, 'pCO2');
        fprintf(fid, delim);
        fprintf(fid, 'pH Total');
        fprintf(fid, delim);
        fprintf(fid, 'pH Free');
        fprintf(fid, '\n');
    elseif(delim_type == '2')
        delim = ' ';
    else
        delim = '\t';
    end
        
        
    
    
    for i=1 : length(CO2_r)
        fprintf(fid, '%i', CO2_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%i', HCO3_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%i', CO3_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%i', DIC_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%i', ALK_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%i', PC02_r(i));
        fprintf(fid, delim);

        if(phflag == 0)
            phtotal = -log10(H_r(i));
            oh = Kw./H_r(i);
            phfree = -log10(H_r(i)/total2free);
            fprintf(fid, '%i', phtotal);
            fprintf(fid, delim);
            fprintf(fid, '%i\n', phfree);
        end
        
        if(phflag == 1)
            phfree  = -log10(H_r(i));
            oh = Kw./H_r(i);   
            phtotal = -log10(H_r(i)*total2free);
            fprintf(fid, '%i', phtotal);
            fprintf(fid, delim);
            fprintf(fid, '%i\n', phfree);
        end
    end
    fclose('all');
end

fprintf('\n\n***Output Complete***');
return;