% ============================================================================
% MarChemSpec Template for Standard Seawater (Calculation Type = 7)
%
% Purpose: Calculate the fractional speciation of one or more GEOTRACES trace 
% metals in solutions of user-specified composition (containing the ions of
% seawater) and temperature. These solutions can also include additions of 
% phosphate and silicate. 
% 
% Note: In the future it will be possible to equilibrate the solution to
% user-specified values of pairs of the quantities pCO2, pH(T), DIC, and
% alkalinity. This feature has not yet been programmed. 
%
% The 'base' seawater composition in MarChemSpec is that of the standard 
% seawater of Millero et al. (Deep-Sea Research I, 55, 50-72 (2008)). 
%
% 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       = 7;                 % defines the type of medium to be modelled
P           = 1;                 % pressure in atmospheres. The first release 
                                 %   of MarChemSpec is limited to a pressure 
                                 %   of 1 atmosphere
nOut        = 40;                % The number of trace metal parameters to be 
                                 %   returned.
iFix        = 0;                 % reserved for future development
ValuesFixed = [0 0; 0 0];        % reserved for future development
nSpecies    = 39;                % row length of nSpecies_in
Species_in  = zeros(3,nSpecies); % ensures that the Species_in array is empty
%-----------------------------------------------------------------------------

% STEP 2. User-defined input variables:                           
                   
tC                = 25;   % temperature in degrees Celsius (allowed range -5 to 45)
iConc             = 2;    % 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 seawater
% components. This example is for standard seawater. When making any
% changes to the composition, it is important to maintain the charge balance

% Cations:
Species_in(1,2:4)  = [0.46896739 0.05281709 0.01028205];    % Na, Mg, Ca
Species_in(1,5:6)  = [0.01020766 0.00009069];               % K, Sr
Species_in(1,7:16) = 10*1e-8;                               % trace metals
% Anions:
Species_in(2,1:2)  = [0.54586956 0.02823522];               % Cl, SO4
Species_in(2,3)    = 0.00230441;                            % alkalinity
Species_in(2,4)    = 0.00084211;                            % Br
Species_in(2,8)    = 0.00006831;                            % F
Species_in(2,9)    = 1e-8;                                  % P
% Neutrals:
Species_in(3,1:2)  = [0.00041517 0.00196633];               % B, DIC  
Species_in(3,4)    = 1e-8;                                  % Si


% STEP 2a. Check for charge balance. Note that this test does not
% include the trace metals which have very low concentrations such
% that their total contribution falls below the tolerance of 1.E-5. 
% (This test, in complete form, is also done within MarChemSpec but  
% the check is added here for convenience.)

Cations = Species_in(1,2) + 2 * (Species_in(1,3)+Species_in(1,4))...
    + Species_in(1,5) + 2 * Species_in(1,6);
Anions = Species_in(2,1) + 2 * Species_in(2,2)...
    + Species_in(2,3) + Species_in(2,4) + Species_in(2,8)...
    + Species_in(2,9);
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 natural water 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, iFix, ValuesFixed, nOut, strFilesDirectory,strCovmxDirectory);

% ============================================================================
%
% MARCHEMSPEC OUTPUTS, AND WHAT THEY MEAN
%
% Calls of MarChemSpec return acid-base and other properties of the natural
% water 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 results summarising the trace metal
%               speciation and the second column their uncertainties.
%
% 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).
%
% ============================================================================
