% _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
%
% 	file: csysUI.m
%
%
%   Richard E. Zeebe* & Dieter A. Wolf-Gladrow**
%
%
%  *School of Ocean and Earth 
%   Science and Technology 
%   Department of Oceanography 
%   University of Hawaii at Manoa
%   1000 Pope Road, MSB 504
%   Honolulu, HI 96822
%   USA
%   email: zeebe@soest.hawaii.edu
%   http://www.soest.hawaii.edu/oceanography/faculty/zeebe.html
%
%               			     
% **Alfred Wegener Institute for       
%	Helmholtz Centre for Polar and Marine Research (AWI)	  
%	P.O. Box 12 01 61		       
%	D-27515 Bremerhaven		
%	Germany	      
%	e-mail: Dieter.Wolf-Gladrow@awi.de  	
%	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:
%
% 04/07/14 DWG & REZ: ocdflag for comparison with CO2SYS
%          KW P-coefficients water -> seawater
% 03/26/14 REZ: equic.m: added p2f (pCO2 -> fCO2)
%          csys.m: output pCO2 and fCO2
%          R: 83.131 -> 83.14510. Kw: 148.96502 -> 148.9652 
%          (use old R to check old values)
% 11/07/12 REZ: clear all, '%' removed from '%h1'
% 03/02/10 REZ: boric acid pressure coefficient corrected
%          (see equic.m)
% 10/24/04 user interface implemented
%          data reading from external file 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 'csysUI'
%
% ---------------------------------------------------------------

clear all;
runUserInterface = 1;
    %  = 0 -> will NOT run user interface
    % >= 1 -> will run user interface
    
% ---------------------------------------------------------------
% 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     = 15;                  % Given Variables (see chart)
    phflag   = 0;                   % Total Scale (0) or Free Scale (1)
    k1k2flag = 1;                   % Roy (0) or Mehrbach (1)
    ocdflag  = 0;                   % csys (0) or CO2SYS V1.1 (1)
    TC = 25;                        % Temperature   (Celcius)
    S =  35;                        % Salinity
    P =  000;                       % Pressure      (bar)
    ph1   = [8.1, 8.2, 8.3];        % pH
    s1    = [9, 10, 11];            % [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 input,   set runUserInterface = 1
%   - If data is read from a file, set 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 CO2 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



% ---------------------------------------------------------------
% ocd Section
%   - ocdflag will be set with a numeric value either 
%     corresponding to the users input or the default value
% ---------------------------------------------------------------
    %----- compare output to CO2SYS? use phflag=0, k1k2flag=1
    %
    % 1: should be close to CO2SYS V1.1 Matlab 09/2011
    % 0: else
    %
    valid = 0;
    while(valid == 0)
        fprintf('\nCompare output to CO2SYS?\n');
        fprintf('0 - No.\t\t[Default]\n');
        fprintf('1 - Yes. Compare to CO2SYS V1.1 Matlab 09/2011 \n');
        
        % Gets user input.  If the user does not enter a selection, the
        % default (0) will be used.
        str_ocdflag = input('Selection: ', 's');
        if(~isempty(str_ocdflag))
            ocdflag = str2num(str_ocdflag);
            % Error checking
            if(isnumeric(ocdflag))
                if(ocdflag == 0 | ocdflag == 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 [0]\n\n');
            ocdflag = 0;
            valid = 1;
        end
    end
    if(ocdflag == 1)
        fprintf('\nTotal pH scale and K1, K2 from ');        
        fprintf('Lueker et al. (2000) will be used.\n');
        phflag = 0;
        k1k2flag = 1;
    end
            
          
if(ocdflag == 0)
% ---------------------------------------------------------------
% 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
    
end % ocdflag == 0


% ---------------------------------------------------------------
% 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 % runUserInterface

% ======================================================   
if(ocdflag == 1 & phflag == 1)
  str1 = '>> Compare to CO2SYS on total scale: ';
  str2 = 'for ocdflag = 1, set phflag = 0';
  error([str1 str2])
end
if(ocdflag == 1 & k1k2flag == 0)
  str1 = '>> Compare to CO2SYS: for ';
  str2 = 'ocdflag = 1, set k1k2flag = 1';
  error([str1 str2])
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*1.e-6;
end

if(flag == 2 | flag == 6 | flag == 10 | flag == 11 | flag == 12)
    hco31 = hco31*1.e-6;
end

if(flag == 3 | flag == 7 | flag == 10 | flag == 13 | flag == 14)
    co31 = co31*1.e-6;
end

if(flag == 4 | flag == 8 | flag == 11 | flag == 13 | flag == 15)
    alk1 = alk1*1.e-6;
end

if(flag == 5 | flag == 9 | flag == 12 | flag == 14 | flag ==  15)
    dic1 = dic1*1.e-6;
end

if(flag == 101)
    pco21 = pco21*1.e-6;
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
if(ocdflag)
bor = 1.*(415.7*(S/35.))*1.e-6;  % (mol/kg), CO2SYS
end
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));
fCO2_r = zeros(size(numOfInputs));
pCO2_r = zeros(size(numOfInputs));
H_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;                        
    fco2 = s/Kh;
    pco2 = fco2/p2f;

    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*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);
    fco2 = pco2*p2f;
    
    s = Kh*fco2;    
    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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;

    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;

    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;

    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;

    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    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;
    fco2 = s/Kh;
    pco2 = fco2/p2f;
    
    % ----------- 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;
    pCO2_r(i)   = pco2*1.e6;
    fCO2_r(i)   = fco2*1.e6;
    H_r(i)      = h;
    
    i = i + 1;
