      PROGRAM DIEWAT
C =====================================================================|

      INTEGER I, J
 
      REAL * 8 ALPHA, BETA
      REAL * 8 DADT, DETADP, DETADT, DETDT2
      REAL * 8 DEDP, DEDT, D2EDT2, DIECAL
      REAL * 8 E(0:4,0:4), ETA, RHO, T, TC
      REAL * 8 G, PRESS, Q, Y, X

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

      ALPHA = 0.0D+0
      BETA  = 0.0D+0
      G     = 0.0D+0
      PRESS = 0.0D+0
      RHO   = 0.0D+0

      DO 20 I = 0, 4
        DO 10 J = 0, 4
          E(I,J) = 0.0D+0
   10   CONTINUE
   20 CONTINUE

      E(0,0) =   4.39109592E+02    
      E(0,1) =  -2.18995148E+02
      E(0,2) =   1.82898246E+01
      E(0,3) =   1.54886800E+02
      E(0,4) =  -6.13542375E+01

      E(1,0) =  -2.33277456E+00 
      E(1,1) =   1.00498361E+00
      E(1,2) =  -2.08896146E-01
      E(1,3) =  -6.13941874E-02

      E(2,0) =   4.61662109E-03 
      E(2,1) =  -1.35650709E-03
      E(2,2) =   1.60491325E-04

      E(3,0) =  -4.03643333E-06
      E(3,1) =   5.94046919E-07

      E(4,0) =   1.31604037E-09

      PRINT *, 'Enter T in C '
      READ  *, TC

      T  = TC + 273.15

      PRINT *, 'Enter density '
      READ  *, RHO

C      PRINT *, 'Enter (d alpha / d T) '
C      READ  *, DADT

      CALL WAPVT (TC, RHO, PRESS, BETA, G, ALPHA)

      ETA    = DIECAL (E, T, RHO)

      DETADP = DEDP (E, T, RHO, BETA)
      DETADT = DEDT (E, T, RHO, ALPHA)
      DETDT2 = D2EDT2 (E, T, RHO, ALPHA, DADT)

      Q = DETADP * 1 / ETA**2
      Y = DETADT * 1 / ETA**2
      X = (DETDT2 - 2 / ETA * DETADT**2) * 1 / ETA**2

      Q = Q *   1.0D+06
      Y = Y * (-1.0D+05)
      X = X * (-1.0D+07)

      PRINT *, 'Dielectric Constant = ', ETA
      PRINT *, 'Q = ', Q
      PRINT *, 'Y = ', Y
      PRINT *, 'X = ', X
      PRINT '(A, F10.5)', ' Pressure = ', PRESS
      PRINT *, 'Alpha = ', ALPHA
      PRINT *, 'Beta = ' , BETA

      END





      REAL * 8 FUNCTION DIECAL (E, T, RHO)
C =====================================================================|

      INTEGER I, J
      REAL * 8 DD, E(0:4,0:4), RHO, T

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

      DD = 0.0 

      DO 20 I = 0, 4
        DO 10 J = 0, 4-I
          DD = DD + E(I,J) * T**I * RHO**J
   10   CONTINUE
   20 CONTINUE

      DIECAL = DD

      RETURN 
      END





      REAL * 8 FUNCTION DEDP (E, T, RHO, BETA)
C =====================================================================|

      INTEGER I, J
      REAL * 8 BETA, DD, E(0:4,0:4), RHO, T
C =====================================================================|

      DD = 0.0 

      DO 20 I = 0, 4
        DO 10 J = 0, 4-I
          DD = DD + J * (E(I,J) * T**I * RHO**J)
   10   CONTINUE
   20 CONTINUE

      DEDP = DD * BETA

      RETURN 
      END






      REAL * 8 FUNCTION DEDT (E, T, RHO, ALPHA)
C =====================================================================|

      INTEGER I, J
      REAL * 8 ALPHA, DD, E(0:4,0:4), RHO, T
C =====================================================================|

      DD = 0.0 

      DO 20 I = 0, 4
        DO 10 J = 0, 4-I
          DD = DD + (E(I,J) * T**I * RHO**J) * (I/T - J * ALPHA)
   10   CONTINUE
   20 CONTINUE

      DEDT = DD

      RETURN 
      END






      REAL * 8 FUNCTION D2EDT2 (E, T, RHO, ALPHA, DADT)
C =====================================================================|

      INTEGER I, J
      REAL * 8 ALPHA, DADT, DA, DB, DC, DD, E(0:4,0:4), RHO, T
