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