      PROGRAM WATER
C =====================================================================|

      INTEGER ID, IH, II(40), IOPT, IP, IT, JJ(40), NC

      REAL * 8 AA, AD, BASE
      REAL * 8 C, CJTH, CJTT, CP, CPD, CV, CVD
      REAL * 8 D, DD, DGSS, DLL, DPDD, DPDT, DPDT1
      REAL * 8 DQ, DVDT, DVV, DZ
      REAL * 8 FD, FH, FP, FT
      REAL * 8 G(40), GASCON, GD, H, HD
      REAL * 8 P, PRES, PSAT
      REAL * 8 Q0, Q5, RT
      REAL * 8 S, SD, SREF, T, TT, TTT, TZ
      REAL * 8 U, UD, UREF, VL, WM, X, Y, Z, ZDUM

      CHARACTER * 1 NT
      CHARACTER * 3 NS, NSS1, NSS2
      CHARACTER * 8 ND, NH, NP

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF
      COMMON /FCTS  / AD, GD, SD, UD, HD, CVD, CPD, DPDT, DVDT, CJTT, 
     &                CJTH
      COMMON /NCONST/ G, II, JJ, NC
      COMMON /QQQQ  / Q0, Q5
      COMMON /UNITS / IT, ID, IP, IH, NT, ND, NP, NH, FT, FD, FP, FH

      DATA NSS1 /' M '/
      DATA NSS2 /'HFT'/

      EXTERNAL BASE, BB, DFIND, PCORR, QQ, THERM, TTT, UNIT

C =====================================================================|

      CALL UNIT

      NS = NSS1
      IF (ID .EQ. 4) NS = NSS2

  100 CONTINUE
      WRITE (6,10)
      WRITE (6, '(/)')

      READ (5, *, END = 200) IOPT, X, TT
      IF (IOPT .LE. 0) GOTO 200

      T = TTT(TT)
      RT = GASCON * T
      CALL BB(T)

      IF (IOPT .EQ. 1) THEN
        DD = X
        D  = DD * FD

        CALL QQ(T,D)
        ZDUM = BASE(D,T)
        PRES = FP * (RT * D * Z + Q0)
        DQ   = RT * (Z + Y * DZ) + Q5
      ELSE IF (IOPT .EQ. 2) THEN
        PRES = X
        P = PRES / FP
        DGSS = P / T / .4D0
        PSAT = 2.0D4
        DLL = 0.0D0
        DVV = 0.0D0
        IF (T .LE. TZ) CALL PCORR (T, PSAT, DLL, DVV)
        IF (P .GT. PSAT) DGSS = DLL
        CALL DFIND (D, P, DGSS, T, DQ)
        DD = D / FD
      ELSE 
        GOTO 100
      END IF

      CALL THERM (D, T)
      U = UD * RT * FH
      C = DSQRT(DABS(CPD * DQ * 1.0D3 / CVD))
      IF (ID .EQ. 4) C = C * 3.280833D0
      H  = HD * RT * FH
      S  = SD * GASCON * FH * FT
      CP = CPD * GASCON * FH * FT
      CV = CVD * GASCON * FH * FT
      VL = 1.0D0 / D
      DPDD = DQ * FD * FP
      DPDT1 = DPDT * FP * FT

      WRITE (6, 15) TT, NT
      WRITE (6, 20) PRES, NP
      WRITE (6, 25) DD, ND
      WRITE (6, 30) DPDT1, DPDD
      WRITE (6, 35) CV, NH, NT 
      WRITE (6, 40) CP, S
      WRITE (6, 45) H, NH, U
      WRITE (6, 50) C, NS
      WRITE (6, 55) CJTT, CJTH
      WRITE (6, 60) DVDT

      GOTO 100
  200 CONTINUE
      STOP

   10 FORMAT (' Enter Option, X, and T, Where for option = 1, ', 
     &        'X = density and '/
     &        '                               for option = 2, ',
     &        'X = pressure. (Enter 0 to quit)')
   15 FORMAT (' T = ', F12.4, ' DEG ', A1)
   20 FORMAT (' P = ', F13.6, 1X, A6)
   25 FORMAT (' D = ', F14.10, 1X, A8)
   30 FORMAT (' DP/DT = ', F16.9, 6X, 'DP/DD = ', F16.5)
   35 FORMAT (' CV = ', F12.6, 1X, A6, A1, 5X)
   40 FORMAT (' CP = ', F12.6, 6X, 'S = ', F12.6)
   45 FORMAT (' H = ', F14.6, 1X, A6, 5X, 'U = ', F14.6)
   50 FORMAT (' VEL SND = ', F14.6, A2, '/SEC')
   55 FORMAT (' JT(T) = ', F11.5, 5X, 'JT(H) = ', F11.5)
   60 FORMAT (' DV/DT = ', F12.6)

      END




      SUBROUTINE BB (T)
