function [hp,co2,co3,hco3,hs,h2s,h3po4,h2po4,hpo4,po4,nh3,nh4,boh3,boh4] = ph_eqcharge(k,ssum,tco2,th2s,tboh,tnh4,tpo4,mnii,feii,ca,so4,no3)
% DESCRIPTION:
%   This function calculates the H+ concentration and the concentration
%   of the dissociated species from the INPUT variables, listed later
%   on the bases of the charge balance equation.
%   To calculate the pH from H+ use: ph=-log10(hp/1.0E+03)
%
% INPUT:
%  k    = equilibrium constants from ph_thermo
%  ssum = sum concentration of the conservative species in seawater
%         from ph_cons
%  tco2 = Total carbonate concentration [mmol/l]
%  th2s = Total suffate concentration [mmol/l]
%  tboh = Total borate concentration [mmol/l]
%  tnh4 = Total amonium concentration [mmol/l]
%  tpo4 = Total phospate concentration [mmol/l]
%  mnii = Mn2+ concentration [mmol/l]
%  feii = Fe2+ concentration [mmol/l]
%  ca   = concentration of total calcite in seawater
%  so4  = concentration of total sulfate in seawater
%  no3  = Nitrate concentration [mmol/l] 
%
% OUTPUT:
%   hp   = Equilibrium concentration H+ [mmol/l]
%   co2  = Equilibrium concentration of CO2 [mmol/l]
%   co3  = Equilibrium concentration of CO3[mmol/l]
%   hco3 = Equilibrium concentration of HCO3 [mmol/l]
%   hs   = Equilibrium concentration of HS [mmol/l]
%   h2s  = Equilibrium concentration of H2S [mmol/l]
%   h3po4 = Equilibrium concentration of H3PO4 [mmol/l]
%   h2po4 = Equilibrium concentration of H2PO4 [mmol/l]
%   hpo4 = Equilibrium concentration of HPO4 [mmol/l]
%   nh3 = Equilibrium concentration of NH3 [mmol/l]
%   nh4 = Equilibrium concentration of NH4 [mmol/l]
%   po4 = Equilibrium concentration of PO4 [mmol/l]
%   boh3 = Equilibrium concentration of BOH3 [mmol/l]
%   boh4 = Equilibrium concentration of BOH4 [mmol/l]
%
% AUTHOR: 
%   Bernard Boudreau(1),  Roger Luff(2) 
%   (1) Dept. Oceanography 
%       Dalhousie University
%       Halifax, NS B3H 4J1,Canada 
%       bernie@boudreau.ocean.dal.ca
%
%   (2) GEOMAR Research Centre for Marine Geosciences
%       Wischhofstr. 1-3, D-24148 Kiel, Germany
%       visumod@technologist.com
%
% DISCLAIMER:
%   This software is provided "as it is" without warranty of any kind.
%
% REFERENCES:
%   Boudreau B. P. (1996) A method-of-lines code for carbon and 
%   nutrient diagenesis in aquatic sediments.Computers & Geosciences
%   22(5), 479-496.
%   Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN
%   and MATLAB libraries to calculate pH distributions in
%   marine systems", Computers & Geosciences.
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% CHECK INPUT ARGUMENTS
%-------------------------------------------------------------------
if nargin ~=12
   error('ph_eqcharge.m: Must pass 12 parameters')
end
%-------------------------------------------------------------------
% SET A FIRST GUESS TO THE FINAL EQUILIBRIUM VALUES BASED ON THE 
% SUM CONCENTRATIONS
%-------------------------------------------------------------------
co3  = 0.11*tco2;
co2  = 0.09*tco2;
hco3 = 0.80*tco2;
hs = th2s/2;
h2s=th2s/2;
boh3=tboh/2;
boh4=tboh/2;
h3po4= tpo4*0.01 ;
h2po4= tpo4*0.45 ;
hpo4 = tpo4*0.45 ;
po4  = tpo4*0.09; 
nh3  = tnh4*0.5;
nh4  = tnh4*0.5;
%-------------------------------------------------------------------
% CALCULATE THE H+ CONCENTRATION ON THE BASE OF THE FIRTS GUESS
%-------------------------------------------------------------------
if co2<0.1 
  hp   = hco3*k(3)/co3;
else
  hp   = co2*k(2)/hco3;
end
oh   = k(1)/hp;
%-------------------------------------------------------------------
% SET UP LOCAL VARIABLES
%-------------------------------------------------------------------
maxiter = 1000;
iter = 0;
atol = 100;
tol = 1e-11;
%-------------------------------------------------------------------
% ITERATION LOOP STARTS HERE 
%-------------------------------------------------------------------
while atol>tol & iter<maxiter    
  iter = iter + 1;
%-------------------------------------------------------------------
% SET UP THE JACOBI
%-------------------------------------------------------------------
  matrix=zeros(15);
% line 1
  matrix(1,1) = -1.0;
  matrix(1,7) = k(2);