C =====================================================================|

      DA = 0.0
      DB = 0.0
      DC = 0.0
      DD = 0.0 

      DO 20 I = 0, 4
        DO 10 J = 0, 4-I
          DA = E(I,J) * RHO**J
          DB = J * ALPHA * (J * ALPHA * T**I - 2 * I * T**(I-1))
          DC = I * (I - 1) * T**(I-2) - J * T**I * DADT
          DD = DD + DA * (DB + DC)
   10   CONTINUE
   20 CONTINUE

      D2EDT2 = DD

      RETURN 
      END




      SUBROUTINE WAPVT (TC, RHO, PSAT, BETA, G, ALPHA)
C =====================================================================|
C     Routine to calculate the properties of pure water                |
C     from the steam tables of keenan et al (1969)                     |
C     See Helgeson and Kirkham (1974)                                  |
C ---------------------------------------------------------------------|
C     A         DBL   Fit coefficents A(i,j)                           |
C     AF        DBL   Helmholtz free energy as a function of rho and T |
C     AF0       DBL   First term of Helmholtz function                 |
C     ALPHA     DBL   Coefficient of thermal expansion                 |
C     BETA      DBL   Coefficient of compressibility                   |
C     C         DBL   Fit coefficients C(i)                            |
C     DIDX      DBL   First derivative of SI with respect to density   |
C     DIDX2     DBL   Second derivative of SI with respect to density  |
C     DKDX      DBL   First derivative of SK with respect to density   |
C     DKDX2     DBL   Second derivative of SK with respect to density  |
C     DPDT      DBL   Derivative of pressure with respect to temp.     |
C     DPDTR     DBL   Derivative of Q by tau with constant rho         |
C     DPDX      DBL   Derivative of pressure by density at constant T  |
C     DQDTR     DBL   Intermediate storage variable                    |
C     DQDX      DBL   First derivative of Q by density at constant T   |
C     DQDX2     DBL   Second derivative of Q by density at constant T  |
C     DTDTR     DBL   First derivative of TJ with respect to Tau       |
C     DTRDT     DBL   Derivative of TR with respect to temperature     |
C     E         DBL   Fit constant                                     |
C     FK        DBL   Intermediate storage variable                    |
C     G         DBL   Gibbs Free Energy in joules per mole             |
C     I         INT   Loop variable                                    |
C     J         INT   Loop variable                                    |
C     MW        DBL   Molecular weight of water                        |
C     PRESS     DBL   Pressure (bar)                                   |
C     Q         DBL   Intermediate storage variable                    |
C     RC        DBL   R in j/ g K                                      |
C     RP        DBL   R in bar cm(3)/g K                               |
C     RHO       DBL   Density (g/cm3)                                  |
C     RHOAJ     DBL   Density(a,j)                                     |
C     RHODIF    DBL   Difference between density and density (a,j)     |
C     SC        DBL   Intermediate storage variable                    |
C     SI        DBL   Sum of A(i,j) * TI(i), i=1,8  in Q               |
C     SK        DBL   Sum of A(i,j), i=9,10  in Q                      |
C     T         DBL   Temperature in deg K                             |
C     TC        DBL   Temperature in deg C                             |
C     TAUCRT    DBL   1000 / critical temperature of water in kelvin   |
C     TF        DBL   Delta Tau - Tau (critical)                       |
C     TI        DBL   Array of RhoDifs to the (i-1)th power            |
C     TJ        DBL   Intermediate storage variable                    |
C     TAUAJ     DBL   Tau(a,j)                                         |
C     TAU       DBL   1000 / current temperature of water in kelvin    |
C     TAUDIF    DBL   Difference between tau and tau (a,j)             |
C ---------------------------------------------------------------------|

      DOUBLE PRECISION MW
      PARAMETER (MW = 18.0153D+00)

      INTEGER I, J

      DOUBLE PRECISION A(10,7), AF, AF0, ALPHA, BETA, C(8)
      DOUBLE PRECISION DIDX, DIDX2, DKDX, DKDX2, DPDT, DPDTR, DPDX
      DOUBLE PRECISION DQ, DQQ, DQDRHO, DQDTR, DQDX, DQDX2
      DOUBLE PRECISION DADT
      DOUBLE PRECISION DTDTR, DTRDT
      DOUBLE PRECISION E, F(8), FK, G, PCRT, PRESS, PSAT
      DOUBLE PRECISION Q, Q1, QQ, QQQ
      DOUBLE PRECISION RC, RP, SC, SI, SK, SUMF
      DOUBLE PRECISION RHO, RHOAJ, RHODIF
      DOUBLE PRECISION T, TC, TCRT, TF, TI(8), TJ
      DOUBLE PRECISION TAU, TAUAJ, TAUCRT, TAUDIF

      INTRINSIC DEXP, DLOG

      DATA TAUCRT, E, RP/ 1.544912D+0, 4.8D+0, 4.6151D+0/
      DATA TCRT, PCRT /374.136, 220.88/
      DATA A / 2.9492937D+1, -1.3213917D+2, 2.7464632D+2, -3.6093828D+2,
     .         3.4218431D+2, -2.4450042D+2, 1.5518535D+2,  5.9728487D+0,
     .        -4.1030848D+2, -4.1605860D+2,-5.1985860D+0,  7.7779182D+0,
     .        -3.3301902D+1, -1.6254622D+1,-1.7731074D+2,  1.2748742D+2,
     .         1.3746153D+2,  1.5597836D+2, 3.3731180D+2, -2.0988866D+2,
     .         6.8335354D+0, -2.6149751D+1, 6.5326396D+1, -2.6181978D+1,
     .         0.0000000D+0,  0.0000000D+0, 0.0000000D+0,  0.0000000D+0,
     .        -1.3746618D+2, -7.3396848D+2,-1.5641040D-1, -7.2546108D-1,
     .        -9.2734289D+0,  4.3125840D+0, 0.0000000D+0,  0.0000000D+0,
     .         0.0000000D+0,  0.0000000D+0, 6.7874983D+0,  1.0401717D+1,
     .        -6.3972405D+0,  2.6409282D+1,-4.7740374D+1,  5.6323130D+1,
     .         0.0000000D+0,  0.0000000D+0, 0.0000000D+0,  0.0000000D+0,
     .         1.3687317D+2,  6.4581880D+2,-3.9661401D+0,  1.5453061D+1,
     .        -2.9142470D+1,  2.9568796D+1, 0.0000000D+0,  0.0000000D+0,
     .         0.0000000D+0,  0.0000000D+0, 7.9847970D+1,  3.9917570D+2,
     .        -6.9048554D-1,  2.7407416D+0,-5.1028070D+0,  3.9636085D+0,
     .         0.0000000D+0,  0.0000000D+0, 0.0000000D+0,  0.0000000D+0,
     .         1.3041253D+1,  7.1531353D+1/
      DATA C / 1.8570650D+3,  3.2291200D+3,-4.1946500D+2,  3.6664900D+1,
     .        -2.0551600D+1,  4.8523300D+0, 4.6000000D+1, -1.0112490D+3/
      DATA F /-7.4192420D+2, -2.9721000D+1,-1.1552860D+1, -0.8685635D+0,
     .         0.1094098D+0,  0.4399930D+0, 0.2520658D+0,  0.0521868D+4/
      DATA TI / 8*0.0D+0/

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

      RC = RP / 10.0D+0
      T  = TC + 273.15D+00
      TAU = 1.0D+3 / T
      TF = TAU - TAUCRT
      FK = DEXP(-E * RHO)