end





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

% Prints to the screen
if(runUserInterface ~= 2)
    fprintf('\nResults:');
    
    % CO2
    fprintf('\n\tCO2  = ');
    for i=1 : length(CO2_r)
        fprintf('%e\t', CO2_r(i));
    end
    
    % HCO3
    fprintf('\n\tHCO3 = ');
    for i=1 : length(HCO3_r)
        fprintf('%e\t', HCO3_r(i));
    end
    
    % CO3
    fprintf('\n\tCO3  = ');
    for i=1 : length(CO3_r)
        fprintf('%e\t', CO3_r(i));
    end
    
    % DIC
    fprintf('\n\tDIC  = ');
    for i=1 : length(DIC_r)
        fprintf('%e\t', DIC_r(i));
    end
    
    % ALK
    fprintf('\n\tALK  = ');
    for i=1 : length(ALK_r)
        fprintf('%e\t', ALK_r(i));
    end
    
    % pCO2
    fprintf('\n\tpCO2 = ');
    for i=1 : length(pCO2_r)
        fprintf('%e\t', pCO2_r(i));
    end
    
    % fCO2
    fprintf('\n\tfCO2 = ');
    for i=1 : length(fCO2_r)
        fprintf('%e\t', fCO2_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('%e\t', phtotal(i));
            i = i + 1;
        end
        
        i = 1;
        fprintf('\n\tpH Free  = ');
        while (i <= length(phfree))
            fprintf('%e\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('%e\t', phtotal(i));
        end
        
        fprintf('\n\tpH Free  = ');
        for i = 1 : length(phfree)
            fprintf('%e\t', phfree(i));
        end
    end
    

% Prints to file with user-given filename.
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, 'CO2');
        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, 'fCO2');
        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, '%e', CO2_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%e', HCO3_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%e', CO3_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%e', DIC_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%e', ALK_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%e', pCO2_r(i));
        fprintf(fid, delim);
        fprintf(fid, '%e', fCO2_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, '%e', phtotal);
            fprintf(fid, delim);
            fprintf(fid, '%e\n', phfree);
        end
        
        if(phflag == 1)
            phfree  = -log10(H_r(i));
            oh = Kw./H_r(i);   
            phtotal = -log10(H_r(i)*total2free);
            fprintf(fid, '%e', phtotal);
            fprintf(fid, delim);
            fprintf(fid, '%e\n', phfree);
        end
    end
    fclose('all');
end

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

% ======================================================
if(TC < 0.0 | TC > 35.)
  str1 = '>> Temperature out of valid range for K1 and K2. ';
  str2 = 'Limits: ~0-35 deg C';
  error([str1 str2])
end
if(S < 20. | S > 40.)
  str1 = '>> Salinity out of valid range for K1 and K2. ';
  str2 = 'Limits: ~20-40';
  error([str1 str2])
end

return;