% line 2
  matrix(2,1) = -1.0;
  matrix(2,7) = -hp;
  matrix(2,8) = k(3);
  matrix(2,15) = -1.0;
% line 3
  matrix(3,1) = -1.0;
  matrix(3,8) = -hp;
  matrix(3,15) = -2.0;
% line 4
  matrix(4,2) = -1.0;
  matrix(4,14) = k(10);
% line 5
  matrix(5,2) = -1.0;
  matrix(5,14) = -hp;
  matrix(5,15) = -1.0;
% line 6
  matrix(6,3) = -1.0;
  matrix(6,10) = k(5);
% line 7
  matrix(7,3) = -1.0;
  matrix(7,10) = -hp;
  matrix(7,11) = k(6);
  matrix(7,15) = -1.0;
% line 8
  matrix(8,3) = -1.0;
  matrix(8,11) = -hp;
  matrix(8,12) = k(7) ;
  matrix(8,15) = -2.0;
% line 9
  matrix(9,3) = -1.0;
  matrix(9,12) = -hp;
  matrix(9,15) = -3.0;
% line 10
  matrix(10,4) = -1.0;
  matrix(10,13) = -hp;
% line 11
  matrix(11,4) = -1.0;
  matrix(11,13) = k(9);
  matrix(11,15) = 1.0;
% line 12
  matrix(12,5) = -1.0;
  matrix(12,9) = k(4);
% line 13
  matrix(13,5) = -1.0;
  matrix(13,9) = -hp;
  matrix(13,15) = -1.0;
% line 14
  matrix(14,6) = -oh;
  matrix(14,7) = -hco3;
  matrix(14,8) = -co3;
  matrix(14,9) = -boh4;
  matrix(14,10) = -h2po4;
  matrix(14,11) = -hpo4;
  matrix(14,12) = -po4;
  matrix(14,13) = -nh3;
  matrix(14,14) = -hs;
  matrix(14,15) = 1.0 ;
% line 15
  matrix(15,6) = -hp;
  matrix(15,15) = -1.0;
%-------------------------------------------------------------------
% INVERT THE JACOBI
%-------------------------------------------------------------------
  matrix=inv(matrix);
%-------------------------------------------------------------------
% SET UP THE FUNCTIONS WHICH HAVE TO BE MINIMIZED
%-------------------------------------------------------------------
  f(1) = tco2 - (co2+co3+hco3);
  f(2) = th2s - (h2s+hs);
  f(3) = tpo4 - (h3po4 + h2po4 + hpo4 + po4);
  f(4) = tnh4 - (nh3 + nh4);
  f(5) = tboh - (boh3 + boh4);
  f(6) = k(1) - hp*oh;
  f(7) = k(2)*co2 - hp*hco3;
  f(8) = k(3)*hco3 - hp*co3;
  f(9) = k(4)*boh3 - hp*boh4;
  f(10)=  k(5)*h3po4 - hp*h2po4;
  f(11) = k(6)*h2po4 - hp*hpo4;
  f(12) = k(7)*hpo4  - hp*po4;
  f(13) = k(9)* nh4  - hp*nh3;
  f(14) = k(10)*h2s  - hp*hs;
  f(15) = ssum+2*ca+hp+nh4+2*(mnii+feii)-hco3-h2po4-hs-boh4-2*so4-2*co3-2*hpo4-3*po4-oh-no3;
  dd(1:15)=0;
%-------------------------------------------------------------------
% CALCULATE THE DISTANCE TO THE NEXT INTERATION STEP
%-------------------------------------------------------------------
  for kk=1:15
    for j=1:15
      dd(kk) =dd(kk) +matrix(j,kk)*f(j);
    end
  end;
%-------------------------------------------------------------------
% CALCULATE THE NEW DISTRIBUTION
%-------------------------------------------------------------------
  co2=co2-dd(1);
  hco3=hco3-dd(2);
  co3=co3-dd(3);
  h2s=h2s-dd(4);
  hs=hs-dd(5);
  h3po4=h3po4-dd(6);
  h2po4=h2po4-dd(7);
  hpo4=hpo4-dd(8);
  po4=po4-dd(9);
  nh3=nh3-dd(10);
  nh4=nh4-dd(11);
  boh3=boh3-dd(12);
  boh4=boh4-dd(13);
  hp =hp-dd(14);
  oh =oh-dd(15);
%-------------------------------------------------------------------
% CHECK THE CONVERGENCE
%-------------------------------------------------------------------
  atol = sum(abs(dd));
end;
%-------------------------------------------------------------------
% IN THE CASE THAT THE FUNCTION DOES NOT CONVERGENTE, END WITH AND
% ERROR MESSAGE OR WITH A H+ CONCENTRATION OF NaN
%-------------------------------------------------------------------
if (hp <= 0)
  hp =NaN;
  warning ('Solver ends with an error, bad choices of input parameters')
end;
if (iter==maxiter)
  error('ph_eqcharge: too much iterations');
end;