C =====================================================================|
C     This subroutine calculates the b's of eqs. 3 & 4 using           |
C     coefficients from blockdata.  It also calculates the first and   |
C     second derivatives with respest to temperature.  The b's         |
C     calculated are in cm3/g                                          |
C ---------------------------------------------------------------------|

      INTEGER I

      REAL * 8 AA
      REAL * 8 B1, B1T, B1TT, B2, B2T, B2TT, BP(10), BQ(10)
      REAL * 8 DZ, G1, G2, GASCON, GF
      REAL * 8 SREF, T, TZ, UREF, V(10), WM, Y, Z

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF
      COMMON /BCONST/ BP, BQ
      COMMON /ELLCON/ G1, G2, GF, B1, B2, B1T, B2T, B1TT, B2TT

C =====================================================================|

      V(1) = 1.0D0
      DO 10 I = 2, 10
        V(I) = V(I-1) * TZ / T
   10 CONTINUE

      B1 = BP(1) + BP(2) * DLOG (1.0D0 / V(2))
      B1T = BP(2) * V(2) / TZ
      B1TT = 0.0D0

      B2 = BQ(1)
      B2T = 0.0D0
      B2TT = 0.0D0

      DO 20 I = 3, 10
        B1 = B1 + BP(I) * V(I-1)
        B1T = B1T - (I-2) * BP(I) * V(I-1) /T
        B1TT = B1TT + BP(I) * (I-2)**2 * V(I-1) / T / T

        B2   = B2 + BQ(I) * V(I-1)
        B2T  = B2T - (I-2) * BQ(I) * V(I-1) / T
        B2TT = B2TT + BQ(I) * (I-2)**2 * V(I-1) / T / T
   20 CONTINUE

      B1TT = B1TT - B1T / T
      B2TT = B2TT - B2T / T

      RETURN
      END





      SUBROUTINE CORR (T, P, DL, DV, DELG)
C =====================================================================|
C     Subroutine corr will calculate, for an input T and P at or       |
C     near the vapor pressure, the correspoinding liquid and vapor     |
C     densities and also del g = (gl -gv) / rt for use in calculating  |
C     the correction to the vapor pressure for del g = 0.              |
C ---------------------------------------------------------------------|

      REAL * 8 AA, AD, BASE, CJTH, CJTT, CPD, CVD
      REAL * 8 DELG, DL, DLIQ, DPDT, DQ, DV, DVAP, DVDT, DZ
      REAL * 8 GASCON, GD, GL, GV, HD
      REAL * 8 P, Q0, Q5, RT, SD, SREF, T, TAU, TZ
      REAL * 8 UD, UREF, WM, Y, Z, ZDUM

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF
      COMMON /FCTS  / AD, GD, SD, UD, HD, CVD, CPD, DPDT, DVDT, CJTT, 
     &                CJTH
      COMMON /QQQQ  / Q0, Q5

C =====================================================================|

      RT = GASCON * T
      IF (T .GT. 646.3D0) GO TO 101
      DLIQ = DL 
      IF (DL .LE. 0.0D0) DLIQ = 1.11D0 - .0004D0 * T
      CALL BB(T)
      CALL DFIND (DL, P, DLIQ, T, DQ)
      CALL THERM (DL, T)
      GL = GD
      DVAP = DV
      IF (DV .LE. 0.0D0) DVAP = P / RT
      CALL DFIND (DV, P, DVAP, T, DQ)
      IF (DV .LT. 5.0D-7) DV = 5.0D-7
      CALL THERM (DV, T)
      GV = GD
      DELG = GL - GV
      RETURN

  101 CONTINUE
      P = 0.0D0
      IF (T .GT. 647.126D0) RETURN
      DELG = 0.0D0
      CALL BB(T)
      TAU = 0.657128D0 * (1.0D0 - T / 647.126D0)**0.325D0
      DL = 0.322D0 + TAU
      DV = 0.322D0 - TAU
      ZDUM = BASE(DV, T)
      CALL QQ(T, DV)
      P = RT * DV * Z + Q0

      RETURN
      END





      SUBROUTINE DFIND (DOUT, P, D, T, DPD)
C =====================================================================|
C     Routine to find density corresponding to input pressure p(mpa),  |
C     and temperature t(k), using initial guess denstiy d(g/cm3).      |
C     The output density is in g/cm3, also, as a byproduct, DP/DRH0 is |
C     calculated ("DPD", MPA cm3/g)                                    |
C ---------------------------------------------------------------------|

      INTEGER L

      REAL * 8 AA, BASE, D, DD, DOUT, DP, DPD, DPDX, DZ, GASCON
      REAL * 8 P, PP, Q0, Q5, RT, SREF, T, TZ, UREF, WM, X, Y, Z

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF
      COMMON /QQQQ  / Q0, Q5