C ---------------------------------------------------------------------|
C     Calculate psi sub 0                                              |
C ---------------------------------------------------------------------|

      SC = 0.0D+0

      DO 10 I = 1, 6
        SC = SC + C(I) / TAU**(I-1)
   10 CONTINUE

      AF0   = SC + C(7) * DLOG(T) + C(8) * DLOG(T / TAU)

C ---------------------------------------------------------------------|
C     Calculate Q term by term                                         |
C ---------------------------------------------------------------------|

      Q     = 0.0D+0
      DQDX  = 0.0D+0
      DQDX2 = 0.0D+0
      DPDTR = 0.0D+0
      DQDTR = 0.0D+0

      DO 50 J = 1, 7
        IF (J .GT. 1) THEN
          TAUAJ = 2.5D+0
          RHOAJ = 1.0D+0
        ELSE
          TAUAJ = TAUCRT
          RHOAJ = 0.634D+0
        END IF

   30   CONTINUE
        RHODIF = RHO - RHOAJ
        TAUDIF = TAU - TAUAJ
        TI(1) = 1.0D+0

        DO 40 I = 2, 8
          TI(I) = TI(I-1) * RHODIF
   40   CONTINUE

        SI = A(1,J) + A(2,J) * TI(2) + A(3,J) * TI(3) + A(4,J) * TI(4) +
     .                A(5,J) * TI(5) + A(6,J) * TI(6) + A(7,J) * TI(7) +
     .                A(8,J) * TI(8)
        DIDX = A(2,J) + 2.0*A(3,J) * TI(2) + 3.0*A(4,J) * TI(3) +
     .                  4.0*A(5,J) * TI(4) + 5.0*A(6,J) * TI(5) +
     .                  6.0*A(7,J) * TI(6) + 7.0*A(8,J) * TI(7)
        DIDX2 = 2.0*A(3,J) + 6.0*A(4,J) * TI(2) + 12.0*A(5,J) * TI(3) +
     .                      20.0*A(6,J) * TI(4) + 30.0*A(7,J) * TI(5) +
     .                      42.0*A(8,J) * TI(6)

        SK = FK * (A(9,J) + A(10,J) * RHO)
        DKDX = FK * A(10,J) - E * SK
        DKDX2 = -E * FK * A(10,J) - E * DKDX

        TJ = TAUDIF**(J-2)
        DTDTR = (J-2) * TAUDIF**(J-3)

        Q = Q + TJ * (SI + SK)
        DQDX  = DQDX  + TJ * (DIDX  + DKDX)
        DQDX2 = DQDX2 + TJ * (DIDX2 + DKDX2)
        DPDTR = DPDTR + DTDTR * (SI + SK)
        DQDTR = DQDTR + DTDTR * (DIDX + DKDX)
   50 CONTINUE


      Q1 = 0.0D+0
      DQDRHO = 0.0D+0

      DO 52 J = 1, 7
        IF (J .EQ. 1) THEN
          TAUAJ = TAUCRT
          RHOAJ = 0.634D+0
        ELSE
          TAUAJ = 2.5D+0
          RHOAJ = 1.0D+0
        END IF

        TAUDIF = TAU - TAUAJ
        RHODIF = RHO - RHOAJ

        QQ = 0.0D+0
        DQ = 0.0D+0
        DO 54 I = 1, 8
          QQ = QQ + A(I,J) * RHODIF**(I-1)
          DQ = DQ + (I-1) * A(I,J) * RHODIF**(I-2)
   54   CONTINUE

        QQQ = 0.0D+0
        DQQ = 0.0D+0

        DO 56 I = 9, 10
          QQQ = QQQ + DEXP(-E * RHO) * A(I,J) * RHO**(I-9)
          DQQ = DQQ - E * DEXP(-E * RHO) * A(I,J) * RHO**(I-9)
   56   CONTINUE

        QQ = QQ + QQQ
        DQQ = DQQ + A(10,J) * DEXP(-E * RHO)
        DQ = DQ + DQQ

        Q1 = Q1 + TAUDIF**(J-2) * QQ
        DQDRHO = DQDRHO +  TAUDIF**(J-2) * DQ
   52 CONTINUE

      Q1     = Q1 * (TAU - TAUCRT)
      DQDRHO = DQDRHO * (TAU - TAUCRT)

      DPDTR = Q + TF * DPDTR
      DQDTR = DQDX + TF * DQDTR
      Q = TF * Q
      DQDX  = TF * DQDX
      DQDX2 = TF * DQDX2
      DTRDT = -1.0D+3 / T**2
      DPDT  = DTRDT * (DPDTR + DQDTR * RHO)

      AF = AF0 + RC * T * (DLOG(RHO) + RHO * Q)

      PRESS = RHO * RP * T * (1.0D+0 + RHO * Q1 + RHO**2 * DQDRHO)

      SUMF = 0.0D+0

      DO 60 I = 1, 8
        SUMF = SUMF + F(I) * (0.65 - 0.01 * TC)**(I-1)
   60 CONTINUE

      SUMF = SUMF * (TCRT - TC) * 1.0D-05 * TAU
      PSAT = PCRT * DEXP(SUMF) 

      DPDX = PSAT / RHO + RHO * RP * T * 
     &      (Q + 3.0 * RHO * DQDX + RHO**2 * DQDX2)

      BETA = 1.0 / (RHO * DPDX)

      ALPHA = BETA * (PSAT / T + T * RP * DPDT * RHO**2)

      DADT = ALPHA / BETA * 

      G  = (AF + (PSAT / RHO) * (RC / RP)) * MW

      RETURN
      END






