function rsw= ph_density(t,s,pp)
%-------------------------------------------------------------------
% PH_DENSITY $Revision: 1.0 $ $Date: 1998/12/07 $
%          Copyright (C) GEOMAR, Roger Luff 1998.
%
% USAGE: rsw = ph_density(t,s,pp)
%
% DESCRIPTION:
%   density calculates the density of seawater using the unesco
%   (1983) formulae, which are valid for a temperature range 
%   of 0-40 deg c and 0.5-43 ppt
%
% INPUT:
%  t     : temperature of seawater at the specific location [deg C]
%  s     : salinity of seawater at the specific location  [PSU]
%  pp    : pressure at the specific location  [atm, (with 1 atm
%          at sea surface)]
%
% OUTPUT:
%  rsw   : density of a "seawater" [kg/m^3]
%
% ORIGINAL AUTHOR: 
%   Bernard Boudreau 
%   Dept. Oceanography 
%   Dalhousie University
%   Halifax, NS B3H 4J1,Canada 
%   bernie@boudreau.ocean.dal.ca
% MATLAB VERSION
%   Roger Luff 98-12-07
%   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.
%  Fofonoff P. and Millard R. C. (1983) "Algorithms for computation 
%   of fundamental properties of seawater." UNESCO Technical Paper 
%   in Marine Science 44, 53.
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% CHECK INPUT ARGUMENTS
%-------------------------------------------------------------------
if nargin ~=3
   error('ph_density.m: Must pass 3 parameters')
end
%-------------------------------------------------------------------
%  convert pressure to bars of sea pressure only
%-------------------------------------------------------------------
p = (pp-1.0)*1.013;
%-------------------------------------------------------------------
%  density of pure reference water
%-------------------------------------------------------------------
a0 = 999.842594;
a1 = 6.793952d-02;
a2 = - 9.095290d-03;
a3 = 1.001685d-04;
a4 = - 1.120083d-06;
a5 = 6.536332d-09;
rw = ((((a5*t + a4)*t + a3)*t + a2)*t + a1)*t + a0;
%-------------------------------------------------------------------
%  density at t, s and 1 atm
%-------------------------------------------------------------------
b0 = 8.24493d-01;
b1 = - 4.0899d-03;
b2 = 7.6438d-05;
b3 = - 8.2467d-07;
b4 = 5.3875d-09;
c0 = - 5.72466d-03;
c1 = 1.0227d-04;
c2 = - 1.6546d-06;
d0 = 4.8314d-04;
rts = rw + ((((b4*t + b3)*t + b2)*t +b1)*t + b0)*s + ((c2*t + c1)*t + c0)*sqrt(s*s*s) + d0*s*s;
if p==0.0
   rsw = rts;
else
%-------------------------------------------------------------------
% density at t, s and p
%-------------------------------------------------------------------
  e0 = 19652.21;
  e1 = 148.4206;
  e2 = - 2.327105;
  e3 = 1.360477d-02;
  e4 = - 5.155288d-05;
  kw = (((e4*t + e3)*t + e2)*t + e1)*t + e0;
  h0 = 3.239908;
  h1 = 1.43713d-03;
  h2 = 1.16092d-04;
  h3 = - 5.77905d-07;
  aw = ((h3*t + h2)*t + h1)*t + h0;
  k0 = 8.50935d-05;
  k1 = - 6.12293d-06;
  k2 = 5.2728d-08;
  bw = (k2*t + k1)*t + k0;
  m0 = - 9.9348d-07;
  m1 = 2.0816d-08;
  m2 = 9.1697d-010;
  b = bw + ((m2*t + m1)*t + m0)*s;
  i0 = 2.2838d-03;
  i1 = - 1.0981d-05;
  i2 = - 1.6078d-06;
  j0 = 1.91075d-04;
  a = aw + ((i2*t + i1)*t + i0)*s + j0*sqrt(s*s*s);
  f0 = 54.6746;
  f1 = - 0.603459;
  f2 = 1.09987d-02;
  f3 = - 6.1670d-05;
  g0 = 7.944d-02;
  g1 = 1.6483d-02;
  g2 = - 5.3009d-04;
  kst0 = kw + (((f3*t + f2)*t + f1)*t + f0)*s + ((g2*t + g1)*t + g0)*sqrt(s*s*s);
  kstp = kst0 + a*p + b*p*p;
  rsw =  rts/(1.0 - p/kstp);
end;