C =====================================================================|

      DD = D
      RT = GASCON * T
      IF (DD .LE. 0.0D0) DD = 1.0D-08
      IF (DD .GT. 1.9D0) DD = 1.9D+00

      L = 0
    9 CONTINUE
      L = L + 1
   11 CONTINUE
      IF (DD .LE. 0.0D0) DD = 1.0D-08
      IF (DD .GT. 1.9D0) DD = 1.9D+00
      CALL QQ (T, DD)
      PP = RT * DD * BASE (DD, T) + Q0
      DPD = RT * (Z + Y * DZ) + Q5

C ---------------------------------------------------------------------|
C     The following 3 lines check for negative DP/DRH0, and if so      |
C     assume guess to be in 2-phase region, and adjust guess           |
C     accordingly.                                                     |
C ---------------------------------------------------------------------|

      IF (DPD .GT. 0.0000D0) GO TO 13
      IF (D   .GE. 0.2967D0) DD = DD * 1.02D0
      IF (D   .LT. 0.2967D0) DD = DD * 0.98D0
      IF (L .LE. 10) GO TO 9

   13 CONTINUE
      DPDX = DPD * 1.1D0
      IF (DPDX .LT. 0.1D0) DPDX = 0.1D0
      DP = DABS(1.0D0 - PP/P)
      IF (DP .LT. 1.0D-08) GO TO 20
      IF (D .GT. 0.3D0 .AND. DP .LT. 1.0D-07) GO TO 20
      IF (D .GT. 0.7D0 .AND. DP .LT. 1.0D-06) GO TO 20
      X = (P - PP) / DPDX
      IF (DABS(X) .GT. 0.1D0) X = X * 0.1D0 / DABS(X)
      DD = DD + X
      IF (DD .LE. 0.0D0) DD = 1.0D-08
   19 CONTINUE
      IF (L .LE. 30) GO TO 9

   20 CONTINUE
      DOUT = DD

      RETURN
      END




      SUBROUTINE IDEAL (T)
C =====================================================================|
C     This subroutine calculates the thermodynamic properties for      |
C     water in the ideal gas state from function of H.W. Woolley       |
C ---------------------------------------------------------------------|
      
      INTEGER I 

      REAL * 8 AI, C(18), CPI, CVI, GI, HI, SI, T, TL, TT, UI

      COMMON /IDF   / AI, GI, SI, UI, HI, CVI, CPI

      DATA C /19.730271018D0, 20.9662681977D0, -0.483429455355D0, 
     &         6.05743189245D0, 22.56023885D0, -9.87532442D0, 
     &        -4.3135538513D0, 0.458155781D0, -0.47754901883D-1, 
     &         0.41238460633D-2, -0.27929052852D-3, 0.14481695261D-4,
     &        -0.56473658748D-6, 0.16200446D-7, -0.3303822796D-9, 
     &         0.451916067368D-11, -0.370734122708D-13,
     &         0.137546068238D-15/
C =====================================================================|
      
      TT = T / 100.0D0
      TL = DLOG(TT)

      GI  = -(C(1) / TT + C(2)) * TL
      HI  =   C(2) + C(1) * (1.0D0 - TL) / TT
      CPI =   C(2) - C(1) / TT

      DO 10 I = 3, 18
        GI  = GI  - C(I) * TT**(I-6)
        HI  = HI  + C(I) * (I-6) * TT**(I-6)
        CPI = CPI + C(I) * (I-6) * (I-5) * TT**(I-6)
   10 CONTINUE

      AI  = GI  - 1.0D0
      UI  = HI  - 1.0D0
      CVI = CPI - 1.0D0
      SI  = UI - AI

      RETURN 
      END




      
      SUBROUTINE PCORR (T, P, DL, DV)
C =====================================================================|
C     Subroutine pcorr will calculate the vapor pressure P and the     |
C     liquid, and the vapor densities corresponding to the input T,    |
C     corrected such that gl - gv = 0.  The function ps is required    |
C     which will give a reasonably good approximation to the vapor     |
C     pressure to be used as the starting point for the iteration.     |
C ---------------------------------------------------------------------|

      REAL * 8 AA, DELG, DL, DP, DV, DZ, GASCON, P, PS, SREF
      REAL * 8 T, TZ, WM, UREF, Y, Z

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF

