SUBROUTINE ph_tdconst(K,T,S,P)
C PH_eqconst $Revision: 1.0 $ $Date: 1998/12/07 $
C
C USAGE: 
C      double precision k(12),t,s,p
C      call ph_tdconst(k,t,s,p)
C   SEE ph_test.f FOR DETAILS OF USING THIS PROGRAM!!!!!!!!!!!!!!!
C
C DESCRIPTION:
C   ph_tdconst calculates the values of thermodynamic constants
C   k(i), in seawater for given temperature, salinity, pressure. The
C   calculated constants are on the sws ph scale.
C   these constants assume the concentrations are in molality.
C   the equilibrium constants to be calculated are:
C   k(1) = dissociation constant of h2o
C   k(2) = 1st dissociation constant of co2(aq)
C   k(3) = 2nd dissociation constant of co2(aq)
C   k(4) = dissociation constant of b(oh)3
C   k(5) = 1st dissociation constant of h3po4
C   k(6) = 2nd dissociation constant of h3po4
C   k(7) = 3nd dissociation constant of h3po4
C   k(8) = dissociation constant of hno3
C   k(9) = dissociation constant of nh4
C   k(10) = 1st dissociation constant of h2s
C   k(11) = solubility product of caco3 (aragonite)
C   k(12) = solubility product of caco3 (calcite)
C  correction to in situ pressure uses the basic correction formula:
C        ln(kp/k0) = - dv/(r*tk)*pp + 0.5*dk/(r*tk)*pp*pp
C  where dv is the change in partial molal volume and dk is the
C  compressibility. pressure pp is water pressure over atmospheric
C  in bars.
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  p     : pressure at the specific location  [atm, (with 1 atm
C          at sea surface)]
C
C OUTPUT:
C  k(12) : equilibrium constants
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 MINOR CHANGES:
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 it 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  Clegg S. L. and Whitefield M. (1995) "A chemical model of seawater
C   including dissolved ammonia and the stoichiometric dissociation
C   constant of ammonia in estuarine water and seawater from -2 to
C   40C." Geochim. Cosmochim. Acta 59(12), 2403-2421.
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  Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN
C   and MATLAB libraries to calculate pH distributions in 
C   marine systems", Computers & Geosciences.
C  Millero F. J. (1983) "Influence of pressure on chemical processes
C   in the sea." In Chemical Oceanography 2nd ed., Vol. 8 (ed. J. P.
C   Riley and R. Chester), pp. 1-88. Academic Press.
C  Millero F. J. (1995) "Thermodynamics of the carbon dioxide system
C   in the oceans." Geochim. Cosmochim. Acta 59(4), 661-677.
C-------------------------------------------------------------------

      IMPLICIT double precision (A-H,K,L,O-Z)
      double precision K(12)
      TK = T + 273.15D+00
      LNTK = LOG(TK)
      S2 = S*S
      SQS = SQRT(S)
      S23 = SQRT(S*S*S)
      PP = (P-1.0)*1.013D+00
      R = 83.147D+00
C-------------------------------------------------------------------
C       1) WATER:
C-------------------------------------------------------------------
      LNK1 = 148.9802 - 13847.26/TK - 23.6521*LNTK
      LNK1 = LNK1 + (-5.977 + 118.67/TK + 1.0495*LNTK)*SQS - 0.01615*S
      DVW = - 20.02 + 0.1119*T - 1.4096D-03*T*T
      DKW = (-5.13 + 0.0794*T)/1000.0
      LNK1 = - DVW/(R*TK)*PP + 0.5*DKW/(R*TK)*PP*PP + LNK1
      K(1) = EXP(LNK1)
C-------------------------------------------------------------------
C     2) CARBONIC ACID   (using the combined, i.e. COM, equations)
C-------------------------------------------------------------------
      LNKC1 = 2.18867 - 2275.036/TK - 1.468591*LNTK
      LNKC1 = LNKC1 + (-0.138681 - 9.33291/TK)*SQS
      LNKC1 = LNKC1 + 0.0726483*S - 0.00574938*S23
      A0 = 25.50
      A1 = 0.151
      A2 = - 0.1271
      DV = -(A0 + A1*(S-34.8) + A2*T)
      B0 = 3.08
      B1 = 0.578
      B2 = - 0.0877
      DK = - (B0 + B1*(S-34.8) + B2*T)/1000.0
      LNKC1 = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNKC1
      K(2) = EXP(LNKC1)

      LNKC2 = -0.84226 - 3741.1288/TK - 1.437139*LNTK
      LNKC2 = LNKC2 + (-0.128417 - 24.41239/TK)*SQS
      LNKC2 = LNKC2 + 0.1195308*S - 0.0091284*S23
      A0 = 15.82
      A1 = - 0.321
      A2 = 0.0219
      DV = -(A0 + A1*(S-34.8) + A2*T)
      B0 = - 1.13
      B1 = 0.314
      B2 = 0.1475
      DK = - (B0 + B1*(S-34.8) + B2*T)/1000.0
      LNKC2 = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNKC2
      K(3) = EXP(LNKC2)
C-------------------------------------------------------------------
C     3) BORIC ACID
C-------------------------------------------------------------------
      LNKB = (-8966.90 -2890.51*SQS -77.942*S +1.726*S23 -0.0993*S2)/TK
      LNKB = LNKB + (148.0248 + 137.194*SQS + 1.62247*S)
      LNKB = LNKB + (-24.4344 - 25.085*SQS - 0.2474*S)*LNTK
      LNKB = LNKB + 0.053105*SQS*TK
      A0 = 29.48
      A1 = - 0.295
      A2 = - 0.1622
      A3 = 2.608D-03
      DV = -(A0 + A1*(S-34.8) + A2*T + A3*T*T)
      B0 = 2.84
      B1 = - 0.354
      DK = - (B0 + B1*(S-34.8))/1000.0
      LNKB = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNKB
      K(4) = EXP(LNKB)
