*=====PC2H6V91====1994.12.15===========================================* *=====PAC2H6======1997. 3.10===========================================* 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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(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 S99C06(FUN) WTDD=-1.0E+30 RETURN END C C==================================================================== *------------------------------------------------- F1C06 = AIPPT REAL FUNCTION AIPPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C06('AIPPT') AIPPT=-1.0E+30 RETURN END *------------------------------------------------- F2C06 = ALAPP REAL FUNCTION ALAPP(P) REAL P,PI PI=P CALL S99C06('ALAPP') ALAPP=-1.0E+30 RETURN END *------------------------------------------------- F3C06 = ALAPT REAL FUNCTION ALAPT(T) REAL T,TI TI=T CALL S99C06('ALAPT') ALAPT=-1.0E+30 RETURN END *------------------------------------------------- F4C06 = ALHP REAL FUNCTION ALHP(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) ALHP=F4C06(PI) IF(ALHP.EQ.-1.0E+10) THEN CALL S97C06('ALHP') ELSE IF(ALHP.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','ALHP') END IF RETURN END *------------------------------------------------- F5C06 = ALHT REAL FUNCTION ALHT(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) ALHT=F5C06(TI) IF(ALHT.EQ.-1.0E+10) THEN CALL S97C06('ALHT') ELSE IF(ALHT.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','ALHT') END IF RETURN END *------------------------------------------------- F6C06 = ALMPD REAL FUNCTION ALMPD(P) REAL P,PI PI=P CALL S99C06('ALMPD') ALMPD=-1.0E+30 RETURN END *------------------------------------------------- F7C06 = ALMPDD REAL FUNCTION ALMPDD(P) REAL P,PI PI=P CALL S99C06('ALMPDD') ALMPDD=-1.0E+30 RETURN END *------------------------------------------------- F8C06 = ALMPT REAL FUNCTION ALMPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C06('ALMPT') ALMPT=-1.0E+30 RETURN END *------------------------------------------------- F9C06 = ALMTD REAL FUNCTION ALMTD(T) REAL T,TI TI=T CALL S99C06('ALMTD') ALMTD=-1.0E+30 RETURN END *------------------------------------------------- F10C06 = ALMTDD REAL FUNCTION ALMTDD(T) REAL T,TI TI=T CALL S99C06('ALMTDD') ALMTDD=-1.0E+30 RETURN END *------------------------------------------------- F11C06 = AMUPD REAL FUNCTION AMUPD(P) REAL P,PI PI=P CALL S99C06('AMUPD') AMUPD=-1.0E+30 RETURN END *------------------------------------------------- F12C06 = AMUPDD REAL FUNCTION AMUPDD(P) REAL P,PI PI=P CALL S99C06('AMUPDD') AMUPDD=-1.0E+30 RETURN END *------------------------------------------------- F13C06 = AMUPT REAL FUNCTION AMUPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C06('AMUPT') AMUPT=-1.0E+30 RETURN END *------------------------------------------------- F14C06 = AMUTD REAL FUNCTION AMUTD(T) REAL T,TI TI=T CALL S99C06('AMUTD') AMUTD=-1.0E+30 RETURN END *------------------------------------------------- F15C06 = AMUTDD REAL FUNCTION AMUTDD(T) REAL T,TI TI=T CALL S99C06('AMUTDD') AMUTDD=-1.0E+30 RETURN END *------------------------------------------------- F16C06 = CPPD REAL FUNCTION CPPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) CPPD=F16C06(PI) IF(CPPD.EQ.-1.0E+10) THEN CALL S97C06('CPPD') ELSE IF(CPPD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','CPPD') END IF RETURN END *------------------------------------------------- F17C06 = CPPDD REAL FUNCTION CPPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) CPPDD=F17C06(PI) IF(CPPDD.EQ.-1.0E+10) THEN CALL S97C06('CPPDD') ELSE IF(CPPDD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','CPPDD') END IF RETURN END *------------------------------------------------- F18C06 = CPPT REAL FUNCTION CPPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) CPPT=F18C06(PI,TI) IF(CPPT.EQ.-1.0E+10) THEN CALL S97C06('CPPT') ELSE IF(CPPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','CPPT') END IF RETURN END *------------------------------------------------- F19C06 = CPTD REAL FUNCTION CPTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) CPTD=F19C06(TI) IF(CPTD.EQ.-1.0E+10) THEN CALL S97C06('CPTD') ELSE IF(CPTD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','CPTD') END IF RETURN END *------------------------------------------------- F20C06 = CPTDD REAL FUNCTION CPTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) CPTDD=F20C06(TI) IF(CPTDD.EQ.-1.0E+10) THEN CALL S97C06('CPTDD') ELSE IF(CPTDD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','CPTDD') END IF RETURN END *------------------------------------------------- F21C06 = CRP REAL FUNCTION CRP(A) CHARACTER A*1 REAL PBAR,T0K,FF INTEGER KPA,MESS COMMON/UNIT/KPA,MESS 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=F21C06(A) IF(FF.EQ.-1.0E+20) THEN IF (MESS.NE.0) THEN WRITE(6,2000) A 2000 FORMAT(1H ,5X,'**** OUT OF RANGE AT CRP FOR ETHANE', - ' WHEN A =',A,' ****') END IF FF=-1.0E+20 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 *------------------------------------------------- F22C06 = EPSPT REAL FUNCTION EPSPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C06('EPSPT') EPSPT=-1.0E+30 RETURN END C------------------------------------------------- F89 = FC C************************************************ C FUNCTION FOR FUNDDAMENTAL CONSTANTS C PROPATH VER.7.1, MAY 8, 1990 C USAGE: B=FC(A) C A, B : CHARACTER TYPE VALIABLES C B='30.0694' WHEN A='M' C B='276.507' WHEN A='R' C************************************************ REAL FUNCTION FC(A) CHARACTER A*1,MSG*120 COMMON/UNIT/KPA,MESS IF (A.EQ.'M') THEN FC=30.0694 ELSE IF (A.EQ.'R') THEN FC=276.507 ELSE FC=-1.E+20 IF (MESS.NE.0) THEN MSG='**** OUT OF RANGE AT FC FOR ETHANE WHEN A=''' & //A//''' ****' WRITE(6,'(1H ,A)') MSG END IF END IF RETURN END *------------------------------------------------- F23C06 = HPD REAL FUNCTION HPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) HPD=F23C06(PI) IF(HPD.EQ.-1.0E+10) THEN CALL S97C06('HPD') ELSE IF(HPD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','HPD') END IF RETURN END *------------------------------------------------- F24C06 = HPDD REAL FUNCTION HPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) HPDD=F24C06(PI) IF(HPDD.EQ.-1.0E+10) THEN CALL S97C06('HPDD') ELSE IF(HPDD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','HPDD') END IF RETURN END *------------------------------------------------- F25C06 = HPT REAL FUNCTION HPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) HPT=F25C06(PI,TI) IF(HPT.EQ.-1.0E+10) THEN CALL S97C06('HPT') ELSE IF(HPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','HPT') END IF RETURN END *------------------------------------------------- F26C06 = HPX REAL FUNCTION HPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) HPX=F26C06(PI,X) IF(HPX.EQ.-1.0E+10) THEN CALL S97C06('HPX') ELSE IF((HPX.EQ.-1.0E+20).AND.(MESS.NE.0)) THEN CALL S98C06(2,P,X,'P','X','HPX') END IF RETURN END *------------------------------------------------- F27C06 = HTD REAL FUNCTION HTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) HTD=F27C06(TI) IF(HTD.EQ.-1.0E+10) THEN CALL S97C06('HTD') ELSE IF(HTD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','HTD') END IF RETURN END *------------------------------------------------- F28C06 = HTDD REAL FUNCTION HTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) HTDD=F28C06(TI) IF(HTDD.EQ.-1.0E+10) THEN CALL S97C06('HTDD') ELSE IF(HTDD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','HTDD') END IF RETURN END *------------------------------------------------- F29C06 = HTX REAL FUNCTION HTX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) HTX=F29C06(TI,X) IF(HTX.EQ.-1.0E+10) THEN CALL S97C06('HTX') ELSE IF((HTX.EQ.-1.0E+20).AND.(MESS.NE.0)) THEN CALL S98C06(2,T,X,'T','X','HTX') END IF RETURN END *------------------------------------------------- F84C06 = IDENTF C************************************************ C FUNCTION FOR IDENTIFICATION OF SUBSTANCE C PROPATH VER.12.1, MAY 2, 2001 C USAGE: B=IDENTF(A) C A, B : CHARACTER TYPE VALIABLES C B='ETHANE' WHEN A='S' C B='C2H6' 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='ETHANE' ELSE IF (A.EQ.'C') THEN IDENTF='C2H6' ELSE IF (A.EQ.'V') THEN IDENTF='12.1' ELSE IDENTF='????????????????????' IF (MESS.NE.0) THEN MSG='**** OUT OF RANGE AT IDENTF FOR ETHANE WHEN A=''' & //A//''' ****' WRITE(6,'(1H ,A)') MSG END IF END IF RETURN END *------------------------------------------------- F30C06 = PST REAL FUNCTION PST(T) REAL TI,T INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) PST=F30C06(TI) IF(PST.EQ.-1.0E+10) THEN CALL S97C06('PST') RETURN ELSE IF(PST.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','PST') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.2)) RETURN PST=PST*1.0E+05 RETURN END *------------------------------------------------- F31C06 = SIGP REAL FUNCTION SIGP(P) REAL P,PI PI=P CALL S99C06('SIGP') SIGP=-1.0E+30 RETURN END *------------------------------------------------- F32C06 = SIGT REAL FUNCTION SIGT(T) REAL T,TI TI=T CALL S99C06('SIGT') SIGT=-1.0E+30 RETURN END *------------------------------------------------- F33C06 = SPD REAL FUNCTION SPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) SPD=F33C06(PI) IF(SPD.EQ.-1.0E+10) THEN CALL S97C06('SPD') ELSE IF(SPD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','SPD') END IF RETURN END *------------------------------------------------- F34C06 = SPDD REAL FUNCTION SPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) SPDD=F34C06(PI) IF(SPDD.EQ.-1.0E+10) THEN CALL S97C06('SPDD') ELSE IF(SPDD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','SPDD') END IF RETURN END *------------------------------------------------- F35C06 = SPT REAL FUNCTION SPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) SPT=F35C06(PI,TI) IF(SPT.EQ.-1.0E+10) THEN CALL S97C06('SPT') ELSE IF(SPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','SPT') END IF RETURN END *------------------------------------------------- F36C06 = SPX REAL FUNCTION SPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) SPX=F36C06(PI,X) IF(SPX.EQ.-1.0E+10) THEN CALL S97C06('SPX') ELSE IF((SPX.EQ.-1.0E+20).AND.(MESS.NE.0)) THEN CALL S98C06(2,P,X,'P','X','SPX') END IF RETURN END *------------------------------------------------- F37C06 = STD REAL FUNCTION STD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) STD=F37C06(TI) IF(STD.EQ.-1.0E+10) THEN CALL S97C06('STD') ELSE IF(STD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','STD') END IF RETURN END *------------------------------------------------- F38C06 = STDD REAL FUNCTION STDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) STDD=F38C06(TI) IF(STDD.EQ.-1.0E+10) THEN CALL S97C06('STDD') ELSE IF(STDD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','STDD') END IF RETURN END *------------------------------------------------- F39C06 = STX REAL FUNCTION STX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) STX=F39C06(TI,X) IF(STX.EQ.-1.0E+10) THEN CALL S97C06('STX') ELSE IF(STX.EQ.-1.0E+20) THEN CALL S98C06(2,T,X,'T','X','STX') END IF RETURN END *------------------------------------------------- F40C06 = TSP REAL FUNCTION TSP(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TSP=F40C06(PI) IF(TSP.EQ.-1.0E+10) THEN CALL S97C06('TSP') RETURN ELSE IF(TSP.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','TSP') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TSP=TSP+273.15 RETURN END *------------------------------------------------- F41C06 = TRPL REAL FUNCTION TRPL(A) CHARACTER A*1 REAL PBAR,T0K,FF INTEGER KPA,MESS COMMON/UNIT/KPA,MESS 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=F41C06(A) IF(FF.EQ.-1.0E+20) THEN IF (MESS.NE.0) THEN WRITE(6,2000) A 2000 FORMAT(1H ,5X,'**** OUT OF RANGE AT TRPL FOR ETHANE', - ' WHEN A =',A,' ****') END IF FF=-1.0E+20 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 *------------------------------------------------- F42C06 = UPD REAL FUNCTION UPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) UPD=F42C06(PI) IF(UPD.EQ.-1.0E+10) THEN CALL S97C06('UPD') ELSE IF(UPD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','UPD') END IF RETURN END *------------------------------------------------- F43C06 = UPDD REAL FUNCTION UPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) UPDD=F43C06(PI) IF(UPDD.EQ.-1.0E+10) THEN CALL S97C06('UPDD') ELSE IF(UPDD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','UPDD') END IF RETURN END *------------------------------------------------- F44C06 = UPT REAL FUNCTION UPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) UPT=F44C06(PI,TI) IF(UPT.EQ.-1.0E+10) THEN CALL S97C06('UPT') ELSE IF(UPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','UPT') END IF RETURN END *------------------------------------------------- F45C06 = UPX REAL FUNCTION UPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) UPX=F45C06(PI,X) IF(UPX.EQ.-1.0E+10) THEN CALL S97C06('UPX') ELSE IF(UPX.EQ.-1.0E+20) THEN CALL S98C06(2,P,X,'P','X','UPX') END IF RETURN END *------------------------------------------------- F46C06 = UTD REAL FUNCTION UTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) UTD=F46C06(TI) IF(UTD.EQ.-1.0E+10) THEN CALL S97C06('UTD') ELSE IF(UTD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','UTD') END IF RETURN END *------------------------------------------------- F47C06 = UTDD REAL FUNCTION UTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) UTDD=F47C06(TI) IF(UTDD.EQ.-1.0E+10) THEN CALL S97C06('UTDD') ELSE IF(UTDD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','UTDD') END IF RETURN END *------------------------------------------------- F48C06 = UTX REAL FUNCTION UTX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) UTX=F48C06(TI,X) IF(UTX.EQ.-1.0E+10) THEN CALL S97C06('UTX') ELSE IF(UTX.EQ.-1.0E+20) THEN CALL S98C06(2,T,X,'T','X','UTX') END IF RETURN END *------------------------------------------------- F49C06 = VPD REAL FUNCTION VPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) VPD=F49C06(PI) IF(VPD.EQ.-1.0E+10) THEN CALL S97C06('VPD') ELSE IF(VPD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','VPD') END IF RETURN END *------------------------------------------------- F50C06 = VPDD REAL FUNCTION VPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) VPDD=F50C06(PI) IF(VPDD.EQ.-1.0E+10) THEN CALL S97C06('VPDD') ELSE IF(VPDD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','VPDD') END IF RETURN END *------------------------------------------------- F51C06 = VPT REAL FUNCTION VPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) VPT=F51C06(PI,TI) IF(VPT.EQ.-1.0E+10) THEN CALL S97C06('VPT') ELSE IF(VPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','VPT') END IF RETURN END *------------------------------------------------- F52C06 = VPX REAL FUNCTION VPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) VPX=F52C06(PI,X) IF(VPX.EQ.-1.0E+10) THEN CALL S97C06('VPX') ELSE IF(VPX.EQ.-1.0E+20) THEN CALL S98C06(2,P,X,'P','X','VPX') END IF RETURN END *------------------------------------------------- F53C06 = VTD REAL FUNCTION VTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) VTD=F53C06(TI) IF(VTD.EQ.-1.0E+10) THEN CALL S97C06('VTD') ELSE IF(VTD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','VTD') END IF RETURN END *------------------------------------------------- F54C06 = VTDD REAL FUNCTION VTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) VTDD=F54C06(TI) IF(VTDD.EQ.-1.0E+10) THEN CALL S97C06('VTDD') ELSE IF(VTDD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','VTDD') END IF RETURN END *------------------------------------------------- F55C06 = VTX REAL FUNCTION VTX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) VTX=F55C06(TI,X) IF(VTX.EQ.-1.0E+10) THEN CALL S97C06('VTX') ELSE IF(VTX.EQ.-1.0E+20) THEN CALL S98C06(2,T,X,'T','X','VTX') END IF RETURN END *------------------------------------------------- F56C06 = XPH REAL FUNCTION XPH(P,H) REAL P,PI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) XPH=F56C06(PI,H) IF(XPH.EQ.-1.0E+10) THEN CALL S97C06('XPH') ELSE IF(XPH.EQ.-1.0E+20) THEN CALL S98C06(2,P,H,'P','H','XPH') END IF RETURN END *------------------------------------------------- F57C06 = XPS REAL FUNCTION XPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) XPS=F57C06(PI,S) IF(XPS.EQ.-1.0E+10) THEN CALL S97C06('XPS') ELSE IF(XPS.EQ.-1.0E+20) THEN CALL S98C06(2,P,S,'P','S','XPS') END IF RETURN END *------------------------------------------------- F58C06 = XPU REAL FUNCTION XPU(P,U) REAL P,PI,U INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) XPU=F58C06(PI,U) IF(XPU.EQ.-1.0E+10) THEN CALL S97C06('XPU') ELSE IF(XPU.EQ.-1.0E+20) THEN CALL S98C06(2,P,U,'P','U','XPU') END IF RETURN END *------------------------------------------------- F59C06 = XPV REAL FUNCTION XPV(P,V) REAL P,PI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) XPV=F59C06(PI,V) IF(XPV.EQ.-1.0E+10) THEN CALL S97C06('XPV') ELSE IF(XPV.EQ.-1.0E+20) THEN CALL S98C06(2,P,V,'P','V','XPV') END IF RETURN END *------------------------------------------------- F60C06 = XTH REAL FUNCTION XTH(T,H) REAL T,TI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) XTH=F60C06(TI,H) IF(XTH.EQ.-1.0E+10) THEN CALL S97C06('XTH') ELSE IF(XTH.EQ.-1.0E+20) THEN CALL S98C06(2,T,H,'T','H','XTH') END IF RETURN END *------------------------------------------------- F61C06 = XTS REAL FUNCTION XTS(T,S) REAL T,TI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) XTS=F61C06(TI,S) IF(XTS.EQ.-1.0E+10) THEN CALL S97C06('XTS') ELSE IF(XTS.EQ.-1.0E+20) THEN CALL S98C06(2,T,S,'T','S','XTS') END IF RETURN END *------------------------------------------------- F62C06 = XTU REAL FUNCTION XTU(T,U) REAL T,TI,U INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) XTU=F62C06(TI,U) IF(XTU.EQ.-1.0E+10) THEN CALL S97C06('XTU') ELSE IF(XTU.EQ.-1.0E+20) THEN CALL S98C06(2,T,U,'T','U','XTU') END IF RETURN END *------------------------------------------------- F63C06 = XTV REAL FUNCTION XTV(T,V) REAL T,TI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) XTV=F63C06(TI,V) IF(XTV.EQ.-1.0E+10) THEN CALL S97C06('XTV') ELSE IF(XTV.EQ.-1.0E+20) THEN CALL S98C06(2,T,V,'T','V','XTV') END IF RETURN END *------------------------------------------------- F64C06 = TPH REAL FUNCTION TPH(P,H) REAL P,PI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TPH=F64C06(PI,H) IF(TPH.EQ.-1.0E+10) THEN CALL S97C06('TPH') RETURN ELSE IF(TPH.EQ.-1.0E+20) THEN CALL S98C06(2,P,H,'P','H','TPH') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPH=TPH+273.15 RETURN END *------------------------------------------------- F65C06 = TPS REAL FUNCTION TPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TPS=F65C06(PI,S) IF(TPS.EQ.-1.0E+10) THEN CALL S97C06('TPS') RETURN ELSE IF(TPS.EQ.-1.0E+20) THEN CALL S98C06(2,P,S,'P','S','TPS') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPS=TPS+273.15 RETURN END *------------------------------------------------- F66C06 = PLDT REAL FUNCTION PLDT(T) REAL T,TI TI=T CALL S99C06('PLDT') PLDT=-1.0E+30 RETURN END *------------------------------------------------- F67C06 = TLDP REAL FUNCTION TLDP(P) REAL P,PI PI=P CALL S99C06('TLDP') TLDP=-1.0E+30 RETURN END *------------------------------------------------- F68C06 = PMLT REAL FUNCTION PMLT(T) REAL TI,T INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) PMLT=F68C06(TI) IF(PMLT.EQ.-1.0E+10) THEN CALL S97C06('PMLT') RETURN ELSE IF(PMLT.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','PMLT') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.2)) RETURN PMLT=PMLT*1.0E+05 RETURN END *------------------------------------------------- F69C06 = TMLP REAL FUNCTION TMLP(P) REAL PI,P INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TMLP=F69C06(PI) IF(TMLP.EQ.-1.0E+10) THEN CALL S97C06('TMLP') RETURN ELSE IF(TMLP.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','TMLP') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TMLP=TMLP+273.15 RETURN END *------------------------------------------------- F70C06 = TPV REAL FUNCTION TPV(P,V) REAL P,PI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TPV=F70C06(PI,V) IF(TPV.EQ.-1.0E+10) THEN CALL S97C06('TPV') RETURN ELSE IF(TPV.EQ.-1.0E+20) THEN CALL S98C06(2,P,V,'P','V','TPV') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPV=TPV+273.15 RETURN END *------------------------------------------------- F71C06 = HPS REAL FUNCTION HPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) HPS=F71C06(PI,S) IF(HPS.EQ.-1.0E+10) THEN CALL S97C06('HPS') ELSE IF(HPS.EQ.-1.0E+20) THEN CALL S98C06(2,P,S,'P','S','HPS') END IF RETURN END *------------------------------------------------- F72C06 = PSTD REAL FUNCTION PSTD(T) REAL T,TI TI=T CALL S99C06('PSTD') PSTD=-1.0E+30 RETURN END *------------------------------------------------- F73C06 = PSTDD REAL FUNCTION PSTDD(T) REAL T,TI TI=T CALL S99C06('PSTDD') PSTDD=-1.0E+30 RETURN END *------------------------------------------------- F74C06 = TSPD REAL FUNCTION TSPD(P) REAL P,PI PI=P CALL S99C06('TSPD') TSPD=-1.0E+30 RETURN END *------------------------------------------------- F75C06 = TSPDD REAL FUNCTION TSPDD(P) REAL P,PI PI=P CALL S99C06('TSPDD') TSPDD=-1.0E+30 RETURN END *------------------------------------------------- F76C06 = CVPDD REAL FUNCTION CVPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) CVPDD=F76C06(PI) IF(CVPDD.EQ.-1.0E+10) THEN CALL S97C06('CVPDD') ELSE IF(CVPDD.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','CVPDD') END IF RETURN END *------------------------------------------------- F77C06 = CVPT REAL FUNCTION CVPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) CVPT=F77C06(PI,TI) IF(CVPT.EQ.-1.0E+10) THEN CALL S97C06('CVPT') ELSE IF(CVPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','CVPT') END IF RETURN END *------------------------------------------------- F78C06 = CVTDD REAL FUNCTION CVTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C06(KPA,T) CVTDD=F78C06(TI) IF(CVTDD.EQ.-1.0E+10) THEN CALL S97C06('CVTDD') ELSE IF(CVTDD.EQ.-1.0E+20) THEN CALL S98C06(1,T,T,'T','T','CVTDD') END IF RETURN END *------------------------------------------------- F79C06 = UPS REAL FUNCTION UPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) UPS=F79C06(PI,S) IF(UPS.EQ.-1.0E+10) THEN CALL S97C06('UPS') ELSE IF(UPS.EQ.-1.0E+20) THEN CALL S98C06(2,P,S,'P','S','UPS') END IF RETURN END *------------------------------------------------- F80C06 = VPS REAL FUNCTION VPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) VPS=F80C06(PI,S) IF(VPS.EQ.-1.0E+10) THEN CALL S97C06('VPS') ELSE IF(VPS.EQ.-1.0E+20) THEN CALL S98C06(2,P,S,'P','S','VPS') END IF RETURN END *------------------------------------------------- F81C06 = PRPT REAL FUNCTION PRPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C06('PRPT') PRPT=-1.0E+30 RETURN END *------------------------------------------------- F82C06 = AKPT REAL FUNCTION AKPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) AKPT=F82C06(PI,TI) IF(AKPT.EQ.-1.0E+10) THEN CALL S97C06('AKPT') ELSE IF(AKPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','AKPT') END IF RETURN END *------------------------------------------------- F83C06 = WPT REAL FUNCTION WPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) WPT=F83C06(PI,TI) IF(WPT.EQ.-1.0E+10) THEN CALL S97C06('WPT') ELSE IF(WPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','WPT') END IF RETURN END *------------------------------------------------- F85C06 = PRPD REAL FUNCTION PRPD(P) REAL P,PI PI=P CALL S99C06('PRPD') PRPD=-1.0E+30 RETURN END *------------------------------------------------- F86C06 = PRPDD REAL FUNCTION PRPDD(P) REAL P,PI PI=P CALL S99C06('PRPDD') PRPDD=-1.0E+30 RETURN END *------------------------------------------------- F87C06 = PRTD REAL FUNCTION PRTD(T) REAL T,TI TI=T CALL S99C06('PRTD') PRTD=-1.0E+30 RETURN END *------------------------------------------------- F88C06 = PRTDD REAL FUNCTION PRTDD(T) REAL T,TI TI=T CALL S99C06('PRTDD') PRTDD=-1.0E+30 RETURN END *------------------------------------------------- F90C06 = BSPT REAL FUNCTION BSPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) BSPT=F90C06(PI,TI) IF(BSPT.EQ.-1.0E+10) THEN CALL S97C06('BSPT') RETURN ELSE IF(BSPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','BSPT') RETURN END IF IF ((KPA.EQ.1).OR.(KPA.EQ.2)) RETURN BSPT=BSPT*1.0E-05 RETURN END *------------------------------------------------- F91C06 = BTPT REAL FUNCTION BTPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) BTPT=F91C06(PI,TI) IF(BTPT.EQ.-1.0E+10) THEN CALL S97C06('BTPT') RETURN ELSE IF(BTPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','BTPT') RETURN END IF IF ((KPA.EQ.1).OR.(KPA.EQ.2)) RETURN BTPT=BTPT*1.0E-05 RETURN END *------------------------------------------------- F92C06 = BPPT REAL FUNCTION BPPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) BPPT=F92C06(PI,TI) IF(BPPT.EQ.-1.0E+10) THEN CALL S97C06('BPPT') ELSE IF(BPPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','BPPT') END IF RETURN END *------------------------------------------------- F93C06 = BVPT REAL FUNCTION BVPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) BVPT=F93C06(PI,TI) IF(BVPT.EQ.-1.0E+10) THEN CALL S97C06('BVPT') ELSE IF(BVPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','BVPT') END IF RETURN END *------------------------------------------------- F94C06 = AJTPT REAL FUNCTION AJTPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TI=G99C06(KPA,T) AJTPT=F94C06(PI,TI) IF(AJTPT.EQ.-1.0E+10) THEN CALL S97C06('AJTPT') RETURN ELSE IF(AJTPT.EQ.-1.0E+20) THEN CALL S98C06(2,P,T,'P','T','AJTPT') RETURN END IF IF ((KPA.EQ.1).OR.(KPA.EQ.2)) RETURN AJTPT=AJTPT*1.0E-05 RETURN END *------------------------------------------------- F95C06 = GAMPT REAL FUNCTION GAMPT(P,T) REAL P,PI PI=P TI=T CALL S99C06('GAMPT') GAMPT=-1.0E+30 RETURN END *------------------------------------------------- F96C06 = GAMPDD REAL FUNCTION GAMPDD(P) REAL P,PI PI=P CALL S99C06('GAMPDD') GAMPDD=-1.0E+30 RETURN END *------------------------------------------------- F97C06 = GAMTDD REAL FUNCTION GAMTDD(T) REAL T,TI TI=T CALL S99C06('GAMTDD') GAMTDD=-1.0E+30 RETURN END *------------------------------------------------- F98C06 = TPSEUP REAL FUNCTION TPSEUP(P) REAL PI,P INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C06(KPA,P) TPSEUP=F98C06(PI) IF(TPSEUP.EQ.-1.0E+10) THEN CALL S97C06('TPSEUP') RETURN ELSE IF(TPSEUP.EQ.-1.0E+20) THEN CALL S98C06(1,P,P,'P','P','TPSEUP') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPSEUP=TPSEUP+273.15 RETURN END *------------------------------------------------- F99C06 = PSBT REAL FUNCTION PSBT(T) REAL T,TI TI=T CALL S99C06('PSBT') PSBT=-1.0E+30 RETURN END *------------------------------------------------- F100C06 = TSBP REAL FUNCTION TSBP(P) REAL P,PI PI=P CALL S99C06('TSBP') TSBP=-1.0E+30 RETURN END *------------------------------------------------- G98C06 REAL FUNCTION G98C06(KPA,P) INTEGER KPA REAL P,PBAR IF (KPA.EQ.1) THEN PBAR=1.0E+00 ELSE IF (KPA.EQ.2) THEN PBAR=1.0E+00 ELSE IF (KPA.EQ.3) THEN PBAR=1.0E-05 ELSE PBAR=1.0E-05 END IF G98C06=P*PBAR RETURN END *------------------------------------------------- G99C06 REAL FUNCTION G99C06(KPA,T) INTEGER KPA REAL T,T0K IF (KPA.EQ.1) THEN T0K=0.0E+00 ELSE IF (KPA.EQ.2) THEN T0K=273.15E+00 ELSE IF (KPA.EQ.3) THEN T0K=0.0E+00 ELSE T0K=273.15E+00 END IF G99C06=T-T0K RETURN END ***** PROPATH V81:ETHANE FUNCTIONS ****1990.9.28************************ REAL FUNCTION F4C06(FP) ***** FP=INPUT PRESSURE IN BAR ***** IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F4C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F4C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAUS,OMEGD,HD) CALL S9C06(TAUS,OMEGDD,HDD) ALHP=HDD-HD F4C06=REAL(ALHP) RETURN END REAL FUNCTION F5C06(FT) ***** FT=INPUT TEMPERATURE IN C(DEGREE CELSIUS) ***** IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F5C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F5C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAU,OMEGD,HD) CALL S9C06(TAU,OMEGDD,HDD) ALHT=HDD-HD F5C06=REAL(ALHT) RETURN END REAL FUNCTION F16C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F16C06=-1.0E+20 IF ((FP.LT.0.34E-02).OR.(FP.GT.48.715)) RETURN F16C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S12C06(TAUS,OMEGD,CPPD) F16C06=REAL(CPPD) RETURN END REAL FUNCTION F17C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F17C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F17C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S12C06(TAUS,OMEGDD,CPPDD) F17C06=REAL(CPPDD) RETURN END REAL FUNCTION F18C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F18C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F18C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S12C06(TAU,OMEGA,CPPT) F18C06=REAL(CPPT) RETURN END REAL FUNCTION F19C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F19C06=-1.0E+20 IF ((FT.LT.-153.25).OR.(FT.GT.32.181)) RETURN F19C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S12C06(TAU,OMEGD,CPTD) F19C06=REAL(CPTD) RETURN END REAL FUNCTION F20C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F20C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F20C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S12C06(TAU,OMEGDD,CPTDD) F20C06=REAL(CPTDD) RETURN END REAL FUNCTION F21C06(A) CHARACTER*1 A,B(1:5) DOUBLE PRECISION CRP(1:5),HCR,SCR DATA B(1)/'H'/,B(2)/'P'/,B(3)/'S'/,B(4)/'T'/,B(5)/'V'/ CALL S9C06(1.0D0,1.0D0,HCR) CALL S8C06(1.0D0,1.0D0,SCR) CRP(1)=HCR CRP(2)=4.8714D+01 CRP(3)=SCR CRP(4)=32.18D+00 CRP(5)=1.0D0/2.0446D+02 DO 10 I=1,5 F21C06=REAL(CRP(I)) IF(A.EQ.B(I)) RETURN 10 CONTINUE F21C06=-1.0E+20 RETURN END REAL FUNCTION F23C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F23C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F23C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAUS,OMEGD,HPD) F23C06=REAL(HPD) RETURN END REAL FUNCTION F24C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F24C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F24C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S9C06(TAUS,OMEGDD,HPDD) F24C06=REAL(HPDD) RETURN END REAL FUNCTION F25C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F25C06=-1.0E+20 CALL S90C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F25C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S9C06(TAU,OMEGA,HPT) F25C06=REAL(HPT) RETURN END REAL FUNCTION F26C06(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F26C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F26C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAUS,OMEGD,HD) CALL S9C06(TAUS,OMEGDD,HDD) HPX=HD+X*(HDD-HD) F26C06=REAL(HPX) RETURN END REAL FUNCTION F27C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F27C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F27C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAU,OMEGD,HTD) F27C06=REAL(HTD) RETURN END REAL FUNCTION F28C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F28C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F28C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S9C06(TAU,OMEGDD,HTDD) F28C06=REAL(HTDD) RETURN END REAL FUNCTION F29C06(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F29C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F29C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR X=DBLE(FX) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAU,OMEGD,HD) CALL S9C06(TAU,OMEGDD,HDD) HTX=HD+X*(HDD-HD) F29C06=REAL(HTX) RETURN END REAL FUNCTION F30C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00,PCR=4.8714D+06) F30C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F30C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR CC TEST FOR EPS IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN PST=PIS*PCR F30C06=REAL(PST*1.0D-05) RETURN END REAL FUNCTION F33C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F33C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F33C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGD,SPD) F33C06=REAL(SPD) RETURN END REAL FUNCTION F34C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F34C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F34C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGDD,SPDD) F34C06=REAL(SPDD) RETURN END REAL FUNCTION F35C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F35C06=-1.0E+20 CALL S90C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F35C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S8C06(TAU,OMEGA,SPT) F35C06=REAL(SPT) RETURN END REAL FUNCTION F36C06(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F36C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F36C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGD,SD) CALL S8C06(TAUS,OMEGDD,SDD) SPX=SD+X*(SDD-SD) F36C06=REAL(SPX) RETURN END REAL FUNCTION F37C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F37C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F37C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAU,OMEGD,STD) F37C06=REAL(STD) RETURN END REAL FUNCTION F38C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F38C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F38C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S8C06(TAU,OMEGDD,STDD) F38C06=REAL(STDD) RETURN END REAL FUNCTION F39C06(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F39C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F39C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR X=DBLE(FX) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAU,OMEGD,SD) CALL S8C06(TAU,OMEGDD,SDD) STX=SD+X*(SDD-SD) F39C06=REAL(STX) RETURN END REAL FUNCTION F40C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F40C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F40C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN TSP=TCR*TAUS F40C06=REAL(TSP-273.15D0) RETURN END REAL FUNCTION F41C06(A) ***** TRIPLE POINT FROM TABLE ON PAGE 93 ***** CHARACTER*1 A,B(1:2) DOUBLE PRECISION TRPL(1:2) DATA B(1)/'P'/,B(2)/'T'/ TRPL(1)=1.13D-05 TRPL(2)=-182.802D+00 DO 10 I=1,2 F41C06=REAL(TRPL(I)) IF(A.EQ.B(I)) RETURN 10 CONTINUE F41C06=-1.0E+20 RETURN END REAL FUNCTION F42C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F42C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F42C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S10C06(TAUS,OMEGD,UPD) F42C06=REAL(UPD) RETURN END REAL FUNCTION F43C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F43C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F43C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S10C06(TAUS,OMEGDD,UPDD) F43C06=REAL(UPDD) RETURN END REAL FUNCTION F44C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F44C06=-1.0E+20 CALL S90C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F44C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S10C06(TAU,OMEGA,UPT) F44C06=REAL(UPT) RETURN END REAL FUNCTION F45C06(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F45C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F45C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S10C06(TAUS,OMEGD,UD) CALL S10C06(TAUS,OMEGDD,UDD) UPX=UD+X*(UDD-UD) F45C06=REAL(UPX) RETURN END REAL FUNCTION F46C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F46C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F46C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S10C06(TAU,OMEGD,UTD) F46C06=REAL(UTD) RETURN END REAL FUNCTION F47C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F47C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F47C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S10C06(TAU,OMEGDD,UTDD) F47C06=REAL(UTDD) RETURN END REAL FUNCTION F48C06(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F48C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F48C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR X=DBLE(FX) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S10C06(TAU,OMEGD,UD) CALL S10C06(TAU,OMEGDD,UDD) UTX=UD+X*(UDD-UD) F48C06=REAL(UTX) RETURN END REAL FUNCTION F49C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,RHOCR=2.044571662D+02) F49C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F49C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN VPD=1.0D0/(OMEGD*RHOCR) F49C06=REAL(VPD) RETURN END REAL FUNCTION F50C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,RHOCR=2.044571662D+02) F50C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F50C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGDD.LE.-1.0D+10) RETURN VPDD=1.0D0/(OMEGDD*RHOCR) F50C06=REAL(VPDD) RETURN END REAL FUNCTION F51C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00,RHOCR=2.044571662D+02) F51C06=-1.0E+20 CALL S90C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F51C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN VPT=1.0D0/(OMEGA*RHOCR) F51C06=REAL(VPT) RETURN END REAL FUNCTION F52C06(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,RHOCR=2.044571662D+02) F52C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F52C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN VD=1.0D0/(OMEGD*RHOCR) VDD=1.0D0/(OMEGDD*RHOCR) VPX=VD+X*(VDD-VD) F52C06=REAL(VPX) RETURN END REAL FUNCTION F53C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00,RHOCR=2.044571662D+02) F53C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F53C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN VTD=1.0D0/(OMEGD*RHOCR) F53C06=REAL(VTD) RETURN END REAL FUNCTION F54C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00,RHOCR=2.044571662D+02) F54C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F54C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGDD.LE.-1.0D+10) RETURN VTDD=1.0D0/(OMEGDD*RHOCR) F54C06=REAL(VTDD) RETURN END REAL FUNCTION F55C06(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00,RHOCR=2.044571662D+02) F55C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F55C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR X=DBLE(FX) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN VD=1.0D0/(OMEGD*RHOCR) VDD=1.0D0/(OMEGDD*RHOCR) VTX=VD+X*(VDD-VD) F55C06=REAL(VTX) RETURN END REAL FUNCTION F56C06(FP,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F56C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.67)) RETURN F56C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR H=DBLE(FH) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAUS,OMEGD,HPD) CALL S9C06(TAUS,OMEGDD,HPDD) HMIN=HPD-ABS(HPD)*1.0D-05 HMAX=HPDD+ABS(HPDD)*1.0D-05 IF ((H.LT.HMIN).OR.(H.GT.HMAX)) THEN XPH=-1.0D+20 ELSE XPH=(H-HPD)/(HPDD-HPD) IF (XPH.LT.0.0D0) XPH=0.0D0 IF (XPH.GT.1.0D0) XPH=1.0D0 END IF F56C06=REAL(XPH) RETURN END REAL FUNCTION F57C06(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F57C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.67)) RETURN F57C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGD,SPD) CALL S8C06(TAUS,OMEGDD,SPDD) SMIN=SPD-ABS(SPD)*1.0D-05 SMAX=SPDD+ABS(SPDD)*1.0D-05 IF ((S.LT.SMIN).OR.(S.GT.SMAX)) THEN XPS=-1.0D+20 ELSE XPS=(S-SPD)/(SPDD-SPD) IF (XPS.LT.0.0D0) XPS=0.0D0 IF (XPS.GT.1.0D0) XPS=1.0D0 END IF F57C06=REAL(XPS) RETURN END REAL FUNCTION F58C06(FP,FU) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F58C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.67)) RETURN F58C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR U=DBLE(FU) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S10C06(TAUS,OMEGD,UPD) CALL S10C06(TAUS,OMEGDD,UPDD) UMIN=UPD-ABS(UPD)*1.0D-05 UMAX=UPDD+ABS(UPDD)*1.0D-05 IF ((U.LT.UMIN).OR.(U.GT.UMAX)) THEN XPU=-1.0D+20 ELSE XPU=(U-UPD)/(UPDD-UPD) IF (XPU.LT.0.0D0) XPU=0.0D0 IF (XPU.GT.1.0D0) XPU=1.0D0 END IF F58C06=REAL(XPU) RETURN END REAL FUNCTION F59C06(FP,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,RHOCR=2.044571662D+02) F59C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.67)) RETURN F59C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR V=DBLE(FV) CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN VPD=1.0D0/(OMEGD*RHOCR) VPDD=1.0D0/(OMEGDD*RHOCR) VMIN=VPD*0.99999D0 VMAX=VPDD*1.00001D0 IF ((V.LT.VMIN).OR.(V.GT.VMAX)) THEN XPV=-1.0D+20 ELSE XPV=(V-VPD)/(VPDD-VPD) IF (XPV.LT.0.0D0) XPV=0.0D0 IF (XPV.GT.1.0D0) XPV=1.0D0 END IF F59C06=REAL(XPV) RETURN END REAL FUNCTION F60C06(FT,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F60C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.15)) RETURN F60C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR H=DBLE(FH) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAU,OMEGD,HTD) CALL S9C06(TAU,OMEGDD,HTDD) HMIN=HTD-ABS(HTD)*1.0D-05 HMAX=HTDD+ABS(HTDD)*1.0D-05 IF ((H.LT.HMIN).OR.(H.GT.HMAX)) THEN XTH=-1.0D+20 ELSE XTH=(H-HTD)/(HTDD-HTD) IF (XTH.LT.0.0D0) XTH=0.0D0 IF (XTH.GT.1.0D0) XTH=1.0D0 END IF F60C06=REAL(XTH) RETURN END REAL FUNCTION F61C06(FT,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F61C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.15)) RETURN F61C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR S=DBLE(FS) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAU,OMEGD,STD) CALL S8C06(TAU,OMEGDD,STDD) SMIN=STD-ABS(STD)*1.0D-05 SMAX=STDD+ABS(STDD)*1.0D-05 IF ((S.LT.SMIN).OR.(S.GT.SMAX)) THEN XTS=-1.0D+20 ELSE XTS=(S-STD)/(STDD-STD) IF (XTS.LT.0.0D0) XTS=0.0D0 IF (XTS.GT.1.0D0) XTS=1.0D0 END IF F61C06=REAL(XTS) RETURN END REAL FUNCTION F62C06(FT,FU) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F62C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.15)) RETURN F62C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR U=DBLE(FU) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S10C06(TAU,OMEGD,UTD) CALL S10C06(TAU,OMEGDD,UTDD) UMIN=UTD-ABS(UTD)*1.0D-05 UMAX=UTDD+ABS(UTDD)*1.0D-05 IF ((U.LT.UMIN).OR.(U.GT.UMAX)) THEN XTU=-1.0D+20 ELSE XTU=(U-UTD)/(UTDD-UTD) IF (XTU.LT.0.0D0) XTU=0.0D0 IF (XTU.GT.1.0D0) XTU=1.0D0 END IF F62C06=REAL(XTU) RETURN END REAL FUNCTION F63C06(FT,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00,RHOCR=2.044571662D+02) F63C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.15)) RETURN F63C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR V=DBLE(FV) CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGD.LE.-1.0D+10) RETURN VTD=1.0D0/(OMEGD*RHOCR) VTDD=1.0D0/(OMEGDD*RHOCR) VMIN=VTD*0.99999D0 VMAX=VTDD*1.00001D0 IF ((V.LT.VMIN).OR.(V.GT.VMAX)) THEN XTV=-1.0D+20 ELSE XTV=(V-VTD)/(VTDD-VTD) IF (XTV.LT.0.0D0) XTV=0.0D0 IF (XTV.GT.1.0D0) XTV=1.0D0 END IF F63C06=REAL(XTV) RETURN END REAL FUNCTION F64C06(FP,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F64C06=-1.0E+20 IF ((FP.GE.0.9999).AND.(FP.LE.40.421)) THEN FHMIN=F25C06(FP,-182.2) ELSE IF ((FP.GT.40.421).AND.(FP.LE.800.1)) THEN FTMIN=F69C06(FP) FHMIN=F25C06(FP,FTMIN) ELSE RETURN END IF FHMAX=F25C06(FP,426.86) IF ((FHMIN.LE.-1.0E+10).OR.(FHMAX.LE.-1.0E+10)) RETURN FHMIN=FHMIN-ABS(FHMIN)*1.0E-03 FHMAX=FHMAX+ABS(FHMAX)*1.0E-03 IF ((FH.LT.FHMIN).OR.(FH.GT.FHMAX)) RETURN F64C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR H=DBLE(FH) CALL S9C06(1.0D0,1.0D0,HCR) IF ((PI.GT.0.9992D0).AND.(PI.LT.1.025D0).AND. - (ABS(H/HCR-1.0D0).LT.0.0325D0)) THEN F64C06=REAL(TCR-273.15D0) RETURN END IF IF (PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S4C06(PI,TAU0,OMEGA0) IF (OMEGA0.LE.-1.0D+10) RETURN CALL S9C06(TAU0,OMEGA0,H0) CALL S12C06(TAU0,OMEGA0,CP0) ELSE CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S9C06(TAUS,OMEGD,HD) CALL S9C06(TAUS,OMEGDD,HDD) IF ((H.GE.HD).AND.(H.LE.HDD)) THEN F64C06=REAL(TCR*TAUS-273.15D0) RETURN END IF IF (H.GT.HDD) THEN CALL S12C06(TAUS,OMEGDD,CPDD) TAU0=TAUS H0=HDD CP0=CPDD ELSE IF (H.LT.HD) THEN CALL S12C06(TAUS,OMEGD,CPD) TAU0=TAUS H0=HD CP0=CPD END IF END IF TAU1=TAU0+(H-H0)/(CP0*TCR) CALL S4C06(PI,TAU1,OMEGA1) IF (OMEGA1.LE.-1.0D+10) RETURN CALL S9C06(TAU1,OMEGA1,H1) TH0=TAU0 TH1=TAU1 CALL S40C06(1,PI,H,TH0,TH1,H0,H1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=THW TPH=TCR*TAUW F64C06=REAL(TPH-273.15D0) RETURN END REAL FUNCTION F65C06(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F65C06=-1.0E+20 IF ((FP.GE.0.9999).AND.(FP.LE.40.421)) THEN FSMIN=F35C06(FP,-182.2) ELSE IF ((FP.GT.40.421).AND.(FP.LE.800.1)) THEN FTMIN=F69C06(FP) FSMIN=F35C06(FP,FTMIN) ELSE RETURN END IF FSMAX=F35C06(FP,426.86) IF ((FSMIN.LE.-1.0E+10).OR.(FSMAX.LE.-1.0E+10)) RETURN FSMIN=FSMIN-ABS(FSMIN)*1.0E-03 FSMAX=FSMAX+ABS(FSMAX)*1.0E-03 IF ((FS.LT.FSMIN).OR.(FS.GT.FSMAX)) RETURN F65C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS) CALL S8C06(1.0D0,1.0D0,SCR) IF ((PI.GT.0.9992D0).AND.(PI.LT.1.025D0).AND. - (ABS(S/SCR-1.0D0).LT.0.0205D0)) THEN F65C06=REAL(TCR-273.15D0) RETURN END IF IF (PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S4C06(PI,TAU0,OMEGA0) IF (OMEGA0.LE.-1.0D+10) RETURN CALL S8C06(TAU0,OMEGA0,S0) CALL S12C06(TAU0,OMEGA0,CP0) IF (CP0.GT.50.0D+03) CP0=50.0D+03 ELSE CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGD,SD) CALL S8C06(TAUS,OMEGDD,SDD) IF ((S.GE.SD).AND.(S.LE.SDD)) THEN F65C06=REAL(TCR*TAUS-273.15D0) RETURN END IF IF (S.GT.SDD) THEN CALL S12C06(TAUS,OMEGDD,CPDD) TAU0=TAUS S0=SDD CP0=CPDD ELSE IF (S.LT.SD) THEN CALL S12C06(TAUS,OMEGD,CPD) TAU0=TAUS S0=SD CP0=CPD END IF END IF TAU1=TAU0*EXP((S-S0)/CP0) CALL S4C06(PI,TAU1,OMEGA1) IF (OMEGA1.LE.-1.0D+10) RETURN CALL S8C06(TAU1,OMEGA1,S1) TH0=TAU0 TH1=TAU1 CALL S40C06(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=THW TPS=TCR*TAUW F65C06=REAL(TPS-273.15D0) RETURN END REAL FUNCTION F68C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F68C06=-1.0E+20 IF ((FT.LT.-182.802).OR.(FT.GT.-170.8)) RETURN T=DBLE(FT)+273.15D0 PMLT=1.1208D-05+255.965D+01*((T/90.348D0)**2.179D0-1.0D0) F68C06=REAL(PMLT) RETURN END REAL FUNCTION F69C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F69C06=-1.0E+20 IF ((FP.LT.39.9).OR.(FP.GT.800.1)) RETURN P=DBLE(FP) TMLP=90.348D0*(1.0D0+(P-1.1208D-05)/255.965D+01) - **(1.0D0/2.179D0)-273.15D0 F69C06=REAL(TMLP) RETURN END REAL FUNCTION F70C06(FP,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,RHOCR=2.044571662D+02,TCR=305.33D+00) F70C06=-1.0E+20 IF ((FP.GE.0.9999).AND.(FP.LE.40.421)) THEN FVMIN=F51C06(FP,-182.2) ELSE IF ((FP.GT.40.421).AND.(FP.LE.800.1)) THEN FTMIN=F69C06(FP) FVMIN=F51C06(FP,FTMIN) ELSE RETURN END IF FVMAX=F51C06(FP,426.86) IF ((FVMIN.LE.-1.0E+10).OR.(FVMAX.LE.-1.0E+10)) RETURN IF ((FV.LT.FVMIN*0.9996).OR.(FV.GT.FVMAX*1.002)) RETURN F70C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR OMEGA=1.0D0/(DBLE(FV)*RHOCR) IF (PI.GE.1.0D0) THEN CALL S5C06(PI,OMEGA,TAU) IF (TAU.LE.-1.0D+10) RETURN TAUW=TAU ELSE CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN IF ((TAUS.EQ.1.0D0).AND.(OMEGA.GE.0.7849D0).AND. - (OMEGA.LE.1.2380D0)) THEN F70C06=REAL(TCR-273.15D0) RETURN END IF IF ((OMEGA.GE.OMEGDD).AND.(OMEGA.LE.OMEGD)) THEN TAUW=TAUS ELSE CALL S5C06(PI,OMEGA,TAU) IF (TAU.LE.-1.0D+10) RETURN TAUW=TAU END IF END IF TPV=TCR*TAUW F70C06=REAL(TPV-273.15D0) RETURN END REAL FUNCTION F71C06(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F71C06=-1.0E+20 IF ((FP.GE.0.9999).AND.(FP.LE.40.421)) THEN FSMIN=F35C06(FP,-182.2) ELSE IF ((FP.GT.40.421).AND.(FP.LE.800.1)) THEN FTMIN=F69C06(FP) FSMIN=F35C06(FP,FTMIN) ELSE RETURN END IF FSMAX=F35C06(FP,426.86) IF ((FSMIN.LE.-1.0E+10).OR.(FSMAX.LE.-1.0E+10)) RETURN FSMIN=FSMIN-ABS(FSMIN)*1.0E-03 FSMAX=FSMAX+ABS(FSMAX)*1.0E-03 IF ((FS.LT.FSMIN).OR.(FS.GT.FSMAX)) RETURN F71C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS) CALL S8C06(1.0D0,1.0D0,SCR) IF ((PI.GT.0.9992D0).AND.(PI.LT.1.025D0).AND. - (ABS(S/SCR-1.0D0).LT.0.0205D0)) THEN CALL S9C06(1.0D0,1.0D0,HCR) F71C06=REAL(HCR) RETURN END IF IF (PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S4C06(PI,TAU0,OMEGA0) IF (OMEGA0.LE.-1.0D+10) RETURN CALL S8C06(TAU0,OMEGA0,S0) CALL S12C06(TAU0,OMEGA0,CP0) ELSE CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGD,SD) CALL S8C06(TAUS,OMEGDD,SDD) IF ((S.GE.SD).AND.(S.LE.SDD)) THEN XPS=(S-SD)/(SDD-SD) CALL S9C06(TAUS,OMEGD,HD) CALL S9C06(TAUS,OMEGDD,HDD) HPS=HD+XPS*(HDD-HD) F71C06=REAL(HPS) RETURN END IF IF (S.GT.SDD) THEN CALL S12C06(TAUS,OMEGDD,CPDD) TAU0=TAUS S0=SDD CP0=CPDD ELSE IF (S.LT.SD) THEN CALL S12C06(TAUS,OMEGD,CPD) TAU0=TAUS S0=SD CP0=CPD END IF END IF TAU1=TAU0*EXP((S-S0)/CP0) CALL S4C06(PI,TAU1,OMEGA1) IF (OMEGA1.LE.-1.0D+10) RETURN CALL S8C06(TAU1,OMEGA1,S1) TH0=TAU0 TH1=TAU1 CALL S40C06(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=THW CALL S4C06(PI,TAUW,OMEGAW) IF (OMEGAW.LE.-1.0D+10) RETURN CALL S9C06(TAUW,OMEGAW,HW) HPS=HW F71C06=REAL(HPS) RETURN END REAL FUNCTION F76C06(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F76C06=-1.0E+20 IF ((FP.LT.0.18E-04).OR.(FP.GT.48.715)) RETURN F76C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S11C06(TAUS,OMEGDD,CVPDD) F76C06=REAL(CVPDD) RETURN END REAL FUNCTION F77C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F77C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F77C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S11C06(TAU,OMEGA,CVPT) F77C06=REAL(CVPT) RETURN END REAL FUNCTION F78C06(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=305.33D+00) F78C06=-1.0E+20 IF ((FT.LT.-181.16).OR.(FT.GT.32.181)) RETURN F78C06=-1.0E+10 TAU=(DBLE(FT)+273.15D0)/TCR IF (ABS(TAU-1.0D0).LT.0.9D-04) THEN TAU=1.0D0 END IF CALL S2C06(TAU,OMEGD,OMEGDD,PIS) IF (OMEGDD.LE.-1.0D+10) RETURN CALL S11C06(TAU,OMEGDD,CVTDD) F78C06=REAL(CVTDD) RETURN END REAL FUNCTION F79C06(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06) F79C06=-1.0E+20 IF ((FP.GE.0.9999).AND.(FP.LE.40.421)) THEN FSMIN=F35C06(FP,-182.2) ELSE IF ((FP.GT.40.421).AND.(FP.LE.800.1)) THEN FTMIN=F69C06(FP) FSMIN=F35C06(FP,FTMIN) ELSE RETURN END IF FSMAX=F35C06(FP,426.86) IF ((FSMIN.LE.-1.0E+10).OR.(FSMAX.LE.-1.0E+10)) RETURN FSMIN=FSMIN-ABS(FSMIN)*1.0E-03 FSMAX=FSMAX+ABS(FSMAX)*1.0E-03 IF ((FS.LT.FSMIN).OR.(FS.GT.FSMAX)) RETURN F79C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS) CALL S8C06(1.0D0,1.0D0,SCR) IF ((PI.GT.0.9992D0).AND.(PI.LT.1.025D0).AND. - (ABS(S/SCR-1.0D0).LT.0.0205D0)) THEN CALL S10C06(1.0D0,1.0D0,UCR) F79C06=REAL(UCR) RETURN END IF IF (PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S4C06(PI,TAU0,OMEGA0) IF (OMEGA0.LE.-1.0D+10) RETURN CALL S8C06(TAU0,OMEGA0,S0) CALL S12C06(TAU0,OMEGA0,CP0) ELSE CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGD,SD) CALL S8C06(TAUS,OMEGDD,SDD) IF ((S.GE.SD).AND.(S.LE.SDD)) THEN XPS=(S-SD)/(SDD-SD) CALL S10C06(TAUS,OMEGD,UD) CALL S10C06(TAUS,OMEGDD,UDD) UPS=UD+XPS*(UDD-UD) F79C06=REAL(UPS) RETURN END IF IF (S.GT.SDD) THEN CALL S12C06(TAUS,OMEGDD,CPDD) TAU0=TAUS S0=SDD CP0=CPDD ELSE IF (S.LT.SD) THEN CALL S12C06(TAUS,OMEGD,CPD) TAU0=TAUS S0=SD CP0=CPD END IF END IF TAU1=TAU0*EXP((S-S0)/CP0) CALL S4C06(PI,TAU1,OMEGA1) IF (OMEGA1.LE.-1.0D+10) RETURN CALL S8C06(TAU1,OMEGA1,S1) TH0=TAU0 TH1=TAU1 CALL S40C06(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=THW CALL S4C06(PI,TAUW,OMEGAW) IF (OMEGAW.LE.-1.0D+10) RETURN CALL S10C06(TAUW,OMEGAW,UW) UPS=UW F79C06=REAL(UPS) RETURN END REAL FUNCTION F80C06(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,RHOCR=204.4571662D+00) F80C06=-1.0E+20 IF ((FP.GE.0.9999).AND.(FP.LE.40.421)) THEN FSMIN=F35C06(FP,-182.2) ELSE IF ((FP.GT.40.421).AND.(FP.LE.800.1)) THEN FTMIN=F69C06(FP) FSMIN=F35C06(FP,FTMIN) ELSE RETURN END IF FSMAX=F35C06(FP,426.86) IF ((FSMIN.LE.-1.0E+10).OR.(FSMAX.LE.-1.0E+10)) RETURN FSMIN=FSMIN-ABS(FSMIN)*1.0E-03 FSMAX=FSMAX+ABS(FSMAX)*1.0E-03 IF ((FS.LT.FSMIN).OR.(FS.GT.FSMAX)) RETURN F80C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS) CALL S8C06(1.0D0,1.0D0,SCR) IF ((PI.GT.0.9992D0).AND.(PI.LT.1.025D0).AND. - (ABS(S/SCR-1.0D0).LT.0.0205D0)) THEN VPS= 1.0D0/RHOCR F80C06=REAL(VPS) RETURN END IF IF (PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S4C06(PI,TAU0,OMEGA0) IF (OMEGA0.LE.-1.0D+10) RETURN CALL S8C06(TAU0,OMEGA0,S0) CALL S12C06(TAU0,OMEGA0,CP0) ELSE CALL S3C06(PI,OMEGD,OMEGDD,TAUS) IF (OMEGD.LE.-1.0D+10) RETURN CALL S8C06(TAUS,OMEGD,SD) CALL S8C06(TAUS,OMEGDD,SDD) IF ((S.GE.SD).AND.(S.LE.SDD)) THEN XPS=(S-SD)/(SDD-SD) VD=1.0D0/(RHOCR*OMEGD) VDD=1.0D0/(RHOCR*OMEGDD) VPS=VD+XPS*(VDD-VD) F80C06=REAL(VPS) RETURN END IF IF (S.GT.SDD) THEN CALL S12C06(TAUS,OMEGDD,CPDD) TAU0=TAUS S0=SDD CP0=CPDD ELSE IF (S.LT.SD) THEN CALL S12C06(TAUS,OMEGD,CPD) TAU0=TAUS S0=SD CP0=CPD END IF END IF TAU1=TAU0*EXP((S-S0)/CP0) CALL S4C06(PI,TAU1,OMEGA1) IF (OMEGA1.LE.-1.0D+10) RETURN CALL S8C06(TAU1,OMEGA1,S1) TH0=TAU0 TH1=TAU1 CALL S40C06(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=THW CALL S4C06(PI,TAUW,OMEGAW) IF (OMEGAW.LE.-1.0D+10) RETURN VPS=1.0D0/(RHOCR*OMEGAW) F80C06=REAL(VPS) RETURN END REAL FUNCTION F82C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F82C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F82C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S14C06(TAU,OMEGA,AKPT) F82C06=REAL(AKPT) RETURN END REAL FUNCTION F83C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F83C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F83C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S13C06(TAU,OMEGA,WPT) F83C06=REAL(WPT) RETURN END REAL FUNCTION F90C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F90C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F90C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S15C06(TAU,OMEGA,BSPT) F90C06=REAL(BSPT*1.0D+05) RETURN END REAL FUNCTION F91C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F91C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F91C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S16C06(TAU,OMEGA,BTPT) F91C06=REAL(BTPT*1.0D+05) RETURN END REAL FUNCTION F92C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F92C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F92C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S17C06(TAU,OMEGA,BPPT) F92C06=REAL(BPPT) RETURN END REAL FUNCTION F93C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F93C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F93C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S18C06(TAU,OMEGA,BVPT) F93C06=REAL(BVPT) RETURN END REAL FUNCTION F94C06(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=4.8714D+06,TCR=305.33D+00) F94C06=-1.0E+20 CALL S91C06(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F94C06=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=(DBLE(FT)+273.15D0)/TCR CALL S4C06(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S19C06(TAU,OMEGA,AJTPT) F94C06=REAL(AJTPT*1.0D+05) RETURN END C F98 + PSEUDO BOILING POINT AT P FUNCTION F98C06(P) DIMENSION T(2),C(2),TL(3),TR(3),CL(3),CR(3) P1=F21C06('P') PP=ABS((P-P1)/P1) T1=F21C06('T') IF (PP.LT.1.0E-5) THEN F98C06=T1 RETURN ENDIF IF (P.LT.P1.OR.P.GT.200.001D00) THEN F98C06=-1.0E+20 RETURN ENDIF T2=0.85*T1 P2=F30C06(T2) TM0=T1+(T1-T2)*(P-P1)/(P1-P2) C-----TM0 KINJICHI TC=T1 T(1)=TM0 150 EPS=1.0E-6 DEL=T1*0.05E00 IREP=0 IREM=5000 KCONT=0 ICONT=0 C(1)=F18C06(P,T(1)) T(2)=T(1)-DEL C(2)=F18C06(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)=F18C06(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=-F18C06(P,TT) 3000 CONV=ABS((CC-C(2))/CC) IF(CONV.LT.EPS) THEN F98C06=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=-F18C06(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=F18C06(P,TA) CB=F18C06(P,TB) 6050 KCONT=KCONT+1 IF (KCONT.GT.IREM) GO TO 8000 DELT=ABS((TA-TB)/TA) IF (DELT.LT.EPS) GO TO 7000 CC=F18C06(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)=F18C06(P,TL(2)) CR(2)=F18C06(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=F18C06(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=F18C06(P,TB) ELSE TB=TR(MR+1) CB=CR(MR+1) ENDIF 6135 TC=TR(MR) CC=CR(MR) ENDIF GO TO 6050 7000 F98C06=TC RETURN 8000 F98C06=-1.0E+10 RETURN END SUBROUTINE S1C06(K,TAU,OMEGA,AK) ***** SUBROUTINE TO CALCULATE A0, A1, A2,A3, A4 AND A5 DEFINED BY ***** EQ.(3.7) ON P.35 AND A6 ***** A6=SIGMA(SIGMA(1/I)*B(I,J)*OMEGAI/TAUJ) ***** INPUT1 : K = 0,1,2,3,4, 5 OR 6 ***** INPUT2 : TAU = T/TCR ***** INPUT3 : OMEGA = RHO/RHOCR ***** OUTPUT : AK = A0, A1, A2, A3, A4, A5, A6 IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION B(1:10,0:7),C(1:10,0:7),JS(1:10) DATA (JS(I),I=1,10)/ 7, 6, 5, 5, 4, 4, 3, 3, 2, 1/ DATA (B(1,J),J=0,7)/+0.6523112D+00,-0.1420959D+01,-0.8281694D+00, - +0.9628378D+00,-0.4873274D+00,-0.1120178D+00,+0.4053669D-01, - +0.6643199D-02/ DATA (B(2,J),J=0,6)/-0.1717300D+00,+0.1342033D+01,-0.5419403D+00, - -0.3585280D+00,+0.3413308D+00,-0.1419773D+00,-0.8327400D-01/ DATA (B(3,J),J=0,5)/+0.1816776D+00,-0.1159004D+01,+0.6856036D-01, - +0.4834712D+00,+0.3294358D+00,+0.2712144D+00/ DATA (B(4,J),J=0,5)/+0.7302986D-01,+0.6713792D+00,-0.4315169D+00, - -0.1305074D+00,-0.2605725D+00,-0.1298954D-01/ DATA (B(5,J),J=0,4)/-0.3324578D-01,+0.8053416D-01,+0.7465193D-01, - +0.5459819D-01,-0.3786991D-01/ DATA (B(6,J),J=0,4)/-0.1392303D+00,-0.2013963D-01,-0.9262326D-01, - -0.3878733D-01,+0.1381212D-01/ DATA (B(7,J),J=0,3)/+0.1066015D+00,-0.2039723D-01,+0.5628173D-01, - +0.7784005D-02/ DATA (B(8,J),J=0,3)/-0.2233251D-01,+0.2384036D-01,+0.2426002D-02, - -0.3414020D-02/ DATA (B(9,J),J=0,2)/-0.1016497D-01,-0.1018997D-01,+0.2835872D-02/ DATA (B(10,J),J=0,1)/+0.4957046D-02,-0.1722518D-02/ * IF (K.EQ.0) THEN DO 10 I=1,10 DO 10 J=0,JS(I) C(I,J)=B(I,J) 10 CONTINUE ELSE IF (K.EQ.1) THEN DO 11 I=1,10 DO 11 J=0,JS(I) C(I,J)=(DBLE(I)+1.0D0)*B(I,J) 11 CONTINUE ELSE IF (K.EQ.2) THEN DO 12 I=1,10 DO 12 J=0,JS(I) C(I,J)=(1.0D0-DBLE(J))*B(I,J) 12 CONTINUE ELSE IF (K.EQ.3) THEN DO 13 I=1,10 DO 13 J=0,JS(I) C(I,J)=((DBLE(I)+DBLE(J))/DBLE(I))*B(I,J) 13 CONTINUE ELSE IF (K.EQ.4) THEN DO 14 I=1,10 DO 14 J=0,JS(I) C(I,J)=((DBLE(J)-1.0D0)/DBLE(I))*B(I,J) 14 CONTINUE ELSE IF (K.EQ.5) THEN DO 15 I=1,10 DO 15 J=0,JS(I) C(I,J)=-((DBLE(J)-1.0D0)*DBLE(J)/DBLE(I))*B(I,J) 15 CONTINUE ELSE IF (K.EQ.6) THEN DO 16 I=1,10 DO 16 J=0,JS(I) C(I,J)=(1.0D0/DBLE(I))*B(I,J) 16 CONTINUE ELSE AK=-1.0D+30 RETURN END IF * SUMA=0.0D0 DO 21 I=10,1,-1 STAU=0.0D0 DO 20 J=JS(I),0,-1 STAU=C(I,J)+STAU/TAU 20 CONTINUE SUMA=STAU+SUMA*OMEGA 21 CONTINUE * AK=SUMA*OMEGA RETURN END SUBROUTINE S2C06(TAU,OMEGD,OMEGDD,PIS) ***** SUBROUTINE TO FIND THE SATURATED LIQUID AND VAPOR DENSITIES ***** CORRESPONDING TO THE INPUT TEMPERATURE BY SOLVING THE MAXWELL ***** RELATION USING NEWTON METHOD AND TO CALCULATE ITS SATURATION ***** PRESSURE FROM THE EQUATION OF STATE, EQ.(3.7) ON P.35 ***** INPUT : TAU = T/TCR ***** OUTPUT1 : OMEGD =RHOD/RHOCR ***** OUTPUT2 : OMEGDD = RHODD/RHOCR ***** OUTPUT3 : PIS = PS/PCR IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, - PCR=4.8714D+06,TCR=305.33D+00,RHOCR=204.4571662D+00) PARAMETER(ITMAX=10000) DATA D1/ 0.7386784D+01/,D2/-0.9743759D+02/,D3/ 0.7440263D+03/, - D4/-0.3221553D+04/,D5/ 0.8705971D+04/,D6/-0.1513969D+05/, - D7/ 0.1687902D+05/,D8/-0.1159862D+05/,D9/ 0.4442097D+04/, - D10/-0.7183749D+03/ DATA DD1/ 0.3294644D+05/,DD2/-0.2078230D+06/, - DD3/ 0.5248058D+06/,DD4/-0.6633872D+06/, - DD5/ 0.4198705D+06/,DD6/-0.1064971D+06/ DATA EE1/-0.6334833D+01/,EE2/ 0.6719200D+02/, - EE3/-0.4169151D+03/,EE4/ 0.1305206D+04/, - EE5/-0.2268328D+04/,EE6/ 0.2061080D+04/, - EE7/-0.7768386D+03/ * OMEGD=-1.0D+20 OMEGDD=-1.0D+20 PIS=-1.0D+20 IF (TAU.LE.0.0D0) RETURN * OMEGD=1.0D0 OMEGDD=1.0D0 PIS=1.0D0 IF (ABS(TAU-1.0D0).LE.1.0D-06) RETURN *** INITIAL GUESSES FOR OMD AND OMDD X=ABS(1.0D0-TAU) QX=X**(1.0D0/3.0D0) OMDW=1.0D0+(D1+(D2+(D3+(D4+(D5+(D6+(D7+(D8+(D9+D10*QX)*QX) - *QX)*QX)*QX)*QX)*QX)*QX)*QX)*QX IF (TAU.LE.0.5568D0) THEN OMDDW=EXP((DD1+(DD2+(DD3+(DD4+(DD5+DD6*QX)*QX)*QX)*QX)*QX)*QX) ELSE IF ((TAU.GT.0.5568D0).AND.(TAU.LE.0.9989D0)) THEN OMDDW=EXP((EE1+(EE2+(EE3+(EE4+(EE5+(EE6+EE7*QX)*QX)*QX)*QX) - *QX)*QX)*QX) END IF IF ((TAU.GT.0.9989D0).AND.(TAU.LT.0.99949D0)) THEN OMDW= 1.0D0/(1.0D0-6.3150852D0*SQRT(X)) OMDDW=1.0D0/(1.0D0+9.5053177D0*SQRT(X)) ELSE IF ((TAU.GE.0.99949D0).AND.(TAU.LT.0.99975D0)) THEN OMDW= 1.0D0/(1.0D0-8.793211D0*SQRT(X)) OMDDW=1.0D0/(1.0D0+12.82828D0*SQRT(X)) ELSE IF ((TAU.GE.0.99975D0).AND.(TAU.LT.0.999846D0)) THEN OMDW= 1.0D0/(1.0D0-12.32222D0*SQRT(X)) OMDDW=1.0D0/(1.0D0+17.70592D0*SQRT(X)) ELSE IF ((TAU.GE.0.999846D0).AND.(TAU.LT.1.0D0)) THEN OMDW= 1.0D0/(1.0D0-15.53562D0*SQRT(X)) OMDDW=1.0D0/(1.0D0+22.19489D0*SQRT(X)) END IF * *** NEWTON METHOD IF (TAU.GT.0.9989D0) THEN EPS=1.0D-10 ELSE EPS=1.0D-14 END IF * DO 1 IT=1,ITMAX CALL S1C06(0,TAU,OMDW,A0D) CALL S1C06(0,TAU,OMDDW,A0DD) CALL S1C06(1,TAU,OMDW,A1D) CALL S1C06(1,TAU,OMDDW,A1DD) CALL S1C06(6,TAU,OMDW,A6D) CALL S1C06(6,TAU,OMDDW,A6DD) F=OMDW*(1.0D0+A0D)-OMDDW*(1.0D0+A0DD) G=LOG(OMDW/OMDDW)+A0D-A0DD+A6D-A6DD DOMDW= (OMDW/(OMDW-OMDDW))*(F-G*OMDDW)/(1.0D0+A1D) DOMDDW=(OMDDW/(OMDW-OMDDW))*(F-G*OMDW)/(1.0D0+A1DD) OMDW=OMDW-DOMDW OMDDW=OMDDW-DOMDDW IF ((ABS(DOMDW).LT.EPS).AND.(ABS(DOMDDW).LT.EPS)) THEN OMEGD=OMDW OMEGDD=OMDDW PIS=((GASCON*RHOCR*TCR/PCR)*TAU) - *(OMDW*OMDDW/(OMDW-OMDDW))*(LOG(OMDW/OMDDW)+A6D-A6DD) RETURN END IF 1 CONTINUE * OMEGD=-1.0D+10 OMEGDD=-1.0D+10 PIS=-1.0D+10 RETURN END SUBROUTINE S3C06(PI,OMEGD,OMEGDD,TAUS) ***** SUBROUTINE TO FIND THE SATURATED LIQUID, VAPOR DENSITIESAND ***** SATURATION TEMPERATURE CORRESPONDING TO THE INPUT PRESSURE ***** USING AN ITERATIVE METHOD ***** INPUT : PI = P/PCR ***** OUTPUT1 : OMEGD =RHOD/RHOCR ***** OUTPUT2 : OMEGDD = RHODD/RHOCR ***** OUTPUT3 : TAUS = TS/TCR IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, TTR=90.348D+00, - PCR=4.8714D+06,TCR=305.33D+00,RHOCR=204.4571662D+00) PARAMETER(EPS1=1.0D-12, IT1MAX=10000, IT2MAX=10000) DATA CA/-11.38996D+00/, CB/18.84523D+00/,CC/-7.635425D+00/, - CD/5.428443D+00/, CE/-1.362327D+00/,CF/0.7692493D+00/, - CEPS/1.3D+00/ * OMEGD=-1.0D+20 OMEGDD=-1.0D+20 TAUS=-1.0D+20 IF (PI.LE.0.0D0) RETURN * OMEGD=1.0D0 OMEGDD=1.0D0 TAUS=1.0D0 IF (ABS(PI-1.0D0).LE.7.0D-04) RETURN *** INITIAL GUESSES FOR TAUS *** TAU0=1.0D0/(1.0D0-(LOG(PI))/6.690605D0) RTAU=TTR/TCR PBAR=PI*PCR*1.0D-05 DO 1 IT1=1,IT1MAX X=ABS((1.0D0-(RTAU/TAU0))/(1.0D0-RTAU)) U=ABS((TAU0-RTAU)/(1.0D0-RTAU)) U1=ABS(1.0D0-U) F=CA+CB*X+(CC+(CD+CE*U)*U)*U+CF*U*(U1**CEPS)-LOG(PBAR) DFDT=(CB*RTAU/TAU0**2+CC+2.0D0*CD*U+3.0D0*CE*U**2 - +CF*(U1**(CEPS-1.0D0))*(1.0D0-(1.0D0+CEPS)*U)) - /(1.0D0-RTAU) TAU0=TAU0-F/DFDT IF (ABS(-F/(DFDT*TAU0)).LT.EPS1) GO TO 10 1 CONTINUE GO TO 1000 10 TAUSW=TAU0 *** ITERATIVE METHOD *** EPS2=1.0D-12 * ZCR=PCR/(RHOCR*GASCON*TCR) DO 2 IT2=1,IT2MAX CALL S2C06(TAUSW,OMDW,OMDDW,PISW) DPI=PI-PISW IF (ABS(DPI/PI).LE.EPS2) THEN OMEGD=OMDW OMEGDD=OMDDW TAUS=TAUSW RETURN END IF CALL S1C06(4,TAUSW,OMDW,A4DW) CALL S1C06(4,TAUSW,OMDDW,A4DDW) DPIDTA=(1.0D0/ZCR)*(OMDW*OMDDW/(OMDW-OMDDW)) - *(LOG(OMDW/OMDDW)+A4DDW-A4DW) DTAUSW=DPI/DPIDTA TAUSW=TAUSW+DTAUSW 2 CONTINUE 1000 OMEGD=-1.0D+10 OMEGDD=-1.0D+10 TAUS=-1.0D+10 RETURN END SUBROUTINE S4C06(PI,TAU,OMEGA) ***** SUBROUTINE TO FIND THE DENSITY CORRESPONDING TO THE INPUT ***** PRESSURE AND TEMPERATURE BY SOLVING THE EQUATION OF STATE, ***** EQ.(3.7) USING NEWTON METHOD ***** INPUT1 : PI = P/PCR ***** INPUT2 : TAU = T/TCR ***** OUTPUT : OMEGA = RHO/RHOCR IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, - PCR=4.8714D+06,TCR=305.33D+00,RHOCR=204.4571662D+00) PARAMETER(EPS=1.0D-12,ITMAX=10000) DATA C1/-0.2091683D+01/,C2/+0.9990607D+00/, - C3/+0.3234132D-03/,C4/-0.7919523D-03/, - C5/-0.3856794D-02/,C6/-0.1581330D-02/, - C7/+0.1213998D-03/,C8/+0.7409514D-04/ OMEGA=1.0D0 IF ((ABS(PI-1.0D0).LE.1.0D-05).AND.(ABS(TAU-1.0D0).LE.1.0D-05)) - RETURN ZCR=PCR/(RHOCR*GASCON*TCR) *** INITIAL GUESS FOR OMEGA *** IF ((PI.GE.1.0D0).AND.(PI.LE.9.443)) THEN PI1=ABS(PI-1.0D0) TAUCR=1.0D0+(0.153668D0+(-0.161860D-02+0.16481D-03*PI1) - *PI1)*PI1 IF (TAU.GT.TAUCR) THEN PIL=LOG(PI) OM=EXP(C1+(C2+(C3+(C4+(C5+(C6+(C7+C8*PIL)*PIL)*PIL)*PIL) - *PIL)*PIL)*PIL) GO TO 5 END IF END IF * IF (PI.GE.1.0D0) THEN IF (TAU.GE.0.3334D0) THEN OM=EXP(0.124832D+01+(-0.462670D-01+(-0.643043D+00 - +(0.512701D+00+(-0.201559D+00+0.310794D-01*TAU)*TAU) - *TAU)*TAU)*TAU) GO TO 5 ELSE OM=3.25D0 GO TO 5 END IF END IF * IF ((PI.LT.1.0D0).AND.(TAU.GT.1.0D0)) THEN OM=ZCR*PI/TAU ELSE CALL S3C06(PI,OMD,OMDD,TAUS) IF (OMD.LE.-1.0D+10) GO TO 1000 IF ((PI.LE.1.0D0).AND.(TAU.GT.TAUS)) THEN IF (ABS(TAUS/TAU-1.0D0).LT.1.0D-05) THEN OMEGA=OMDD RETURN END IF OM=OMDD ELSE IF ((PI.LE.1.0D0).AND.(TAU.LT.TAUS)) THEN IF (ABS(TAU/TAUS-1.0D0).LT.1.0D-05) THEN OMEGA=OMD RETURN END IF IF (TAU.GT.0.720D0) THEN OM=2.5D0 ELSE OM=EXP( 0.4472284D+03+(-0.7547690D+04+( 0.5553781D+05 - +(-0.2332401D+06+( 0.6162415D+06+(-0.1062743D+07 - +( 0.1197017D+07+(-0.8497473D+06+( 0.3452629D+06 - -0.6122880D+05*TAU)*TAU)*TAU)*TAU)*TAU)*TAU)*TAU) - *TAU)*TAU) END IF END IF END IF *** NEWTON METHOD *** 5 ZCRPT=ZCR*PI/TAU DO 10 IT=1,ITMAX CALL S1C06(0,TAU,OM,A0) CALL S1C06(1,TAU,OM,A1) DOM=(-OM*(1.0D0+A0)+ZCRPT)/(1.0D0+A1) IF (ABS(DOM/OM).LT.EPS) THEN OMEGA=OM RETURN END IF OM=OM+DOM 10 CONTINUE 1000 OMEGA=-1.0D+10 RETURN END SUBROUTINE S5C06(PI,OMEGA,TAU) ***** SUBROUTINE TO FIND THE TEMPERATURE CORRESPONDING TO THE INPUT ***** PRESSURE AND DENSITY BY SOLVING THE EQUATION OF STATE, ***** EQ.(3.7) USING NEWTON METHOD ***** INPUT1 : PI = P/PCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : TAU = T/TCR IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, - PCR=4.8714D+06,TCR=305.33D+00,RHOCR=204.4571662D+00) PARAMETER(EPS=1.0D-12,ITMAX=10000) TAU=1.0D0 IF ((ABS(PI-1.0D0).LE.1.0D-05).AND.(ABS(OMEGA-1.0D0).LE.1.0D-05)) - RETURN ZCR=PCR/(RHOCR*GASCON*TCR) *** INITIAL GUESS FOR TAU *** IF (PI.GE.1.0D0) THEN IF (OMEGA.LE.1.0D0) THEN TAUW=700.0D0/TCR ELSE IF ((OMEGA.GT.1.0D0).AND.(OMEGA.LE.2.930)) THEN TAUW=1.0D0+(0.153668D+00+(-0.161860D-02+0.16481D-03*PI) - *PI)*PI ELSE IF ((OMEGA.GT.2.93D0).AND.(OMEGA.LE.3.1D0)) THEN TAUW=140.0D0/TCR ELSE IF (OMEGA.GT.3.1D0) THEN TAUW=120.0D0/TCR END IF ELSE IF (OMEGA.GT.1.0D0) THEN OM=OMEGA TAUW=EXP(0.23471017D0+(-0.60873692D0+(0.50706395D0 - -0.14272629D0*OM)*OM)*OM) ELSE TAUW=ZCR*PI/OMEGA END IF END IF *** NEWTON METHOD *** ZCRPOM=ZCR*PI/OMEGA DO 1 IT=1,ITMAX CALL S1C06(0,TAUW,OMEGA,A0) CALL S1C06(2,TAUW,OMEGA,A2) DTAU=(-TAUW*(1.0D0+A0)+ZCRPOM)/(1.0D0+A2) IF (ABS(DTAU/TAUW).LT.EPS) THEN TAU=TAUW RETURN END IF TAUW=TAUW+DTAU 1 CONTINUE TAU=-1.0D+10 RETURN END SUBROUTINE S6C06(TAU,OMEGA,PI) ***** SUBROUTINE TO CALCULATE THE PRESSURE CORRESPONDING TO THE INPUT ***** TEMPERATURE AND DENSITY FROM THE EQUATION OF STATE, EQ.(3.7) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : PI = P/PCR IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, - PCR=4.8714D+06,TCR=305.33D+00,RHOCR=204.4571662D+00) ZCR=PCR/(RHOCR*GASCON*TCR) CALL S1C06(0,TAU,OMEGA,A0) PI=(1.0D0/ZCR)*OMEGA*TAU*(1.0D0+A0) RETURN END SUBROUTINE S7C06(IDSHCP,TAU,SHCP) ***** SUBROUTINE TO CALCULATE THE ENTROPY, ENTHALPY AND ISOBARIC ***** SPECIFIC HEAT IN THE IDEAL GAS STATE, EQS.(3.1)-(3.5) ***** INPUT1 : IDSHCP = 1 FOR OUTPUT = S(ENTROPY) ***** IDSHCP = 2 FOR OUTPUT = H(ENTHALPY) ***** IDSHCP = 3 FOR OUTPUT = CP(ISOBARIC SPECIFIC HEAT) ***** INPUT2 : TAU = T/TCR ***** OUTPUT : SHCP = S IN J/(KG*K) FOR IDSHCP=1 ***** SHCP = H IN J/KG FOR IDSHCP=2 ***** SHCP = CP IN J/(KG*K) FOR IDSHCP=3 IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, TCR=305.33D+00) PARAMETER(T0=100.0D+00, S0=0.0D0, H0=968.426D+03, - S00=22.1158D+00, H00=4.0670D+00) DIMENSION A(0:6),AS(0:6),AH(0:6),B(1:5),BS(1:5),BH(1:5) DATA (A(I),I=0,6) /+0.68120976D+02,-0.30634058D+02, - +0.95275029D+01,-0.16947102D+01,+0.17630585D+00, - -0.99545402D-02,+0.23536430D-03/ DATA (B(I),I=1,5) /-0.87407084D+02,+0.78481374D+02, - -0.44865859D+02,+0.14654346D+02,-0.20518393D+01/ T=TAU*TCR TH=T/T0 * IF (IDSHCP.EQ.1) THEN SAS=0.0D0 SASTH=0.0D0 DO 110 J=6,1,-1 AS(J)=A(J)/DBLE(J) SAS=SAS+AS(J) SASTH=SASTH*TH+AS(J) 110 CONTINUE SASTH=SASTH*TH SBS=0.0D0 SBSTH=0.0D0 DO 112 J=5,1,-1 BS(J)=B(J)/DBLE(J) SBS=SBS+BS(J) SBSTH=SBSTH/TH+BS(J) 112 CONTINUE SBSTH=SBSTH/TH SID=GASCON*(SASTH-SBSTH+A(0)*LOG(TH)+SBS-SAS+S00)+S0 SHCP=SID * ELSE IF (IDSHCP.EQ.2) THEN SAH=0.0D0 SAHTH=0.0D0 DO 120 J=6,0,-1 AH(J)=A(J)/(DBLE(J)+1.0D0) SAH=SAH+AH(J) SAHTH=SAHTH*TH+AH(J) 120 CONTINUE SBH=0.0D0 SBHTH=0.0D0 DO 122 J=5,2,-1 BH(J)=B(J)/(DBLE(J)-1.0D0) SBH=SBH+BH(J) SBHTH=SBHTH/TH+BH(J) 122 CONTINUE SBHTH=SBHTH/TH**2 HID=GASCON*T*(SAHTH-SBHTH - +(1.0D0/TH)*(B(1)*LOG(TH)+SBH-SAH+H00))+H0 SHCP=HID * ELSE IF (IDSHCP.EQ.3) THEN SACPTH=0.0D0 DO 130 J=6,0,-1 SACPTH=SACPTH*TH+A(J) 130 CONTINUE SBCPTH=0.0D0 DO 132 J=5,1,-1 SBCPTH=SBCPTH/TH+B(J) 132 CONTINUE SBCPTH=SBCPTH/TH CPID=GASCON*(SACPTH+SBCPTH) SHCP=CPID END IF RETURN END SUBROUTINE S8C06(TAU,OMEGA,S) ***** SUBROUTINE TO CALCULATE THE ENTROPY BY EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : S = ENTROPY IN (J/(KG*K)) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, - PST=0.101325D+06,TCR=305.33D+00,RHOCR=204.4571662D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S7C06(1,TAUW,S0) CALL S1C06(4,TAUW,OMEGAW,A4) S=GASCON*(A4-LOG(OMEGAW/(PST/(RHOCR*GASCON*TAUW*TCR))))+S0 RETURN END SUBROUTINE S9C06(TAU,OMEGA,H) ***** SUBROUTINE TO CALCULATE THE ENTHALPY BY EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : H = ENTHALPY IN (J/KG) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, TCR=305.33D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S7C06(2,TAUW,H0) CALL S1C06(3,TAUW,OMEGAW,A3) H=GASCON*(TAUW*TCR)*A3+H0 RETURN END SUBROUTINE S10C06(TAU,OMEGA,U) ***** SUBROUTINE TO CALCULATE THE INTERNAL ENERGY BY EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : U = INTERNAL ENERGY IN (J/KG) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, TCR=305.33D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S9C06(TAUW,OMEGAW,H) CALL S1C06(0,TAUW,OMEGAW,A0) U=H-GASCON*(TAUW*TCR)*(1.0D0+A0) RETURN END SUBROUTINE S11C06(TAU,OMEGA,CV) ***** SUBROUTINE TO CALCULATE THE ISOCHORIC SPECIFIC HEAT, EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : CV = ISOCHORIC SPECIFIC HEAT IN J/(KG*K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S7C06(3,TAUW,CP0) CALL S1C06(5,TAUW,OMEGAW,A5) CV0=CP0-GASCON CV=CV0+GASCON*A5 RETURN END SUBROUTINE S12C06(TAU,OMEGA,CP) ***** SUBROUTINE TO CALCULATE THE ISOBARIC SPECIFIC HEAT, EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : CP = ISOBARIC SPECIFIC HEAT IN J/(KG*K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S11C06(TAUW,OMEGAW,CV) CALL S1C06(1,TAUW,OMEGAW,A1) CALL S1C06(2,TAUW,OMEGAW,A2) CP=CV+GASCON*(1.0D0+A2)**2/ABS(1.0D0+A1) RETURN END SUBROUTINE S13C06(TAU,OMEGA,W) ***** SUBROUTINE TO CALCULATE THE VELOCUTY OF SOUND, EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : W = VELOCITY OF SOUND IN M/S IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, TCR=305.33D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S11C06(TAUW,OMEGAW,CV) CALL S12C06(TAUW,OMEGAW,CP) CALL S1C06(1,TAUW,OMEGAW,A1) W0=SQRT(GASCON*(TAUW*TCR)*(CP/CV)) W=W0*SQRT(ABS(1.0D0+A1)) RETURN END SUBROUTINE S14C06(TAU,OMEGA,AK) ***** SUBROUTINE TO CALCULATE THE ADIABATIC INDEX, EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : AK = ADIABATIC INDEX IN (-) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C PARAMETER(GASCON=276.507D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S11C06(TAUW,OMEGAW,CV) CALL S12C06(TAUW,OMEGAW,CP) CALL S1C06(0,TAUW,OMEGAW,A0) CALL S1C06(1,TAUW,OMEGAW,A1) AK0=CP/CV AK=AK0*ABS((1.0D0+A1)/(1.0D0+A0)) RETURN END SUBROUTINE S15C06(TAU,OMEGA,BS) ***** SUBROUTINE TO CALCULATE THE ADIABATIC COMPRESSIBILITY ***** EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : BS = ADIABATIC COMPRESSIBILITY IN (1/PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, - TCR=305.33D+00, RHOCR=204.4571662D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S11C06(TAUW,OMEGAW,CV) CALL S12C06(TAUW,OMEGAW,CP) CALL S1C06(1,TAUW,OMEGAW,A1) BS=(CV/CP)/((TAUW*TCR)*(OMEGAW*RHOCR)*GASCON*(1.0D0+A1)) RETURN END SUBROUTINE S16C06(TAU,OMEGA,BT) ***** SUBROUTINE TO CALCULATE THE ISOTHERMAL COMPRESSIBILITY ***** EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : BT = ISOTHERMAL COMPRESSIBILITY IN (1/PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=276.507D+00, - TCR=305.33D+00, RHOCR=204.4571662D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S1C06(1,TAUW,OMEGAW,A1) BT=1.0D0/((TAUW*TCR)*(OMEGAW*RHOCR)*GASCON*(1.0D0+A1)) RETURN END SUBROUTINE S17C06(TAU,OMEGA,BP) ***** SUBROUTINE TO CALCULATE THE VOLUMETRIC EXPANSION COEFF. ***** EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : BP = VOLUMETRIC EXPANSION COEFF. IN (1/K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(TCR=305.33D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S1C06(1,TAUW,OMEGAW,A1) CALL S1C06(2,TAUW,OMEGAW,A2) BP=(1.0D0/(TAUW*TCR))*((1.0D0+A2)/(1.0D0+A1)) RETURN END SUBROUTINE S18C06(TAU,OMEGA,BV) ***** SUBROUTINE TO CALCULATE THE PRESSURE COEFF. ***** EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : BV = PRESSURE COEFF. IN (1/K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(TCR=305.33D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S1C06(0,TAUW,OMEGAW,A0) CALL S1C06(2,TAUW,OMEGAW,A2) BV=(1.0D0/(TAUW*TCR))*((1.0D0+A2)/(1.0D0+A0)) RETURN END SUBROUTINE S19C06(TAU,OMEGA,AJT) ***** SUBROUTINE TO CALCULATE THE JOU;LETHOMSON COEFF. ***** EQ.(2.4) ***** INPUT1 : TAU = T/TCR ***** INPUT2 : OMEGA = RHO/RHOCR ***** OUTPUT : AJT = JOULE-THOMSON COEFF IN (K/PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(RHOCR=204.4571662D+00) TAUW=TAU OMEGAW=OMEGA IF ((ABS(TAUW-1.0D0).LE.1.0D-05) - .AND.(ABS(OMEGAW-1.0).LE.1.0D-05)) THEN TAUW=1.0D0 OMEGAW=1.0D0 END IF CALL S1C06(1,TAUW,OMEGAW,A1) CALL S1C06(2,TAUW,OMEGAW,A2) CALL S12C06(TAUW,OMEGAW,CP) AJT=(1.0D0/((OMEGAW*RHOCR)*CP))*((A2-A1)/(1.0D0+A1)) RETURN END SUBROUTINE S40C06(IHS,PI,Z,TH1,TH2,Z1,Z2,TH) *** IHS : IHS=1 FOR ARGUMENTS(FP,FH):F64 *** IHS=2 FOR ARGUMENTS(FP,FS):F65, F71, F79 AND F80 *** TH : OUTPUT=T/TCR=TAU IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(EPS=1.0D-07,ITMAX=10000) DO 1 IT=1,ITMAX DTH=(Z-Z2)*(TH2-TH1)/(Z2-Z1) IF (ABS(DTH/TH1).LT.EPS) THEN TH=TH2 RETURN END IF TH1=TH2 Z1=Z2 TH2=TH2+DTH IF (TH2.LE.0.298D0) THEN ****** LOWER LIMIT OF TAU FOR ETHANE TH2=0.298D0 END IF CALL S4C06(PI,TH2,OMEGA2) IF (OMEGA2.LE.-1.0D+10) GO TO 1000 IF (IHS.EQ.1) THEN CALL S9C06(TH2,OMEGA2,Z2) ELSE IF (IHS.EQ.2) THEN CALL S8C06(TH2,OMEGA2,Z2) END IF 1 CONTINUE 1000 TH=-1.0D+10 RETURN END SUBROUTINE S90C06(FP,FT,ILL) ***** SUBROUTINE TO CHECK IF THE RANGE OF ARGUMENTS:(P,T) IS PROPER***** ***** FOR F25C06:HPT, F35C06:SPT, F44C06:UPT, F51C06:VPT. ****** IMPLICIT REAL(F) ILL=10000 IF ((FP.LT.0.999).OR.(FP.GT.800.1)) RETURN FTMAX=426.95 IF ((FP.GE.0.999).AND.(FP.LE.40.421)) THEN FTMIN=-182.25 ELSE IF ((FP.GT.40.421).AND.(FP.LE.800.1)) THEN FTMIN=F69C06(FP)-0.1 END IF IF ((FT.GT.FTMIN).AND.(FT.LT.FTMAX)) THEN ILL=0 END IF RETURN END SUBROUTINE S91C06(FP,FT,ILL) ***** SUBROUTINE TO CHECK IF THE RANGE OF ARGUMENTS:(P,T) IS PROPER***** ***** FOR F18C06:CPPT, F77C06:CVPT, F82C06:AKPT, F83C06:WPT, ****** ***** F90C06:BSPT, F91C06:BTPT, F92C06:BPPT, F93C06:BVPT, F94C06:AJTPT** IMPLICIT REAL(F) ILL=10000 IF ((FP.LT.0.999).OR.(FP.GT.800.1)) RETURN FTMIN=-153.25 FTMAX=426.95 IF ((FT.GT.FTMIN).AND.(FT.LT.FTMAX)) THEN ILL=0 END IF RETURN END SUBROUTINE S97C06(NFUN) *** LEVEL 1 ERROR MESSAGE *** CHARACTER NFUN*6, MSG*125 INTEGER KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN MSG='**** NO CONVERGENCE AT '//NFUN//' FOR ETAHNE ****' WRITE(6,1000) MSG 1000 FORMAT(1H ,5X,A) END IF RETURN END SUBROUTINE S98C06(IARG,ARG1,ARG2,NARG1,NARG2,NFUN) *** IARG=1 FOR ONE ARGUMENT (SECOND ARUMENT IS DUMMY) *** IARG=2 FOR TW0 ARUMENTS *** LEVEL 2 ERROR MESSAGE *** CHARACTER NFUN*6, NARG1*1,NARG2*1 INTEGER KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN IF (IARG.EQ.1) THEN WRITE(6,2010) NFUN,NARG1,ARG1 ELSE IF (IARG.EQ.2) THEN WRITE(6,2020) NFUN,NARG1,ARG1,NARG2,ARG2 END IF END IF 2010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ETHANE', - ' WHEN ',A1,' =', 1PE14.7,' ****') 2020 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR ETHANE', - ' WHEN ',A1,' =',1PE14.7,' AND ',A1,' =',1PE14.7,' ****') RETURN END SUBROUTINE S99C06(NFUN) *** LEVEL 3 ERROR MESSAGE *** CHARACTER NFUN*6, MSG*125 INTEGER KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN MSG='**** FUNCTION '//NFUN//' UNAVAILABLE FOR ETHANE ****' WRITE(6,3000) MSG 3000 FORMAT(1H ,5X,A) END IF RETURN END REAL FUNCTION T90(T) * Conversion of temperature scal from IPTS-68 to ITS-90 * Range -200C<=T68<=630C anf T68>=1064 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 * Check of argument range 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 REAL FUNCTION T68(T) * Conversion of temperature scal from ITS-90 to IPTS-68 * Range -200C<=T68<=630C anf T68>=1064 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 * Subroutine Program Specifying KPA and MESS * MS-FORTRAN and MS-C Mixed Langage Programing * SUBROUTINE KPAMES(KPAC,MESSC) COMMON/UNIT/KPA,MESS INTEGER KPAC INTEGER MESSC KPA=KPAC MESS=MESSC RETURN END