C =====================================================================|

      P = PS(T)
   10 CONTINUE
        CALL CORR (T, P, DL, DV, DELG)
        DP = DELG * GASCON * T / (1.0D0 / DV - 1.0D0 / DL)
        P = P + DP
        IF (DABS(DELG) .LT. 1.0D-4) RETURN
      GO TO 10
      END





      SUBROUTINE QQ (T, D)
C =====================================================================|
C     This routine calculates, for a given t(k) and d(g/cm3), the      |
C     residual contributions to:  pressure(q), helmholtz fct(ar),      |
C     dp/drh0(q5), and also the gibbs function, entropy, internal      |
C     energy, enthalpy, isochoric heat capacity and dpdt. (eqn. 5)     |
C     Terms 37 thru 39 are the additional terms affecting only the     |
C     immediate vicinity of the critical point, and term 40 is the     |
C     additional term improving the low t, high p region.              |
C ---------------------------------------------------------------------|

      INTEGER I, II(40), J, JJ(40), K, KM, L, N

      REAL * 8 AA, AAD(4), AAT(4), ADZ(4), AR, ATT, ATZ(4), CVR
      REAL * 8 D, D2F, DADT, DD, DDZ, DEL, DEX, DFDT, DPDTR, DPT, DZ
      REAL * 8 E, EX1, EX2, FCT
      REAL * 8 G(40), GASCON, GR, HR
      REAL * 8 Q, Q0, Q2A, Q5, Q5T, Q10, Q20, QM, QP
      REAL * 8 QR(11), QT(10), QZR(0:9), QZT(9)
      REAL * 8 RT, SR, SREF, T, TAU, TEX, TX, TZ
      REAL * 8 UR, UREF, V, WM, Y, Z, ZZ

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF
      COMMON /ADDCON/ ATZ, ADZ, AAT, AAD
      COMMON /NCONST/ G, II, JJ, N
      COMMON /RESF  / AR, GR, SR, UR, HR, CVR, DPDTR
      COMMON /QQQQ  / Q0, Q5

      EQUIVALENCE (QR(2), QZR(0)), (QT(2), QZT(1))

C =====================================================================|

      RT    = GASCON * T
      QR(1) = 0.0D+0
      Q5    = 0.0D+0
      Q     = 0.0D+0
      AR    = 0.0D+0
      DADT  = 0.0D+0
      CVR   = 0.0D+0
      DPDTR = 0.0D+0

      E     = DEXP(-AA * D)
      Q10   = D * D * E
      Q20   = 1.0D+0 - E
      QR(2) = Q10
      V     = TZ / T
      QT(1) = T / TZ

      DO 4 I = 2, 10
        QR(I+1) = QR(I) * Q20
        QT(I)   = QT(I-1) * V
    4 CONTINUE

      DO 10 I = 1, N
        K = II(I) + 1
        L = JJ(I)
        ZZ = K
        QP = G(I) * AA * QZR(K-1) * QZT(L)
        Q = Q + QP
        Q5 = Q5 + AA * (2.0D0 / D -
     &            AA * (1.0D0 - E * (K-1) / Q20)) * QP
        AR = AR + G(I) * QZR(K) * QZT(L) / Q10 / ZZ / RT
        DFDT = Q20**K * (1-L) * QZT(L+1) /TZ / K
        D2F = L * DFDT
        DPT = DFDT * Q10 * AA * K / Q20
        DADT = DADT + G(I) * DFDT
        DPDTR = DPDTR + G(I) * DPT
        CVR = CVR + G(I) * D2F / GASCON
   10 CONTINUE

       QP = 0.0D+0
       Q2A = 0.0D+0

       IF (T .LT. 642 .OR. T .GT. 652) GO TO 21

       DO 20 J = 37, 40
         IF (G(J) .EQ. 0.0D+0) GOTO 20
         K = II(J)
         KM = JJ(J)
         DDZ = ADZ (J-36)
         DEL = D / DDZ - 1.0D+0
         IF (DABS(DEL) .LT. 1.0D-10) DEL = 1.0D-10
         DD = DEL * DEL
         EX1 = -AAD(J-36) * DEL**K
         DEX = DEXP(EX1) * DEL**KM
         ATT = AAT(J-36)
         TX = ATZ(J-36)
         TAU = T / TX - 1.0D+0
         EX2 = -ATT * TAU * TAU
         TEX = DEXP(EX2)
         Q10 = DEX * TEX
         QM = KM / DEL - K * AAD(J-36) * DEL**(K-1)
         FCT = QM * D**2 * Q10 / DDZ
         Q5T = FCT * (2.0D0 / D + QM / DDZ) - (D / DDZ)**2 * Q10 *
     &         (KM / DEL / DEL + K * (K-1) * AAD(J-36) * DEL**(K-2))
         Q5 = Q5 + Q5T * G(J)
         QP = QP + G(J) * FCT
         DADT = DADT - 2.0D+0 * G(J) * ATT * TAU * Q10 / TX
         DPDTR = DPDTR - 2.0D+0 * G(J) * ATT * TAU * FCT /TX
         Q2A = Q2A + T * G(J) * (4.0D+0 * ATT * EX2 + 2.0D+0 * ATT) *
     &         Q10 / TX / TX
         AR = AR + Q10 * G(J) / RT
   20 CONTINUE

   21 CONTINUE
      SR = - DADT / GASCON
      UR = AR + SR
      CVR = CVR + Q2A / GASCON
      Q = Q + QP

      RETURN
      END





      SUBROUTINE THERM (D, T)
