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