% ============================================================================
% MarChemSpec Template for Standard Seawater (Calculation Type = 3)
%
%
% Purpose: Calculate state variables and thermodynamic quantities of 
% solutions containing the ions of artificial seawater, and also Tris buffer, 
% at different practical salinities and temperatures. The buffer components of 
% the solution are entered as H+ and Tris molalities.
%
% The 'base' artificial seawater composition and method of making up the 
% buffer solutions is that of DelValls and Dickson (Deep-Sea Research I, 
% 45, 1541-1554 (1998)). 
%
% Further information can be found in the MarChemSpec manual.
% ============================================================================

% STEP 0: Don't forget to specify the directory that contains the MATLAB  
%         MEX file, using the addPath() command. It can also be done 
%         in this script. (Linux users: *before* starting MATLAB you must
%         make the assignments to LD_LIBRARY_PATH and LD_PRELOAD as described
%         in the manual).


% STEP 1: Assign the two directories that contain MarChemSpec data files, 
%         and the parent directory for those (ASWbuffer, and seawater) 
%         containing files of variances. They must be terminated with 
%         either a '\' (Windows) or '/' (Linux).
%    e.g. Windows: strFilesDirectory='{drive and path}\';
%         Linux:   strFilesDirectory='{path}/';

strFilesDirectory = 
strCovmxDirectory = 

% ----------------------------------------------------------------------------
% The following input variables have fixed values for these calculations.
% (Do not change any of these values!)
iCalc          = 3;                 % defines the type of medium to be modelled
P              = 1.0;               % pressure in atmospheres. The first release 
                                    %   of MarChemSpec is limited to a pressure 
                                    %   of 1 atmosphere
iConc          = 1;
nOut           = 8;                 % this is the number of state variables and 
                                    %   stoichiometric equilibrium constants to 
                                    %   be returned 
nSpecies       = 7;                 %  row length of nSpecies_in
Species_in     = zeros(3,nSpecies);
% ----------------------------------------------------------------------------

% STEP 2. User-defined input variables:                           
                   
tC                = 25;   % temperature in degrees Celsius (allowed range -5 to 45)
iConc             = 1;    % unit for Species_in: 1 = molality (mol/kg water); 
                          %                      2 = amount content (mol/kg solution)

% The Species_in array is used to define the totals of the artifical seawater
% components. This example is for a standard 0.04m equimolal buffer. When 
% making any changes to the composition, it is important to maintain the
% charge balance.

% Cations: 
Species_in(1,1:3) = [0.04 0.4461818 0.05474];  % H, Na, Mg
Species_in(1,4:5) = [0.01075 0.01058];         % Ca, K
% Anions:
Species_in(2,1:2) = [0.569202 0.02927];        % Cl, SO4
% Neutral:
Species_in(3,1)   = 0.08;                      % Tris

% The input molality for Tris is the sum of Tris and TrisH+, the molality
% of TrisH+ is defined by the input molality for H. The Na+ molality
% is reduced by the H molality to maintain charge balance.

% STEP 2a. Check for charge balance (this is also done within 
% MarChemSpec but the check is added here for convenience).

Cations = Species_in(1,1) + Species_in(1,2)...
    + 2 * (Species_in(1,3)+Species_in(1,4))...
    + Species_in(1,5) ;
Anions = Species_in(2,1) + 2 * Species_in(2,2);
Diff = Cations - Anions;
Diff2 = Anions - Cations;
Chargebalance = abs(Cations - Anions)/(0.5 * (Cations + Anions));
if(Chargebalance > 1e-5)
    disp ("Charge not balanced")
    if(Diff > 0) 
        disp ("Cation excess")
        disp (Diff)
    else
        disp ("Anion excess")
        disp (Diff2)
    end
    disp("Correct and run again")
    return
end


% STEP 3: Call the code to obtain the calculated properties of the artificial 
%         seawater you have specified. The contents of the output arrays are 
%         summarised further below. Return to step 2 for more calculations, 
%         or to correct errors. 

[Outputs, NamesOutputs, cSpecies_out, mSpecies_out,  ActCoeffs, NamesSpecies,iFail] = ...
  MarChemSpec(iCalc, tC, P, iConc, Species_in, nOut, strFilesDirectory,strCovmxDirectory);

% ============================================================================
%
% MARCHEMSPEC OUTPUTS, AND WHAT THEY MEAN
%
% Calls of MarChemSpec return acid-base and other properties of the
% artificial seawater being modelled, and the outputs are described in the 
% MarChemSpec documentation. Here is a summary:
%
% Outputs:      Double precision array, with two columns. The first column  
%               contains values of the variables that characterise the system 
%               modelled, and the second column their uncertainties. The
%               missing value code, for quantities not calculated, is 99.
%
% NamesOutputs: Column of shortform names for the variables in Outputs.
%
% cSpecies_out: Double precision array, with three rows, containing:
%               First row: cation amount contents (moles per kg of solution).
%               Second row: anion amount contents
%               Third row: neutral (uncharged species) amount contents.
%
% mSpecies_out: Same as cSpecies_out, but containing molalities (moles per
%               kg of pure water).
%
% Actcoeffs:    Same as cSpecies_out, but containing molality-based activity
%               coefficients.
%
% NameSpecies:  Character array with three rows, containing species names 
%               corresponding to the values in the three arrays above. 
%
% iFail:        On output, this contains error codes (a value of zero 
%               indicates complete success).
%
% ============================================================================