C =====================================================================|
C     This subroutine calculates the thermodynamic functions in        |
C     dimensionless units (AD = A/RT, GD = G / RT, etc.                |
C ---------------------------------------------------------------------|

      REAL * 8 AA, AB, AD, AR, AI
      REAL * 8 CJTH, CJTT, CPD, CPI, CVB, CVD, CVI, CVR
      REAL * 8 D, DPDD, DPDTB, DPDT, DPDTR, DVDT, DZB
      REAL * 8 GASCON, GB, GD, GI, GR, HB, HD, HI, HR
      REAL * 8 Q0, Q5, RT
      REAL * 8 SB, SD, SI, SR, SREF
      REAL * 8 T, TZ, UB, UD, UI, UR, UREF, WM
      REAL * 8 Y, Z, ZB

      COMMON /ACONST/ WM, GASCON, TZ, AA, ZB, DZB, Y, UREF, SREF
      COMMON /BASEF / AB, GB, SB, UB, HB, CVB, DPDTB
      COMMON /FCTS  / AD, GD, SD, UD, HD, CVD, CPD, DPDT, DVDT, CJTT, 
     &                CJTH
      COMMON /IDF   / AI, GI, SI, UI, HI, CVI, CPI
      COMMON /QQQQ  / Q0, Q5
      COMMON /RESF  / AR, GR, SR, UR, HR, CVR, DPDTR

C =====================================================================|

      CALL IDEAL (T)

      RT = GASCON * T
      Z = ZB + Q0 / RT / D
      DPDD = RT * (ZB + Y * DZB) + Q5
      AD = AB + AR + AI - UREF / T + SREF
      GD = AD + Z
      UD = UB + UR + UI - UREF / T
      DPDT = RT * D * DPDTB + DPDTR
      CVD = CVB + CVR + CVI
      CPD = CVD + T * DPDT**2 / (D * D * DPDD * GASCON)
      HD = UD + Z
      SD = SB + SR + SI - SREF
      DVDT = DPDT / DPDD / D / D
      CJTT = 1.0D0 / D - T * DVDT
      CJTH = CJTT / CPD / GASCON


      RETURN
      END




      SUBROUTINE UNIT
C =====================================================================|
C     This subroutine queries the user as to his choice of units and   |
C     sets internal parameters appropriately.  The internal units of   |
C     the program are temperatures in K, densities in g/cm3, all other |
C     quantities are calculated in dimensionless units and dimensioned |
C     at output time.                                                  |
C ---------------------------------------------------------------------|

      INTEGER ID, IH, IP, IT
 
      REAL * 8 FD, FFD(4), FFH(6), FFP(5), FH, FP, FT

      CHARACTER * 1 NT, NNT(4)
      CHARACTER * 8 A1, A2, A3, A4
      CHARACTER * 8 ND, NH, NND(4), NNH(6), NNP(5), NP

      COMMON /UNITS / IT, ID, IP, IH, NT, ND, NP, NH, FT, FD, FP, FH

      DATA A1 /'Temperat'/
      DATA A2 /'Density '/
      DATA A3 /'Pressure'/
      DATA A4 /'Energy  '/
      DATA FFD /1.0D-3, 1.0D+0, 1.80152D-2, 1.6018D-2/
      DATA FFH /1.0D+0, 1.0D+0, 1.80152D+1, 0.2388459D+0, 
     &          4.30285666D+0, 0.4299226D0/
      DATA FFP /1.0D+0, 1.0D+1, 9.869232667D+0, 145.038D+0, 10.1971D+0/
      DATA NND /' Kg/m3  ', ' g/cm3  ', ' mol/L  ', ' lb/ft3 '/
      DATA NNH /' kJ/kg  ', '  J/g   ', ' J/mol  ', ' cal/g  ',
     &          'cal/mol ', ' BTU/lb '/ 
      DATA NNP /'  MPa   ', '  Bar   ', '  Atm   ', '  PSI   ', 
     &          ' kg/cm2 '/
      DATA NNT /'K', 'C', 'R', 'F'/

C =====================================================================|

      WRITE (6, 11) A1
   30 CONTINUE
      WRITE (6, 12)
      READ (5, *, END = 99) IT
      IF (IT .EQ. 0) STOP
      IF (IT .GT. 4) GOTO 30
      NT = NNT(IT)

      WRITE (6, 11) A2
   31 CONTINUE 
      WRITE (6, 13)
      READ (5, *, END = 99) ID
      IF (ID .GT. 4 .OR. ID .LT. 1) GOTO 31
      ND = NND(ID)
      FD = FFD(ID)

      WRITE (6, 11) A3
   32 CONTINUE 
      WRITE (6, 14)
      READ (5, *, END = 99) IP
      IF (IP .GT. 5 .OR. IP .LT. 1) GOTO 32
      NP = NNP(IP)
      FP = FFP(IP)

      WRITE (6, 11) A4
   33 CONTINUE 
      WRITE (6, 15)
      READ (5, *, END = 99) IH
      IF (IH .GT. 6 .OR. IH .LT. 1) GOTO 33
      NH = NNH(IP)
      FH = FFH(IP)

      RETURN
   99 CONTINUE
      STOP 

   11 FORMAT (' Enter units chosen for ', A8)
   12 FORMAT (' Choose from 1 = Deg K, 2 = Deg C, 3 = Deg R, 4 = Deg F')
   13 FORMAT (' Choose from 1 = Kg/m3, 2 = g/cm3, 3 = mol/L,',
     &        ' 4 = lb/ft3')
   14 FORMAT (' Choose from 1 = MPa, 2 = Bar, 3 = atm, 4 = PSIA,', 
     &        ' 5 = Kg/cm2')
   15 FORMAT (' Choose from 1 = KJ/Kg, 2 = J/g, 3 = J/mol, 4 = Cal/g,'
     &        ' 5 = Cal/mol, 6 = BTU/lb')

      END





      REAL * 8 FUNCTION BASE (D, T)
C =====================================================================|
C     This function calculates Z (=pbase/(DRT) of eq. q)               |
C     (called BASE),  and also Abase, Gbase, Sbase, Ubase, Hbase,      |
C     CVbase, and 1/(DRT) * Dp / Dt for the base fct (called DPDTB)    |
C     The AB, GB, SB, UB, HB, and CVB are calculated in dimensionless  |
C     units:  AB/RT, GB/RT, SB/R, etc.                                 |
C ---------------------------------------------------------------------|

      REAL * 8 AA, AB, B1, B2, B1T, B2T, B1TT, B2TT, BB2TT, CVB
      REAL * 8 D, DPDTB, DZ, DZ0
      REAL * 8 G1, G2, GASCON, GB, GF, HB
      REAL * 8 SB, SREF, T, TZ, UB, UREF, WM, X, Y, Z, Z0

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF
      COMMON /BASEF / AB, GB, SB, UB, HB, CVB, DPDTB
      COMMON /ELLCON/ G1, G2, GF, B1, B2, B1T, B2T, B1TT, B2TT

C ---------------------------------------------------------------------|
C     G1, G2, and GF are the alpha, beta, and gamma of eq. 2, which    |
C     are supplied by the blockdata routine.  B1 and B2 are the        |
C     "excluded volume" and "2nd virial" (eqs 3 and 4) supplied by the |
C     subroutine bb(t), which also supplies the 1st and 2nd            |
C     derivatives with respect to t                                    |
C ---------------------------------------------------------------------|

C =====================================================================|

      Y   = 0.25D0 * B1 * D
      X   =  1.0D0 - Y
      Z0  = (1.0D0 + G1 * Y + G2 * Y * Y) / X**3
      Z   = Z0 + 4.0D0 * Y * (B2 / B1 - GF)
      DZ0 = (G1 + 2.0D0 * G2 * Y) / X**3 + 
     &      3.0D0 * (1.0D0 + G1 * Y + G2 * Y * Y ) / X**4
      DZ  = DZ0 + 4.0D0 * (B2 / B1 - GF)

      AB = -DLOG(X) - (G2 - 1.0D0) / X + 28.16666667D0 / X / X + 
     &     4.0D0 * Y * (B2 / B1 - GF) + 15.166666667D0 + 
     &      DLOG(D * T * GASCON / 0.101325D0)
      GB = AB + Z

      BASE = Z
      BB2TT = T * T * B2TT
      UB = -T * B1T * (Z - 1.0D0 - D * B2) / B1 - D * T * B2T
      HB = Z + UB
      CVB = 2.0D0 * UB + (Z0 - 1.0D0) * ((T * B1T / B1)**2 - 
     &      T * T * B1TT / B1) - D * (BB2TT - GF * B1TT * T * T) - 
     &      (T * B1T / B1)**2 * Y * DZ0
      DPDTB = BASE / T + BASE * D / Z * (DZ * B1T / 4.0D0 + 
     &        B2T - B2 / B1 * B1T)
      SB = UB - AB

      RETURN
      END





      REAL * 8 FUNCTION PS (T)
C =====================================================================|
C     This function calculates an approximation to the vapor pressure, |
C     PS, as a function of the input temperature.  The vapor pressure  |
C     calculated agrees with the vapor pressure predicted by the       |
C     surface to within 0.02% to within a degree or so of the critical |
C     temperature, and can serve as an initial guess for further       |
C     refinement by imposing the condition that g1 = gv.               |
C ---------------------------------------------------------------------|

      INTEGER I

      REAL * 8 A(8), B, PL, Q, T, V, W, Z

      DATA A /-7.8889166D0, 2.5514255D0, -6.716169D0, 33.239495D0,
     &        -105.38479D0, 174.35319D0, -148.39348D0, 48.631602D0/


C =====================================================================|

      IF (T .GT. 314.0D0) THEN
        V = T / 647.25D0
        W = DABS (1.0D0 - V)
        B = 0.0D0
        DO 10 I = 1, 8
          Z = I
          B = B + A(I) * W**((Z + 1.0D0) / 2.0D0)
   10   CONTINUE
        Q = B / V
        PS = 22.093D0 * DEXP(Q)
      ELSE 
        PL = 6.3573118D0 - 8858.843D0 / T + 607.56335D0 * T**(-0.6D0)
        PS = 0.1D0 * DEXP(PL)
      END IF      

      RETURN
      END





      REAL * 8 FUNCTION TSAT (P)
C =====================================================================|
C     This function calculates the saturation temperature for a given  |
C     pressure, by an iterative process using PS and TDPSDT.           |
C ---------------------------------------------------------------------|

      INTEGER K

      REAL * 8 DP, P, PL, PP, PS, TDPSDT, TG

C =====================================================================|

      TSAT = 0.0D0
      IF (P .GT. 22.05D0) RETURN
      K = 0

C ---------------------------------------------------------------------|
C     Convert from bars to MPa                                         |
C ---------------------------------------------------------------------|

      PL = 2.302585D0 + DLOG(P)

      TG = 372.83D0 + PL * (27.7589D0 + PL * (2.3819D0 + 
     &                PL * (0.24834D0 + PL * 0.0193855D0)))

   10 CONTINUE
      IF (TG .LT. 273.15D0) TG = 273.15D0
      IF (TG .GT. 647.00D0) TG = 647.0D0
      IF (K .LT. 8) GO TO 20
      WRITE (6, *) K, P, PP, TG
      GO TO 30

   20 CONTINUE
      K = K + 1
      PP = PS(TG)
      DP = TDPSDT(TG)
      IF (ABS(1.0D0 - PP / P) .LT. 1.0D-05) GO TO 30
      TG = TG * (1.0D0 + (P - PP) / DP)
      GO TO 10

   30 CONTINUE
      TSAT = TG

      RETURN
      END

      



      REAL * 8 FUNCTION TDPSDT (T)
C =====================================================================|
C     This function calculates t * (dPs / dt), and is used by the      |
C     function tsat                                                    |
C ---------------------------------------------------------------------|

      INTEGER I

      REAL * 8 A(8), B, C, Q, T, V, W, Y, Z

      DATA A /-7.8889166D0, 2.5514255D0, -6.716169D0, 33.239495D0,
     &        -105.38479D0, 174.35319D0, -148.39348D0, 48.631602D0/

C =====================================================================|

      V = T / 647.25D0
      W = 1.0D0 - V
      B = 0.0D0
      C = 0.0D0

      DO 10 I = 1, 8
        Z = I
        Y = A(I) * W**((Z + 1.0D0) / 2.0D0)
        C = C + Y / W * (0.5D0 - 0.5D0 * Z - 1.0D0 / V)
        B = B + Y
   10 CONTINUE

      Q = B / V
      TDPSDT = 22.093D0 * DEXP(Q) * C

      RETURN
      END





      REAL * 8 FUNCTION TTT(T)
C =====================================================================|
C     Function to convert input temperatures in external units         |
C     to deg K                                                         |
C ---------------------------------------------------------------------|

      INTEGER ID, IH, IP, IT

      REAL * 8 FD, FH, FP, FT
      REAL * 8 T

      CHARACTER * 1 NT
      CHARACTER * 8 ND, NH, NP

      COMMON /UNITS / IT, ID, IP, IH, NT, ND, NP, NH, FT, FD, FP, FH

C =====================================================================|

      IF (IT .EQ. 1) THEN
        TTT = T
        FT  = 1.0D+0
        RETURN
      ELSE IF (IT .EQ. 2) THEN 
        TTT = T + 273.15D+0
        FT  = 1.0D+0
        RETURN
      ELSE IF (IT .EQ. 3) THEN
        TTT = T / 1.8D+0
        FT  = 0.5555555555556D+0
        RETURN 
      ELSE IF (IT .EQ. 4) THEN
        TTT = (T + 459.67D+0) / 1.8D+0
        FT  = 0.5555555555556D+0
        RETURN 
      END IF 

      STOP 
      END




      BLOCK DATA
C =====================================================================|
C     This blockdata subroutine supplies most of the fixed parameters  |
C     used in the rest of the routines.  BP is the b(i) of table I,    |
C     BQ is the B(i) of table I, G1, G2, and GF are the alpha, beta,   |
C     and gamma of eq. 2, and G, II, JJ are the g(i), k(i), and l(k)   |
C     of eq. 5                                                         |
C ---------------------------------------------------------------------|

      INTEGER II(40), JJ(40), NC

      REAL * 8 AA, AAD(4), AAT(4), ADZ(4), ATZ(4)
      REAL * 8 B1, B1T, B1TT, B2, B2T, B2TT, BP(10), BQ(10)
      REAL * 8 DZ, G(40), G1, G2, GASCON, GF
      REAL * 8 SREF, TZ, UREF, WM, Y, Z

      COMMON /ACONST/ WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF
      COMMON /ADDCON/ ATZ, ADZ, AAT, AAD
      COMMON /BCONST/ BP, BQ
      COMMON /ELLCON/ G1, G2, GF, B1, B2, B1T, B2T, B1TT, B2TT
      COMMON /NCONST/ G, II, JJ, NC

C ---------------------------------------------------------------------|
C     See table A.2                                                    |
C ---------------------------------------------------------------------|

      DATA AAD /34.0D0, 4.0D1, 3.0D1, 1.05D3/
      DATA AAT /2.0D4, 2.0D4, 4.0D4, 25.0D0/
      DATA ADZ / 0.319D0, 0.319D0, 0.319D0, 1.55D0/
      DATA ATZ /64.0D1,  64.0D1, 641.6D0,  27.0D1/

      DATA AA /1.0D0/
      DATA BP /0.7478629D0, -0.3540782D0, 0.0D0, 0.0D0, 0.7159876D-2,
     &         0.0D0, -0.3528426D-2, 0.0D0, 0.0D0, 0.0D0/
      DATA BQ /1.1278334D0, 0.0D0, -0.5944001D0, -5.010996D0, 0.0D0,
     &         0.63684256D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0/
      DATA G /-530.62968529023D0, 2274.4901424408D0,  787.79333020687D0, 
     &        -69.830527374994D0, 17863.832875422D0, -39514.731563338D0, 
     &         33803.884280753D0, -13855.050202703D0, -256374.3661326D0, 
     &         482125.75981415D0, -341830.1696966D0, 122231.56417448D0, 
     &         1179743.3655832D0, -2173481.0110373D0, 1082995.216862D0, 
     &        -254419.98064049D0, -3137777.4947767D0, 5291191.0757704D0, 
     &        -1380257.7177877D0, -251099.14369001D0, 4656182.6115608D0, 
     &        -7275277.3275387D0, 417742.46148294D0, 1401635.8244616D0, 
     &        -3155523.1392127D0, 4792966.6384584D0, 409126.64781209D0, 
     &        -1362636.9388386D0, 696252.20862664D0, -1083490.0096447D0, 
     &        -227228.27401688D0, 383654.8600066D0, 6883.3257944332D0, 
     &         21757.245522644D0, -2662.794482977D0, -70730.418082074D0, 
     &        -0.225D0, -1.68D0, 0.055D0, -93.0D0/
      DATA G1 /11.0D0/
      DATA G2 /44.333333333333D0/
      DATA GASCON /0.461522D0/
      DATA GF /3.5D0/
      DATA II /0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 
     &         4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 8, 8, 8, 8, 
     &         2, 2, 0, 4, 2, 2, 2, 4/
      DATA JJ /2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 
     &         2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7,
     &         1, 4, 4, 4, 0, 2, 0, 0/
      DATA NC   /36/
      DATA SREF /7.6180802D0/
      DATA TZ   /647.073D0/
      DATA UREF /-4328.455039D0/
      DATA WM   /18.0152D0/

      END






