SUBROUTINE ph_density(RSW,S,T,PP)
c-------------------------------------------------------------------
c ph_density $Revision: 1.0 $ $Date: 1998/12/07 $
c            Copyright (C) GEOMAR, Roger Luff 1998.
c
c USAGE:
C      double precision RSW,S,T,PP
C      call ph_density(RSW,S,T,PP)
C   SEE ph_test.f FOR DETAILS OF USING THIS PROGRAM!!!!!!!!!!!!!!!
c
c DESCRIPTION:
c   density calculates the density of seawater using the unesco
c   (1983) formulae, which are valid for a temperature range 
c   of 0-40 deg c and 0.5-43 ppt
c
c INPUT:
c  t     : temperature of seawater at the specific location [deg C]
c  s     : salinity of seawater at the specific location  [PSU]
c  pp    : pressure at the specific location  [atm, (with 1 atm
c          at sea surface)]
c
c OUTPUT:
c  rsw   : density of a "seawater" [kg/m^3]
c
c ORIGINAL AUTHOR: 
c   Bernard Boudreau 
c   Dept. Oceanography 
c   Dalhousie University
c   Halifax, NS B3H 4J1,Canada
c   bernie@boudreau.ocean.dal.ca
c REVISION BY:
c   Roger Luff and Matthias Haeckel 98-12-07
c   GEOMAR Research Centre for Marine Geosciences
c   Wischhofstr. 1-3, D-24148 Kiel, Germany
c   visumod@technologist.com
c
c DISCLAIMER:
c   This software is provided "as is" without warranty of any kind.
c
c REFERENCES:
c  Boudreau B. P. (1996) "A method-of-lines code for carbon and
c   nutrient diagenesis in aquatic sediments". Computers & Geosciences
c   22(5), 479-496.
c  Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN
c   and MATLAB libraries to calculate pH distributions in marine
c   systems", Computers & Geosciences.
c  Fofonoff P. and Millard R. C. (1983) "Algorithms for computation
c   of fundamental properties of seawater." UNESCO Technical Paper
c   in Marine Science 44, 53.
c-------------------------------------------------------------------

      implicit none
      double precision A0,A1,A2,A3,A4,A5,RW
      double precision B0,B1,B2,B3,B4
      double precision C0,C1,C2,sqsss
      double precision D0,RTS,RSW,P,PP,S,T
      double precision E0,E1,E2,E3,E4,KW
      double precision H0,H1,H2,H3,AW
      double precision K0,K1,K2,BW
      double precision M0,M1,M2,B,KSTP
      double precision I0,I1,I2,J0,A
      double precision F0,F1,F2,F3,G0,G1,G2,KST0

c-------------------------------------------------------------------
C  Convert pressure to bars of sea pressure only
c-------------------------------------------------------------------
      P = (PP-1.0)*1.013
      sqsss = SQRT(S*S*S)
c-------------------------------------------------------------------
C  Density of pure reference water
c-------------------------------------------------------------------
      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
c-------------------------------------------------------------------
C  Density at t, S and 1 atm
c-------------------------------------------------------------------
      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)*sqsss + D0*S*S
      IF(P.EQ.0.0) THEN
         RSW = RTS
         RETURN
      ENDIF
c-------------------------------------------------------------------
C  Density at t, S and P
c-------------------------------------------------------------------
      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*sqsss
      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)*sqsss
      KSTP = KST0 + A*P + B*P*P
      RSW =  RTS/(1.0 - P/KSTP)
      RETURN
      END