C-------------------------------------------------------------------
C     4) PHOSPHORIC ACID
C-------------------------------------------------------------------
      LNKP1 = 115.54 - 4576.752/TK - 18.453*LNTK
      LNKP1 = LNKP1 + (0.69171 - 106.736/TK)*SQS
      LNKP1 = LNKP1 + (-0.01844 - 0.65643/TK)*S
      DV = - 14.51 + 0.1211*T - 0.321D-03*T*T
      DK = (-2.67 + 0.0427*T)/1000.0
      LNKP1 = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNKP1
      K(5) = EXP(LNKP1)

      LNKP2 = 172.1033 - 8814.715/TK - 27.927*LNTK
      LNKP2 = LNKP2 + (1.3566 - 160.34/TK)*SQS
      LNKP2 = LNKP2 + (-0.05778 + 0.37335/TK)*S
      DV = - 23.12 + 0.1758*T - 2.647D-03*T*T
      DK = (-5.15 + 0.0900*T)/1000.0
      LNKP2 = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNKP2
      K(6) = EXP(LNKP2)

      LNKP3 = -18.126 - 3070.75/TK
      LNKP3 = LNKP3 + (2.81197 + 17.27039/TK)*SQS
      LNKP3 = LNKP3 + (-0.09984 - 44.99486/TK)*S
      DV = - 26.57 + 0.2020*T - 3.042D-03*T*T
      DK = (-4.08 + 0.0714*T)/1000.0
      LNKP3 = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNKP3
      K(7) = EXP(LNKP3)
C-------------------------------------------------------------------
C     5) NITRIC ACID
C-------------------------------------------------------------------
      K(8) = 23.44
C-------------------------------------------------------------------
C     6) AMMONIA
C-------------------------------------------------------------------
      LNK9 = - 6285.33/TK + 0.0001635*TK - 0.25444
     &       + (0.46532 - 123.7184/TK)*SQS

     &       + (-0.01992 + 3.17556/TK)*S
C-------------------------------------------------------------------
C     note: Millero (1983) gives the DV and DK for the reaction
C            NH3 + H2O = NH4 + OH
C     while the value calculated by this program is for the reaction
C                  NH4 = NH3 + H
C     To get the desired value, the DV for water is subtracted from
C     the Millero DV and the the negative of this resulting number is
C     used.
C-------------------------------------------------------------------
      DV = - 26.43 + 0.0889*T - 0.905-03*T*T
      DK = (-5.03 + 0.0814*T)/1000.0
      DV = - (DV - DVW)
      DK = - (DK - DKW)
      LNK9 = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNK9
      K(9) = EXP(LNK9)
C-------------------------------------------------------------------
C     7) HYDROGEN SULFIDE
C-------------------------------------------------------------------
      LNKS = 225.838 - 13275/TK - 34.6435*LNTK + 0.3449*SQS - 0.0274*S
      DV = - 11.07 + 0.009*T - 0.942-03*T*T
      DK = (-2.89 + 0.054*T)/1000.0
      LNKS = - DV/(R*TK)*PP + 0.5*DK/(R*TK)*PP*PP + LNKS
      K(10) = EXP(LNKS)
C-------------------------------------------------------------------
C     8) CALCITE & ARAGONITE
C-------------------------------------------------------------------
      LOGK12 = - 171.9065 - 0.077993*TK + 2839.319/TK + 71.595*LOG10(TK)
      LOGK12 = LOGK12 + (-0.77712 + 2.8426D-03*TK + 178.34/TK)*SQS
     &         - 0.07711*S + 4.1249D-03*S23
      K(12) = 10**(LOGK12)
      LNK12 = LOG(K(12))
      DVC = - (48.76 - 0.5304*T)
      DKC = - (11.76 - 0.3692*T)/1000.0
      LNK12 = - DVC/(R*TK)*PP + 0.5*DKC/(R*TK)*PP*PP + LNK12
      K(12) = EXP(LNK12)
      LOGK11 = - 171.945 - 0.077993*TK + 2903.293/TK + 71.595*LOG10(TK)
      LOGK11 = LOGK11 + (-0.068393 + 1.7276D-03*TK + 88.135/TK)*SQS
     &         - 0.10018*S + 5.9415D-03*S23
      K(11) = 10**(LOGK11)
      DVA = DVC + 2.8
      DKA = DKC
      LNK11 = LOG(K(11))
      LNK11 = - DVC/(R*TK)*PP + 0.5*DKC/(R*TK)*PP*PP + LNK11
      K(11) = EXP(LNK11)
C
      K(1) = K(1)*1.0D+06
      K(2) = K(2)*1.0D+03
      K(3) = K(3)*1.0D+03
      K(4) = K(4)*1.0D+03
      K(5) = K(5)*1.0D+03
      K(6) = K(6)*1.0D+03
      K(7) = K(7)*1.0D+03
      K(8) = K(8)*1.0D+03
      K(9) = K(9)*1.0D+03
      K(10) = K(10)*1.0D+03
      K(11) = K(11)*1.0D+06
      K(12) = K(12)*1.0D+06
      RETURN
      END