C ---------------------------------------------------------- C ** Thermophysical Properties of Methane ** C ** Literature Reference : ** C ** Daniel G. Friend, James F. Ely, and Hepburn Ingham ** C ** Thermophysical Properties of Methane, J.Phys.Chem. ** C ** Ref.Data, Vol.18, No.2, 1989, pp583--638 ** C ** Programmed by Cheng Wenlong ** C ** Dept. of Thermal Science and Engergy Engineering ** C ** University of Science and Technology of China ** C ---------------------------------------------------------- C==================================================================== C May 2, 2001: modified by Ryo Akasaka for ver.12.1 C The following 18 functions were added but not implemented. C AKPD(8A), AKPDD(8B), AKTD(8C), AKTDD(8D), CVPD(7A), CVTD(7B), C EPSPD(2A), EPSPDD(2B), EPSTD(2C), EPSTDD(2D), GAMPD(9A), C GAMTD(9B), TPH2(6H), TPS2(6S), WPD(8E), WPDD(8F), WTD(8G), C WTDD(8H) C C------------------------------------------------- F8A = AKPD REAL FUNCTION AKPD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='AKPD' IF (MESS.NE.0) CALL S99CH4(FUN) AKPD=-1.0E+30 RETURN END C------------------------------------------------- F8B = AKPDD REAL FUNCTION AKPDD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='AKPDD' IF (MESS.NE.0) CALL S99CH4(FUN) AKPDD=-1.0E+30 RETURN END C------------------------------------------------- F8C = AKTD REAL FUNCTION AKTD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='AKTD' IF (MESS.NE.0) CALL S99CH4(FUN) AKTD=-1.0E+30 RETURN END C------------------------------------------------- F8D = AKTDD REAL FUNCTION AKTDD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='AKTDD' IF (MESS.NE.0) CALL S99CH4(FUN) AKTDD=-1.0E+30 RETURN END C------------------------------------------------- F7A = CVPD REAL FUNCTION CVPD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='CVPD' IF (MESS.NE.0) CALL S99CH4(FUN) CVPD=-1.0E+30 RETURN END C------------------------------------------------- F7B = CVTD REAL FUNCTION CVTD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='CVTD' IF (MESS.NE.0) CALL S99CH4(FUN) CVTD=-1.0E+30 RETURN END C------------------------------------------------- F2A = EPSPD REAL FUNCTION EPSPD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='EPSPD' IF (MESS.NE.0) CALL S99CH4(FUN) EPSPD=-1.0E+30 RETURN END C------------------------------------------------- F2B = EPSPDD REAL FUNCTION EPSPDD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='EPSPDD' IF (MESS.NE.0) CALL S99CH4(FUN) EPSPDD=-1.0E+30 RETURN END C------------------------------------------------- F2C = EPSTD REAL FUNCTION EPSTD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='EPSTD' IF (MESS.NE.0) CALL S99CH4(FUN) EPSTD=-1.0E+30 RETURN END C------------------------------------------------- F2D = EPSTDD REAL FUNCTION EPSTDD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='EPSTDD' IF (MESS.NE.0) CALL S99CH4(FUN) EPSTDD=-1.0E+30 RETURN END C------------------------------------------------- F9A = GAMPD REAL FUNCTION GAMPD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='GAMPD' IF (MESS.NE.0) CALL S99CH4(FUN) GAMPD=-1.0E+30 RETURN END C------------------------------------------------- F9B = GAMTD REAL FUNCTION GAMTD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='GAMTD' IF (MESS.NE.0) CALL S99CH4(FUN) GAMTD=-1.0E+30 RETURN END C------------------------------------------------- F6H = TPH2 REAL FUNCTION TPH2(P,H) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='TPH2' IF (MESS.NE.0) CALL S99CH4(FUN) TPH2=-1.0E+30 RETURN END C------------------------------------------------- F6S = TPS2 REAL FUNCTION TPS2(P,S) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='TPS2' IF (MESS.NE.0) CALL S99CH4(FUN) TPS2=-1.0E+30 RETURN END C------------------------------------------------- F8E = WPD REAL FUNCTION WPD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='WPD' IF (MESS.NE.0) CALL S99CH4(FUN) WPD=-1.0E+30 RETURN END C------------------------------------------------- F8F = WPDD REAL FUNCTION WPDD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='WPDD' IF (MESS.NE.0) CALL S99CH4(FUN) WPDD=-1.0E+30 RETURN END C------------------------------------------------- F8G = WTD REAL FUNCTION WTD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='WTD' IF (MESS.NE.0) CALL S99CH4(FUN) WTD=-1.0E+30 RETURN END C------------------------------------------------- F8H = WTDD REAL FUNCTION WTDD(T) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='WTDD' IF (MESS.NE.0) CALL S99CH4(FUN) WTDD=-1.0E+30 RETURN END C C==================================================================== C C ------------------------------------------------------- C ** F8 = ALMPT(P,T) [W/mK] ** C ** Thermal Conductivity of Methane ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [K],[C] ** C ------------------------------------------------------- REAL FUNCTION ALMPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F8CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'ALMPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F8CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF ALMPT=FF RETURN END C ------------------------------------------------------- C ** F8CH4(P,T) , ALMPT(P,T) : Thermal Conductivity ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F8CH4(P,T) IMPLICIT DOUBLE PRECISION(A-H,L-Z) TM=T+273.15D00 TP=F69CH4(P) IF((T.GT.TP.AND.TM.GE.91).AND.TM.LE.623.AND.P.LE.620) GOTO 100 GOTO 900 100 CONTINUE RO=F51CH4(P,T) LO=1.0D0/RO IF(RO.EQ.-1.0E+20) GO TO 900 F8CH4=THECON(LO,TM)*0.001D0 RETURN 900 F8CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** Thermal Conductivity [mW/(mK] ** C ** INPUT : ** C ** LO : Density [kg/m^3] ** C ** T : Temperature [K] ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION THECON(LO,T) DOUBLE PRECISION LO,T,LABDA0,LABDAEX,LABDACR THECON=LABDA0(LO,T)+LABDAEX(LO,T)+LABDACR(LO,T) RETURN END C ----------------------------------------------------------- C ** The Modified Eucken Model ** C ** for the thermal conductivity of the dilute gas ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION LABDA0(LO,T) DOUBLE PRECISION T,LO DOUBLE PRECISION FIA0,FINT,F1,F2,TI,K,E,I4CH4 DATA F1/1.458850D0/, F2/-0.4377162D0/ DATA K/1.380658D-23/,E/240.234492D-23/ TI=K*T/E FINT=F1+F2/TI LABDA0=0.51826*FIA0(T)*(3.75D0-FINT*(I4CH4(LO,T)+1.5D0)) RETURN END C ------------------------------------------------------- C ** Excess Thermal Conductivity ** C ** LABDAEX(LO,T) [mW/(mK)] ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION LABDAEX(LO,T) DOUBLE PRECISION LO,T,TM DOUBLE PRECISION F54CH4,LV,LC,TC,DLTAV DOUBLE PRECISION R,S,J DIMENSION R(7),S(7),J(7) DATA TC/190.551D0/,LC/162.66/ DATA(R(I),I=1,7)/1.0D0, 3.0D0, 4.0D0, 4.0D0, 5.0D0, 5.0D0, 2.0D0/ DATA(S(I),I=1,7)/0.0D0, 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.0D0, 0.0D0/ DATA(J(I),I=1,7)/2.4149207D0, 0.55166331D0, -0.52837734D0, & 0.073809553D0, 0.24465507D0, -0.047613626D0, 1.5554612D0/ TM=T-273.15D0 TAO=TC/T DLTA=LO/LC IF(T.LT.TC.AND.LO.LT.LC) THEN LV=1.D0/F54CH4(TM) DLTAV=LV/LC ELSE DLTAV=1.D0 ENDIF LABDAEX=J(7)*DLTA*DLTA/DLTAV DO 111 I=1,6 LABDAEX=LABDAEX+J(I)*DLTA**R(I)*TAO**S(I) 111 CONTINUE LABDAEX=6.29638*LABDAEX RETURN END C ------------------------------------------------------- C ** Critical Enhancement Thermal Conductivity ** C ** LABDACR(LO,T) [mW/(mK)] ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION LABDACR(LO,T) IMPLICIT DOUBLE PRECISION(A-Z) DATA TC/190.551D0/,LC/162.66D0/ DATA FT/2.646D0/,FP/2.678D0/,FA/-0.637D0/, & BLTA/0.355D0/,A/3.352D0/,B/0.732D0/, & E/0.287D0/,R/0.535D0/,Q/0.1133D0/,S/-6.098D0/,W/-1.401D0/ TAO=TC/T DLTA=LO/LC TAOS=1.D0-1.D0/TAO DLTAS=1-DLTA F=DEXP(-(FT*(DABS(TAOS))**0.5D0+FP*DLTAS*DLTAS+FA*DLTAS)) IF(DABS(DLTAS).LT.1.0D-7.AND.DABS(TAOS).LT.0.03D0) THEN XT=0.0801D0*(DABS(TAOS))**(-1.190D0) ELSE IF(DABS(TAOS).LT.0.03D0.AND.DABS(DLTAS).LT.0.25D0) THEN TP=-(DABS(DLTAS))**(1.D0/BLTA)/S IF(TAOS.LT.TP) THEN SITA=1.D0+E*(1.D0+S*TAOS*DABS(DLTAS)**(-1.D0/BLTA)) & **(2.D0*BLTA) ELSE SITA=1.D0 ENDIF OMIGA=W*TAOS*(DABS(DLTAS))**(-1.0/BLTA) XT=Q*(DABS(DLTAS))**(-A)*SITA**B/(SITA+OMIGA*(SITA+R)) ELSE XT=0.28631*DLTA*TAO/(1.D0+2.D0*R1CH4(LO,T)+R3CH4(LO,T)) ENDIF LABDACR=91.855/(COEVIS(LO,T)*TAO*TAO)*F* & (1.D0+R1CH4(LO,T)-R5CH4(LO,T))**2.D0*XT**0.4681D0 RETURN END C ------------------------------------------------------- C ** F14 = AMUTD(T) [Pas] ** C ** Coefficient of Viscosity of Saturated Liquid ** C ------------------------------------------------------- REAL FUNCTION AMUTD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F14CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'AMUTD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F14CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF AMUTD=FF RETURN END C ------------------------------------------------------- C ** F15 = AMUTDD(T) [Pas] ** C ** Coefficient of Viscosity of Saturated Vapor ** C ------------------------------------------------------- REAL FUNCTION AMUTDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F15CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'AMUTDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F15CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF AMUTDD=FF RETURN END C ------------------------------------------------------- C ** F14CH4(T) , AMUTD(T) : ** C ** COEFFICIENT OF VISCOSITY OF SATURATED LIQUID ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F14CH4(T) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA TC/190.551D00/,TL/90.6854D00/ TM=T+273.15D00 IF(TM.LT.TL) GO TO 600 IF(TM.GT.TC) GO TO 600 F1=0. RO=F53CH4(T) LO=1.0D0/RO IF(RO.EQ.-1.0E+20) GO TO 600 F14CH4=1.0D-06*COEVIS(LO,TM) RETURN 600 F14CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F15CH4(T) , AMUTDD(T) : ** C ** COEFFICIENT OF VISCOSITY OF SATURATED VAPOR ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F15CH4(T) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA TC/190.551D0/,TL/90.6854D0/ TM=T+273.15D00 IF(TM.LT.TL) GO TO 600 IF(TM.GT.TC) GO TO 600 F1=0. RO=F54CH4(T) LO=1.0D0/RO IF(RO.EQ.-1.0E+20) GO TO 600 F15CH4=1.0D-06*COEVIS(LO,TM) RETURN 600 F15CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F13 = AMUPT(P,T) [Pas] ** C ** Coefficient of Viscosity ** C ----------------------------------------------------------- REAL FUNCTION AMUPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F13CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'AMUPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F13CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF AMUPT=FF RETURN END C ------------------------------------------------------------ C ** F13CH4(P,T) , AMUPT(P,T) : COEFFICIENT OF VISCOSITY ** C ------------------------------------------------------------ DOUBLE PRECISION FUNCTION F13CH4(P,T) IMPLICIT DOUBLE PRECISION(A-H,L-Z) TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 TM=T+273.15D00 IF((T.GT.TP.AND.TM.GE.91).AND.TM.LE.623.AND.P.LE.1000) GOTO 100 GO TO 900 100 CONTINUE RO=F51CH4(P,T) LO=1.0D0/RO IF(RO.EQ.-1.0E+20) GO TO 900 F13CH4=1.0D-6*COEVIS(LO,TM) RETURN 900 F13CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F12 = AMUPDD(P) [Pas] ** C ** Coefficient of Viscosity of Saturated Vapor ** C ------------------------------------------------------- REAL FUNCTION AMUPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F12CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'AMUPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F12CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF AMUPDD=FF RETURN END C ------------------------------------------------------- C ** F11CH4(P) , AMUPDD ** C ** COEFFICIENT OF VISCOSITY OF SATURATED VAPOR ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F12CH4(P) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF(P.LT.0.198790D00.OR.P.GT.45.9952D00) GO TO 600 T=F40CH4(P) IF(T.EQ.-1.0E+20) GO TO 600 TM=T+273.15D00 F1=0. RO=F54CH4(T) LO=1.0D0/RO IF(RO.EQ.-1.0E+20) GO TO 600 F12CH4=1.0D-06*COEVIS(LO,TM) RETURN 600 F12CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F11 = AMUPD ** C ** Coefficient of Viscosity of Saturated Liquid ** C ------------------------------------------------------- REAL FUNCTION AMUPD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F11CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'AMUPD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F11CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF AMUPD=FF RETURN END C ------------------------------------------------------- C ** F11CH4(P) , AMUPD ** C ** COEFFICIENT OF VISCOSITY OF SATURATED LIQUID ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F11CH4(P) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF(P.LT.0.198790D00.OR.P.GT.45.9921D00) GO TO 600 T=F40CH4(P) IF(T.EQ.-1.0E+20) GO TO 600 TM=T+273.15D00 F1=0. RO=F53CH4(T) LO=1.0D0/RO IF(RO.EQ.-1.0E+20) GO TO 600 F11CH4=COEVIS(LO,TM)*1.0D-6 RETURN 600 F11CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** Coefficient of Viscosoty [uPas] ** C ** INPUT : ** C ** LO : Density [kg/m^3] ** C ** T : Temperature [K] ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION COEVIS(LO,T) DOUBLE PRECISION LO,T,FIA0,FAIEX COEVIS=FIA0(T)+FAIEX(LO,T) RETURN END C ----------------------------------------------------------- C ** The Chapman-Enskog theory for ** C ** the dilute gas viscosity ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION FIA0(T) DOUBLE PRECISION T DOUBLE PRECISION TI,C,OMIGA,K,E DIMENSION C(9) DATA K/1.380658D-23/,E/240.234492D-23/ DATA (C(J),J=1,9) /-3.0328138281D0,16.918880086D0,-37.189364917D0, & 41.288861858D0,-24.615921140D0,8.9488430959D0,-1.8739245042D0, & 0.20966101390D0,-9.6570437074D-3/ TI=K*T/E OMIGA=0.0D0 DO 111 J=1,9 OMIGA=OMIGA+C(J)*TI**((J-1.D0)/3.D0-1.D0) 111 CONTINUE OMIGA=1.D0/OMIGA FIA0=10.50D0*DSQRT(TI)/OMIGA RETURN END C --------------------------------------------------- C ** Excess Viscosity FAIEX(LO,T) [uPas] ** C --------------------------------------------------- DOUBLE PRECISION FUNCTION FAIEX(LO,T) DOUBLE PRECISION LO,T DOUBLE PRECISION LC,TC,R,S,G,TAO,DLTA,FF1,FF2 DIMENSION R(11),S(11),G(11) DATA LC/162.66D0/,TC/190.551D0/ DATA (R(J),J=1,11) /1.D0, 1.D0, 2.D0, 2.D0, 2.D0, 3.D0, & 3.D0, 4.D0, 4.D0, 1.D0, 1.D0/ DATA (S(J),J=1,11) /0.D0, 1.D0, 0.D0, 1.D0, 1.5D0, 0.D0, & 2.D0, 0.D0, 1.D0, 0.D0, 1.D0/ DATA (G(J),J=1,11) /0.41250137D0, -0.14390912D0, 0.10366993D0, & 0.40287464D0, -0.24903524D0, -0.12953131D0, 0.06575776D0, & 0.02566628D0, -0.03716526D0, -0.38798341D0, 0.03533815D0/ TAO=TC/T DLTA=LO/LC FF1=0.D0 FF2=0.D0 DO 111 J=1,9 FF1=FF1+G(J)*DLTA**R(J)*TAO**S(J) 111 CONTINUE DO 222 J=10,11 FF2=FF2+G(J)*DLTA**R(J)*TAO**S(J) 222 CONTINUE FAIEX=12.149*FF1/(1.D0+FF2) RETURN END C ----------------------------------------------------------- C ** F95 = GAMPT(P,T) ** C ** Ratio of Specific Heats of Saturated Vapor [-] ** C ** INPUT : ** C ** T : Temperature [C],[K] ** C ** P : Pressure [Pa],[bar] ** C ----------------------------------------------------------- REAL FUNCTION GAMPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F95CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'GAMPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F95CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF GAMPT=FF RETURN END C ---------------------------------------------------------- C ** F95CH4(P,T) , GAMPT(P,T) : RATIO OF SPECIFIC HEATS ** C ---------------------------------------------------------- DOUBLE PRECISION FUNCTION F95CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 CP=F18CH4(P,T) IF(CP.EQ.-1.0E+20) GO TO 900 CV=F77CH4(P,T) IF(CV.EQ.-1.0E+20) GO TO 900 F95CH4=CP/CV RETURN 900 F95CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F97 = GAMTDD(T) ** C ** Ratio of Specific Heats of Saturated Vapor [-] ** C ** INPUT : ** C ** T : Temperature [C],[K] ** C ----------------------------------------------------------- REAL FUNCTION GAMTDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F97CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'GAMTDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F97CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF GAMTDD=FF RETURN END C ------------------------------------------------------------- C ** GAMTDD(T) RATIO OF SPECIFIC HEATS OF SATURATED VAPOR ** C ------------------------------------------------------------- DOUBLE PRECISION FUNCTION F97CH4(T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) CPDD=F20CH4(T) IF(CPDD.EQ.-1.0E+20) GO TO 900 CVDD=F78CH4(T) IF(CVDD.EQ.-1.0E+20) GO TO 900 F97CH4=CPDD/CVDD RETURN 900 F97CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F96 = GAMPDD(P) ** C ** Ratio of Specific Heats of Saturated Vapor [-] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ----------------------------------------------------------- REAL FUNCTION GAMPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F96CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'GAMPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F96CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF GAMPDD=FF RETURN END C ------------------------------------------------------- C ** F96CH4(P) , GAMPDD : RATIO OF SPECIFIC HEATS ** C ** OF SATURATED ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F96CH4(P) IMPLICIT DOUBLE PRECISION (A-H,L-Z) IF(P.LT.0.11696D00.OR.P.GT.45.992D00) GO TO 900 CPDD=F17CH4(P) IF(CPDD.EQ.-1.0E+20) GO TO 900 CVDD=F76CH4(P) IF(CVDD.EQ.-1.0E+20) GO TO 900 F96CH4=CPDD/CVDD RETURN 900 F96CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F4= ALHP(P) : Latent Heat of Vaporization [J/kg] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ----------------------------------------------------------- REAL FUNCTION ALHP(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F4CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'ALHP'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F4CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF ALHP=FF RETURN END C ----------------------------------------------------------- C ** F5= ALHT(T) : Latent Heat of Vaporization [J/kg] ** C ** INPUT : ** C ** T : Temperature [C],[K] ** C ----------------------------------------------------------- REAL FUNCTION ALHT(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F5CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'ALHT'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F5CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF ALHT=FF RETURN END C ------------------------------------------------------- C ** F4CH4(P) , ALHP : LATENT HEAT OF VAPORIZATION ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F4CH4(P) IMPLICIT DOUBLE PRECISION(A-H,O-Z) IF(P.LT.0.11696D00.OR.P.GT.45.992D00) GO TO 900 HL=F23CH4(P) HV=F24CH4(P) IF(HL.EQ.-1.0E+20.OR.HV.EQ.-1.0E+20) GO TO 900 F4CH4=HV-HL RETURN 900 F4CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F5CH4(T) , ALHT : LATENT HEAT OF VAPORIZATION ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F5CH4(T) IMPLICIT DOUBLE PRECISION(A-H,L-Z) HL=F27CH4(T) HV=F28CH4(T) IF(HL.EQ.-1.0E+20.OR.HV.EQ.-1.0E+20) GO TO 900 F5CH4=HV-HL RETURN 900 F5CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F82= AKPT(P,T) : Adiabatic Exponent [-] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [C],[K] ** C ----------------------------------------------------------- REAL FUNCTION AKPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F82CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'AKPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F82CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF AKPT=FF RETURN END C --------------------------------------------------- C ** F82CH4(P,T) , AKPT(P,T) : ADIABATIC EXPONENT ** C --------------------------------------------------- DOUBLE PRECISION FUNCTION F82CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DATA R/518.264D00/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 LOO=F51CH4(P,T) TM=T+273.15D0 IF(LOO.EQ.-1.0E+20) GO TO 900 LO=1.D00/LOO CP=F18CH4(P,T) IF(CP.EQ.-1.0E+20) GO TO 900 CV=F77CH4(P,T) IF(CV.EQ.-1.0E+20) GO TO 900 DPDL=-R*TM*LO*LO*(1.0D0+2.0D0*R1CH4(LO,TM)+R3CH4(LO,TM)) VPT=F51CH4(P,T) IF(VPT.EQ.-1.0E+20) GO TO 900 F82CH4=-CP/CV*VPT*DPDL/(1.0D05*P) RETURN 900 F82CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F93= BVPT(P,T) ** C ** Pressure Coefficient [1/K] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [C],[K] ** C ----------------------------------------------------------- REAL FUNCTION BVPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F93CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'BVPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F93CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF BVPT=FF RETURN END C ----------------------------------------------------- C ** F93CH4(P,T) , BVPT(P,T) : PRESSURE COEFFICIENT ** C ----------------------------------------------------- DOUBLE PRECISION FUNCTION F93CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DATA R/518.264D0/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 VV=F51CH4(P,T) IF(VV.EQ.-1.0E+20) GO TO 900 LO=1.0D00/VV TM=T+273.15D0 DPDT=-LO*R*(-1-R1CH4(LO,TM)+R5CH4(LO,TM)) F93CH4=DPDT/(P*1.0D05) RETURN 900 F93CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F90= BSPT(P,T) ** C ** Adiabatic Compressibity [1/Pa] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [C],[K] ** C ----------------------------------------------------------- REAL FUNCTION BSPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F90CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'BSPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F90CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF BSPT=FF RETURN END C ----------------------------------------------------- C ** F90CH4(P,T) , BSPT : ADIABATIC COMPRESSIBILITY ** C ----------------------------------------------------- DOUBLE PRECISION FUNCTION F90CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 VV=F51CH4(P,T) IF(VV.EQ.-1.0E+20) GO TO 900 WC=F83CH4(P,T) IF(WC.EQ.-1.0E+20) GO TO 900 BBS=VV/(WC**2) F90CH4=BBS RETURN 900 F90CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F91= BTPT(P,T) ** C ** Isothermal Compressibility [1/Pa] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [C],[K] ** C ----------------------------------------------------------- REAL FUNCTION BTPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F91CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'BTPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F91CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF BTPT=FF RETURN END C ------------------------------------------------------- C ** F91CH4(P,T) , BTPT : ISOTHERMAL COMPRESSIBILITY ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F91CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 CP=F18CH4(P,T) IF(CP.EQ.-1.0E+20) GO TO 900 CV=F77CH4(P,T) IF(CV.EQ.-1.0E+20) GO TO 900 BS=F90CH4(P,T) IF(BS.EQ.-1.0E+20) GO TO 900 F91CH4=(CP/CV)*BS RETURN 900 F91CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F92= BPPT(P,T) ** C ** Volumetric Coefficient of Expansion [1/K] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [C],[K] ** C ----------------------------------------------------------- REAL FUNCTION BPPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F92CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'BPPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F92CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF BPPT=FF RETURN END C ----------------------------------------------------------- C ** F92CH4(P,T) , BPPT ** C ** VOLUMETRIC COEFFICIENT OF EXPANSION ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F92CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) C DATA RR/8.314510D00/,R/518.264D00/ DATA R/518.264D0/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 VV=F51CH4(P,T) TM=T+273.15D0 IF(VV.EQ.-1.0E+20) GO TO 900 LO=1.0D00/VV CP=F18CH4(P,T) IF(CP.EQ.-1.0E+20) GO TO 900 CV=F77CH4(P,T) IF(CV.EQ.-1.0E+20) GO TO 900 TM=T+273.15D00 DPDT=-LO*R*(-1-R1CH4(LO,TM)+R5CH4(LO,TM)) DVDT=(CP-CV)/(TM*DPDT) F92CH4=(1.0D00/VV)*DVDT RETURN 900 F92CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F94 = AJTPT(P,T) : Joule-Thomson Coefficient [K/Pa] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [K],[C] ** C ---------------------------------------------------------- REAL FUNCTION AJTPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F94CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'AJTPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F94CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF AJTPT=FF RETURN END C ----------------------------------------------------------- C ** F94CH4(P,T) , AJTPT(P,T) : JOULE-THOMSON COEFFICIENT ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F94CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 VV=F51CH4(P,T) IF(VV.EQ.-1.0E+20) GO TO 900 BP=F92CH4(P,T) IF(BP.EQ.-1.0E+20) GO TO 900 VT=BP*VV TM=T+273.15D00 CP=F18CH4(P,T) IF(CP.EQ.-1.0E+20) GO TO 900 F94CH4=(TM*VT-VV)/CP RETURN 900 F94CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F31 = SIGP(P) : Surface Tension [N/m] ** C ** INPUT : P Pressure [Pa],[bar] ** C ----------------------------------------------------------- REAL FUNCTION SIGP(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F31CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'SIGP'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F31CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF SIGP=FF RETURN END C ------------------------------------------------ C ** F31CH4(P) , SIGP : SURFACE TENSION *** C ------------------------------------------------ DOUBLE PRECISION FUNCTION F31CH4(P) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DATA TL/90.68D00/,TC/190.55D00/,T1/104.99D00/,SIGMA1/15.026D00/ T=F40CH4(P) IF(T.EQ.-1.0E+20) GO TO 900 TM=T+273.15D00 IF(TM.LT.TL.OR.TM.GT.TC) GO TO 900 F31CH4=(SIGMA1*((TC-TM)/(TC-T1))**1.3941D00)/1000.0 RETURN 900 F31CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F32 = SIGT(T) : Surface Tension [N/m] ** C ** INPUT : T Temperature [K],[C] ** C ----------------------------------------------------------- REAL FUNCTION SIGT(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F32CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'SIGT'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F32CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF SIGT=FF RETURN END C ------------------------------------------------ C ** F32CH4(T) , SIGT : SURFACE TENSION *** C ------------------------------------------------ DOUBLE PRECISION FUNCTION F32CH4(T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DATA TL/90.6854D00/,TC/190.551D00/,T1/104.99D00/,SIGMA1/15.026D00/ TM=T+273.15D00 IF(TM.LT.TL.OR.TM.GT.TC) GO TO 900 F32CH4=(SIGMA1*((TC-TM)/(TC-T1))**1.3941D00)/1000.0 RETURN 900 F32CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F98 = TPSEUP(P) : Pseudo Boiling Point [K],[C] ** C ------------------------------------------------------- FUNCTION TPSEUP(P) CHARACTER FUN*6 REAL PBAR,T0K,P,FF INTEGER KPA DOUBLE PRECISION F98CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'TPSEUP'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF PI=P*PBAR DBP=DBLE(PI) FF = F98CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 TPSEUP=FF+T0K RETURN END C ------------------------------------------------------------ C ** F98CH4(P) , TPSEUP : PSEUDO BOILING POINT AT P ** C ------------------------------------------------------------ DOUBLE PRECISION FUNCTION F98CH4(P) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION T(2),C(2),TL(3),TR(3),CL(3),CR(3) P1=F21CH4('P') PP=DABS((P-P1)/P1) T1=F21CH4('T') IF (PP.LT.1.0D-5) THEN F98CH4=T1 RETURN ENDIF IF (P.LT.P1.OR.P.GT.500.001D00) THEN F98CH4=-1.0E+20 RETURN ENDIF T2=0.85*T1 P2=F30CH4(T2) TM0=T1+(T1-T2)*(P-P1)/(P1-P2) TC=T1 T(1)=TM0 WRITE(6,901) TM0 901 FORMAT(1H ,' +++++ TM0 : ',G15.7) IF (P.GT.400.0) T(1)=-30.0 150 EPS=1.0D-7 DEL=T1*0.05 IREP=0 IREM=5000 KCONT=0 ICONT=0 C(1)=F18CH4(P,T(1)) T(2)=T(1)-DEL C(2)=F18CH4(P,T(2)) 1000 RINC=C(2)-C(1) IF(RINC.GT.0.3)THEN GOTO 1500 ELSE T(2)=T(1) C(2)=C(1) T(1)=T(1)-DEL C(1)=F18CH4(P,T(1)) GOTO 1000 ENDIF 1500 C(1)=-C(1) C(2)=-C(2) 2000 IREP=IREP+1 IF(IREP.GT.IREM) GO TO 8000 TT=T(2)+1.3*(T(2)-T(1)) CC=-F18CH4(P,TT) 3000 CONV=DABS((CC-C(2))/CC) IF(CONV.LT.EPS) THEN F98CH4=TT RETURN ENDIF 4000 DEC=C(1)-CC ICONT=ICONT+1 IF(DEC.GE.0.0) THEN T(1)=TT C(1)=CC ELSE IF (ICONT.GT.5) GO TO 6000 TT=(TT+T(2))*0.5 CC=-F18CH4(P,TT) GOTO 4000 ENDIF 5000 DEC=C(1)-C(2) IF(DEC.LT.0.0) THEN TT=T(1) CC=C(1) T(1)=T(2) C(1)=C(2) T(2)=TT C(2)=CC ENDIF GOTO 2000 6000 TA=T(1) TB=TT IF (TB.LT.TA) THEN TA=TT TB=T(1) ENDIF TC=TA+0.5*(TB-TA) CA=F18CH4(P,TA) CB=F18CH4(P,TB) 6050 KCONT=KCONT+1 IF (KCONT.GT.IREM) GO TO 8000 DELT=DABS((TA-TB)/TA) IF (DELT.LT.EPS) GO TO 7000 CC=F18CH4(P,TC) DTA=(TC-TA)*0.3 TL(1)=TA TL(2)=TC-DTA TL(3)=TC DTB=(TB-TC)*0.3 TR(1)=TC TR(2)=TC+DTB TR(3)=TB CL(1)=CA CL(3)=CC CR(1)=CC CR(3)=CB CL(2)=F18CH4(P,TL(2)) CR(2)=F18CH4(P,TR(2)) CMXL=CL(1) ML=1 DO 6120 I=2,3 IF(CL(I).GT.CMXL) THEN CMXL=CL(I) ML=I ENDIF 6120 CONTINUE CMXR=CR(1) MR=1 DO 6130 I=2,3 IF(CR(I).GT.CMXR) THEN CMXR=CR(I) MR=I ENDIF 6130 CONTINUE IF(CMXL.GT.CMXR) THEN IF(ML.EQ.1) THEN TA=TL(1)-DTA CA=F18CH4(P,TA) ELSE TA=TL(ML-1) CA=CL(ML-1) ENDIF IF(ML.EQ.3) THEN TB=TR(2) CB=CR(2) ELSE TB=TL(ML+1) CB=CL(ML+1) ENDIF TC=TL(ML) CC=CL(ML) ELSE IF(MR.EQ.1) THEN TA=TL(2) CA=CL(2) ELSE TA=TR(MR-1) CA=CR(MR-1) ENDIF IF(MR.EQ.3) THEN TB=TR(3)+DTB CB=F18CH4(P,TB) ELSE TB=TR(MR+1) CB=CR(MR+1) ENDIF 6135 TC=TR(MR) CC=CR(MR) ENDIF GO TO 6050 7000 F98CH4=TC RETURN 8000 F98CH4=-1.0E+10 RETURN END C ----------------------------------------------- C ** F64 = TPH(P,H) : TEMPERATURE [K],[C] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** H : Specific Enthalpy [J/kg] ** C ----------------------------------------------- REAL FUNCTION TPH(P,H) CHARACTER FUN*6,FLUID*8 REAL T0K,PBAR,P,H,FF INTEGER KPA DOUBLE PRECISION F64CH4,DBP,HH COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'TPH'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF PI=P*PBAR DBP=DBLE(PI) HH=DBLE(H) FF = F64CH4(DBP,HH) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,H 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND H =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 TPH=FF+T0K RETURN END C ----------------------------------------------- C ** F64CH4(P,H) , TPH : TEMPERATURE ** C ----------------------------------------------- DOUBLE PRECISION FUNCTION F64CH4(P,H) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA PC/45.95001D00/,ERR/1.D-07/ IF(P.GE.0.11719D00.AND.P.LE.PC) GO TO 50 40 IF(P.GE.0.11719D00.AND.P.LT.450.D00) GO TO 100 IF(P.GE.450.D00.AND.P.LE.10000.D00) GO TO 120 GO TO 900 50 TS=F40CH4(P) H0=F23CH4(P) H1=F24CH4(P) IF(H.GE.H0.AND.H.LE.H1) GO TO 800 GO TO 40 100 TL=F69CH4(P) TCC=620.D00-273.15D00 GO TO 150 120 TL=F69CH4(P) TCC=470.D00-273.15D00 150 H0=F25CH4(P,TL ) H1=F25CH4(P,TCC) IF(H.LT.H0.OR.H.GT.H1) GO TO 900 T0=TL G0=H0-H IF(H0.EQ.-1.0E+20) GO TO 900 T1=TCC T3=T1 G1=H1-H 200 T2=(T0+T1)/2.D00 T=T2 H2=F25CH4(P,T) G2=H2-H IF(DABS(T2-T3).LT.ERR) GO TO 350 T3=T2 IF((G0*G2).LT.0.) GO TO 300 G0=G2 T0=T2 GO TO 200 300 G1=G2 T1=T2 GO TO 200 350 F64CH4=T2 RETURN 800 F64CH4=TS RETURN 900 F64CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F65 = TPS(P,S) : Temperature ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** S : Specific Entropy [J/(kg.K)] ** C ------------------------------------------------------- REAL FUNCTION TPS(P,S) CHARACTER FUN*6,FLUID*8 REAL T0K,PBAR,P,S,FF INTEGER KPA DOUBLE PRECISION F65CH4,DBP,SS COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'TPS'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF PI=P*PBAR DBP=DBLE(PI) SS=DBLE(S) FF = F65CH4(DBP,SS) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,S 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND S =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 TPS=FF+T0K RETURN END C ----------------------------------------------- C ** F65CH4(P,S) : TPS(P,S) TEMPERATURE ** C ----------------------------------------------- DOUBLE PRECISION FUNCTION F65CH4(P,S) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA PC/45.992D00/,ERR/1.D-07/ IF(P.GE.0.11696D00.AND.P.LE.PC) GO TO 50 40 IF(P.GE.0.11696D00.AND.P.LT.450.D00) GO TO 100 IF(P.GE.450.D00.AND.P.LE.10000.D00) GO TO 120 GO TO 900 50 TS=F40CH4(P) S0=F33CH4(P) S1=F34CH4(P) IF(S.GE.S0.AND.S.LE.S1) GO TO 800 GO TO 40 100 TL=F69CH4(P) TCC=620.D00-273.15D00 GO TO 150 120 TL=F69CH4(P) TCC=470.D00-273.15D00 150 S0=F35CH4(P,TL) S1=F35CH4(P,TCC) IF(S.LT.S0.OR.S.GT.S1) GO TO 900 T0=TL G0=S0-S IF(S0.EQ.-1.0E+20) GO TO 900 T1=TCC T3=T1 G1=S1-S 200 T2=(T0+T1)/2.D00 T=T2 S2=F35CH4(P,T) G2=S2-S IF(DABS(T2-T3).LT.ERR) GO TO 350 T3=T2 IF((G0*G2).LT.0.) GO TO 300 G0=G2 T0=T2 GO TO 200 300 G1=G2 T1=T2 GO TO 200 350 F65CH4=T2 RETURN 800 F65CH4=TS RETURN 900 F65CH4=-1.0E+20 RETURN END C --------------------------------------------------- C ** F71 = HPS(P,S) : Specific Enthalpy [J/kg] ** C ** INPUT ** C ** P: PRESSURE [Pa],[bar] ** C ** S: Specific Entropy ** C --------------------------------------------------- REAL FUNCTION HPS(P,S) CHARACTER FUN*6,FLUID*8 REAL P,S,FF INTEGER KPA DOUBLE PRECISION F71CH4,DBP,SS COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'HPS'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) SS=DBLE(S) FF = F71CH4(DBP,SS) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,S 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND S =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF HPS=FF RETURN END C ------------------------------------------------------- C ** F71CH4(P,S) , HPS(P,S) : SPECIFIC ENTHALPY ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F71CH4(P,S) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA PC/45.992D00/ IF(P.GE.0.11696D00.AND.P.LE.PC) GO TO 50 40 IF(P.GE.0.11696D00.AND.P.LT.450.D00) GO TO 100 IF(P.GE.450.D00.AND.P.LE.10000.D00) GO TO 120 GO TO 900 50 S0=F33CH4(P) S1=F34CH4(P) IF(S.LT.S0.OR.S.GT.S1) GO TO 40 X=(S-S0)/(S1-S0) H0=F23CH4(P) H1=F24CH4(P) F71CH4=H0+X*(H1-H0) RETURN 100 TL=F69CH4(P) TCC=620.D00-273.15D00 GO TO 150 120 TL=F69CH4(P) TCC=470.D00-273.15D00 150 SMIN=F35CH4(P,TL) SMAX=F35CH4(P,TCC) IF(S.LT.SMIN.OR.S.GT.SMAX) GO TO 900 T=F65CH4(P,S) IF(T.EQ.-1.0E+20.OR.S.EQ.-1.0E+20) GO TO 900 F71CH4=F25CH4(P,T) RETURN 900 F71CH4=-1.0E+20 RETURN END C ------------------------------------------------- C ** F83= WPT(P,T) : Velocity of Sound [m/s] ** C ------------------------------------------------- REAL FUNCTION WPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F83CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'WPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F83CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF WPT=FF RETURN END C ------------------------------------------- C ** F83CH4(P,T) , WPT Velocity of Sound ** C ------------------------------------------- DOUBLE PRECISION FUNCTION F83CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DATA U/1.6605402D-27/,MR/16.043D00/,NA/6.0221367D023/ DATA RM/8.314510D0/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 LOO=F51CH4(P,T) TM=T+273.15D0 IF(LOO.EQ.-1.0E+20) GO TO 900 LO=1.D00/LOO CP=F18CH4(P,T) IF(CP.EQ.-1.0E+20) GO TO 900 CV=F77CH4(P,T) IF(CV.EQ.-1.0E+20) GO TO 900 W2=RM*TM/(U*NA*MR)*CP/CV*(1+2.D0*R1CH4(LO,TM)+R3CH4(LO,TM)) F83CH4=DSQRT(W2) RETURN 900 F83CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F16 = CPPD(P) : Isobaric Specific Heat of ** C ** Saturated Liquid [J/(kg.K)] ** C ----------------------------------------------------------- REAL FUNCTION CPPD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F16CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'CPPD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F16CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF CPPD=FF RETURN END C ------------------------------------------------------- C ** F16CH4(P) , CPPD : SPECIFIC HEAT CAPACITY OF ** C ** SATURATED LIQUID ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F16CH4(P) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ IF(P.LT.0.11696D00.OR.P.GT.45.992D00) GO TO 900 VL=F49CH4(P) T=F40CH4(P) TM=T+273.15 IF(T.EQ.-1.0E+20.OR.VL.EQ.-1.0E+20) GO TO 900 LO=1.D00/VL F16CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) F16CH4=F16CH4+R*(1.D0+R1CH4(LO,TM)-R5CH4(LO,TM))**2.D0/ & (1.D0+2.0D0*R1CH4(LO,TM)+R3CH4(LO,TM)) RETURN 900 F16CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F78 = CVTDD(T) : Isochoric Specific Heat of ** C ** Saturated Vapor [J/(kg.K)] ** C ------------------------------------------------------- REAL FUNCTION CVTDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F78CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'CVTDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F78CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF CVTDD=FF RETURN END C ------------------------------------------------------- C ** F78CH4(T) , CVTDD : SPECIFIC HEAT CAPACITY OF ** C ** SATURATED VAPOR ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F78CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TM=T+273.15D00 VV=F54CH4(T) IF(VV.EQ.-1.0E+20) GO TO 900 LO=1.D00/VV F78CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) RETURN 900 F78CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F76 = CVPDD(P) : Isochoric Specific Heat of ** C ** Saturated Vapor [J/(kg.K)] ** C ------------------------------------------------------- REAL FUNCTION CVPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F76CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'CVPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F76CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF CVPDD=FF RETURN END C ------------------------------------------------------------------- C ** F76CH4(P) , CVPDD : SPECIFIC HEAT CAPACITY OF SATURATED VAPOR ** C ------------------------------------------------------------------- DOUBLE PRECISION FUNCTION F76CH4(P) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ IF(P.LT.0.11696D00.OR.P.GT.45.992D00) GO TO 900 VV=F50CH4(P) IF(VV.EQ.-1.0E+20) GO TO 900 LO=1.D00/VV T=F40CH4(P) IF(T.EQ.-1.0E+20) GO TO 900 TM=T+273.15D0 F76CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) RETURN 900 F76CH4=-1.0E+20 RETURN END C ------------------------------------------------------------ C ** F77 = CVPT(P,T) : Isochoric Specific Heat [J/(kg.K)] ** C ------------------------------------------------------------ REAL FUNCTION CVPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F77CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'CVPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F77CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF CVPT=FF RETURN END C ------------------------------------------------------- C ** F77CH4(P,T) , CVPT : ISOCHORIC SPECIFIC HEAT ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F77CH4(P,T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LE.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.197.86D00)) GO TO 100 GO TO 900 100 TM=T+273.15D00 LOO=F51CH4(P,T) LO=1.0D0/LOO F77CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) RETURN 900 F77CH4=-1.0E+20 RETURN END C --------------------------------------------------- C ** F20 = CPTDD(T) : Isobaric Specific Heat of ** C ** Saturated Vapor [J/(kg.K)] ** C --------------------------------------------------- REAL FUNCTION CPTDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F20CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'CPTDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F20CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF CPTDD=FF RETURN END C ------------------------------------------------------- C ** F20CH4(T) , CPTDD : SPECIFIC HEAT CAPACITY OF ** C ** SATURATED VAPOR ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F20CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.25095D00/ TM=T+273.15D00 VV=F54CH4(T) IF(VV.EQ.-1.0E+20) GO TO 900 LO=1.D00/VV F20CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) F20CH4=F20CH4+R*(1.D0+R1CH4(LO,TM)-R5CH4(LO,TM))**2.0D0/ & (1.0D0+2.0D0*R1CH4(LO,TM)+R3CH4(LO,TM)) RETURN 900 F20CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F19 = CPTD(T) : Isobaric Specific Heat of ** C ** Saturated Liquid [J/(kg.K)] ** C ---------------------------------------------------------- REAL FUNCTION CPTD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F19CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'CPTD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F19CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF CPTD=FF RETURN END C ------------------------------------------------------- C ** F19CH4(T) , CPTD : SPECIFIC HEAT CAPACITY OF ** C ** SATURATED LIQUID ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F19CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TM=T+273.15D00 VL=F53CH4(T) IF(VL.EQ.-1.0E+20) GO TO 900 LO=1.D00/VL F19CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) F19CH4=F19CH4+R*(1.D0+R1CH4(LO,TM)-R5CH4(LO,TM))**2.0D0/ & (1.0D0+2.0D0*R1CH4(LO,TM)+R3CH4(LO,TM)) RETURN 900 F19CH4=-1.0E+20 RETURN END C ----------------------------------------------- C ** F18 = CPPT(P,T) Isobaric Specific Heat ** C ----------------------------------------------- REAL FUNCTION CPPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F18CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'CPPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F18CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF CPPT=FF RETURN END C ------------------------------------------------------------ C ** F18CH4(P,T), CPPT : SPECIFIC HEAT CAPACITY [J/(kg.k)] ** C ------------------------------------------------------------ DOUBLE PRECISION FUNCTION F18CH4(P,T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 LOO=F51CH4(P,T) IF(LOO.EQ.-1.0E+20) GO TO 900 LO=1.D00/LOO TM=T+273.15D0 F18CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) F18CH4=F18CH4+R*(1.D0+R1CH4(LO,TM)-R5CH4(LO,TM))**2.0D0/ & (1.0D0+2.0D0*R1CH4(LO,TM)+R3CH4(LO,TM)) RETURN 900 F18CH4=-1.0E+20 RETURN END C --------------------------------------------------- C ** F17 = CPPDD(P) : Isobaric Specific Heat of ** C ** Saturated Vapor [J/(kg.K)] ** C --------------------------------------------------- REAL FUNCTION CPPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F17CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'CPPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F17CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF CPPDD=FF RETURN END C ------------------------------------------------------- C ** F17CH4(P) , CPPD : SPECIFIC HEAT CAPACITY OF ** C ** SATURATED VAPOR [J/(kg.K)] ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F17CH4(P) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ IF(P.LT.0.11696D00.OR.P.GT.45.95001D00) GO TO 900 VV=F50CH4(P) IF(VV.EQ.-1.0E+20) GO TO 900 LO=1.D00/VV T=F40CH4(P) TM=T+273.15D0 IF(T.EQ.-1.0E+20) GO TO 900 F17CH4=-R*(I4CH4(LO,TM)+R4CH4(LO,TM)) F17CH4=F17CH4+R*(1.D0+R1CH4(LO,TM)-R5CH4(LO,TM))**2.0D0/ & (1.0D0+2.0D0*R1CH4(LO,TM)+R3CH4(LO,TM)) RETURN 900 F17CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F61 = XTS(T,S) : Dryness Fraction ** C ------------------------------------------------------- REAL FUNCTION XTS(T,S) CHARACTER FUN*6,FLUID*8 REAL T,S,FF INTEGER KPA DOUBLE PRECISION F61CH4,DBT,SS COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XTS'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) SS=DBLE(S) FF = F61CH4(DBT,SS) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,S 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND S =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XTS=FF RETURN END C --------------------------------------------------- C ** F62 = XTU(T,U) Dryness Fraction ** C --------------------------------------------------- REAL FUNCTION XTU(T,U) CHARACTER FUN*6,FLUID*8 REAL T,U,FF INTEGER KPA DOUBLE PRECISION F62CH4,DBT,UU COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XTU'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) UU=DBLE(U) FF = F62CH4(DBT,UU) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,U 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND U =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XTU=FF RETURN END C ----------------------------------------------- C ** F61CH4(T,S) , XTS(T,S) : DRYNESS FRACTION ** C ----------------------------------------------- DOUBLE PRECISION FUNCTION F61CH4(T,S) IMPLICIT DOUBLE PRECISION(A-H,L-Z) SL=F37CH4(T) SV=F38CH4(T) IF(S.LT.SL.OR.S.GT.SV) GO TO 900 IF(SL.EQ.-1.0E+20.OR.SV.EQ.-1.0E+20) GO TO 900 F61CH4=(S-SL)/(SV-SL) RETURN 900 F61CH4=-1.0E+20 RETURN END C ------------------------------------------------- C ** F62CH4(T,U) , XTU(T,U) : DRYNESS FRACTION ** C ------------------------------------------------- DOUBLE PRECISION FUNCTION F62CH4(T,U) IMPLICIT DOUBLE PRECISION(A-H,L-Z) UL=F46CH4(T) UV=F47CH4(T) IF(UL.EQ.-1.0E+20.OR.UV.EQ.-1.0E+20) GO TO 900 IF(U.LT.UL.OR.U.GT.UV) GO TO 900 F62CH4=(U-UL)/(UV-UL) RETURN 900 F62CH4=-1.0E+20 RETURN END C ----------------------------------------------- C ** F60 = XTH(T,H) : Dryness Fraction [-] ** C ----------------------------------------------- REAL FUNCTION XTH(T,H) CHARACTER FUN*6,FLUID*8 REAL T,H,FF INTEGER KPA DOUBLE PRECISION F60CH4,DBT,HH COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XTH'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) HH=DBLE(H) FF = F60CH4(DBT,HH) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,H 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND H =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XTH=FF RETURN END C --------------------------------------------- C ** F60CH4(T,H), XTH DRYNESS FRACTION ** C --------------------------------------------- DOUBLE PRECISION FUNCTION F60CH4(T,H) IMPLICIT DOUBLE PRECISION(A-H,L-Z) HL=F27CH4(T) HV=F28CH4(T) IF(HL.EQ.-1.0E+20.OR.HV.EQ.-1.0E+20) GO TO 900 IF(H.LT.HL.OR.H.GT.HV) GO TO 900 F60CH4=(H-HL)/(HV-HL) RETURN 900 F60CH4=-1.0E+20 RETURN END C ----------------------------------------------- C ** F58 = XPU(P,U) : Dryness Fraction ** C ----------------------------------------------- REAL FUNCTION XPU(P,U) CHARACTER FUN*6,FLUID*8 REAL P,U,FF INTEGER KPA DOUBLE PRECISION F58CH4,DBP,UU COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XPU'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) UU=DBLE(U) FF = F58CH4(DBP,UU) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,U 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND U =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XPU=FF RETURN END C ----------------------------------------------- C ** F58CH4(P,U) , XPU : DRYNESS FRACTION ** C ----------------------------------------------- DOUBLE PRECISION FUNCTION F58CH4(P,U) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA PC/45.992D00/,PL/0.11696D00/ IF(P.LT.PL.OR.P.GT.PC) GO TO 900 UL=F42CH4(P) UV=F43CH4(P) IF(U.LT.UL.OR.U.GT.UV) GO TO 900 IF(UL.EQ.-1.0E+20.OR.UV.EQ.-1.0E+20) GO TO 900 F58CH4=(U-UL)/(UV-UL) RETURN 900 F58CH4=-1.0E+20 RETURN END C --------------------------------------------------- C ** F57 = XPS(P,S) : Dryness Fraction ** C --------------------------------------------------- REAL FUNCTION XPS(P,S) CHARACTER FUN*6,FLUID*8 REAL P,S,FF INTEGER KPA DOUBLE PRECISION F57CH4,DBP,SS COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XPS'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) SS=DBLE(S) FF = F57CH4(DBP,SS) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,S 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND S =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XPS=FF RETURN END C --------------------------------------------------- C ** F57CH4(P,S) : XPS(P,S) DRYNESS FRACTION ** C --------------------------------------------------- DOUBLE PRECISION FUNCTION F57CH4(P,S) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA PC/45.992D00/,PL/0.11696D00/ IF(P.LT.PL.OR.P.GT.PC) GO TO 900 SL=F33CH4(P) SV=F34CH4(P) IF(S.LT.SL.OR.S.GT.SV) GO TO 900 IF(SL.EQ.-1.0E+20.OR.SV.EQ.-1.0E+20) GO TO 900 F57CH4=(S-SL)/(SV-SL) RETURN 900 F57CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F56 = XPH : Dryness Fraction [-] ** C ------------------------------------------------------- REAL FUNCTION XPH(P,H) CHARACTER FUN*6,FLUID*8 REAL P,H,FF INTEGER KPA DOUBLE PRECISION F56CH4,DBP,HH COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XPH'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) HH=DBLE(H) FF = F56CH4(DBP,HH) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,H 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND H =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XPH=FF RETURN END C ------------------------------------------------------- C ** F56CH4(P,H),XPH : DRYNESS FRACTION ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F56CH4(P,H) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA PC/45.992D00/,PL/0.11696D00/ IF(P.LT.PL.OR.P.GT.PC) GO TO 900 HL=F23CH4(P) HV=F24CH4(P) IF(H.LT.HL.OR.H.GT.HV) GO TO 900 IF(HL.EQ.-1.0E+20.OR.HV.EQ.-1.0E+20) GO TO 900 F56CH4=(H-HL)/(HV-HL) RETURN 900 F56CH4=-1.0E+20 RETURN END C ------------------------------------------------------------------ C ** F45 = UPX(P,X) : Specific Internal Energy of Mixture [J/kg] ** C ------------------------------------------------------------------ REAL FUNCTION UPX(P,X) CHARACTER FUN*6,FLUID*8 REAL P,X,FF INTEGER KPA DOUBLE PRECISION F45CH4,DBP,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'UPX'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) XX=DBLE(X) FF = F45CH4(DBP,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF UPX=FF RETURN END C ------------------------------------------------------------------ C ** F48 = UTX(P,X) : Specific Internal Energy of Mixture [J/kg] ** C ------------------------------------------------------------------ REAL FUNCTION UTX(T,X) CHARACTER FUN*6,FLUID*8 REAL T,FF INTEGER KPA DOUBLE PRECISION F48CH4,DBT,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'UTX'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) XX=DBLE(X) FF = F48CH4(DBT,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF UTX=FF RETURN END C ---------------------------------------------------------------- C ** F45CH4(P,X),UPX(P,X) : SPECIFIC INTERNAL ENERGY OF MIXTURE ** C ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION F45CH4(P,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF((P.LT.0.11696D00.OR.P.GT.45.992D00) & .OR.(X.LT.0.0.OR.X.GT.1.0)) GO TO 900 UL=F42CH4(P) UV=F43CH4(P) IF(UL.EQ.-1.0E+20.OR.UV.EQ.-1.0E+20) GO TO 900 F45CH4=UL+X*(UV-UL) RETURN 900 F45CH4=-1.0E+20 RETURN END C ---------------------------------------------------------------- C ** F48CH4(P,X),UTX(P,X) : SPECIFIC INTERNAL ENERGY OF MIXTURE ** C ---------------------------------------------------------------- DOUBLE PRECISION FUNCTION F48CH4(T,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF(X.LT.0.0.OR.X.GT.1.0) GO TO 900 UL=F46CH4(T) UV=F47CH4(T) IF(UL.EQ.-1.0E+20.OR.UV.EQ.-1.0E+20) GO TO 900 F48CH4=UL+X*(UV-UL) RETURN 900 F48CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F46 = UTD(T) Specific Internal Energy of ** C ** Saturated Liquid [J/KG] ** C ----------------------------------------------------------- REAL FUNCTION UTD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F46CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'UTD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F46CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF UTD=FF RETURN END C ----------------------------------------------------------- C ** F47 = UTDD(T) Specific Internal Energy of ** C ** Saturated VAPOR [J/kg] ** C ----------------------------------------------------------- REAL FUNCTION UTDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F47CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'UTDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F47CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF UTDD=FF RETURN END C ------------------------------------------------------- C ** F46CH4(T) , UTD : Specific Internal Energy of ** C ** SATURATED LIQUID ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F46CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TM=T+273.15D00 VL=F53CH4(T) IF(VL.EQ.-1.0E+20) GO TO 900 LO=1.D00/VL F46CH4=R*TM*(I2CH4(LO,TM)+R2CH4(LO,TM)) RETURN 900 F46CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F47CH4(T) , UTDD : Specific Internal Energy of ** C ** SATURATED VAPOR ** C ------------------------------------------------------- C *** UTDD(T) SPECIFIC INTERNAL ENERGY OF SATURATED VAPOR DOUBLE PRECISION FUNCTION F47CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TM=T+273.15D00 VV=F54CH4(T) IF(VV.EQ.-1.0E+20) GO TO 900 LO=1.D00/VV F47CH4=R*TM*(I2CH4(LO,TM)+R2CH4(LO,TM)) RETURN 900 F47CH4=-1.0E+20 RETURN END C ------------------------------------------------------------- C ** F42=UPD(P) Specific Internal Energy of Saturated LIQUID ** C ------------------------------------------------------------- REAL FUNCTION UPD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F42CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'UPD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F42CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF UPD=FF RETURN END C ------------------------------------------------------------- C ** F43=UPDD(P) Specific Internal Energy of Saturated VAPOR ** C ------------------------------------------------------------- REAL FUNCTION UPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F43CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'UPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F43CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF UPDD=FF RETURN END C ------------------------------------------------------ C ** F42CH4(P) , UPD(P) : SPECIFIC INTERNAL ENERGY OF ** C ** SATURATED LIQUID ** C ------------------------------------------------------ DOUBLE PRECISION FUNCTION F42CH4(P) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ IF(P.LT.0.11719D00.OR.P.GT.45.95001D00) GO TO 900 T=F40CH4(P) TM=T+273.15D00 VL=F53CH4(T) IF(T.EQ.-1.0E+20.OR.VL.EQ.-1.0E+20) GO TO 900 LO=1.D0/VL F42CH4=R*TM*(I2CH4(LO,TM)+R2CH4(LO,TM)) RETURN 900 F42CH4=-1.0E+20 RETURN END C ------------------------------------------------------ C ** F42CH4(P) , UPDD(P): SPECIFIC INTERNAL ENERGY OF ** C ** SATURATED VAPOR ** C ------------------------------------------------------ DOUBLE PRECISION FUNCTION F43CH4(P) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ T=F40CH4(P) VL=F54CH4(T) IF(T.EQ.-1.0E+20.OR.VL.EQ.-1.0E+20) GO TO 900 TM=T+273.15D00 LO=1.D0/VL F43CH4=R*TM*(I2CH4(LO,TM)+R2CH4(LO,TM)) RETURN 900 F43CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F44 = UPT(P,T) : Specific Internal Energy [J/kg] ** C ----------------------------------------------------------- REAL FUNCTION UPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F44CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'UPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F44CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF UPT=FF RETURN END C --------------------------------------------------- C ** F44CH4(P,T), UPT : SPECIFIC INTERNAL ENERGY ** C --------------------------------------------------- DOUBLE PRECISION FUNCTION F44CH4(P,T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 TM=T+273.15D00 VLM=F51CH4(P,T) IF(VLM.EQ.-1.0E+20) GO TO 900 LO=1.D00/VLM C CALL S03CH4(T,HID) C HID=HID*1000.D00/16.043D00 C CALL S12CH4(T,LO,ULL) C UID=HID-R*TM C F44CH4=UID+ULL F44CH4=R*TM*(I2CH4(LO,TM)+R2CH4(LO,TM)) RETURN 900 F44CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F52 = VPX: Specific Volume of Mixture [m**3/kg] ** C ---------------------------------------------------------- REAL FUNCTION VPX(P,X) CHARACTER FUN*6,FLUID*8 REAL P,X,FF INTEGER KPA DOUBLE PRECISION F52CH4,DBP,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'VPX'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) XX=DBLE(X) FF = F52CH4(DBP,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF VPX=FF RETURN END C ------------------------------------------------------- C ** F52CH4(P,X) : VPX SPECIFIC VOLUME OF MIXTURE ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F52CH4(P,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF(X.LT.0.0.OR.X.GT.1.0) GO TO 900 VL=F49CH4(P) VV=F50CH4(P) IF(VL.EQ.-1.0E+20.OR.VV.EQ.-1.0E+20) GO TO 900 F52CH4=VL+X*(VV-VL) RETURN 900 F52CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F55 = VTX(T,X) Specific Volume of Mixture [M^3/kg] ** C ----------------------------------------------------------- REAL FUNCTION VTX(T,X) CHARACTER FUN*6,FLUID*8 REAL T,FF INTEGER KPA DOUBLE PRECISION F55CH4,DBT,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'VTX'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) XX=DBLE(X) FF = F55CH4(DBT,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF VTX=FF RETURN END C ------------------------------------------------------- C ** F55CH4 : VTX(T,X) SPECIFIC VOLUME OF MIXTURE ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F55CH4(T,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) TM=T+273.15D00 IF(X.LT.0.0.OR.X.GT.1.0) GO TO 900 VL=F53CH4(T) VV=F54CH4(T) IF(VL.EQ.-1.0E+20.OR.VV.EQ.-1.0E+20) GO TO 900 F55CH4=VL+X*(VV-VL) RETURN 900 F55CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F36 = SPX : Specific Entropy of Mixture [J/(kg.K)] ** C ---------------------------------------------------------- REAL FUNCTION SPX(P,X) CHARACTER FUN*6,FLUID*8 REAL P,X,FF INTEGER KPA DOUBLE PRECISION F36CH4,DBP,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'SPX'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) XX=DBLE(X) FF = F36CH4(DBP,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF SPX=FF RETURN END C ------------------------------------------------------- C ** F36CH4 : SPX(P,X) SPECIFIC ENTROPY OF MIXTURE ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F36CH4(P,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF((P.LT.0.11696D00.OR.P.GT.45.992D00) & .OR.(X.LT.0.0.OR.X.GT.1.0)) GO TO 900 SL=F33CH4(P) SV=F34CH4(P) IF(SL.EQ.-1.0E+20.OR.SV.EQ.-1.0E+20) GO TO 900 F36CH4=SL+X*(SV-SL) RETURN 900 F36CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F23=HPD(P) Specific Enthalpy of Saturated Liquid ** C ---------------------------------------------------------- REAL FUNCTION HPD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F23CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'HPD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F23CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF HPD=FF RETURN END C ---------------------------------------------------------- C ** F24 = HPDD(P) : Specific Enthalpy of Saturated Vapor ** C ---------------------------------------------------------- REAL FUNCTION HPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F24CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'HPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F24CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF HPDD=FF RETURN END C ----------------------------------------------------------- C ** F23CH4 HPD(P) SPECIFIC ENTHALPY OF SATURATED LIQUID ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F23CH4(P) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ VL=F49CH4(P) IF (VL.EQ.-1.0E+20) GO TO 600 LO=1.D00/VL T=F40CH4(P) TM=T+273.15D00 F23CH4=R*TM*(1.D0+I2CH4(LO,TM)+R2CH4(LO,TM)+R1CH4(LO,TM)) RETURN 600 F23CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F24CH4 HPDD(P) SPECIFIC ENTHALPY OF SATURATED VAPOR ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F24CH4(P) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ LV=F50CH4(P) IF(LV.EQ.-1.0E+20) GO TO 600 LO=1.D00/LV T=F40CH4(P) TM=T+273.15D00 F24CH4=R*TM*(1.D0+I2CH4(LO,TM)+R2CH4(LO,TM)+R1CH4(LO,TM)) RETURN 600 F24CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F26=HPX(P,X) : Specific Enthalpy of Mixture [J/kg] ** C ---------------------------------------------------------- REAL FUNCTION HPX(P,X) CHARACTER FUN*6,FLUID*8 REAL P,X,FF INTEGER KPA DOUBLE PRECISION F26CH4,DBP,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'HPX'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) XX=DBLE(X) FF = F26CH4(DBP,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF HPX=FF RETURN END C ------------------------------------------------------------ C ** F26CH4 : HPX(P,X) SPECIFIC ENTHALPY OF MIXTURE [J/kg] ** C ------------------------------------------------------------ DOUBLE PRECISION FUNCTION F26CH4(P,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF((P.LT.0.11696D00.OR.P.GT.45.992D00) & .OR.(X.LT.0.0.OR.X.GT.1.0)) GO TO 900 HL=F23CH4(P) HV=F24CH4(P) IF(HL.EQ.-1.0E+20.OR.HV.EQ.-1.0E+20) GO TO 900 F26CH4=HL+X*(HV-HL) RETURN 900 F26CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F39=STX(T,X) Specific Entropy of Mixture [J/(kg.K)] ** C ---------------------------------------------------------- REAL FUNCTION STX(T,X) CHARACTER FUN*6,FLUID*8 REAL T,FF INTEGER KPA DOUBLE PRECISION F39CH4,DBT,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'STX'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) XX=DBLE(X) FF = F39CH4(DBT,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF STX=FF RETURN END C ----------------------------------------------------- C ** F39CH4 = STX(T,X) : SPECIFIC ENTROPY OF MIXTURE ** C ----------------------------------------------------- DOUBLE PRECISION FUNCTION F39CH4(T,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF(X.LT.0.0.OR.X.GT.1.0) GO TO 900 SL=F37CH4(T) SV=F38CH4(T) IF(SL.EQ.-1.0E+20.OR.SV.EQ.-1.0E+20) GO TO 900 F39CH4=SL+X*(SV-SL) RETURN 900 F39CH4=-1.0E+20 RETURN END C --------------------------------------------------- C ** F29=HTX(T,X) : Specific Enthalpy of Mixture ** C --------------------------------------------------- REAL FUNCTION HTX(T,X) CHARACTER FUN*6,FLUID*8 REAL T,FF INTEGER KPA DOUBLE PRECISION F29CH4,DBT,XX COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'HTX'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) XX=DBLE(X) FF = F29CH4(DBT,XX) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,X 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND X =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF HTX=FF RETURN END C ------------------------------------------------------ C ** F29CH4, HTX(T,X) SPECIFIC ENTHALPY OF MIXTURE ** C ------------------------------------------------------ DOUBLE PRECISION FUNCTION F29CH4(T,X) IMPLICIT DOUBLE PRECISION(A-H,L-Z) TM=T+273.15D00 IF(X.LT.0.0.OR.X.GT.1.0) GO TO 900 HL=F27CH4(T) HV=F28CH4(T) IF(HL.EQ.-1.0E+20.OR.HV.EQ.-1.0E+20) GO TO 900 F29CH4=HL+X*(HV-HL) RETURN 900 F29CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F37 , STD(T) : Specific Entropy of Saturated Liquid ** C ---------------------------------------------------------- REAL FUNCTION STD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F37CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'STD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F37CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF STD=FF RETURN END C ---------------------------------------------------------- C ** F38 , STDD(T) : Specific Entropy of Saturated Vapor ** C ---------------------------------------------------------- REAL FUNCTION STDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F38CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'STDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F38CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF STDD=FF RETURN END C --------------------------------------------------------- C ** F37CH4(T) : SPECIFIC ENTROPY OF SATURATED LIQUID ** C --------------------------------------------------------- DOUBLE PRECISION FUNCTION F37CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA DLT/1.D-05/, & R/518.2645D00/,TL/90.6854D00/,TC/190.551D00/, & SL/4234.869D00/,SC/7511.064D00/ TM=T+273.15D00 IF (DABS((TM-TL)/TL).LT.DLT) GO TO 300 IF (DABS((TM-TC)/TC).LT.DLT) GO TO 350 IF(TM.LT.TL) GO TO 900 IF(TM.GT.TC) GO TO 900 LO=1.D00/F53CH4(T) F37CH4=-R*(I0CH4(LO,TM)+R0CH4(LO,TM)-I2CH4(LO,TM)-R2CH4(LO,TM)) RETURN 300 F37CH4=SL RETURN 350 F37CH4=SC RETURN 900 F37CH4=-1.0E+20 RETURN END C ------------------------------------------------------------ C ** F38CH4 , STDD(T) : SPECIFIC ENTROPY OF SATURATED VAPOR ** C ------------------------------------------------------------ DOUBLE PRECISION FUNCTION F38CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA DLT/1.D-05/,R/518.264D00/,TC/190.551D00/, & SL/10222.527D00/,SC/7511.064D00/, & TL/90.6854D00/ TM=T+273.15D00 IF(DABS((TM-TC)/TC).LT.DLT) GO TO 300 IF(DABS((TM-TL)/TL).LT.DLT) GO TO 350 IF(TM.LT.TL) GO TO 600 IF(TM.GT.TC) GO TO 600 LO=1.D00/F54CH4(T) F38CH4=-R*(I0CH4(LO,TM)+R0CH4(LO,TM)-I2CH4(LO,TM)-R2CH4(LO,TM)) RETURN 300 F38CH4=SC RETURN 350 F38CH4=SL RETURN 600 F38CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F33 , SPD(P): Specific Entropy of Saturated Liquid ** C ---------------------------------------------------------- REAL FUNCTION SPD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F33CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'SPD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F33CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF SPD=FF RETURN END C ---------------------------------------------------------- C ** F34 , SPDD(P) : Specific Entropy of Saturared Vapor ** C ---------------------------------------------------------- REAL FUNCTION SPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F34CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'SPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F34CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF SPDD=FF RETURN END C ----------------------------------------------------------- C ** F33CH4 : SPD(P) SPECIFIC ENTROPY OF SATURATED LIQUID ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F33CH4(P) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DOUBLE PRECISION I0CH4,I2CH4 DATA PC/45.992D00/,PL/0.11696D00/,DLT/1.D-05/,R/518.264D00/, & SL/4234.8688D00/,SC/7511.064D00/ IF(DABS((P-PC)/PC).LT.DLT) GO TO 300 IF(DABS((P-PL)/PL).LT.DLT) GO TO 350 IF(P.LT.PL) GO TO 600 IF(P.GT.PC) GO TO 600 T=F40CH4(P) LO=1.D00/F53CH4(T) IF(T.EQ.-1.0E+20.OR.LO.EQ.-1.0E+20) GO TO 600 TM=T+273.15D00 F33CH4=-R*(I0CH4(LO,TM)+R0CH4(LO,TM)-I2CH4(LO,TM)-R2CH4(LO,TM)) RETURN 300 F33CH4=SC RETURN 350 F33CH4=SL RETURN 600 F33CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F34, SPDD(P) : SPECIFIC ENTROPY OF SATURATED VAPOR ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F34CH4(P) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DOUBLE PRECISION I0CH4,I2CH4 DATA PC/45.992D00/,PL/0.11696D00/,DLT/1.D-05/,R/518.264D00/, & SL/10222.527/,SC/7155.064D00/ IF(DABS((P-PC)/PC).LT.DLT) GO TO 300 IF(DABS((P-PL)/PL).LT.DLT) GO TO 350 IF(P.LT.PL) GO TO 900 IF(P.GT.PC) GO TO 900 T=F40CH4(P) TM=T+273.15D00 LO=1.D00/F54CH4(T) IF(T.EQ.-1.0E+20.OR.LO.EQ.-1.0E+20) GO TO 900 F34CH4=-R*(I0CH4(LO,TM)+R0CH4(LO,TM)-I2CH4(LO,TM)-R2CH4(LO,TM)) RETURN 300 F34CH4=SC RETURN 350 F34CH4=SL RETURN 900 F34CH4=-1.0E+20 RETURN END C --------------------------------------------------- C ** F35 , SPT : Specific Entropy [J/(kg.k)] ** C --------------------------------------------------- REAL FUNCTION SPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F35CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'SPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F35CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF SPT=FF RETURN END C ------------------------------------------------------- C ** F35CH4 , SPT(P,T) : SPECIFIC ENTROPY [J/(KG.K)] ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F35CH4(P,T) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DOUBLE PRECISION I0CH4,I2CH4 DATA ERR/1.D-7/,R/518.264D00/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 TM=T+273.15D00 LOO=F51CH4(P,T) IF(LOO.EQ.ERR) GO TO 900 LO=1.D00/LOO F35CH4=-R*(I0CH4(LO,TM)+R0CH4(LO,TM)-I2CH4(LO,TM)-R2CH4(LO,TM)) RETURN 900 F35CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F25 : HPT(P,T) Specific Enthalpy [J/kg] ** C ---------------------------------------------------------- REAL FUNCTION HPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F25CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'HPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F25CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF HPT=FF RETURN END C --------------------------------------------------- C ** F25CH4(P,T): HPT(P,T) SPECIFIC ENTHALPY ** C --------------------------------------------------- DOUBLE PRECISION FUNCTION F25CH4(P,T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LT.346.86D00)) GO TO 100 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.196.86D00)) GO TO 100 GO TO 900 100 LOO=F51CH4(P,T) IF(LOO.EQ.-1.0E+20) GO TO 900 TM=T+273.15D00 LO=1.D00/LOO F25CH4=R*TM*(1.D0+I2CH4(LO,TM)+R2CH4(LO,TM)+R1CH4(LO,TM)) RETURN 900 F25CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F27 : HTD(T) Specific Enthalpy of Saturated Liquid ** C ----------------------------------------------------------- REAL FUNCTION HTD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F27CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'HTD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F27CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF HTD=FF RETURN END C ----------------------------------------------------------- C ** F28 : HTDD(T) Specific Enthalpy of Saturated Vapor ** C ----------------------------------------------------------- REAL FUNCTION HTDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F28CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'HTDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F28CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF HTDD=FF RETURN END C ----------------------------------------------------------- C ** F27CH4: HTD(T) SPECIFIC ENTHALPY OF SATURATED LIQUID ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F27CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TM=T+273.15D00 P=F30CH4(T) IF(P.EQ.-1.0E+20) GO TO 900 VL=F53CH4(T) IF(VL.EQ.-1.0E+20) GO TO 900 LO=1.D0/VL F27CH4=R*TM*(1+I2CH4(LO,TM)+R2CH4(LO,TM)+R1CH4(LO,TM)) RETURN 900 F27CH4=-1.0E+20 RETURN END C ----------------------------------------------------------- C ** F28CH4: HTDD(T) SPECIFIC ENTHALPY OF SATURATED VAPOR ** C ----------------------------------------------------------- DOUBLE PRECISION FUNCTION F28CH4(T) IMPLICIT DOUBLE PRECISION(A-I,L-Z) DATA R/518.264D00/ TM=T+273.15D00 P=F30CH4(T) VL=F54CH4(T) IF(VL.EQ.-1.0E+20.OR.P.EQ.-1.0E+20) GO TO 900 LO=1.D0/VL F28CH4=R*TM*(1+I2CH4(LO,TM)+R2CH4(LO,TM)+R1CH4(LO,TM)) RETURN 900 F28CH4=-1.0E+20 RETURN END C ----------------------------------------------- C ** F63 = XTV ** C ** Dryness Fraction [-] ** C ----------------------------------------------- REAL FUNCTION XTV(T,V) CHARACTER FUN*6,FLUID*8 REAL T,V,FF INTEGER KPA DOUBLE PRECISION F63CH4,DBT,VV COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XTV'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) VV=DBLE(V) FF = F63CH4(DBT,VV) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,T,V 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN T =',1PE14.7,' AND V =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XTV=FF RETURN END C --------------------------------------------------- C ** F63CH4 : XTV(T,V) DRYNESS FRACTION ** C --------------------------------------------------- DOUBLE PRECISION FUNCTION F63CH4(T,V) IMPLICIT DOUBLE PRECISION(A-H,L-Z) VL=F53CH4(T) VV=F54CH4(T) IF(V.LT.VL.OR.V.GT.VV) GO TO 900 IF(VL.EQ.-1.0E+20.OR.VV.EQ.-1.0E+20) GO TO 900 F63CH4=(V-VL)/(VV-VL) RETURN 900 F63CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F59 = XPV ** C ** Dryness Fraction [-] ** C ---------------------------------------------------------- REAL FUNCTION XPV(P,V) CHARACTER FUN*6,FLUID*8 REAL P,V,FF INTEGER KPA DOUBLE PRECISION F59CH4,DBP,VV COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'XPV'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) VV=DBLE(V) FF = F59CH4(DBP,VV) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,V 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND V =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF XPV=FF RETURN END C ------------------------------------------------------- C ** F59CH4: XPV(P,V) DRYNESS FRACTION ** C ------------------------------------------------------- DOUBLE PRECISION FUNCTION F59CH4(P,V) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA PC/45.95001D00/,PL/0.11719D00/ IF(P.LT.PL.OR.P.GT.PC) GO TO 900 VL=F49CH4(P) VV=F50CH4(P) IF(V.LT.VL.OR.V.GT.VV) GO TO 900 IF(VL.EQ.-1.0E+20.OR.VV.EQ.-1.0E+20) GO TO 900 F59CH4=(V-VL)/(VV-VL) RETURN 900 F59CH4=-1.0E+20 RETURN END C ------------------------------------------------------- C ** F70 = TPV(P,V) ** C ** INPUT: Pressure [Pa],[bar], Specific Volume ** C ** OUTPUT: Temperature [K],[C] ** C ------------------------------------------------------- REAL FUNCTION TPV(P,V) CHARACTER FUN*6,FLUID*8 REAL T0K,PBAR,P,V,FF INTEGER KPA DOUBLE PRECISION F70CH4,DBP,VV COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'TPV'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF PI=P*PBAR DBP=DBLE(PI) VV=DBLE(V) FF = F70CH4(DBP,VV) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,P,V 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN P =',1PE14.7,' AND V =',1PE14.7,' ****') FF=-1.0E+20 END IF END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 TPV=FF+T0K RETURN END C --------------------------------------------------- C ** FUNCTION F70CH4 :TPV(P,V) TEMPERATURE ** C --------------------------------------------------- DOUBLE PRECISION FUNCTION F70CH4(P,V) IMPLICIT DOUBLE PRECISION(A-H,L-Z) DATA R/518.264D00/,ERR/1.D-07/,PC/45.992D00/ IF(P.GE.0.11719D00.AND.P.LE.PC) GO TO 50 40 IF(P.GE.0.11719D00.AND.P.LT.450.D00) GO TO 100 IF(P.GE.450.D00.AND.P.LE.10000.D00) GO TO 120 GO TO 900 50 TS=F40CH4(P) V0=F49CH4(P) V1=F50CH4(P) IF(V.LT.V0.OR.V.GT.V1) GO TO 40 F70CH4=TS RETURN 100 TL=F69CH4(P) TCC=620.D00-273.15D00 X0=TL+273.15D00 X1=TCC+273.15D00 X3=X1 GO TO 150 120 TL=F69CH4(P) TCC=470.D00-273.15D00 X0=TL+273.15D00 X1=TCC+273.15D00 X3=X1 150 V0=F51CH4(P,TL ) V1=F51CH4(P,TCC) IF(V.LT.V0.OR.V.GT.V1) GO TO 900 LO=1.0D0/V F0=1.0D5*P-LO*R*X0*(1+R1CH4(LO,X0)) F1=1.0D5*P-LO*R*X1*(1+R1CH4(LO,X1)) 250 X2=(X0+X1)/2.D00 F2=1.0D5*P-LO*R*X2*(1+R1CH4(LO,X2)) IF(DABS(X2-X3).LT.ERR) GO TO 400 X3=X2 IF((F0*F2).LT.0.) GO TO 350 F0=F2 X0=X2 GO TO 250 350 F1=F2 X1=X2 GO TO 250 400 F70CH4=X2-273.15D00 RETURN 900 F70CH4=-1.0E+20 RETURN END C ---------------------------------------------------------- C ** F51 = VPT ** C ** INPUT : Pressure [PA],[BAR], Temperature [K],[C] ** C ** OUTPUT : Specific Volume [m**3/kg] ** C ---------------------------------------------------------- REAL FUNCTION VPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F51CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'VPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F51CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF VPT=FF RETURN END C ------------------------------------------------ C *** F51CH4(P,T) VPT: SPECIFIC VOLUME *** C ------------------------------------------------ DOUBLE PRECISION FUNCTION F51CH4(P,T) DOUBLE PRECISION P,TP,T,TC,R,ERR,PC,TM,T1,T2 DOUBLE PRECISION X0,X1,X2,F0,F1,F2 DOUBLE PRECISION R1CH4,F53CH4,F54CH4,F40CH4,F69CH4 DATA TC/190.551D00/,R/518.264D00/,ERR/1.D-07/,PC/45.992D00/ TP=F69CH4(P) IF(TP.EQ.-1.0E+20) GO TO 900 IF((P.GE.0.11696D00.AND.P.LT.450.D00) & .AND.(T.GE.TP.AND.T.LE.346.86D00)) GO TO 40 IF((P.GE.450.D00.AND.P.LE.10000.D00) & .AND.(T.GE.TP.AND.T.LE.197.86D00)) GO TO 40 GO TO 900 40 TM=T+273.15D00 IF(TM.GT.TC.OR.P.GT.PC) GO TO 90 T1=TP T2=F40CH4(P) IF(T.GE.T1.AND.T.LE.T2) GO TO 60 50 X0=0. X1=1.D00/F54CH4(T) X3=X1 GO TO 100 60 X0=1.D00/F53CH4(T)-0.5D00 C 60 X0=1.D00/F53CH4(T) X1=600. X3=X1 GO TO 100 C 90 X0=0. 90 IF(TM.GT.TC) THEN X0=0.0D0 ELSE X0=1.D00/F53CH4(T) ENDIF X1=600. X3=X1 100 F0=1.D5*P-X0*R*TM*(1.D0+R1CH4(X0,TM)) F1=1.D5*P-X1*R*TM*(1.D0+R1CH4(X1,TM)) 70 X2=0.5*(X0+X1) F2=1.D5*P-X2*R*TM*(1.D0+R1CH4(X2,TM)) IF(DABS(F2/(1.0D5*P)).LT.ERR) GOTO 400 X3=X2 IF((F0*F2).LT.0.) GO TO 410 F0=F2 X0=X2 GO TO 70 410 F1=F2 X1=X2 GO TO 70 400 F51CH4=1.D00/X2 RETURN 900 F51CH4=-1.0E+20 RETURN END C -------------------------------------------------------------- C ** Exponents and Coefficients for the Residual Free Energy ** C -------------------------------------------------------------- SUBROUTINE COEFF(R,S,N) DOUBLE PRECISION N REAL R,S DIMENSION R(32),S(32),N(32) R(1)=1.E0 R(2)=1.E0 R(3)=1.E0 R(4)=2.E0 R(5)=2.E0 R(6)=2.E0 R(7)=3.E0 R(8)=3.E0 R(9)=3.E0 R(10)=6.E0 R(11)=7.E0 R(12)=7.E0 R(13)=8.E0 R(14)=1.E0 R(15)=1.E0 R(16)=2.E0 R(17)=2.E0 R(18)=3.E0 R(19)=3.E0 R(20)=5.E0 R(21)=6.E0 R(22)=7.E0 R(23)=8.E0 R(24)=10.E0 R(25)=2.E0 R(26)=3.E0 R(27)=3.E0 R(28)=4.E0 R(29)=4.E0 R(30)=5.E0 R(31)=5.E0 R(32)=5.E0 S(1)=0.E0 S(2)=1.5E0 S(3)=2.5E0 S(4)=-0.5E0 S(5)=1.5E0 S(6)=2.E0 S(7)=0.E0 S(8)=1.E0 S(9)=2.5E0 S(10)=0.E0 S(11)=2.E0 S(12)=5.E0 S(13)=2.E0 S(14)=5.E0 S(15)=6.E0 S(16)=3.5E0 S(17)=5.5E0 S(18)=3.E0 S(19)=7.E0 S(20)=6.E0 S(21)=8.5E0 S(22)=4.E0 S(23)=6.5E0 S(24)=5.5E0 S(25)=22.E0 S(26)=11.E0 S(27)=18.E0 S(28)=11.E0 S(29)=23.E0 S(30)=17.E0 S(31)=18.E0 S(32)=23.E0 N(1)=0.38443609966D00 N(2)=-0.17969259880D01 N(3)=0.32944494737D00 N(4)=0.22631272844D-1 N(5)=0.75923676880D-1 N(6)=0.69375844726D-1 N(7)=0.24116326395D-1 N(8)=0.10700992085D-1 N(9)=-0.38093327516D-1 N(10)=0.47153756114D-3 N(11)=0.55660767881D-3 N(12)=0.54875934653D-6 N(13)=-0.99963269997D-4 N(14)=-0.12808797928D00 N(15)=0.38019887338D-1 N(16)=0.13922665055D00 N(17)=-0.87499634886D-1 N(18)=-0.33489416576D-2 N(19)=-0.51757629712D-1 N(20)=0.25283517912D-1 N(21)=0.51870320595D-3 N(22)=-0.16677059452D-2 N(23)=-0.60740192739D-3 N(24)=-0.97291535999D-4 N(25)=-0.29884401046D-4 N(26)=-0.13094011124D-1 N(27)=0.19817583380D-1 N(28)=0.20846576233D-1 N(29)=-0.35802505263D-1 N(30)=-0.20348685174D00 N(31)=0.21596475509D00 N(32)=-0.42934062825D-2 RETURN END C ---------------------------------------------------- C ** the Dimensionless Residual Term of Free Energy ** C ---------------------------------------------------- DOUBLE PRECISION FUNCTION R0CH4(V,T) DOUBLE PRECISION V,T,N,FA1,FA2,FA3,DLTA,TAO REAL S,R DIMENSION N(32),S(32),R(32) DATA TC/190.551D0/,VC/162.66D0/ FA1=0.D0 FA2=0.D0 FA3=0.D0 DLTA=V/VC TAO=TC/T CALL COEFF(R,S,N) DO 111 I=1,13 FA1=FA1+N(I)*DLTA**R(I)*TAO**S(I) 111 CONTINUE DO 222 I=14,24 FA2=FA2+N(I)*DLTA**R(I)*TAO**S(I) 222 CONTINUE FA2=DEXP(-DLTA*DLTA)*FA2 DO 333 I=25,32 FA3=FA3+N(I)*DLTA**R(I)*TAO**S(I) 333 CONTINUE FA3=DEXP(-DLTA**4.)*FA3 R0CH4=FA1+FA2+FA3 R0CH4=R0CH4 RETURN END C -------------------------------------------------- C ** the Derivatives No.1 ** C ** of The Residual Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION R1CH4(V,T) DOUBLE PRECISION V,T,N,FA1,FA2,FA3,DLTA,TAO REAL S,R DIMENSION N(32),S(32),R(32) DATA TC/190.551D0/,VC/162.66D0/ FA1=0.D0 FA2=0.D0 FA3=0.D0 DLTA=V/VC TAO=TC/T CALL COEFF(R,S,N) DO 111 I=1,13 FA1=FA1+R(I)*N(I)*DLTA**R(I)*TAO**S(I) 111 CONTINUE DO 222 I=14,24 FA2=FA2+(R(I)-2.*DLTA*DLTA)*N(I)*DLTA**R(I)*TAO**S(I) 222 CONTINUE FA2=DEXP(-DLTA*DLTA)*FA2 DO 333 I=25,32 FA3=FA3+(R(I)-4.*DLTA**4)*N(I)*DLTA**R(I)*TAO**S(I) 333 CONTINUE FA3=DEXP(-DLTA**4.)*FA3 R1CH4=FA1+FA2+FA3 RETURN END C -------------------------------------------------- C ** the Derivatives No.2 ** C ** of the Residual Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION R2CH4(V,T) DOUBLE PRECISION V,T,N,FA1,FA2,FA3,DLTA,TAO REAL S,R DIMENSION N(32),S(32),R(32) DATA TC/190.551D0/,VC/162.66D0/ FA1=0.D0 FA2=0.D0 FA3=0.D0 DLTA=V/VC TAO=TC/T CALL COEFF(R,S,N) DO 111 I=1,13 FA1=FA1+S(I)*N(I)*DLTA**R(I)*TAO**S(I) 111 CONTINUE DO 222 I=14,24 FA2=FA2+S(I)*N(I)*DLTA**R(I)*TAO**S(I) 222 CONTINUE FA2=DEXP(-DLTA*DLTA)*FA2 DO 333 I=25,32 FA3=FA3+S(I)*N(I)*DLTA**R(I)*TAO**S(I) 333 CONTINUE FA3=DEXP(-DLTA**4.)*FA3 R2CH4=FA1+FA2+FA3 R2CH4=R2CH4 RETURN END C -------------------------------------------------- C ** the Derivatives No.3 ** C ** of the Residual Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION R3CH4(V,T) DOUBLE PRECISION V,T,N,FA1,FA2,FA3,DLTA,TAO REAL S,R DIMENSION N(32),S(32),R(32) DATA TC/190.551D0/,VC/162.66D0/ FA1=0.D0 FA2=0.D0 FA3=0.D0 DLTA=V/VC TAO=TC/T CALL COEFF(R,S,N) DO 111 I=1,13 FA1=FA1+(R(I)*(R(I)-1))*N(I)*DLTA**R(I)*TAO**S(I) 111 CONTINUE DO 222 I=14,24 FA2=FA2+(R(I)*(R(I)-1)-2.*(2.*R(I)+1.)*DLTA*DLTA+4*DLTA**4.) & *N(I)*DLTA**R(I)*TAO**S(I) 222 CONTINUE FA2=DEXP(-DLTA*DLTA)*FA2 DO 333 I=25,32 FA3=FA3+(R(I)*(R(I)-1)-4.D0*(2*R(I)+3)*DLTA**4.+16.*DLTA**8.0) & *N(I)*DLTA**R(I)*TAO**S(I) 333 CONTINUE FA3=DEXP(-DLTA**4)*FA3 R3CH4=FA1+FA2+FA3 R3CH4=R3CH4 RETURN END C -------------------------------------------------- C ** the Derivatives No.4 ** C ** of the Residual Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION R4CH4(V,T) DOUBLE PRECISION V,T,N,FA1,FA2,FA3,DLTA,TAO REAL S,R DIMENSION N(32),S(32),R(32) DATA TC/190.551D0/,VC/162.66D0/ FA1=0.D0 FA2=0.D0 FA3=0.D0 DLTA=V/VC TAO=TC/T CALL COEFF(R,S,N) DO 111 I=1,13 FA1=FA1+S(I)*(S(I)-1.)*N(I)*DLTA**R(I)*TAO**S(I) 111 CONTINUE DO 222 I=14,24 FA2=FA2+S(I)*(S(I)-1.)*N(I)*DLTA**R(I)*TAO**S(I) 222 CONTINUE FA2=DEXP(-DLTA*DLTA)*FA2 DO 333 I=25,32 FA3=FA3+S(I)*(S(I)-1.D0)*N(I)*DLTA**R(I)*TAO**S(I) 333 CONTINUE FA3=DEXP(-DLTA**4)*FA3 R4CH4=FA1+FA2+FA3 R4CH4=R4CH4 RETURN END C -------------------------------------------------- C ** the Derivatives No.5 ** C ** of the Residual Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION R5CH4(V,T) DOUBLE PRECISION V,T,N,FA1,FA2,FA3,DLTA,TAO REAL S,R DIMENSION N(32),S(32),R(32) DATA TC/190.551D0/,VC/162.66D0/ FA1=0.D0 FA2=0.D0 FA3=0.D0 DLTA=V/VC TAO=TC/T CALL COEFF(R,S,N) DO 111 I=1,13 FA1=FA1+R(I)*S(I)*N(I)*DLTA**R(I)*TAO**S(I) 111 CONTINUE DO 222 I=14,24 FA2=FA2+S(I)*(R(I)-2.D0*DLTA*DLTA)* & N(I)*DLTA**R(I)*TAO**S(I) 222 CONTINUE FA2=DEXP(-DLTA*DLTA)*FA2 DO 333 I=25,32 FA3=FA3+S(I)*(R(I)-4.D0*DLTA**4)* & N(I)*DLTA**R(I)*TAO**S(I) 333 CONTINUE FA3=DEXP(-DLTA**4)*FA3 R5CH4=FA1+FA2+FA3 R5CH4=R5CH4 RETURN END C ---------------------------------------------------- C ** the Dimensionless Ideal Gas Term of Free Energy ** C ---------------------------------------------------- DOUBLE PRECISION FUNCTION I0CH4(V,T) DOUBLE PRECISION V,T,DLTA,TAO,Q DIMENSION Q(7) DATA (Q(I),I=1,7) /-10.413865D0, 2.5998324D0, & -3.3854083D0, 1.6900979D0, -0.3911541D0, 4.7206715D0, & -10.543907D0/ DATA TC/190.551D0/,VC/162.66D0/ DLTA=V/VC TAO=TC/T I0CH4=Q(1)+DLOG(DLTA)+Q(2)*DLOG(TAO)+Q(3)*TAO**(-1.D0/3.D0)+ & Q(4)*TAO**(-2.D0/3.D0)+Q(5)/TAO+Q(6)* & DLOG(1.-DEXP(Q(7)*TAO)) RETURN END C -------------------------------------------------- C ** the Derivatives No.1 ** C ** of the Ideal Gas Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION I1CH4(V,T) DOUBLE PRECISION V,T V=V T=T I1CH4=1.D0 RETURN END C -------------------------------------------------- C ** the Derivatives No.2 ** C ** of the Ideal Gas Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION I2CH4(V,T) DOUBLE PRECISION V,T,DLTA,TAO,Q DIMENSION Q(7) DATA (Q(I),I=1,7) /-10.413865D0, 2.5998324D0, & -3.3854083D0, 1.6900979D0, -0.3911541D0, 4.7206715D0, & -10.543907D0/ DATA TC/190.551D0/,VC/162.66D0/ DLTA=V/VC TAO=TC/T I2CH4=Q(2)- & Q(3)/3.D0*TAO**(-1.D0/3.D0)-2.D0/3.D0*Q(4)*TAO**(-2.D0/3.D0) & -Q(5)/TAO-Q(6)*Q(7)*TAO/(DEXP(-Q(7)*TAO)-1) RETURN END C -------------------------------------------------- C ** the Derivatives No.3 ** C ** of the Ideal Gas Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION I3CH4(V,T) DOUBLE PRECISION V,T V=V T=T I3CH4=-1.D0 RETURN END C -------------------------------------------------- C ** the Derivatives No.4 ** C ** of the Ideal Gas Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION I4CH4(V,T) DOUBLE PRECISION V,T,DLTA,TAO,Q DIMENSION Q(7) DATA (Q(I),I=1,7) /-10.413865D0, 2.5998324D0, & -3.3854083D0, 1.6900979D0, -0.3911541D0, 4.7206715D0, & -10.543907D0/ DATA TC/190.551D0/,VC/162.66D0/ DLTA=V/VC TAO=TC/T I4CH4=-Q(2)+4.D0*Q(3)/9.D0*TAO**(-1.D0/3.D0)+ & 10.D0*Q(4)/9.D0*TAO**(-2.D0/3.D0) & +2.D0*Q(5)/TAO-Q(6)*Q(7)*Q(7)*TAO*TAO* & DEXP(Q(7)*TAO)*(DEXP(Q(7)*TAO)-1)**(-2.) RETURN END C -------------------------------------------------- C ** the Derivatives No.5 ** C ** of the Ideal Gas Free Energy ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION I5CH4(V,T) DOUBLE PRECISION V,T V=V T=T I5CH4=0.D0 RETURN END C ---------------------------------------------------------- C ** F6 = ALMPD(P) ** C ---------------------------------------------------------- FUNCTION ALMPD(P) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'ALMPD'/ CALL S99CH4(FUN) P=P ALMPD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F7 = ALMPDD(P) ** C ---------------------------------------------------------- FUNCTION ALMPDD(P) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'ALMPDD'/ CALL S99CH4(FUN) P=P ALMPDD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F22 = EPSPT(P,T) ** C ---------------------------------------------------------- FUNCTION EPSPT(P,T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'EPSPT'/ CALL S99CH4(FUN) P=P T=T EPSPT=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F67 = TLDP(P) ** C ---------------------------------------------------------- FUNCTION TLDP(P) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'TLDP'/ CALL S99CH4(FUN) P=P TLDP=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F66 = PLDT(T) ** C ---------------------------------------------------------- FUNCTION PLDT(T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PLDT'/ CALL S99CH4(FUN) T=T PLDT=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F99 = PSBT(T) ** C ---------------------------------------------------------- FUNCTION PSBT(T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PSBT'/ CALL S99CH4(FUN) T=T PSBT=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F87 = PRTD(T) ** C ---------------------------------------------------------- REAL FUNCTION PRTD(T) REAL T,TI INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PRTD'/ TI=T CALL S99CH4(FUN) PRTD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F88 = PRTDD(T) ** C ---------------------------------------------------------- REAL FUNCTION PRTDD(T) REAL T,TI INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PRTDD'/ TI=T CALL S99CH4(FUN) PRTDD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F85 = PRPD(P) ** C ---------------------------------------------------------- REAL FUNCTION PRPD(P) REAL P,PI INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PRPD'/ PI=P CALL S99CH4(FUN) PRPD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F86 = PRPDD(P) ** C ---------------------------------------------------------- REAL FUNCTION PRPDD(P) REAL P,PI INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PRPDD'/ PI=P CALL S99CH4(FUN) PRPDD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F72 = PSTD(T) ** C ---------------------------------------------------------- FUNCTION PSTD(T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PSTD'/ CALL S99CH4(FUN) T=T PSTD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F73 = PSTDD(T) ** C ---------------------------------------------------------- FUNCTION PSTDD(T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'PSTDD'/ CALL S99CH4(FUN) T=T PSTDD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F74 = TSPD(P) ** C ---------------------------------------------------------- FUNCTION TSPD(P) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'TSPD'/ CALL S99CH4(FUN) P=P TSPD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F75 = TSPDD(P) ** C ---------------------------------------------------------- FUNCTION TSPDD(P) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'TSPDD'/ CALL S99CH4(FUN) P=P TSPDD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F80 = VPS(P,S) ** C ---------------------------------------------------------- FUNCTION VPS(P,S) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'VPS'/ CALL S99CH4(FUN) P=P S=S VPS=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F79 = UPS(P,S) ** C ---------------------------------------------------------- FUNCTION UPS(P,S) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'UPS'/ CALL S99CH4(FUN) P=P S=S UPS=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F9 = ALMTD(T) ** C ---------------------------------------------------------- FUNCTION ALMTD(T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'ALMTD'/ CALL S99CH4(FUN) T=T ALMTD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F10 = ALMTDD(T) ** C ---------------------------------------------------------- FUNCTION ALMTDD(T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'ALMTDD'/ CALL S99CH4(FUN) T=T ALMTDD=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F1 = AIPPT(P,T) ** C ---------------------------------------------------------- FUNCTION AIPPT(P,T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'AIPPT'/ CALL S99CH4(FUN) P=P T=T AIPPT=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F2 = ALAPP(P) ** C ---------------------------------------------------------- FUNCTION ALAPP(P) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'ALAPP'/ P=P CALL S99CH4(FUN) ALAPP=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** F3 = ALAPP(T) ** C ---------------------------------------------------------- FUNCTION ALAPT(T) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'ALAPT'/ T=T CALL S99CH4(FUN) ALAPT=-1.0E+30 RETURN END C ----------------------------------------------- C ** F100 = TSBP ** C ----------------------------------------------- FUNCTION TSBP(P) INTEGER KPA COMMON/UNIT/KPA,MESS CHARACTER FUN*6 DATA FUN/'TSBP'/ CALL S99CH4(FUN) P=P TSBP=-1.0E+30 RETURN END C ---------------------------------------------------------- C ** Conversion of Temperature Scale ** C ** from IPTS-68 to ITS-90 ** C ** Range : -200C<=T68<=630C and T68>=1064C ** C ---------------------------------------------------------- REAL FUNCTION T90(T) REAL T68,T0K REAL*8 CA(8) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS DATA CA/-0.14542, -0.26722, 1.06471, 1.131286, -4.25835, & -1.59924, 7.29176, -3.52573/ IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T0K=0.0 ELSE T0K=273.15 END IF T68=T-T0K IF((-200.0.LE.T68).AND.(T68.LE.630.0)) THEN RT68=T68/630.0 TEMP=T68+RT68*(CA(1)+RT68*(CA(2)+RT68*(CA(3)+RT68*(CA(4) & +RT68*(CA(5)+RT68*(CA(6)+RT68*(CA(7)+RT68*CA(8)))))))) ELSE IF((1064.0.LE.T68).AND.(T68.LE.3000.0)) THEN TEMP=T68-0.25*(T68+273.15)*(T68+273.15)/1337.33/1337.33 ELSE IF(MESS.NE.0) WRITE(6,6020) T 6020 FORMAT(1H ,5X,'**** OUT OF RANGE AT T90', & ' WHEN T68 =',1PE14.7,' ****') T90=-1.0E+20 RETURN END IF T90=TEMP+T0K RETURN END C ---------------------------------------------------------- C ** Conversion of Temperature Scale ** C ** from ITS-90 to IPTS-68 ** C ** Range : -200C<=T68<=630C and T68>=1064C ** C ---------------------------------------------------------- REAL FUNCTION T68(T) T1=T90(T) T2=T ITER=0 DT=T-T1 1000 CONTINUE ITER=ITER+1 IF(ITER.GE.10000) THEN WRITE(6,*) ' ITERATION EXCEEDED MAXIMUM ' T68=-1.0E+10 RETURN END IF T2=T2+DT T1=T90(T2) DT=T-T1 IF(ABS(DT).GE.1.0E-4) GOTO 1000 T68=T2 RETURN END C ---------------------------------------------------------- C ** SUBROUTINE PROGRAM SPECIFYING KPA AND MESS ** C ** MS-FORTRAN and MS-C Mixed Langage Programing ** C ---------------------------------------------------------- SUBROUTINE KPAMES(KPAC,MESSC) COMMON/UNIT/KPA,MESS INTEGER KPAC INTEGER MESSC KPA=KPAC MESS=MESSC RETURN END C ------------------------------------------------------------- C * * FUNCTION FOR SETTING UNITS (PRESSURE) * * C ------------------------------------------------------------- REAL FUNCTION G98CH4(KPA,P) REAL P,PBAR IF(KPA.EQ.1) THEN PBAR=1.0E0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0E0 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 ELSE PBAR=1.0E-05 END IF G98CH4=P*PBAR RETURN END C ------------------------------------------------------------- C * * FUNCTION FOR SETTING UNITS (PRESSURE) * * C ------------------------------------------------------------- REAL FUNCTION G99CH4(KPA,T) REAL T,T0K IF(KPA.EQ.1) THEN T0K=0.0E0 ELSE IF(KPA.EQ.2) THEN T0K=273.15E0 ELSE IF(KPA.EQ.3) THEN T0K=0.0E0 ELSE T0K=273.15E0 END IF G99CH4=T-T0K RETURN END C ------------------------------------------------------------- C ** SUBROUTINE FOR ERROR MESSAGE * C ------------------------------------------------------------- C ** LEVEL 1 ERROR MESSAGE ** C ------------------------------------------------------------- SUBROUTINE S97CH4(FUN) COMMON/UNIT/KPA,MESS CHARACTER FUN*6,MSG*125 IF(MESS.NE.0) THEN MSG='**** NO CONVERGENCE AT '//FUN//' FOR METHANE ****' WRITE(6,*) MSG END IF RETURN END C ------------------------------------------------------------- C ** LEVEL 2 ERROR MESSAGE ** C ------------------------------------------------------------- SUBROUTINE S98CH4(IPT,P,T,N1,N2,FUN) COMMON/UNIT/KPA,MESS CHARACTER FUN*6,N1*1,N2*1 IF(MESS.NE.0) THEN IF (IPT.EQ.1) THEN WRITE(6,6010) FUN,P 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR METHANE', & ' WHEN P =',1PE14.7,' ****') ELSE IF (IPT.EQ.2) THEN WRITE(6,6020) FUN,T 6020 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6, & ' FOR METHANE',' WHEN T =',1PE14.7,' ****') ELSE WRITE(6,6030) FUN,N1,P,N2,T 6030 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6, & ' FOR METHANE',' WHEN ',A1,' =',1PE14.7, & ' AND ',A1,' =',1PE14.7,' ****') END IF END IF RETURN END C ------------------------------------------------------------- C ** LEVEL 3 ERROR MESSAGE ** C ------------------------------------------------------------- SUBROUTINE S99CH4(FUN) COMMON/UNIT/KPA,MESS CHARACTER FUN*6,MSG*125 IF(MESS.NE.0) THEN MSG='**** FUNCTION '//FUN//' UNAVAILABLE FOR METHANE ****' WRITE(6,100) MSG 100 FORMAT(1H ,5X,A) END IF RETURN END C ---------------------------------------------------------- C ** F21 FUNCTION CRP('A') ** C ** CRP: Critical Constants of Methane ** C ---------------------------------------------------------- REAL FUNCTION CRP(A) CHARACTER FUN*6,FLUID*8,A*1 REAL T0K,PBAR,FF INTEGER KPA DOUBLE PRECISION F21CH4 COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'CRP'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF FF = F21CH4(A) IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,A 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN A =',A,' ****') FF=-1.0E+20 END IF END IF IF(A.EQ.'T') THEN IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 FF=FF+T0K ELSE IF(A.EQ.'P') THEN IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) PBAR=1.0 FF=FF/PBAR END IF CRP=FF RETURN END C ------------------------------------------------ C ** CRP(A) QUANTITIES AT THE CRITICAL POINT ** C ------------------------------------------------ DOUBLE PRECISION FUNCTION F21CH4(A) CHARACTER*1 A,B(5) DATA B/'H','P','S','T','V'/ IF((A.EQ.B(1)).OR.(A.EQ.B(2)).OR.(A.EQ.B(3)).OR. & (A.EQ.B(4)).OR.(A.EQ.B(5))) GO TO 5 GO TO 900 5 IF(A.NE.B(1)) GO TO 10 F21CH4=132.3319D3 RETURN 10 IF(A.NE.B(2)) GO TO 20 F21CH4=45.9920D00 RETURN 20 IF(A.NE.B(3)) GO TO 30 F21CH4=7.511064D03 RETURN 30 IF(A.NE.B(4)) GO TO 40 F21CH4=-82.599D00 RETURN 40 IF(A.NE.B(5)) GO TO 900 F21CH4=1.D00/162.66D00 RETURN 900 F21CH4=-1.0E+20 RETURN END C -------------------------------------------------- C ** F89 FUNCTION FC('A') ** C ** FUNCTION FOR FUNDDAMENTAL CONSTANTS ** C -------------------------------------------------- REAL FUNCTION FC(A) CHARACTER A*1,MSG*120 COMMON/UNIT/KPA,MESS IF(A.EQ.'M'.OR.A.EQ.'m') THEN FC=16.043 ELSE IF(A.EQ.'R'.OR.A.EQ.'r') THEN FC=518.26 ELSE FC=-1.0E+20 IF(MESS.NE.0) THEN MSG='**** OUT OF RANGE AT FC FOR METHANE & WHEN A=''' //A//''' ****' WRITE(6,'(1H ,A)') MSG ENDIF ENDIF RETURN END C -------------------------------------------------- C ** FUNCTION FOR IDENTIFICATION OF SUBSTANCE ** C ** F84 USAGE: IDENTF(A) ** C ** A, B : CHARACTER TYPE VALIABLES ** C ** B='METHANE' WHEN A='S' ** C ** B='CH4' WHEN A='C' ** C ** B='12.1' WHEN A='V' ** C -------------------------------------------------- CHARACTER*20 FUNCTION IDENTF(A) CHARACTER A*1,MSG*120 COMMON/UNIT/KPA,MESS IF (A.EQ.'S') THEN IDENTF='METHANE' ELSE IF (A.EQ.'C') THEN IDENTF='CH4' ELSE IF (A.EQ.'V') THEN IDENTF='12.1' ELSE IDENTF='????????????????????' IF (MESS.NE.0) THEN MSG='**** OUT OF RANGE AT IDENTF FOR METHANE WHEN A=''' & //A//''' ****' WRITE(6,'(1H ,A)') MSG END IF END IF RETURN END C ------------------------------------------- C ** F41 FUNCTION TRPL('A') ** C ** FUNCTION OF TRIPLE POINT OF METHANE ** C ** A=P : TRIPLE PRESSURE ** C ** A=T : TRIPLE TEMPERATURE ** C ------------------------------------------- REAL FUNCTION TRPL(A) CHARACTER FUN*6,FLUID*8,A*1 REAL T0K,PBAR,FF INTEGER KPA DOUBLE PRECISION F41CH4 COMMON/UNIT/KPA,MESS DATA FLUID/'METHANE'/, FUN/'TRPL'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF FF = F41CH4(A) IF(FF.EQ.-1.0E+20) THEN IF(MESS.NE.0) THEN WRITE(6,6010) FUN,FLUID,A 6010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ',A, & ' WHEN A =',A,' ****') FF=-1.0E+20 END IF END IF IF(A.EQ.'T') THEN IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 FF=FF+T0K ELSE IF(A.EQ.'P') THEN IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) PBAR=1.0 FF=FF/PBAR END IF TRPL=FF RETURN END C ----------------------------------------------- C ** TRPL(A) QUANTITIES AT THE TRIPLE POINT ** C ----------------------------------------------- DOUBLE PRECISION FUNCTION F41CH4(A) CHARACTER*1 A,B(2) DATA B/'P','T'/ IF (A.EQ.B(1).OR.A.EQ.B(2)) GO TO 5 GO TO 900 5 IF(A.NE.B(1)) GO TO 10 F41CH4=0.11696D00 RETURN 10 IF(A.NE.B(2)) GO TO 900 F41CH4=182.4646D00 RETURN 900 F41CH4=-1.0E+20 RETURN END C -------------------------------------------------- C ** F68 = PMLT ** C ** Pressure on Melting Curve [Pa],[bar] ** C ** Temperature [K],[C] ** C -------------------------------------------------- REAL FUNCTION PMLT(T) CHARACTER FUN*6 REAL PBAR,T0K,T,FF INTEGER KPA DOUBLE PRECISION F68CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'PMLT'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF TI=T-T0K DBT=DBLE(TI) FF = F68CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) PBAR=1.0 PMLT=FF/PBAR RETURN END C -------------------------------------------------- C ** F68 = PMLT ** C ** Pressure on Melting Curve [Pa],[bar] ** C ** Temperature [K],[C] ** C -------------------------------------------------- REAL FUNCTION TMLP(P) CHARACTER FUN*6 REAL PBAR,T0K,P,FF INTEGER KPA DOUBLE PRECISION F69CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'TMLP'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF PI=P*PBAR DBP=DBLE(PI) FF = F69CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 TMLP=FF+T0K RETURN END C ----------------------------------------------- C *** F68CH4 PRESSURE ON THE MELTING LINE *** C ----------------------------------------------- DOUBLE PRECISION FUNCTION F68CH4(T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DIMENSION D(5) DATA TC/260.D00/,TL/90.6854D00/,PC/10416.D00/,PL/0.11696D00/ DATA D/-32.2988083D00,163.1572171D00, & -238.0692255D00,156.5940578D00,-38.78155184D00/,DLT/1.D-05/ TM=T+273.15D00 IF(DABS((TM-TC)/TC).LT.DLT) GO TO 300 IF(DABS((TM-TL)/TL).LT.DLT) GO TO 350 IF(TM.LT.TL) GO TO 600 IF(TM.GT.TC) GO TO 600 P1=D(1)*(TM/TL-1.)**0.1+D(2)*(TM/TL-1.)**0.2 P2=D(3)*(TM/TL-1.)**0.3+ & D(4)*(TM/TL-1.)**0.4+D(5)*(TM/TL-1.)**0.5 P=DLOG(PL)+(P1+P2) F68CH4=DEXP(P) RETURN 300 F68CH4=PC RETURN 350 F68CH4=PL RETURN 600 F68CH4=-1.0E+20 RETURN END C -------------------------------------------------- C *** F69CH4 TEMPERATURE ON THE MELTING LINE *** C -------------------------------------------------- DOUBLE PRECISION FUNCTION F69CH4(P) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DIMENSION D(5) DATA TC/260.D00/,TL/90.6854D00/,PC/10416.D00/,PL/0.11696D00/, & EV/-1.0D+20/,D/-32.2988083D00,163.157217D00, & -238.0692255D00,156.5940578D00,-38.78155184D00/,DLT/1.D-05/ IF(P.LT.PL) GO TO 600 IF(P.GT.PC) GO TO 600 IF(DABS((P-PC)/PC).LT.DLT) GO TO 300 IF(DABS((P-PL)/PL).LT.DLT) GO TO 350 ERR=1.E-7 X0=90.6854D00 X1=260.D00 X3=X1 IC=0 F3=0. F4=0. DO 50 K=1,5 MR=0.1D00*DBLE(K) F3=F3+D(K)*(X0/TL-1.D00)**MR 50 F4=F4+D(K)*(X1/TL-1.D00)**MR F0=DLOG(P/PL)-F3 F1=DLOG(P/PL)-F4 100 IC=IC+1 X2=(X0+X1)/2.D00 F5=0. DO 110 K=1,5 MR=0.1D00*DBLE(K) 110 F5=F5+D(K)*(X2/TL-1.D00)**MR F2=DLOG(P/PL)-F5 IF(DABS(X2-X3).LE.ERR) GO TO 400 X3=X2 IF(F0*F2.LT.0.) GO TO 200 F0=F2 X0=X2 IF(IC.GT.200) GO TO 900 GO TO 100 200 F1=F2 X1=X2 IF(IC.GT.200) GO TO 900 GO TO 100 300 F69CH4=TC-273.15D00 RETURN 350 F69CH4=TL-273.15D00 RETURN 400 F69CH4=X2-273.15D00 RETURN 600 F69CH4=-1.0E+20 RETURN 900 F69CH4=EV RETURN END C -------------------------------------------------- C ** F30 = PST(T) ** C ** Saturation Pressure [Pa], [Bar] ** C ** INPUT : Temperature [K], [C] ** C -------------------------------------------------- REAL FUNCTION PST(T) CHARACTER FUN*6 REAL PBAR,T0K,T,FF INTEGER KPA DOUBLE PRECISION F30CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'PST'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0e-5 T0K=0.0 ELSE PBAR=1.0e-5 T0K=273.15 END IF TI=T-T0K DBT=DBLE(TI) FF = F30CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) PBAR=1.0 PST=FF/PBAR RETURN END C -------------------------------------------------------- C ** F30CH4 SATURATION PRESSURE ** C -------------------------------------------------------- DOUBLE PRECISION FUNCTION F30CH4(T) DOUBLE PRECISION TC,PC,TL,PL,TM,TT,H(5),E,F,DLT,P,T DATA TC/190.551/,PC/45.993D0/,TL/90.6854/,PL/0.11696D0/ DATA DLT/1.0D-5/ DATA (H(I),I=1,5) /-6.589879, 0.6355175, 11.31028, -10.38720, & 3.393075/, E/1.90/ TM=T+273.15D00 IF(DABS((TM-TC)/TC).LT.DLT) GO TO 300 IF(DABS((TM-TL)/TL).LT.DLT) GO TO 350 IF(TM.LT.TL) GO TO 600 IF(TM.GT.TC) GO TO 600 TT=1-TM/TC F=H(1)*TT/(1-TT)+H(2)*TT+H(3)*TT**E+H(4)*TT*TT & +H(5)*TT*TT*TT P=PC*DEXP(F) F30CH4=P RETURN 300 F30CH4=PC RETURN 350 F30CH4=PL RETURN 600 F30CH4=-1.0E+20 RETURN END C -------------------------------------------------- C ** F40 = TSP(P) ** C ** Saturation Temperature [K], [C] ** C ** INPUT : PRESSURE [Pa], [bar] ** C -------------------------------------------------- REAL FUNCTION TSP(P) CHARACTER FUN*6 REAL PBAR,T0K,P,FF INTEGER KPA DOUBLE PRECISION F40CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'TSP'/ IF(KPA.EQ.1) THEN PBAR=1.0 T0K=0.0 ELSE IF(KPA.EQ.2) THEN PBAR=1.0 T0K=273.15 ELSE IF(KPA.EQ.3) THEN PBAR=1.0E-05 T0K=0.0 ELSE PBAR=1.0E-05 T0K=273.15 END IF PI=P*PBAR DBP=DBLE(PI) FF = F40CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF IF((FF.EQ.-1.0E+10).OR.(FF.EQ.-1.0E+20)) T0K=0.0 TSP=FF+T0K RETURN END C -------------------------------------------------- C ** F40CH4(P) SATURATION TEMPERATURE ** C -------------------------------------------------- DOUBLE PRECISION FUNCTION F40CH4(P) IMPLICIT DOUBLE PRECISION (A-H,L-Z) DIMENSION H(5) DATA TC/190.551/,PC/45.993D0/,TL/90.6854/,PL/0.11696D0/ DATA DLT/1.0D-5/,ERR/1.0E-5/ DATA (H(I),I=1,5) /-6.589879, 0.6355175, 11.31028, -10.38720, & 3.393075/, E/1.90/,EV/-1.0E+20/ AA(X)=P/PC- & DEXP(H(1)*S/(1-S)+H(2)*S+H(3)*S**E+H(4)*S*S+H(5)*S*S*S) IF(DABS((P-PC)/PC).LT.DLT) GO TO 300 IF(DABS((P-PL)/PL).LT.DLT) GO TO 350 IF(P.LT.PL) GO TO 600 IF(P.GT.PC) GO TO 600 X0=TL X1=TC X3=X1 IC=0 S=1.D00-X0/TC F0=AA(X0) F1=AA(X1) 100 IC=IC+1 X2=(X0+X1)/2.D00 S=1.D00-X2/TC F2=AA(X2) IF(DABS(X2-X3).LE.ERR) GO TO 400 X3=X2 IF(F0*F2.LT.0.) GO TO 200 F0=F2 X0=X2 IF(IC.GT.200) GO TO 900 GO TO 100 200 F1=F2 X1=X2 IF(IC.GT.200) GO TO 900 GO TO 100 300 F40CH4=TC-273.15D00 RETURN 350 F40CH4=TL-273.15D00 RETURN 400 F40CH4=X2-273.15D00 RETURN 600 F40CH4=-1.0E+20 RETURN 900 F40CH4=EV RETURN END C ------------------------------------------------------ C ** F49 = VPD ** C ** Specific Volume of Saturation Liquid [m**3/kg] ** C ** INPUT Pressure [Pa],[bar] ** C ------------------------------------------------------ REAL FUNCTION VPD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F49CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'VPD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F49CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF VPD=FF RETURN END C ------------------------------------------------------ C ** F50 = VPDD ** C ** Specific Volume of Saturation Vapor [m**3/kg] ** C ** INPUT Pressure [Pa],[bar] ** C ------------------------------------------------------ REAL FUNCTION VPDD(P) CHARACTER FUN*6 REAL P,FF INTEGER KPA DOUBLE PRECISION F50CH4,DBP COMMON/UNIT/KPA,MESS DATA FUN/'VPDD'/ PI=G98CH4(KPA,P) DBP=DBLE(PI) FF = F50CH4(DBP) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(1,P,T,'P','T',FUN) FF=-1.0E+20 END IF VPDD=FF RETURN END C ---------------------------------------------------- C ** F49CH4(P) SPECIFIC VOLUME OF SATURATED LIQUID ** C ---------------------------------------------------- DOUBLE PRECISION FUNCTION F49CH4(P) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF(P.LT.0.11719D00.OR.P.GT.45.95001D00) GO TO 900 T=F40CH4(P) IF(T.EQ.-1.0E+20) GO TO 900 F49CH4=F53CH4(T) RETURN 900 F49CH4=-1.0E+20 RETURN END C ---------------------------------------------------- C ** F50CH4(P) SPECIFIC VOLUME OF SATURATED VAPOR ** C ---------------------------------------------------- DOUBLE PRECISION FUNCTION F50CH4(P) IMPLICIT DOUBLE PRECISION(A-H,L-Z) IF(P.LT.0.11719D00.OR.P.GT.45.95001D00) GO TO 900 T=F40CH4(P) IF(T.EQ.-1.0E+20) GO TO 900 F50CH4=F54CH4(T) RETURN 900 F50CH4=-1.0E+20 RETURN END C ------------------------------------------------------ C ** F53CH4(T) SPECIFIC VOLUME OF SATURATED LIQUID ** C ------------------------------------------------------ DOUBLE PRECISION FUNCTION F53CH4(T) IMPLICIT DOUBLE PRECISION(A-H,J-Z) DIMENSION G(4) DATA TC/190.551D00/,LC/162.66D00/,TL/90.6854D00/, & BELTA/0.355/,LF/451.53/,DLT/1.0D-5/ DATA (G(I),I=1,4)/1.838982,-0.7727452,0.5592446,-0.3807793/ TM=T+273.15D00 IF(DABS((TM-TC)/TC).LT.DLT) GO TO 300 IF(DABS((TM-TL)/TL).LT.DLT) GO TO 350 IF(TM.LT.TL) GO TO 600 IF(TM.GT.TC) GO TO 600 P=F30CH4(T) TT=1-TM/TC VL=LC*(1+(G(1)*TT**BELTA+G(2)*TT*TT+G(3)*TT*TT*TT) & /(1+G(4)*TT**(1-BELTA))) F53CH4=1/VL RETURN 300 F53CH4=1.D00/LC RETURN 350 F53CH4=1.D00/LF RETURN 600 F53CH4=-1.0E+20 RETURN END C ----------------------------------------------------- C ** F54CH4(T) SPECIFIC VOLUME OF SATURATED VAPOR ** C ----------------------------------------------------- DOUBLE PRECISION FUNCTION F54CH4(T) IMPLICIT DOUBLE PRECISION(A-H,J-Z) DIMENSION J(5) DATA(J(I),I=1,5) /-0.7377483,-1.241532,-1.649972,2.281949, & 1.439570/,BELTA/0.355/,DLT/1.0D-5/,PC/45.992/ DATA ZC/0.28631/,TC/190.551D0/,LC/162.66/,LL/251.23/,TL/90.685/ TM=T+273.15D00 IF(DABS((TM-TL)/TL).LT.DLT) GO TO 300 IF(DABS((TM-TC)/TC).LT.DLT) GO TO 350 IF(TM.LT.TL) GO TO 600 IF(TM.GT.TC) GO TO 600 P=F30CH4(T) PP=P/PC TT=1-TM/TC VV=J(1)*TT**BELTA+J(2)*TT**(2*BELTA)+J(3)*(TT+TT**4)+J(4)*TT*TT VV=(1.-1./ZC)*VV/(1+J(5)*TT)+1 VV=VV-1/ZC*(1-(1-TT)**8/PP) VV=LC*(1-TT)**7/VV F54CH4=1.D0/VV RETURN 300 F54CH4=1.D00/LL RETURN 350 F54CH4=1.D00/LC RETURN 600 F54CH4=-1.0E+20 RETURN END C ----------------------------------------------------- C ** F53 = VTD ** C ** Specific Volume of Saturated Liquid [m**3/Kg] ** C ** INPUT T : Temperature [K],[C] ** C ----------------------------------------------------- REAL FUNCTION VTD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F53CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'VTD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F53CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF VTD=FF RETURN END C ----------------------------------------------------- C ** F54 = VTDD ** C ** Specific Volume of Saturated Vapor [m**3/Kg] ** C ** INPUT T : Temperature [K],[C] ** C ----------------------------------------------------- REAL FUNCTION VTDD(T) CHARACTER FUN*6 REAL T,FF INTEGER KPA DOUBLE PRECISION F54CH4,DBT COMMON/UNIT/KPA,MESS DATA FUN/'VTDD'/ TI=G99CH4(KPA,T) DBT=DBLE(TI) FF = F54CH4(DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(2,P,T,'P','T',FUN) FF=-1.0E+20 END IF VTDD=FF RETURN END C ------------------------------------------------------- C ** F81= PRPT(P,T) : Prandtl Number [-] ** C ** INPUT : ** C ** P : Pressure [Pa],[bar] ** C ** T : Temperature [K],[C] ** C ------------------------------------------------------- REAL FUNCTION PRPT(P,T) CHARACTER FUN*6 REAL P,T,FF INTEGER KPA DOUBLE PRECISION F81CH4,DBP,DBT COMMON/UNIT/KPA,MESS DATA FUN/'PRPT'/ PI=G98CH4(KPA,P) TI=G99CH4(KPA,T) DBP=DBLE(PI) DBT=DBLE(TI) FF = F81CH4(DBP,DBT) IF(FF.EQ.-1.0E+10) THEN CALL S97CH4(FUN) FF=-1.0E+10 ELSE IF(FF.EQ.-1.0E+20) THEN CALL S98CH4(3,P,T,'P','T',FUN) FF=-1.0E+20 END IF PRPT=FF RETURN END C ------------------------------------------------- C ** F81CH4(P,T) , PRPT : PRANDTL NUMBER ** C ------------------------------------------------- DOUBLE PRECISION FUNCTION F81CH4(P,T) IMPLICIT DOUBLE PRECISION (A-H,L-Z) TA=T+273.15D00 IF(TA.LT.298.15D00.OR.TA.GT.423.15D00) GO TO 900 IF(TA.GE.298.15D00.AND.TA.LE.423.15D00 & .AND.P.GE.1.D00.AND.P.LE.500.D00) GO TO 100 GO TO 900 100 CP=F18CH4(P,T) IF(CP.EQ.-1.0E+20) GO TO 900 MU=F13CH4(P,T) IF(MU.EQ.-1.0E+20) GO TO 900 LM=F8CH4(P,T) IF(LM.EQ.-1.0E+20) GO TO 900 F81CH4=(CP*MU)/LM RETURN 900 F81CH4=-1.0E+20 RETURN END