*====+==p10f2.f ===FLUORINE===(1996.07.15 by T.SHIGECHI)==========* 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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(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 S99C09(FUN) WTDD=-1.0E+30 RETURN END C C==================================================================== *------------------------------------------------- F1C09 = AIPPT REAL FUNCTION AIPPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C09('AIPPT') AIPPT=-1.0E+30 RETURN END *------------------------------------------------- F2C09 = ALAPP REAL FUNCTION ALAPP(P) REAL P,PI PI=P CALL S99C09('ALAPP') ALAPP=-1.0E+30 RETURN END *------------------------------------------------- F3C09 = ALAPT REAL FUNCTION ALAPT(T) REAL T,TI TI=T CALL S99C09('ALAPT') ALAPT=-1.0E+30 RETURN END *------------------------------------------------- F4C09 = ALHP REAL FUNCTION ALHP(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) ALHP=F4C09(PI) IF(ALHP.EQ.-1.0E+10) THEN CALL S97C09('ALHP') ELSE IF(ALHP.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','ALHP') END IF RETURN END *------------------------------------------------- F5C09 = ALHT REAL FUNCTION ALHT(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) ALHT=F5C09(TI) IF(ALHT.EQ.-1.0E+10) THEN CALL S97C09('ALHT') ELSE IF(ALHT.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','ALHT') END IF RETURN END *------------------------------------------------- F6C09 = ALMPD REAL FUNCTION ALMPD(P) REAL P,PI PI=P CALL S99C09('ALMPD') ALMPD=-1.0E+30 RETURN END *------------------------------------------------- F7C09 = ALMPDD REAL FUNCTION ALMPDD(P) REAL P,PI PI=P CALL S99C09('ALMPDD') ALMPDD=-1.0E+30 RETURN END *------------------------------------------------- F8C09 = ALMPT REAL FUNCTION ALMPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C09('ALMPT') ALMPT=-1.0E+30 RETURN END *------------------------------------------------- F9C09 = ALMTD REAL FUNCTION ALMTD(T) REAL T,TI TI=T CALL S99C09('ALMTD') ALMTD=-1.0E+30 RETURN END *------------------------------------------------- F10C09 = ALMTDD REAL FUNCTION ALMTDD(T) REAL T,TI TI=T CALL S99C09('ALMTDD') ALMTDD=-1.0E+30 RETURN END *------------------------------------------------- F11C09 = AMUPD REAL FUNCTION AMUPD(P) REAL P,PI PI=P CALL S99C09('AMUPD') AMUPD=-1.0E+30 RETURN END *------------------------------------------------- F12C09 = AMUPDD REAL FUNCTION AMUPDD(P) REAL P,PI PI=P CALL S99C09('AMUPDD') AMUPDD=-1.0E+30 RETURN END *------------------------------------------------- F13C09 = AMUPT REAL FUNCTION AMUPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C09('AMUPT') AMUPT=-1.0E+30 RETURN END *------------------------------------------------- F14C09 = AMUTD REAL FUNCTION AMUTD(T) REAL T,TI TI=T CALL S99C09('AMUTD') AMUTD=-1.0E+30 RETURN END *------------------------------------------------- F15C09 = AMUTDD REAL FUNCTION AMUTDD(T) REAL T,TI TI=T CALL S99C09('AMUTDD') AMUTDD=-1.0E+30 RETURN END *------------------------------------------------- F16C09 = CPPD REAL FUNCTION CPPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) CPPD=F16C09(PI) IF(CPPD.EQ.-1.0E+10) THEN CALL S97C09('CPPD') ELSE IF(CPPD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','CPPD') END IF RETURN END *------------------------------------------------- F17C09 = CPPDD REAL FUNCTION CPPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) CPPDD=F17C09(PI) IF(CPPDD.EQ.-1.0E+10) THEN CALL S97C09('CPPDD') ELSE IF(CPPDD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','CPPDD') END IF RETURN END *------------------------------------------------- F18C09 = CPPT REAL FUNCTION CPPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) CPPT=F18C09(PI,TI) IF(CPPT.EQ.-1.0E+10) THEN CALL S97C09('CPPT') ELSE IF(CPPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','CPPT') END IF RETURN END *------------------------------------------------- F19C09 = CPTD REAL FUNCTION CPTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) CPTD=F19C09(TI) IF(CPTD.EQ.-1.0E+10) THEN CALL S97C09('CPTD') ELSE IF(CPTD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','CPTD') END IF RETURN END *------------------------------------------------- F20C09 = CPTDD REAL FUNCTION CPTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) CPTDD=F20C09(TI) IF(CPTDD.EQ.-1.0E+10) THEN CALL S97C09('CPTDD') ELSE IF(CPTDD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','CPTDD') END IF RETURN END *------------------------------------------------- F21C09 = 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=F21C09(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 FLUORINE', - ' 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 *------------------------------------------------- F22C09 = EPSPT REAL FUNCTION EPSPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C09('EPSPT') EPSPT=-1.0E+30 RETURN END *----------------------------------------------------- F89 = FC *************************************************** * FUNCTION FOR FUNDAMENTAL CONSTANTS * PROPATH VER.10.1, JUNE 7, 1996. * USAGE: B=FC(A) * A, B : CHARACTER TYPE VARIABLES * B='37.99681' WHEN A='M' * B='218.8205' WHEN A='R' *************************************************** REAL FUNCTION FC(A) CHARACTER A*1, MSG*120 COMMON /UNIT/KPA,MESS IF (A.EQ.'M') THEN FC=37.99681 ELSE IF (A.EQ.'R') THEN FC=218.8205 ELSE FC=-1.0E+20 IF (MESS.NE.0) THEN MSG='**** OUT OF RANGE AT FC FOR FLUORINE WHEN A=''' - //A//''' ****' WRITE(6,'(1H, A)') MSG END IF END IF RETURN END *------------------------------------------------- F23C09 = HPD REAL FUNCTION HPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) HPD=F23C09(PI) IF(HPD.EQ.-1.0E+10) THEN CALL S97C09('HPD') ELSE IF(HPD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','HPD') END IF RETURN END *------------------------------------------------- F24C09 = HPDD REAL FUNCTION HPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) HPDD=F24C09(PI) IF(HPDD.EQ.-1.0E+10) THEN CALL S97C09('HPDD') ELSE IF(HPDD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','HPDD') END IF RETURN END *------------------------------------------------- F25C09 = HPT REAL FUNCTION HPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) HPT=F25C09(PI,TI) IF(HPT.EQ.-1.0E+10) THEN CALL S97C09('HPT') ELSE IF(HPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','HPT') END IF RETURN END *------------------------------------------------- F26C09 = HPX REAL FUNCTION HPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) HPX=F26C09(PI,X) IF(HPX.EQ.-1.0E+10) THEN CALL S97C09('HPX') ELSE IF(HPX.EQ.-1.0E+20) THEN CALL S98C09(2,P,X,'P','X','HPX') END IF RETURN END *------------------------------------------------- F27C09 = HTD REAL FUNCTION HTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) HTD=F27C09(TI) IF(HTD.EQ.-1.0E+10) THEN CALL S97C09('HTD') ELSE IF(HTD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','HTD') END IF RETURN END *------------------------------------------------- F28C09 = HTDD REAL FUNCTION HTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) HTDD=F28C09(TI) IF(HTDD.EQ.-1.0E+10) THEN CALL S97C09('HTDD') ELSE IF(HTDD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','HTDD') END IF RETURN END *------------------------------------------------- F29C09 = HTX REAL FUNCTION HTX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) HTX=F29C09(TI,X) IF(HTX.EQ.-1.0E+10) THEN CALL S97C09('HTX') ELSE IF(HTX.EQ.-1.0E+20) THEN CALL S98C09(2,T,X,'T','X','HTX') END IF RETURN END *------------------------------------------------- F30C09 = PST REAL FUNCTION PST(T) REAL TI,T INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) PST=F30C09(TI) IF(PST.EQ.-1.0E+10) THEN CALL S97C09('PST') RETURN ELSE IF(PST.EQ.-1.0E+20) THEN CALL S98C09(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 *------------------------------------------------- F31C09 = SIGP REAL FUNCTION SIGP(P) REAL P,PI PI=P CALL S99C09('SIGP') SIGP=-1.0E+30 RETURN END *------------------------------------------------- F32C09 = SIGT REAL FUNCTION SIGT(T) REAL T,TI TI=T CALL S99C09('SIGT') SIGT=-1.0E+30 RETURN END *------------------------------------------------- F33C09 = SPD REAL FUNCTION SPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) SPD=F33C09(PI) IF(SPD.EQ.-1.0E+10) THEN CALL S97C09('SPD') ELSE IF(SPD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','SPD') END IF RETURN END *------------------------------------------------- F34C09 = SPDD REAL FUNCTION SPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) SPDD=F34C09(PI) IF(SPDD.EQ.-1.0E+10) THEN CALL S97C09('SPDD') ELSE IF(SPDD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','SPDD') END IF RETURN END *------------------------------------------------- F35C09 = SPT REAL FUNCTION SPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) SPT=F35C09(PI,TI) IF(SPT.EQ.-1.0E+10) THEN CALL S97C09('SPT') ELSE IF(SPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','SPT') END IF RETURN END *------------------------------------------------- F36C09 = SPX REAL FUNCTION SPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) SPX=F36C09(PI,X) IF(SPX.EQ.-1.0E+10) THEN CALL S97C09('SPX') ELSE IF(SPX.EQ.-1.0E+20) THEN CALL S98C09(2,P,X,'P','X','SPX') END IF RETURN END *------------------------------------------------- F37C09 = STD REAL FUNCTION STD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) STD=F37C09(TI) IF(STD.EQ.-1.0E+10) THEN CALL S97C09('STD') ELSE IF(STD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','STD') END IF RETURN END *------------------------------------------------- F38C09 = STDD REAL FUNCTION STDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) STDD=F38C09(TI) IF(STDD.EQ.-1.0E+10) THEN CALL S97C09('STDD') ELSE IF(STDD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','STDD') END IF RETURN END *------------------------------------------------- F39C09 = STX REAL FUNCTION STX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) STX=F39C09(TI,X) IF(STX.EQ.-1.0E+10) THEN CALL S97C09('STX') ELSE IF(STX.EQ.-1.0E+20) THEN CALL S98C09(2,T,X,'T','X','STX') END IF RETURN END *------------------------------------------------- F40C09 = TSP REAL FUNCTION TSP(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TSP=F40C09(PI) IF(TSP.EQ.-1.0E+10) THEN CALL S97C09('TSP') RETURN ELSE IF(TSP.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','TSP') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TSP=TSP+273.15 RETURN END *------------------------------------------------- F41C09 = 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=F41C09(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 FLUORINE', - ' 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 *------------------------------------------------- F42C09 = UPD REAL FUNCTION UPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) UPD=F42C09(PI) IF(UPD.EQ.-1.0E+10) THEN CALL S97C09('UPD') ELSE IF(UPD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','UPD') END IF RETURN END *------------------------------------------------- F43C09 = UPDD REAL FUNCTION UPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) UPDD=F43C09(PI) IF(UPDD.EQ.-1.0E+10) THEN CALL S97C09('UPDD') ELSE IF(UPDD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','UPDD') END IF RETURN END *------------------------------------------------- F44C09 = UPT REAL FUNCTION UPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) UPT=F44C09(PI,TI) IF(UPT.EQ.-1.0E+10) THEN CALL S97C09('UPT') ELSE IF(UPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','UPT') END IF RETURN END *------------------------------------------------- F45C09 = UPX REAL FUNCTION UPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) UPX=F45C09(PI,X) IF(UPX.EQ.-1.0E+10) THEN CALL S97C09('UPX') ELSE IF(UPX.EQ.-1.0E+20) THEN CALL S98C09(2,P,X,'P','X','UPX') END IF RETURN END *------------------------------------------------- F46C09 = UTD REAL FUNCTION UTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) UTD=F46C09(TI) IF(UTD.EQ.-1.0E+10) THEN CALL S97C09('UTD') ELSE IF(UTD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','UTD') END IF RETURN END *------------------------------------------------- F47C09 = UTDD REAL FUNCTION UTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) UTDD=F47C09(TI) IF(UTDD.EQ.-1.0E+10) THEN CALL S97C09('UTDD') ELSE IF(UTDD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','UTDD') END IF RETURN END *------------------------------------------------- F48C09 = UTX REAL FUNCTION UTX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) UTX=F48C09(TI,X) IF(UTX.EQ.-1.0E+10) THEN CALL S97C09('UTX') ELSE IF(UTX.EQ.-1.0E+20) THEN CALL S98C09(2,T,X,'T','X','UTX') END IF RETURN END *------------------------------------------------- F49C09 = VPD REAL FUNCTION VPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) VPD=F49C09(PI) IF(VPD.EQ.-1.0E+10) THEN CALL S97C09('VPD') ELSE IF(VPD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','VPD') END IF RETURN END *------------------------------------------------- F50C09 = VPDD REAL FUNCTION VPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) VPDD=F50C09(PI) IF(VPDD.EQ.-1.0E+10) THEN CALL S97C09('VPDD') ELSE IF(VPDD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','VPDD') END IF RETURN END *------------------------------------------------- F51C09 = VPT REAL FUNCTION VPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) VPT=F51C09(PI,TI) IF(VPT.EQ.-1.0E+10) THEN CALL S97C09('VPT') ELSE IF(VPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','VPT') END IF RETURN END *------------------------------------------------- F52C09 = VPX REAL FUNCTION VPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) VPX=F52C09(PI,X) IF(VPX.EQ.-1.0E+10) THEN CALL S97C09('VPX') ELSE IF(VPX.EQ.-1.0E+20) THEN CALL S98C09(2,P,X,'P','X','VPX') END IF RETURN END *------------------------------------------------- F53C09 = VTD REAL FUNCTION VTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) VTD=F53C09(TI) IF(VTD.EQ.-1.0E+10) THEN CALL S97C09('VTD') ELSE IF(VTD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','VTD') END IF RETURN END *------------------------------------------------- F54C09 = VTDD REAL FUNCTION VTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) VTDD=F54C09(TI) IF(VTDD.EQ.-1.0E+10) THEN CALL S97C09('VTDD') ELSE IF(VTDD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','VTDD') END IF RETURN END *------------------------------------------------- F55C09 = VTX REAL FUNCTION VTX(T,X) REAL T,TI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) VTX=F55C09(TI,X) IF(VTX.EQ.-1.0E+10) THEN CALL S97C09('VTX') ELSE IF(VTX.EQ.-1.0E+20) THEN CALL S98C09(2,T,X,'T','X','VTX') END IF RETURN END *------------------------------------------------- F56C09 = XPH REAL FUNCTION XPH(P,H) REAL P,PI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) XPH=F56C09(PI,H) IF(XPH.EQ.-1.0E+10) THEN CALL S97C09('XPH') ELSE IF(XPH.EQ.-1.0E+20) THEN CALL S98C09(2,P,H,'P','H','XPH') END IF RETURN END *------------------------------------------------- F57C09 = XPS REAL FUNCTION XPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) XPS=F57C09(PI,S) IF(XPS.EQ.-1.0E+10) THEN CALL S97C09('XPS') ELSE IF(XPS.EQ.-1.0E+20) THEN CALL S98C09(2,P,S,'P','S','XPS') END IF RETURN END *------------------------------------------------- F58C09 = XPU REAL FUNCTION XPU(P,U) REAL P,PI,U INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) XPU=F58C09(PI,U) IF(XPU.EQ.-1.0E+10) THEN CALL S97C09('XPU') ELSE IF(XPU.EQ.-1.0E+20) THEN CALL S98C09(2,P,U,'P','U','XPU') END IF RETURN END *------------------------------------------------- F59C09 = XPV REAL FUNCTION XPV(P,V) REAL P,PI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) XPV=F59C09(PI,V) IF(XPV.EQ.-1.0E+10) THEN CALL S97C09('XPV') ELSE IF(XPV.EQ.-1.0E+20) THEN CALL S98C09(2,P,V,'P','V','XPV') END IF RETURN END *------------------------------------------------- F60C09 = XTH REAL FUNCTION XTH(T,H) REAL T,TI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) XTH=F60C09(TI,H) IF(XTH.EQ.-1.0E+10) THEN CALL S97C09('XTH') ELSE IF(XTH.EQ.-1.0E+20) THEN CALL S98C09(2,T,H,'T','H','XTH') END IF RETURN END *------------------------------------------------- F61C09 = XTS REAL FUNCTION XTS(T,S) REAL T,TI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) XTS=F61C09(TI,S) IF(XTS.EQ.-1.0E+10) THEN CALL S97C09('XTS') ELSE IF(XTS.EQ.-1.0E+20) THEN CALL S98C09(2,T,S,'T','S','XTS') END IF RETURN END *------------------------------------------------- F62C09 = XTU REAL FUNCTION XTU(T,U) REAL T,TI,U INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) XTU=F62C09(TI,U) IF(XTU.EQ.-1.0E+10) THEN CALL S97C09('XTU') ELSE IF(XTU.EQ.-1.0E+20) THEN CALL S98C09(2,T,U,'T','U','XTU') END IF RETURN END *------------------------------------------------- F63C09 = XTV REAL FUNCTION XTV(T,V) REAL T,TI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) XTV=F63C09(TI,V) IF(XTV.EQ.-1.0E+10) THEN CALL S97C09('XTV') ELSE IF(XTV.EQ.-1.0E+20) THEN CALL S98C09(2,T,V,'T','V','XTV') END IF RETURN END *------------------------------------------------- F64C09 = TPH REAL FUNCTION TPH(P,H) REAL P,PI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TPH=F64C09(PI,H) IF(TPH.EQ.-1.0E+10) THEN CALL S97C09('TPH') RETURN ELSE IF(TPH.EQ.-1.0E+20) THEN CALL S98C09(2,P,H,'P','H','TPH') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPH=TPH+273.15 RETURN END *------------------------------------------------- F65C09 = TPS REAL FUNCTION TPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TPS=F65C09(PI,S) IF(TPS.EQ.-1.0E+10) THEN CALL S97C09('TPS') RETURN ELSE IF(TPS.EQ.-1.0E+20) THEN CALL S98C09(2,P,S,'P','S','TPS') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPS=TPS+273.15 RETURN END *------------------------------------------------- F66C09 = PLDT REAL FUNCTION PLDT(T) REAL T,TI TI=T CALL S99C09('PLDT') PLDT=-1.0E+30 RETURN END *------------------------------------------------- F67C09 = TLDP REAL FUNCTION TLDP(P) REAL P,PI PI=P CALL S99C09('TLDP') TLDP=-1.0E+30 RETURN END *------------------------------------------------- F68C09 = PMLT REAL FUNCTION PMLT(T) REAL TI,T INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) PMLT=F68C09(TI) IF(PMLT.EQ.-1.0E+10) THEN CALL S97C09('PMLT') RETURN ELSE IF(PMLT.EQ.-1.0E+20) THEN CALL S98C09(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 *------------------------------------------------- F69C09 = TMLP REAL FUNCTION TMLP(P) REAL PI,P INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TMLP=F69C09(PI) IF(TMLP.EQ.-1.0E+10) THEN CALL S97C09('TMLP') RETURN ELSE IF(TMLP.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','TMLP') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TMLP=TMLP+273.15 RETURN END *------------------------------------------------- F70C09 = TPV REAL FUNCTION TPV(P,V) REAL P,PI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TPV=F70C09(PI,V) IF(TPV.EQ.-1.0E+10) THEN CALL S97C09('TPV') RETURN ELSE IF(TPV.EQ.-1.0E+20) THEN CALL S98C09(2,P,V,'P','V','TPV') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPV=TPV+273.15 RETURN END *------------------------------------------------- F71C09 = HPS REAL FUNCTION HPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) HPS=F71C09(PI,S) IF(HPS.EQ.-1.0E+10) THEN CALL S97C09('HPS') ELSE IF(HPS.EQ.-1.0E+20) THEN CALL S98C09(2,P,S,'P','S','HPS') END IF RETURN END *------------------------------------------------- F72C09 = PSTD REAL FUNCTION PSTD(T) REAL T,TI TI=T CALL S99C09('PSTD') PSTD=-1.0E+30 RETURN END *------------------------------------------------- F73C09 = PSTDD REAL FUNCTION PSTDD(T) REAL T,TI TI=T CALL S99C09('PSTDD') PSTDD=-1.0E+30 RETURN END *------------------------------------------------- F74C09 = TSPD REAL FUNCTION TSPD(P) REAL P,PI PI=P CALL S99C09('TSPD') TSPD=-1.0E+30 RETURN END *------------------------------------------------- F75C09 = TSPDD REAL FUNCTION TSPDD(P) REAL P,PI PI=P CALL S99C09('TSPDD') TSPDD=-1.0E+30 RETURN END *------------------------------------------------- F76C09 = CVPDD REAL FUNCTION CVPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) CVPDD=F76C09(PI) IF(CVPDD.EQ.-1.0E+10) THEN CALL S97C09('CVPDD') ELSE IF(CVPDD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','CVPDD') END IF RETURN END *------------------------------------------------- F77C09 = CVPT REAL FUNCTION CVPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) CVPT=F77C09(PI,TI) IF(CVPT.EQ.-1.0E+10) THEN CALL S97C09('CVPT') ELSE IF(CVPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','CVPT') END IF RETURN END *------------------------------------------------- F78C09 = CVTDD REAL FUNCTION CVTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) CVTDD=F78C09(TI) IF(CVTDD.EQ.-1.0E+10) THEN CALL S97C09('CVTDD') ELSE IF(CVTDD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','CVTDD') END IF RETURN END *------------------------------------------------- F79C09 = UPS REAL FUNCTION UPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) UPS=F79C09(PI,S) IF(UPS.EQ.-1.0E+10) THEN CALL S97C09('UPS') ELSE IF(UPS.EQ.-1.0E+20) THEN CALL S98C09(2,P,S,'P','S','UPS') END IF RETURN END *------------------------------------------------- F80C09 = VPS REAL FUNCTION VPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) VPS=F80C09(PI,S) IF(VPS.EQ.-1.0E+10) THEN CALL S97C09('VPS') ELSE IF(VPS.EQ.-1.0E+20) THEN CALL S98C09(2,P,S,'P','S','VPS') END IF RETURN END *------------------------------------------------- F81C09 = PRPT REAL FUNCTION PRPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99C09('PRPT') PRPT=-1.0E+30 RETURN END *------------------------------------------------- F82C09 = AKPT REAL FUNCTION AKPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) AKPT=F82C09(PI,TI) IF(AKPT.EQ.-1.0E+10) THEN CALL S97C09('AKPT') ELSE IF(AKPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','AKPT') END IF RETURN END *------------------------------------------------- F83C09 = WPT REAL FUNCTION WPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) WPT=F83C09(PI,TI) IF(WPT.EQ.-1.0E+10) THEN CALL S97C09('WPT') ELSE IF(WPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','WPT') END IF RETURN END *------------------------------------------------- F84C09 = 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 VARIABLES C B='FLUORINE' WHEN A='S' C B='F2' 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='FLUORINE' ELSE IF (A.EQ.'C') THEN IDENTF='F2' ELSE IF (A.EQ.'V') THEN IDENTF='12.1' ELSE IDENTF='????????????????????' IF (MESS.NE.0) THEN MSG=' *** OUT OF RANGE AT IDENTF FOR FLUORINE WHEN A=''' & //A//''' ****' WRITE(6,'(1H,A)') MSG END IF END IF RETURN END *------------------------------------------------- F85C09 = PRPD REAL FUNCTION PRPD(P) REAL P,PI PI=P CALL S99C09('PRPD') PRPD=-1.0E+30 RETURN END *------------------------------------------------- F86C09 = PRPDD REAL FUNCTION PRPDD(P) REAL P,PI PI=P CALL S99C09('PRPDD') PRPDD=-1.0E+30 RETURN END *------------------------------------------------- F87C09 = PRTD REAL FUNCTION PRTD(T) REAL T,TI TI=T CALL S99C09('PRTD') PRTD=-1.0E+30 RETURN END *------------------------------------------------- F88C09 = PRTDD REAL FUNCTION PRTDD(T) REAL T,TI TI=T CALL S99C09('PRTDD') PRTDD=-1.0E+30 RETURN END *------------------------------------------------- F90C09 = BSPT REAL FUNCTION BSPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) BSPT=F90C09(PI,TI) IF(BSPT.EQ.-1.0E+10) THEN CALL S97C09('BSPT') RETURN ELSE IF(BSPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','BSPT') RETURN END IF RETURN END *------------------------------------------------- F91C09 = BTPT REAL FUNCTION BTPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) BTPT=F91C09(PI,TI) IF(BTPT.EQ.-1.0E+10) THEN CALL S97C09('BTPT') RETURN ELSE IF(BTPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','BTPT') RETURN END IF RETURN END *------------------------------------------------- F92C09 = BPPT REAL FUNCTION BPPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) BPPT=F92C09(PI,TI) IF(BPPT.EQ.-1.0E+10) THEN CALL S97C09('BPPT') ELSE IF(BPPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','BPPT') END IF RETURN END *------------------------------------------------- F93C09 = BVPT REAL FUNCTION BVPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) BVPT=F93C09(PI,TI) IF(BVPT.EQ.-1.0E+10) THEN CALL S97C09('BVPT') ELSE IF(BVPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','BVPT') END IF RETURN END *------------------------------------------------- F94C09 = AJTPT REAL FUNCTION AJTPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) AJTPT=F94C09(PI,TI) IF(AJTPT.EQ.-1.0E+10) THEN CALL S97C09('AJTPT') RETURN ELSE IF(AJTPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','AJTPT') RETURN END IF RETURN END *------------------------------------------------- F95C09 = GAMPT REAL FUNCTION GAMPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) TI=G99C09(KPA,T) GAMPT=F95C09(PI,TI) IF(GAMPT.EQ.-1.0E+10) THEN CALL S97C09('GAMPT') ELSE IF(GAMPT.EQ.-1.0E+20) THEN CALL S98C09(2,P,T,'P','T','GAMPT') END IF RETURN END *------------------------------------------------- F96C09 = GAMPDD REAL FUNCTION GAMPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98C09(KPA,P) GAMPDD=F96C09(PI) IF(GAMPDD.EQ.-1.0E+10) THEN CALL S97C09('GAMPDD') ELSE IF(GAMPDD.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','GAMPDD') END IF RETURN END *------------------------------------------------- F97C09 = GAMTDD REAL FUNCTION GAMTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99C09(KPA,T) GAMTDD=F97C09(TI) IF(GAMTDD.EQ.-1.0E+10) THEN CALL S97C09('GAMTDD') ELSE IF(GAMTDD.EQ.-1.0E+20) THEN CALL S98C09(1,T,T,'T','T','GAMTDD') END IF RETURN END *-------------------------------------------------- F98C09 = TPSEUP REAL FUNCTION TPSEUP(P) ***** TPSEUP = PSEUDO BOILING POINT REAL P,PI INTEGER KPA,MESS COMMON /UNIT/KPA,MESS PI=G98C09(KPA,P) TPSEUP=F98C09(PI) IF(TPSEUP.EQ.-1.0E+10) THEN CALL S97C09('TPSEUP') ELSE IF(TPSEUP.EQ.-1.0E+20) THEN CALL S98C09(1,P,P,'P','P','TPSEUP') END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPSEUP=TPSEUP+273.15 RETURN END *------------------------------------------------- F99C09 = PSBT REAL FUNCTION PSBT(T) REAL T,TI TI=T CALL S99C09('PSBT') PSBT=-1.0E+30 RETURN END *------------------------------------------------- F100C09 = TSBP REAL FUNCTION TSBP(P) REAL P,PI PI=P CALL S99C09('TSBP') TSBP=-1.0E+30 RETURN END *------------------------------------------------- G98C09 REAL FUNCTION G98C09(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 G98C09=P*PBAR RETURN END *------------------------------------------------- G99C09 REAL FUNCTION G99C09(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 G99C09=T-T0K RETURN END ***** PROPATH V10:FLUORINE FUNCTIONS ***1996/06/10:T.SHIGECHI*********** ***** (SUBSTANCE NUMBER = C09 IN NAGASAKI UNIV.) REAL FUNCTION F4C09(FP) ***** FP=INPUT PRESSURE IN BAR ***** IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F4C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F4C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAUS,HL) CALL S10C09(OMEGAG,TAUS,HG) ALHP=(HG-HL)/AKGMOL F4C09=REAL(ALHP) RETURN END REAL FUNCTION F5C09(FT) ***** FT=INPUT TEMPERATURE IN C(DEGREE CELSIUS) ***** IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F5C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F5C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAU,HL) CALL S10C09(OMEGAG,TAU,HG) ALHT=(HG-HL)/AKGMOL F5C09=REAL(ALHT) RETURN END REAL FUNCTION F16C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F16C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F16C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S12C09(OMEGAL,TAUS,CPL) CPPD=CPL/AKGMOL F16C09=REAL(CPPD) RETURN END REAL FUNCTION F17C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F17C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F17C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S12C09(OMEGAG,TAUS,CPG) CPPDD=CPG/AKGMOL F17C09=REAL(CPPDD) RETURN END REAL FUNCTION F18C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,AKGMOL=3.799681D-02) F18C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F18C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S12C09(OMEGA,TAU,CP) CPPT=CP/AKGMOL F18C09=REAL(CPPT) RETURN END REAL FUNCTION F19C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F19C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F19C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S12C09(OMEGAL,TAU,CPL) CPTD=CPL/AKGMOL F19C09=REAL(CPTD) RETURN END REAL FUNCTION F20C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F20C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F20C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S12C09(OMEGAG,TAU,CPG) CPTDD=CPG/AKGMOL F20C09=REAL(CPTDD) RETURN END REAL FUNCTION F21C09(A) CHARACTER*1 A,B(1:5) DOUBLE PRECISION CRP(1:5),HCR,SCR PARAMETER(AKGMOL=3.799681D-02) DATA B(1)/'H'/,B(2)/'P'/,B(3)/'S'/,B(4)/'T'/,B(5)/'V'/ CALL S10C09(1.0D0,1.0D0,HCR) CALL S8C09(1.0D0,1.0D0,SCR) CRP(1)=HCR/AKGMOL CRP(2)=5.23952D+01 CRP(3)=SCR/AKGMOL CRP(4)=144.414D+00 CRP(5)=(1.0D0/1.56030D+04)/AKGMOL DO 10 I=1,5 F21C09=REAL(CRP(I)) IF(A.EQ.B(I)) RETURN 10 CONTINUE F21C09=-1.0E+20 RETURN END REAL FUNCTION F23C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F23C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F23C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAUS,HL) HPD=HL/AKGMOL F23C09=REAL(HPD) RETURN END REAL FUNCTION F24C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F24C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F24C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S10C09(OMEGAG,TAUS,HG) HPDD=HG/AKGMOL F24C09=REAL(HPDD) RETURN END REAL FUNCTION F25C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,AKGMOL=3.799681D-02) F25C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F25C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S10C09(OMEGA,TAU,H) HPT=H/AKGMOL F25C09=REAL(HPT) RETURN END REAL FUNCTION F26C09(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F26C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F26C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAUS,HL) CALL S10C09(OMEGAG,TAUS,HG) HPX=(HL+X*(HG-HL))/AKGMOL F26C09=REAL(HPX) RETURN END REAL FUNCTION F27C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F27C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F27C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAU,HL) HTD=HL/AKGMOL F27C09=REAL(HTD) RETURN END REAL FUNCTION F28C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F28C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F28C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S10C09(OMEGAG,TAU,HG) HTDD=HG/AKGMOL F28C09=REAL(HTDD) RETURN END REAL FUNCTION F29C09(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F29C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F29C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) X=DBLE(FX) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAU,HL) CALL S10C09(OMEGAG,TAU,HG) HTX=(HL+X*(HG-HL))/AKGMOL F29C09=REAL(HTX) RETURN END REAL FUNCTION F30C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,PCR=5.23952D+06) F30C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F30C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN PST=PIS*PCR F30C09=REAL(PST*1.0D-05) RETURN END REAL FUNCTION F33C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F33C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F33C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAUS,SL) SPD=SL/AKGMOL F33C09=REAL(SPD) RETURN END REAL FUNCTION F34C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F34C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F34C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S8C09(OMEGAG,TAUS,SG) SPDD=SG/AKGMOL F34C09=REAL(SPDD) RETURN END REAL FUNCTION F35C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,AKGMOL=3.799681D-02) F35C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F35C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S8C09(OMEGA,TAU,S) SPT=S/AKGMOL F35C09=REAL(SPT) RETURN END REAL FUNCTION F36C09(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F36C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F36C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAUS,SL) CALL S8C09(OMEGAG,TAUS,SG) SPX=(SL+X*(SG-SL))/AKGMOL F36C09=REAL(SPX) RETURN END REAL FUNCTION F37C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F37C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F37C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAU,SL) STD=SL/AKGMOL F37C09=REAL(STD) RETURN END REAL FUNCTION F38C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F38C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F38C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S8C09(OMEGAG,TAU,SG) STDD=SG/AKGMOL F38C09=REAL(STDD) RETURN END REAL FUNCTION F39C09(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F39C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F39C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) X=DBLE(FX) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAU,SL) CALL S8C09(OMEGAG,TAU,SG) STX=(SL+X*(SG-SL))/AKGMOL F39C09=REAL(STX) RETURN END REAL FUNCTION F40C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F40C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F40C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN TSP=TCR/TAUS F40C09=REAL(TSP-273.15D0) RETURN END REAL FUNCTION F41C09(A) ***** TRIPLE POINT FROM TABLE 6 & 7 ON PAGES 188 AND 189. ***** CHARACTER*1 A,B(1:2) DOUBLE PRECISION TRPL(1:2) DATA B(1)/'P'/,B(2)/'T'/ TRPL(1)=2.39D-03 TRPL(2)=-219.669D+00 DO 10 I=1,2 F41C09=REAL(TRPL(I)) IF(A.EQ.B(I)) RETURN 10 CONTINUE F41C09=-1.0E+20 RETURN END REAL FUNCTION F42C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F42C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F42C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S9C09(OMEGAL,TAUS,UL) UPD=UL/AKGMOL F42C09=REAL(UPD) RETURN END REAL FUNCTION F43C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F43C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F43C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S9C09(OMEGAG,TAUS,UG) UPDD=UG/AKGMOL F43C09=REAL(UPDD) RETURN END REAL FUNCTION F44C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,AKGMOL=3.799681D-02) F44C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F44C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S9C09(OMEGA,TAU,U) UPT=U/AKGMOL F44C09=REAL(UPT) RETURN END REAL FUNCTION F45C09(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F45C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F45C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S9C09(OMEGAL,TAUS,UL) CALL S9C09(OMEGAG,TAUS,UG) UPX=(UL+X*(UG-UL))/AKGMOL F45C09=REAL(UPX) RETURN END REAL FUNCTION F46C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F46C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F46C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S9C09(OMEGAL,TAU,UL) UTD=UL/AKGMOL F46C09=REAL(UTD) RETURN END REAL FUNCTION F47C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F47C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F47C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S9C09(OMEGAG,TAU,UG) UTDD=UG/AKGMOL F47C09=REAL(UTDD) RETURN END REAL FUNCTION F48C09(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F48C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F48C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) X=DBLE(FX) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S9C09(OMEGAL,TAU,UL) CALL S9C09(OMEGAG,TAU,UG) UTX=(UL+X*(UG-UL))/AKGMOL F48C09=REAL(UTX) RETURN END REAL FUNCTION F49C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F49C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F49C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN VPD=(1.0D0/(OMEGAL*RHOCR))/AKGMOL F49C09=REAL(VPD) RETURN END REAL FUNCTION F50C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F50C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F50C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAG.LE.-1.0D+10) RETURN VPDD=(1.0D0/(OMEGAG*RHOCR))/AKGMOL F50C09=REAL(VPDD) RETURN END REAL FUNCTION F51C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,RHOCR=1.56030D+04, - AKGMOL=3.799681D-02) F51C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F51C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN VPT=(1.0D0/(OMEGA*RHOCR))/AKGMOL F51C09=REAL(VPT) RETURN END REAL FUNCTION F52C09(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F52C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F52C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR X=DBLE(FX) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN VL=1.0D0/(OMEGAL*RHOCR) VG=1.0D0/(OMEGAG*RHOCR) VPX=(VL+X*(VG-VL))/AKGMOL F52C09=REAL(VPX) RETURN END REAL FUNCTION F53C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F53C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F53C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN VTD=(1.0D0/(OMEGAL*RHOCR))/AKGMOL F53C09=REAL(VTD) RETURN END REAL FUNCTION F54C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F54C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F54C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAG.LE.-1.0D+10) RETURN VTDD=(1.0D0/(OMEGAG*RHOCR))/AKGMOL F54C09=REAL(VTDD) RETURN END REAL FUNCTION F55C09(FT,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F55C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F55C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) X=DBLE(FX) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN VL=1.0D0/(OMEGAL*RHOCR) VG=1.0D0/(OMEGAG*RHOCR) VTX=(VL+X*(VG-VL))/AKGMOL F55C09=REAL(VTX) RETURN END REAL FUNCTION F56C09(FP,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F56C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F56C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR H=DBLE(FH) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAUS,HL) CALL S10C09(OMEGAG,TAUS,HG) HPD=HL/AKGMOL HPDD=HG/AKGMOL 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 F56C09=REAL(XPH) RETURN END REAL FUNCTION F57C09(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F57C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F57C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAUS,SL) CALL S8C09(OMEGAG,TAUS,SG) SPD=SL/AKGMOL SPDD=SG/AKGMOL 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 F57C09=REAL(XPS) RETURN END REAL FUNCTION F58C09(FP,FU) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F58C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F58C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR U=DBLE(FU) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S9C09(OMEGAL,TAUS,UL) CALL S9C09(OMEGAG,TAUS,UG) UPD=UL/AKGMOL UPDD=UG/AKGMOL 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 F58C09=REAL(XPU) RETURN END REAL FUNCTION F59C09(FP,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F59C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F59C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR V=DBLE(FV) CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN VPD=(1.0D0/(OMEGAL*RHOCR))/AKGMOL VPDD=(1.0D0/(OMEGAG*RHOCR))/AKGMOL 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 F59C09=REAL(XPV) RETURN END REAL FUNCTION F60C09(FT,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F60C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F60C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) H=DBLE(FH) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAU,HL) CALL S10C09(OMEGAG,TAU,HG) HTD=HL/AKGMOL HTDD=HG/AKGMOL 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 F60C09=REAL(XTH) RETURN END REAL FUNCTION F61C09(FT,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F61C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F61C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) S=DBLE(FS) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAU,SL) CALL S8C09(OMEGAG,TAU,SG) STD=SL/AKGMOL STDD=SG/AKGMOL 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 F61C09=REAL(XTS) RETURN END REAL FUNCTION F62C09(FT,FU) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F62C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F62C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) U=DBLE(FU) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S9C09(OMEGAL,TAU,UL) CALL S9C09(OMEGAG,TAU,UG) UTD=UL/AKGMOL UTDD=UG/AKGMOL 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 F62C09=REAL(XTU) RETURN END REAL FUNCTION F63C09(FT,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F63C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F63C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) V=DBLE(FV) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAL.LE.-1.0D+10) RETURN VTD=(1.0D0/(OMEGAL*RHOCR))/AKGMOL VTDD=(1.0D0/(OMEGAG*RHOCR))/AKGMOL 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 F63C09=REAL(XTV) RETURN END REAL FUNCTION F64C09(FP,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,AKGMOL=3.799681D-02) F64C09=-1.0E+20 IF ((FP.GE.2.3899E-03).AND.(FP.LE.200.01E+00))THEN FTMIN=F69C09(FP) FHMIN=F25C09(FP,FTMIN) FHMAX=F25C09(FP,26.851) ELSE RETURN END IF 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 F64C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR H=DBLE(FH)*AKGMOL CALL S10C09(1.0D0,1.0D0,HCR) IF ((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(H/HCR-1.0D0).LT.1.0D-05)) - THEN F64C09=REAL(TCR-273.15D0) RETURN END IF IF (PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S5C09(PI,TAU0,OMEGA0) IF (OMEGA0.LE.-1.0D+10) RETURN CALL S10C09(OMEGA0,TAU0,H0) CALL S10C09(OMEGA0,TAU0,CP0) ELSE CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S10C09(OMEGAL,TAUS,HL) CALL S10C09(OMEGAG,TAUS,HG) IF ((H.GE.HL).AND.(H.LE.HG)) THEN F64C09=REAL(TCR/TAUS-273.15D0) RETURN END IF IF (H.GT.HG) THEN CALL S12C09(OMEGAG,TAUS,CPG) TAU0=TAUS H0=HG CP0=CPG ELSE IF (H.LT.HL) THEN CALL S12C09(OMEGAL,TAUS,CPL) TAU0=TAUS H0=HL CP0=CPL END IF END IF TAU1=1.0D0/(1.0D0/TAU0+(H-H0)/(CP0*TCR)) CALL S5C09(PI,TAU1,OMEGA1) IF (OMEGA1.LE.-1.0D+10) RETURN CALL S10C09(OMEGA1,TAU1,H1) TH0=1.0D0/TAU0 TH1=1.0D0/TAU1 CALL S40C09(1,PI,H,TH0,TH1,H0,H1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=1.0D0/THW TPH=TCR/TAUW F64C09=REAL(TPH-273.15D0) RETURN END REAL FUNCTION F65C09(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,AKGMOL=3.799681D-02) F65C09=-1.0E+20 IF ((FP.GE.2.3899E-03).AND.(FP.LE.200.01E+00))THEN FTMIN=F69C09(FP) FTMAX=26.85 FSMIN=F35C09(FP,FTMIN) FSMAX=F35C09(FP,FTMAX) ELSE RETURN END IF 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 F65C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS)*AKGMOL CALL S8C09(1.0D0,1.0D0,SCR) IF ((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(S/SCR-1.0D0).LT.1.0D-05)) - THEN F65C09=REAL(TCR-273.15D0) RETURN END IF IF (PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S5C09(PI,TAU0,OMEGA0) IF (OMEGA0.LE.-1.0D+10) RETURN CALL S8C09(OMEGA0,TAU0,S0) CALL S12C09(OMEGA0,TAU0,CP0) ELSE CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAUS,SL) CALL S8C09(OMEGAG,TAUS,SG) IF ((S.GE.SL).AND.(S.LE.SG)) THEN F65C09=REAL(TCR/TAUS-273.15D0) RETURN END IF IF (S.GT.SG) THEN CALL S12C09(OMEGAG,TAUS,CPG) TAU0=TAUS S0=SG CP0=CPG ELSE IF (S.LT.SL) THEN CALL S12C09(OMEGAL,TAUS,CPL) TAU0=TAUS S0=SL CP0=CPL END IF END IF TAU1=TAU0/EXP((S-S0)/CP0) CALL S5C09(PI,TAU1,OMEGA1) IF (OMEGA1.LE.-1.0D+10) RETURN CALL S8C09(OMEGA1,TAU1,S1) TH0=1.0D0/TAU0 TH1=1.0D0/TAU1 CALL S40C09(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=1.0D0/THW TPS=TCR/TAUW F65C09=REAL(TPS-273.15D0) RETURN END REAL FUNCTION F68C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F68C09=-1.0E+20 IF ((FT.LT.-219.6691).OR.(FT.GT.-217.750)) RETURN T=DBLE(FT)+273.15D0 CALL S1C09(T,PMELT) PMLT=PMELT*1.0D-05 F68C09=REAL(PMLT) RETURN END REAL FUNCTION F69C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F69C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.200.01E+00)) RETURN F69C09=-1.0E+10 P=DBLE(FP)*1.0D+05 CALL S2C09(P,TMELT) IF (TMELT.LE.-1.0D+10) RETURN TMLP=TMELT F69C09=REAL(TMLP-273.15D0) RETURN END REAL FUNCTION F70C09(FP,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,RHOCR=1.56030D+04,TCR=144.414D+00, - AKGMOL=3.799681D-02) DATA FMIN/0.9998/,FMAX/1.002/ F70C09=-1.0E+20 IF ((FP.GE.2.3899E-03).AND.(FP.LE.200.01E+00))THEN FTMIN=F69C09(FP) FVMIN=F51C09(FP,FTMIN) FVMAX=F51C09(FP,26.851) ELSE RETURN END IF IF ((FVMIN.LE.-1.0E+10).OR.(FVMAX.LE.-1.0E+10)) RETURN IF ((FV.LT.FVMIN*FMIN).OR.(FV.GT.FVMAX*FMAX)) RETURN F70C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR OMEGA=1.0D0/(DBLE(FV)*AKGMOL*RHOCR) IF (PI.GE.1.0D0) THEN CALL S6C09(PI,OMEGA,TAU) IF (TAU.LE.-1.0D+10) RETURN TAUW=TAU ELSE CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAL.LE.-1.0D+10) RETURN IF ((OMEGA.GE.OMEGAG).AND.(OMEGA.LE.OMEGAL)) THEN TAUW=TAUS ELSE CALL S6C09(PI,OMEGA,TAU) IF (TAU.LE.-1.0D+10) RETURN TAUW=TAU END IF END IF TPV=TCR/TAUW F70C09=REAL(TPV-273.15D0) RETURN END REAL FUNCTION F71C09(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F71C09=-1.0E+20 IF ((FP.GE.2.3899E-03).AND.(FP.LE.200.01E+00))THEN FTMIN=F69C09(FP) FTMAX=26.85 FSMIN=F35C09(FP,FTMIN) FSMAX=F35C09(FP,FTMAX) ELSE RETURN END IF 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 F71C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS)*AKGMOL CALL S8C09(1.0D0,1.0D0,SCR) IF ((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(S/SCR-1.0D0).LT.1.0D-05)) - THEN CALL S10C09(1.0D0,1.0D0,HCR) F71C09=REAL(HCR/AKGMOL) RETURN END IF IF(PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S5C09(PI,TAU0,OMEGA0) IF(OMEGA0.LE.-1.0D+10) RETURN CALL S8C09(OMEGA0,TAU0,S0) CALL S12C09(OMEGA0,TAU0,CP0) ELSE CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF(OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAUS,SL) CALL S8C09(OMEGAG,TAUS,SG) IF((S.GE.SL).AND.(S.LE.SG)) THEN XPS=(S-SL)/(SG-SL) CALL S10C09(OMEGAL,TAUS,HL) CALL S10C09(OMEGAG,TAUS,HG) HPS=HL+XPS*(HG-HL) F71C09=REAL(HPS/AKGMOL) RETURN END IF IF(S.GT.SG) THEN CALL S12C09(OMEGAG,TAUS,CPG) TAU0=TAUS S0=SG CP0=CPG ELSE IF(S.LT.SL) THEN CALL S12C09(OMEGAL,TAUS,CPL) TAU0=TAUS S0=SL CP0=CPL END IF END IF TAU1=TAU0/EXP((S-S0)/CP0) CALL S5C09(PI,TAU1,OMEGA1) IF(OMEGA1.LE.-1.0D+10) RETURN CALL S8C09(OMEGA1,TAU1,S1) TH0=1.0D0/TAU0 TH1=1.0D0/TAU1 CALL S40C09(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=1.0D0/THW CALL S5C09(PI,TAUW,OMEGAW) IF (OMEGAW.LE.-1.0D+10) RETURN CALL S10C09(OMEGAW,TAUW,HW) HPS=HW F71C09=REAL(HPS/AKGMOL) RETURN END REAL FUNCTION F76C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F76C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F76C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S11C09(OMEGAG,TAUS,CVG) CVPDD=CVG/AKGMOL F76C09=REAL(CVPDD) RETURN END REAL FUNCTION F77C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00,AKGMOL=3.799681D-02) F77C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F77C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S11C09(OMEGA,TAU,CV) CVPT=CV/AKGMOL F77C09=REAL(CVPT) RETURN END REAL FUNCTION F78C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00,AKGMOL=3.799681D-02) F78C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F78C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S11C09(OMEGAG,TAU,CVG) CVTDD=CVG/AKGMOL F78C09=REAL(CVTDD) RETURN END * REAL FUNCTION F79C09(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,AKGMOL=3.799681D-02) F79C09=-1.0E+20 IF ((FP.GE.2.3899E-03).AND.(FP.LE.200.01E+00))THEN FTMIN=F69C09(FP) FTMAX=26.85 FSMIN=F35C09(FP,FTMIN) FSMAX=F35C09(FP,FTMAX) ELSE RETURN END IF 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 F79C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS)*AKGMOL CALL S8C09(1.0D0,1.0D0,SCR) IF ((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(S/SCR-1.0D0).LT.1.0D-05)) - THEN CALL S9C09(1.0D0,1.0D0,UCR) F79C09=REAL(UCR/AKGMOL) RETURN END IF IF(PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S5C09(PI,TAU0,OMEGA0) IF(OMEGA0.LE.-1.0D+10) RETURN CALL S8C09(OMEGA0,TAU0,S0) CALL S12C09(OMEGA0,TAU0,CP0) ELSE CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF(OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAUS,SL) CALL S8C09(OMEGAG,TAUS,SG) IF((S.GE.SL).AND.(S.LE.SG)) THEN XPS=(S-SL)/(SG-SL) CALL S9C09(OMEGAL,TAUS,UL) CALL S9C09(OMEGAG,TAUS,UG) UPS=UL+XPS*(UG-UL) F79C09=REAL(UPS/AKGMOL) RETURN END IF IF(S.GT.SG) THEN CALL S12C09(OMEGAG,TAUS,CPG) TAU0=TAUS S0=SG CP0=CPG ELSE IF(S.LT.SL) THEN CALL S12C09(OMEGAL,TAUS,CPL) TAU0=TAUS S0=SL CP0=CPL END IF END IF TAU1=TAU0/EXP((S-S0)/CP0) CALL S5C09(PI,TAU1,OMEGA1) IF(OMEGA1.LE.-1.0D+10) RETURN CALL S8C09(OMEGA1,TAU1,S1) TH0=1.0D0/TAU0 TH1=1.0D0/TAU1 CALL S40C09(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=1.0D0/THW CALL S5C09(PI,TAUW,OMEGAW) IF (OMEGAW.LE.-1.0D+10) RETURN CALL S9C09(OMEGAW,TAUW,UW) UPS=UW F79C09=REAL(UPS/AKGMOL) RETURN END REAL FUNCTION F80C09(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,RHOCR=1.56030D+04,AKGMOL=3.799681D-02) F80C09=-1.0E+20 IF ((FP.GE.2.3899E-03).AND.(FP.LE.200.01E+00))THEN FTMIN=F69C09(FP) FTMAX=26.85 FSMIN=F35C09(FP,FTMIN) FSMAX=F35C09(FP,FTMAX) ELSE RETURN END IF 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 F80C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR S=DBLE(FS)*AKGMOL CALL S8C09(1.0D0,1.0D0,SCR) IF ((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(S/SCR-1.0D0).LT.1.0D-05)) - THEN VCR=1.0D0/RHOCR F80C09=REAL(VCR/AKGMOL) RETURN END IF IF(PI.GE.1.0D0) THEN TAU0=1.0D0 CALL S5C09(PI,TAU0,OMEGA0) IF(OMEGA0.LE.-1.0D+10) RETURN CALL S8C09(OMEGA0,TAU0,S0) CALL S12C09(OMEGA0,TAU0,CP0) ELSE CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF(OMEGAL.LE.-1.0D+10) RETURN CALL S8C09(OMEGAL,TAUS,SL) CALL S8C09(OMEGAG,TAUS,SG) IF((S.GE.SL).AND.(S.LE.SG)) THEN XPS=(S-SL)/(SG-SL) VL=1.0D0/(OMEGAL*RHOCR) VG=1.0D0/(OMEGAG*RHOCR) VPS=VL+XPS*(VG-VL) F80C09=REAL(VPS/AKGMOL) RETURN END IF IF(S.GT.SG) THEN CALL S12C09(OMEGAG,TAUS,CPG) TAU0=TAUS S0=SG CP0=CPG ELSE IF(S.LT.SL) THEN CALL S12C09(OMEGAL,TAUS,CPL) TAU0=TAUS S0=SL CP0=CPL END IF END IF TAU1=TAU0/EXP((S-S0)/CP0) CALL S5C09(PI,TAU1,OMEGA1) IF(OMEGA1.LE.-1.0D+10) RETURN CALL S8C09(OMEGA1,TAU1,S1) TH0=1.0D0/TAU0 TH1=1.0D0/TAU1 CALL S40C09(2,PI,S,TH0,TH1,S0,S1,THW) IF (THW.LE.-1.0D+10) RETURN TAUW=1.0D0/THW CALL S5C09(PI,TAUW,OMEGAW) IF (OMEGAW.LE.-1.0D+10) RETURN VPS=1.0D0/(OMEGAW*RHOCR) F80C09=REAL(VPS/AKGMOL) RETURN END REAL FUNCTION F82C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F82C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F82C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S13C09(OMEGA,TAU,AKPT) F82C09=REAL(AKPT) RETURN END REAL FUNCTION F83C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F83C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F83C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S14C09(OMEGA,TAU,WPT) F83C09=REAL(WPT) RETURN END REAL FUNCTION F90C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F90C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F90C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S15C09(OMEGA,TAU,BSPT) F90C09=REAL(BSPT) RETURN END REAL FUNCTION F91C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F91C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F91C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S16C09(OMEGA,TAU,BTPT) F91C09=REAL(BTPT) RETURN END REAL FUNCTION F92C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F92C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F92C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S17C09(OMEGA,TAU,BPPT) F92C09=REAL(BPPT) RETURN END REAL FUNCTION F93C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F93C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F93C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S18C09(OMEGA,TAU,BVPT) F93C09=REAL(BVPT) RETURN END REAL FUNCTION F94C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F94C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF (ILL90.NE.0) RETURN F94C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) IF (OMEGA.LE.-1.0D+10) RETURN CALL S19C09(OMEGA,TAU,AJTPT) F94C09=REAL(AJTPT) RETURN END REAL FUNCTION F95C09(FP,FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06,TCR=144.414D+00) F95C09=-1.0E+20 CALL S90C09(FP,FT,ILL90) IF(ILL90.NE.0) RETURN F95C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR TAU=TCR/(DBLE(FT)+273.15D0) CALL S5C09(PI,TAU,OMEGA) CALL S12C09(OMEGA,TAU,CPPT) CALL S11C09(OMEGA,TAU,CVPT) GAMPT=CPPT/CVPT F95C09=REAL(GAMPT) RETURN END REAL FUNCTION F96C09(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PCR=5.23952D+06) F96C09=-1.0E+20 IF ((FP.LT.2.389E-03).OR.(FP.GT.52.3953E+00)) RETURN F96C09=-1.0E+10 PI=DBLE(FP)*1.0D+05/PCR CALL S4C09(PI,OMEGAL,OMEGAG,TAUS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S12C09(OMEGAG,TAUS,CPG) CALL S11C09(OMEGAG,TAUS,CVG) GAMPDD=CPG/CVG F96C09=REAL(GAMPDD) RETURN END REAL FUNCTION F97C09(FT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TCR=144.414D+00) F97C09=-1.0E+20 IF ((FT.LT.-219.670).OR.(FT.GT.-128.735)) RETURN F97C09=-1.0E+10 TAU=TCR/(DBLE(FT)+273.15D0) CALL S3C09(TAU,OMEGAL,OMEGAG,PIS) IF (OMEGAG.LE.-1.0D+10) RETURN CALL S12C09(OMEGAG,TAU,CPG) CALL S11C09(OMEGAG,TAU,CVG) GAMTDD=CPG/CVG F97C09=REAL(GAMTDD) RETURN END C F98 + PSEUDO BOILING POINT AT P FUNCTION F98C09(P) DIMENSION T(2),C(2),TL(3),TR(3),CL(3),CR(3) P1=F21C09('P') PP=ABS((P-P1)/P1) T1=F21C09('T') IF (PP.LT.1.0E-5) THEN F98C09=T1 RETURN ENDIF IF (P.LT.P1.OR.P.GT.200.001D00) THEN F98C09=-1.0E+20 RETURN ENDIF T2=0.85*T1 P2=F30C09(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.05D00 IREP=0 IREM=5000 KCONT=0 ICONT=0 C(1)=F18C09(P,T(1)) T(2)=T(1)-DEL C(2)=F18C09(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)=F18C09(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=-F18C09(P,TT) 3000 CONV=ABS((CC-C(2))/CC) IF(CONV.LT.EPS) THEN F98C09=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=-F18C09(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=F18C09(P,TA) CB=F18C09(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=F18C09(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)=F18C09(P,TL(2)) CR(2)=F18C09(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=F18C09(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=F18C09(P,TB) ELSE TB=TR(MR+1) CB=CR(MR+1) ENDIF 6135 TC=TR(MR) CC=CR(MR) ENDIF GO TO 6050 7000 F98C09=TC RETURN 8000 F98C09=-1.0E+10 RETURN END *======================================================================* SUBROUTINE S1C09(T,PMELT) ***** SUBROUTINE TO CALCULATE THE MELTING PRESSURE CORRESPONDING ***** ***** TO THE INPUT TEMPERATURE. (BY EQ.(4.26)) ***** **** INPUT : T=TEMPERATURE IN (K) **** OUTPUT : PMELT=MELTING PRESSURE IN (PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(P0=249.975D+6, PT=2.52D+02, - TT=53.4811D+00, C=2.1845D+00) IF (ABS(T/TT-1.0).LT.2.0D-06) THEN PMELT=2.390D+02 RETURN END IF PMELT=PT+P0*((T/TT)**C-1.0D0) RETURN END SUBROUTINE S2C09(P,TMELT) ***** SUBROUTINE TO FIND THE MELTING TEMPERATURE CORRESPONDING TO ***** ***** THE INPUT PRESSURE BY NEWTON METHOD. (BY EQ.(4.26)) ***** **** INPUT : P=PRESSURE IN (PA) **** OUTPUT : TMELT=MELTING TEMPERATURE IN (K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(P0=249.975D+6, PT=2.52D+02, - TT=53.4811D+00, C=2.1845D+00) TMELT=TT*(1+((P-PT)/P0))**(1.0D0/C) RETURN END C======================================================================= SUBROUTINE S3C09(TAU,OMEGAL,OMEGAG,PIS) ***** SUBROUTINE TO FIND THE LIQUID AND VAPOR DENSITIES ***** ***** CORRESPONDING TO THE INPUT TEMPERATURE BY SOLVING THE MAXWELL***** ***** RELATION (EQ.(4.25)) USING NEWTON METHOD AND CALCULATE THE ***** ***** SATURATED PRESSURE. ***** **** INPUT : TAU=TCR/T **** OUTPUT1 : OMEGAL=RHOL/RHOCR **** OUTPUT2 : OMEGAG=RHOG/RHOCR **** OUTPUT3 : PIS=PS/PCR IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00, - PCR=5.23952D+06, TCR=144.414D+00, RHOCR=1.56030D+04) PARAMETER(EPS=1.0D-12,ITMAX=10000) ***** NUMERICAL CONSTANTS OF EQ.(3.2)*********************************** DATA C1/ 1.7036719D+00/, C2/ 2.2817780D+00/, - C3/-4.2038623D+00/, C4/ 2.5128703D+00/, - C5/-5.1162174D-01/ ***** NUMERICAL CONSTANTS OF EQ.(3.3)*********************************** DATA D1/-1.8797889D+00/, D2/ 1.7137650D-01/, - D3/-6.1655536D+00/, D4/ 2.9709048D+00/, - D5/-6.3193094D+00/, D6/ 7.6558470D+00/, - D7/-1.9018882D+01/, D8/ 2.2233488D+01/ - D9/-8.8726855D+00/, D10/ 2.9633564D+00/ ************************************************************************ OMEGAL=-1.0D+20 OMEGAG=-1.0D+20 PIS=-1.0D+20 IF (TAU.LE.0.0D0) RETURN OMEGAL=1.0D0 OMEGAG=1.0D0 PIS=1.0D0 IF (ABS(TAU-1.0D0).LE.2.0D-06) RETURN *** INITIAL GUESSES FOR OML AND OMG ********************************** ***** NEW CRITICAL PARAMETERS TO USE EQS.(3.1), (3.2) AND (3.3) C TC=144.121D+00 C PC=5.1724D+06 C RHOC=0.015027D+00 * TH=ABS(TAU-1.0D0) THQ=TH**(1.0D0/3.0D0) OMLW=1.0D0+(C1+(C2+(C3+(C4+C5*THQ)*THQ)*THQ)*THQ**2)*THQ OMGW=EXP((D1+(D2+(D3+(D4+(D5+(D6+(D7+(D8+(D9+D10*THQ)*THQ**2) - *THQ)*THQ**5)*THQ**2)*THQ**5)*THQ)*THQ)*THQ)*THQ) * *** NEWTON METHOD *** DO 1 IT=1,ITMAX CALL S20C09(OMLW,TAU,ALPHAL) CALL S20C09(OMGW,TAU,ALPHAG) CALL S21C09(OMLW,TAU,DADOL) CALL S21C09(OMGW,TAU,DADOG) CALL S22C09(OMLW,TAU,D2ADOL) CALL S22C09(OMGW,TAU,D2ADOG) * F=LOG(OMLW/OMGW)+(ALPHAL-ALPHAG)+(OMLW*DADOL-OMGW*DADOG) G=(OMLW-OMGW)+(OMLW*OMLW*DADOL-OMGW*OMGW*DADOG) DOMLW=(OMLW/(OMLW-OMGW))*(G-OMGW*F)/ - (1.0D0+2.0D0*OMLW*DADOL+OMLW*OMLW*D2ADOL) DOMGW=(OMGW/(OMLW-OMGW))*(G-OMLW*F)/ - (1.0D0+2.0D0*OMGW*DADOG+OMGW*OMGW*D2ADOG) OMLW=OMLW-DOMLW OMGW=OMGW-DOMGW IF ((ABS(DOMLW).LT.EPS).AND.(ABS(DOMGW).LT.EPS)) THEN OMEGAL=OMLW OMEGAG=OMGW PIS=((GASCON*RHOCR*TCR/PCR)/TAU) - *(OMLW*OMGW/(OMLW-OMGW))*(LOG(OMLW/OMGW)+ALPHAL-ALPHAG) RETURN END IF 1 CONTINUE OMEGAL=-1.0D+10 OMEGAG=-1.0D+10 PIS=-1.0D+10 RETURN END SUBROUTINE S4C09(PI,OMEGAL,OMEGAG,TAUS) ***** SUBROUTINE TO FIND THE SATURATED TEMPERATURE, LIQUID DENSITY ***** ***** AND VAPOR DENSITY CORRESPONDING TO THE INPUT PRESSURE BY ***** ***** SOLVING THE MAXWELL RELATION (EQ.(4.25)) USING AN ITERATIVE ***** ***** METHOD. ***** **** INPUT : PI=P/PCR **** OUTPUT1 : OMEGAL=RHOL/RHOCR **** OUTPUT2 : OMEGAG=RHOG/RHOCR **** OUTPUT3 : TAUS=TCR/TS IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00, - PCR=5.23952D+06, TCR=144.414D+00, RHOCR=1.56030D+04) PARAMETER(EPS1=1.0D-12,EPS2=1.0D-12,IT1MAX=10000,IT2MAX=10000) ***** NUMERICAL CONSTANS OF EQ.(3.1)************************************ DATA B1/-6.1843496D+00/, B2/ 1.2095388D+00/, - B3/-6.1708665D+00/, B4/-3.0963258D-01/, - B5/ 6.8976007D-02/ ************************************************************************ OMEGAL=-1.0D+20 OMEGAG=-1.0D+20 TAUS=-1.0D+20 IF (PI.LE.0.0D0) RETURN OMEGAL=1.0D0 OMEGAG=1.0D0 TAUS=1.0D0 IF (ABS(PI-1.0D0).LE.1.6D-05) RETURN *** INITIAL GUESS FOR TAU *** IF (PI.GT.0.02D0) THEN TH=LOG(PI)/(B1-LOG(PI)) ELSE IF ((PI.GT.2.84D-04).AND.(PI.LT.0.02D0)) THEN TH=0.5D0 ELSE TH=1.5D0 GO TO 10 END IF * DO 1 IT1=1,IT1MAX F=(B1-LOG(PI))*TH+B2*TH**1.5D0+B3*TH**2+B4*TH**3.5D0 - +B5*TH**4.5D0-LOG(PI) DF=B1-LOG(PI)+1.5D0*B2*SQRT(TH)+2.0D0*B3*TH+3.5D0*B4*TH**2.5D0 - +4.5D0*B5*TH**3.5D0 TH=TH-F/DF IF(ABS(-F/(DF*TH)).LT.EPS1) GO TO 10 1 CONTINUE GO TO 1000 10 TAUSW=1.0D0+TH ***** ITERATIVE METHOD *** DO 2 IT2=1,IT2MAX CALL S3C09(TAUSW,OMLW,OMGW,PISW) DP=(PI-PISW)*PCR IF (ABS(DP/(PI*PCR)).LE.EPS2) THEN OMEGAL=OMLW OMEGAG=OMGW TAUS=TAUSW RETURN END IF CALL S20C09(OMLW,TAUSW,ALPHAL) CALL S20C09(OMGW,TAUSW,ALPHAG) CALL S23C09(OMLW,TAUSW,DADTL) CALL S23C09(OMGW,TAUSW,DADTG) * DPSDT=GASCON*RHOCR* - (OMLW*OMGW/(OMLW-OMGW))*(LOG(OMLW/OMGW)+(ALPHAL-ALPHAG) - -TAUSW*(DADTL-DADTG)) DTSW=DP/DPSDT TAUSW=1.0D0/(1.0D0/TAUSW+DTSW/TCR) 2 CONTINUE 1000 OMEGAG=-1.0D+10 OMEGAL=-1.0D+10 TAUS=-1.0D+10 RETURN END SUBROUTINE S5C09(PI,TAU,OMEGA) ***** SUBROUTINE TO FIND THE DENSITY CORRESPONDING TO THE INPUT ***** ***** PRESSURE AND TEMPERATURE BY SOLVING THE EQUATION OF STATE, ***** ***** EQ.(4.4) USING NEWTON METHOD. ***** **** INPUT1 : PI= P/PCR **** INPUT2 : TAU=TCR/T **** OUTPUT : OMEGA=RHO/RHOCR IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00, - PCR=5.23952D+06, TCR=144.414D+00, RHOCR=1.56030D+04) PARAMETER(EPS=1.0D-12,ITMAX=10000) * OMEGA=1.0D0 IF ((ABS(PI-1.0D0).LE.1.0D-05).AND.(ABS(TAU-1.0D0).LE.1.0D-05)) - RETURN *** INITIAL GUESS FOR OMEGA *** IF ((PI.GE.0.998D0).AND.(PI.LE.1.0001D0).AND. - (TAU.GE.0.996D0).AND.(TAU.LE.1.0005D0)) THEN OM=0.95D0 GO TO 10 END IF * ZCR=PCR/(RHOCR*GASCON*TCR) IF ((PI.GT.1.0D0).AND.(PI.LE.3.82D0).AND.(TAU.LE.1.0D0)) THEN PILOG=LOG(PI) OM=0.2828D0*PILOG+0.1443D0 ELSE IF ((PI.LE.1.0D0).AND.(TAU.LE.1.0D0)) THEN OM=ZCR*PI*TAU ELSE IF ((PI.GE.1.0D0).AND.(TAU.GT.1.0D0)) THEN OM=3.0D0 ELSE CALL S4C09(PI,OML,OMG,TAUS) IF (OML.LE.-1.0D+10) GO TO 1000 IF ((PI.LT.1.0D0).AND.(TAU.LE.TAUS)) THEN IF (ABS(TAUS/TAU-1.0D0).LT.1.0D-05) THEN OMEGA=OMG RETURN END IF OM=OMG ELSE IF ((PI.LT.1.0D0).AND.(TAU.GT.TAUS)) THEN IF (ABS(TAU/TAUS-1.0D0).LT.1.0D-05) THEN OMEGA=OML RETURN END IF OM=3.0D0 END IF END IF * 10 CONTINUE * ZCRPT=ZCR*PI*TAU *** NEWTON METHOD *** DO 1 IT=1,ITMAX CALL S21C09(OM,TAU,DADO) CALL S22C09(OM,TAU,D2ADO2) Z=(1.0D0+OM*DADO)*OM DOMZ=1.0D0+(2.0D0*DADO+OM*D2ADO2)*OM F=ZCRPT-Z DF=-DOMZ OM=OM-F/DF IF (ABS(-F/(DF*OM)).LT.EPS) THEN OMEGA=OM RETURN END IF 1 CONTINUE 1000 OMEGA=-1.0D+10 RETURN END C======================================================================= SUBROUTINE S6C09(PI,OMEGA,TAU) ***** SUBROUTINE TO FIND THE TEMPERATURE CORRESPONDING TO THE INPUT***** ***** PRESSURE AND DENSITY BY SOLVING THE EQUATION OF STATE, ***** ***** EQ.(4.4) USING NEWTON METHOD. ***** **** INPUT1 : PI=P/PCR **** INPUT2 : OMEGA=RHO/RHOCR **** OUTPUT : TAU=TCR/T IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00, - PCR=5.23952D+06, TCR=144.414D+00, RHOCR=1.56030D+04) PARAMETER(EPS=1.0D-12,ITMAX=10000) PARAMETER(DELTA=1.0383171848D+00) DIMENSION COA1(1:31),COA2(1:31) DIMENSION CN(1:31),CR(1:31),CS(1:31),CQ(14:31) * TAU=1.0D0 IF ((ABS(PI-1.0D0).LE.1.0D-05).AND. - (ABS(OMEGA-1.0D0).LE.1.0D-05)) RETURN * **** CALCULATION OF THE COEFFICIENTS CALL S26C09(CN,CR,CS,CQ) DOM2=(DELTA*OMEGA)*(DELTA*OMEGA) DO 1 I=1,13 COA1(I)=CN(I)*CR(I)*(OMEGA**(CR(I)-1.0D0)) 1 CONTINUE DO 2 I=14,31 COA2(I)=CN(I)*(OMEGA**(CR(I)-1.0D0)) - *(CR(I)-2.0D0*CQ(I)*DOM2)*EXP(-CQ(I)*DOM2) 2 CONTINUE * ZCR=PCR/(RHOCR*GASCON*TCR) ZCRPOM=ZCR*PI/OMEGA *** INITIAL GUESS FOR TAU *** IF (PI.GT.1.0D+00) THEN TAUW=TCR/300.0D0 ELSE IF (OMEGA.GT.1.0D0) THEN CALL S4C09(PI,OML,OMG,TAUS) TAUW=TAUS ELSE TAUW=OMEGA/(ZCR*PI) END IF END IF * DO 10 IT=1,ITMAX * SZA1=0.0D0 SZB1=0.0D0 DO 101 I=1,13 SZA1=SZA1+COA1(I)*(TAUW**CS(I)) SZB1=SZB1+COA1(I)*(1.0D0-CS(I))*(TAUW**CS(I)) 101 CONTINUE SZA2=0.0D0 SZB2=0.0D0 DO 102 I=14,31 SZA2=SZA2+COA2(I)*(TAUW**CS(I)) SZB2=SZB2+COA2(I)*(1.0D0-CS(I))*(TAUW**CS(I)) 102 CONTINUE DADO=SZA1+SZA2 DZIN=SZB1+SZB2 * Z=(1.0D0+OMEGA*DADO)/TAUW DZTAUW=-(1.0D0+OMEGA*DZIN)/(TAUW*TAUW) F=ZCRPOM-Z DF=-DZTAUW TAUW=TAUW-F/DF IF (ABS(-F/(DF*TAUW)).LT.EPS) THEN TAU=TAUW RETURN END IF 10 CONTINUE TAU=-1.0D+10 RETURN END C======================================================================= SUBROUTINE S7C09(OMEGA,TAU,PI) ***** SUBROUTINE TO CALCULATE THE PRESSURE CORRESPONDING TO THE ***** ***** INPUT TEMPERATURE AND DENSITY FROM THE EQUATION OF STATE, ***** ***** EQ.(4.4). ***** **** INPUT1 : OMEGA=RHO/RHOCR RHO AND RHOCR IN (MOL/M**3) **** INPUT2 : TAU=TCR/T T AND TCR IN (K) **** OUTPUT : PI=P/PCR P AND OCR IN (PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00, - PCR=5.23952D+06, TCR=144.414D+00, RHOCR=1.56030D+04) CALL S21C09(OMEGA,TAU,DADO) P=(OMEGA*RHOCR)*GASCON*(TCR/TAU)*(1.0D0+OMEGA*DADO) PI=P/PCR RETURN END SUBROUTINE S8C09(OMEGA,TAU,S) ***** SUBROUTINE TO CALCULATE THE ENTROPY. (EQ.(4.7)) ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : S=ENTROPY IN (J/(K*MOL)) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00,OMEGAA=0.0025853657D+00) CALL S27C09(TAU,ALPID) CALL S28C09(TAU,DATID) CALL S20C09(OMEGA,TAU,ALPHA) CALL S23C09(OMEGA,TAU,DADT) S=GASCON*(TAU*(DATID+DADT)-ALPHA-LOG(OMEGA/OMEGAA)-ALPID) RETURN END SUBROUTINE S9C09(OMEGA,TAU,U) ***** SUBROUTINE TO CALCULATE THE INTERNAL ENERGY. (EQ.(4.8)) ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : U=INTERNAL ENERGY IN (J/MOL) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00,TCR=144.414D+00) CALL S28C09(TAU,DATID) CALL S23C09(OMEGA,TAU,DADT) U=GASCON*(TCR/TAU)*(TAU*(DATID+DADT)) RETURN END SUBROUTINE S10C09(OMEGA,TAU,H) ***** SUBROUTINE TO CALCULATE THE ENTHALPY. (EQ.(4.8)) ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : H=ENTHALPY IN (J/MOL) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00,TCR=144.414D+00) CALL S28C09(TAU,DATID) CALL S21C09(OMEGA,TAU,DADO) CALL S23C09(OMEGA,TAU,DADT) H=GASCON*(TCR/TAU)*(1.0D0+TAU*(DATID+DADT)+OMEGA*DADO) RETURN END SUBROUTINE S11C09(OMEGA,TAU,CV) ***** SUBROUTINE TO CALCULATE THE ISOCHORIC HEAT CAPACITY.(EQ.(4.11))*** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : CV=ISOCHORIC HEAT CAPACITY IN (J/(K*MOL)) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00) CALL S29C09(TAU,DAT2ID) CALL S24C09(OMEGA,TAU,D2ADT2) CV=GASCON*(-TAU*TAU*(DAT2ID+D2ADT2)) RETURN END SUBROUTINE S12C09(OMEGA,TAU,CP) ***** SUBROUTINE TO CALCULATE THE ISOBARIC HEAT CAPACITY. (EQ.(4.13))*** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : CP=ISOBARIC HEAT CAPACITY IN (J/(K*MOL)) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00) CALL S11C09(OMEGA,TAU,CV) CALL S21C09(OMEGA,TAU,DADO) CALL S22C09(OMEGA,TAU,D2ADO2) CALL S25C09(OMEGA,TAU,D2ADOT) CP=CV+GASCON*(1.0D0+OMEGA*DADO-OMEGA*TAU*D2ADOT)**2 - /(1.0D0+2.0D0*OMEGA*DADO+OMEGA*OMEGA*D2ADO2) RETURN END * SUBROUTINE S13C09(OMEGA,TAU,AK) ***** SUBROUTINE TO CALCULATE THE ISENTROPIC EXPONENT. *** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : AK=ISENTROPIC EXPONENT (DIMENSIONLESS) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CALL S11C09(OMEGA,TAU,CV) CALL S12C09(OMEGA,TAU,CP) CALL S21C09(OMEGA,TAU,DADO) CALL S22C09(OMEGA,TAU,D2ADO2) AK=(CP/CV)*(1.0D0+2.0D0*OMEGA*DADO+OMEGA*OMEGA*D2ADO2) - /(1.0D0+OMEGA*DADO) RETURN END SUBROUTINE S14C09(OMEGA,TAU,W) ***** SUBROUTINE TO CALCULATE THE VELOCITY OF SOUND. (EQ.(4.17))*** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : W=VELOCITY OF SOUND IN (METER/SECOND) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASCON=8.31448D+00,TCR=144.414D+00) PARAMETER(AKGMOL=0.03799681D+00) CALL S21C09(OMEGA,TAU,DADO) CALL S22C09(OMEGA,TAU,D2ADO2) CALL S25C09(OMEGA,TAU,D2ADOT) CALL S24C09(OMEGA,TAU,D2ADT2) CALL S29C09(TAU,DAT2ID) W2RTM=(1.0D0+2.0D0*OMEGA*DADO+OMEGA*OMEGA*D2ADO2) - -(1.0D0+OMEGA*DADO-OMEGA*TAU*D2ADOT)**2 - /(TAU*TAU*(DAT2ID+D2ADT2)) W=SQRT((GASCON*(TCR/TAU)/AKGMOL)*W2RTM) RETURN END SUBROUTINE S15C09(OMEGA,TAU,BS) ***** SUBROUTINE TO CALCULATE THE ISENTROPIC COMPRESSIBILITY ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : BS=ISENTROPIC COMPRESSIBILITY IN (1/PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(RHOCR=1.56030D+04) PARAMETER(AKGMOL=0.03799681D+00) CALL S14C09(OMEGA,TAU,W) BS=1.0D0/((OMEGA*RHOCR*AKGMOL)*W**2) RETURN END SUBROUTINE S16C09(OMEGA,TAU,BT) ***** SUBROUTINE TO CALCULATE THE ISOTHERMAL COMPRESSIBILITY ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : BT=ISOTHERMAL COMPRESSIBILITY IN (1/PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CALL S11C09(OMEGA,TAU,CV) CALL S12C09(OMEGA,TAU,CP) CALL S15C09(OMEGA,TAU,BS) BT=(CP/CV)*BS RETURN END SUBROUTINE S17C09(OMEGA,TAU,BP) ***** SUBROUTINE TO CALCULATE THE VOLUMETRIC COEFFICIENT OF EXPANSION **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : BP= VOLUMETRIC COEFFICIENT OF EXPANSION IN (1/K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(TCR=144.414D+00) CALL S21C09(OMEGA,TAU,DADO) CALL S22C09(OMEGA,TAU,D2ADO2) CALL S25C09(OMEGA,TAU,D2ADOT) BP=(TAU/TCR)*(1.0D0+OMEGA*DADO-OMEGA*TAU*D2ADOT) - /(1.0D0+2.0D0*OMEGA*DADO+OMEGA*OMEGA*D2ADO2) RETURN END SUBROUTINE S18C09(OMEGA,TAU,BV) ***** SUBROUTINE TO CALCULATE THE PRESSURE COEFFICIENT ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : BV= PRESSURE COEFFICIENT IN (1/K) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(TCR=144.414D+00) CALL S21C09(OMEGA,TAU,DADO) CALL S25C09(OMEGA,TAU,D2ADOT) BV=(TAU/TCR)*(1.0D0+OMEGA*DADO-OMEGA*TAU*D2ADOT) - /(1.0D0+OMEGA*DADO) RETURN END SUBROUTINE S19C09(OMEGA,TAU,AJT) ***** SUBROUTINE TO CALCULATE THE JOULE-THOMSON COEEFICIENT ********* **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : AJT=JOULE-THOMSON COEEFICIENT IN (K/PA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(RHOCR=1.56030D+04,TCR=144.414D+00) CALL S12C09(OMEGA,TAU,CP) CALL S17C09(OMEGA,TAU,BP) AJT=(1.0D0/(OMEGA*RHOCR*CP))*((TCR/TAU)*BP-1.0D0) RETURN END SUBROUTINE S20C09(OMEGA,TAU,ALPHA) ***** SUBROUTINE TO CALCULATE ALPHA ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : ALPHA IMPLICIT DOUBLE PRECISION(A-H,N,O-Z) PARAMETER(DELTA=1.0383171848D+00) DIMENSION N(1:31),R(1:31),S(1:31),Q(14:31) CALL S26C09(N,R,S,Q) SOUT1=0.0D0 DO 1 I=1,13 SOUT1=SOUT1+N(I)*(OMEGA**R(I))*(TAU**S(I)) 1 CONTINUE SOUT2=0.0D0 DOM2=(DELTA*OMEGA)*(DELTA*OMEGA) DO 2 I=14,31 SOUT2=SOUT2+N(I)*(OMEGA**R(I))*(TAU**S(I))*EXP(-Q(I)*DOM2) 2 CONTINUE ALPHA=SOUT1+SOUT2 RETURN END SUBROUTINE S21C09(OMEGA,TAU,DADO) ***** SUBROUTINE TO CALCULATE D(ALPHA)/D(OMEGA) ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : DADO=D(ALPHA)/D(OMEGA) IMPLICIT DOUBLE PRECISION(A-H,N,O-Z) PARAMETER(DELTA=1.0383171848D+00) DIMENSION N(1:31),R(1:31),S(1:31),Q(14:31) CALL S26C09(N,R,S,Q) SOUT1=0.0D0 DO 1 I=1,13 SOUT1=SOUT1+N(I)*R(I)*(OMEGA**(R(I)-1.0D0))*(TAU**S(I)) 1 CONTINUE SOUT2=0.0D0 DOM2=(DELTA*OMEGA)*(DELTA*OMEGA) DO 2 I=14,31 SOUT2=SOUT2+N(I)*(OMEGA**(R(I)-1.0D0))*(TAU**S(I)) - *(R(I)-2.0D0*Q(I)*DOM2)*EXP(-Q(I)*DOM2) 2 CONTINUE DADO=SOUT1+SOUT2 RETURN END SUBROUTINE S22C09(OMEGA,TAU,D2ADO2) ***** SUBROUTINE TO CALCULATE D2(ALPHA)/D(OMEGA)2 ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : D2AD02=D2(ALPHA)/D(OMEGA)2 IMPLICIT DOUBLE PRECISION(A-H,N,O-Z) PARAMETER(DELTA=1.0383171848D+00) DIMENSION N(1:31),R(1:31),S(1:31),Q(14:31) CALL S26C09(N,R,S,Q) SOUT1=0.0D0 DO 1 I=1,13 SOUT1=SOUT1+N(I)*R(I)*(R(I)-1.0D0)*(OMEGA**(R(I)-2.0D0)) - *(TAU**S(I)) 1 CONTINUE SOUT2=0.0D0 DOM2=(DELTA*OMEGA)*(DELTA*OMEGA) DO 2 I=14,31 SOUT2=SOUT2+N(I)*(OMEGA**(R(I)-2.0D0))*(TAU**S(I)) - *((R(I)-2.0D0*Q(I)*DOM2)*(R(I)-1.0D0-2.0D0*Q(I)*DOM2) - -4.0D0*Q(I)*DOM2)*EXP(-Q(I)*DOM2) 2 CONTINUE D2ADO2=SOUT1+SOUT2 RETURN END SUBROUTINE S23C09(OMEGA,TAU,DADT) ***** SUBROUTINE TO CALCULATE D(ALPHA)/D(TAU) ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : DADT=D(ALPHA)/D(TAU) IMPLICIT DOUBLE PRECISION(A-H,N,O-Z) PARAMETER(DELTA=1.0383171848D+00) DIMENSION N(1:31),R(1:31),S(1:31),Q(14:31) CALL S26C09(N,R,S,Q) SOUT1=0.0D0 DO 1 I=1,13 SOUT1=SOUT1+N(I)*S(I)*(OMEGA**R(I))*(TAU**(S(I)-1.0D0)) 1 CONTINUE SOUT2=0.0D0 DOM2=(DELTA*OMEGA)*(DELTA*OMEGA) DO 2 I=14,31 SOUT2=SOUT2+N(I)*S(I)*(OMEGA**R(I))*(TAU**(S(I)-1.0D0)) - *EXP(-Q(I)*DOM2) 2 CONTINUE DADT=SOUT1+SOUT2 RETURN END SUBROUTINE S24C09(OMEGA,TAU,D2ADT2) ***** SUBROUTINE TO CALCULATE D2(ALPHA)/D(TAU)2 ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : D2ADT2=D2(ALPHA)/D(TAU)2 IMPLICIT DOUBLE PRECISION(A-H,N,O-Z) PARAMETER(DELTA=1.0383171848D+00) DIMENSION N(1:31),R(1:31),S(1:31),Q(14:31) CALL S26C09(N,R,S,Q) SOUT1=0.0D0 DO 1 I=1,13 SOUT1=SOUT1+N(I)*S(I)*(S(I)-1.0D0)*(OMEGA**R(I)) - *(TAU**(S(I)-2.0D0)) 1 CONTINUE SOUT2=0.0D0 DOM2=(DELTA*OMEGA)*(DELTA*OMEGA) DO 2 I=14,31 SOUT2=SOUT2+N(I)*S(I)*(S(I)-1.0D0) - *(OMEGA**R(I))*(TAU**(S(I)-2.0D0)) - *EXP(-Q(I)*DOM2) 2 CONTINUE D2ADT2=SOUT1+SOUT2 RETURN END SUBROUTINE S25C09(OMEGA,TAU,D2ADOT) ***** SUBROUTINE TO CALCULATE D2(ALPHA)/D(OMEGA)D(TAU) ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT1 : OMEGA=RHO/RHOCR **** INPUT2 : TAU=TCR/T **** OUTPUT : D2ADOT=D2(ALPHA)/D(OMEGA)D(TAU) IMPLICIT DOUBLE PRECISION(A-H,N,O-Z) PARAMETER(DELTA=1.0383171848D+00) DIMENSION N(1:31),R(1:31),S(1:31),Q(14:31) CALL S26C09(N,R,S,Q) SOUT1=0.0D0 DO 1 I=1,13 SOUT1=SOUT1+N(I)*R(I)*S(I)*(OMEGA**(R(I)-1.0D0)) - *(TAU**(S(I)-1.0D0)) 1 CONTINUE SOUT2=0.0D0 DOM2=(DELTA*OMEGA)*(DELTA*OMEGA) DO 2 I=14,31 SOUT2=SOUT2+N(I)*S(I)*(OMEGA**(R(I)-1.0D0))*(TAU**(S(I)-1.0D0)) - *(R(I)-2.0D0*Q(I)*DOM2)*EXP(-Q(I)*DOM2) 2 CONTINUE D2ADOT=SOUT1+SOUT2 RETURN END SUBROUTINE S26C09(AN,AR,AS,AQ) ***** SUBROUTINE TO SUPPLY THE COEFFICIENTS N(I), R(I), S(I), Q(I) ***** ***** LISTED IN TABLE 4.B ON PAGE 38. ***** **** OUTPUT : AN(I) ( I= 1...31 ) FOR N **** OUTPUT : AR(I) ( I= 1...31 ) FOR R **** OUTPUT : AS(I) ( I= 1...31 ) FOR S **** OUTPUT : AQ(I) ( I=14...31 ) FOR Q IMPLICIT DOUBLE PRECISION(A,N,R,S) DIMENSION AN(1:31),N(1:31) DIMENSION AR(1:31),R(1:31),AS(1:31),S(1:31),AQ(14:31) DATA N( 1)/ 1.51144749736D+00/, N( 2)/-2.98666288409D+00/, - N( 3)/ 3.29644905098D+00/, N( 4)/-2.98458624201D+00/, - N( 5)/-2.28688966459D+00/, N( 6)/-1.09492193400D+00/, - N( 7)/ 3.04775277572D+00/, N( 8)/ 1.15689564208D-01/, - N( 9)/-1.16100171627D+00/, N(10)/ 2.95656394476D-01/, - N(11)/ 7.11482542928D-02/, N(12)/-1.71363832155D-03/, - N(13)/ 6.65317955515D-04/, N(14)/ 5.06026676251D+00/, - N(15)/-6.29268435440D+00/, N(16)/ 6.17784808739D+00/, - N(17)/-1.55366191788D+00/, N(18)/-2.87170687343D+00/, - N(19)/ 3.17214480494D+00/, N(20)/-2.67969025215D+00/, - N(21)/ 2.71865479252D+00/, N(22)/-1.07191065039D+00/, - N(23)/ 1.26597342291D+00/, N(24)/-7.06244695489D-01/, - N(25)/ 2.68707888826D-01/, N(26)/ 5.27251190274D-02/, - N(27)/ 5.44411481926D-02/, N(28)/ 2.28949994105D-04/, - N(29)/-5.47908264304D-10/, N(30)/-9.64273224950D-02/, - N(31)/ 3.68084486225D-04/ DATA R( 1)/1.0D0/, R( 2)/1.0D0/, R( 3)/1.0D0/, R( 4)/1.0D0/, - R( 5)/2.0D0/, R( 6)/2.0D0/, R( 7)/3.0D0/, R( 8)/3.0D0/, - R( 9)/4.0D0/, R(10)/4.0D0/, R(11)/5.0D0/, R(12)/8.0D0/, - R(13)/9.0D0/, R(14)/2.0D0/, R(15)/2.0D0/, R(16)/2.0D0/, - R(17)/2.0D0/, R(18)/3.0D0/, R(19)/3.0D0/, R(20)/3.0D0/, - R(21)/4.0D0/, R(22)/4.0D0/, R(23)/4.0D0/, R(24)/5.0D0/, - R(25)/6.0D0/, R(26)/7.0D0/, R(27)/8.0D0/, R(28)/12.0D0/, - R(29)/4.0D0/, R(30)/6.0D0/, R(31)/6.0D0/ DATA S( 1)/0.0D0/, S( 2)/0.5D0/, S( 3)/1.5D0/, S( 4)/2.0D0/, - S( 5)/0.5D0/, S( 6)/1.0D0/, S( 7)/0.5D0/, S( 8)/2.0D0/, - S( 9)/0.5D0/, S(10)/1.0D0/, S(11)/0.0D0/, S(12)/0.5D0/, - S(13)/0.0D0/, S(14)/1.0D0/, S(15)/3.0D0/, S(16)/4.0D0/, - S(17)/5.0D0/, S(18)/1.0D0/, S(19)/4.0D0/, S(20)/5.0D0/, - S(21)/1.0D0/, S(22)/3.0D0/, S(23)/5.0D0/, S(24)/4.0D0/, - S(25)/4.0D0/, S(26)/1.0D0/, S(27)/1.0D0/, S(28)/5.0D0/, - S(29)/30.0D0/, S(30)/20.0D0/, S(31)/25.0D0/ * DO 1 I=1,31 AN(I)=N(I) AR(I)=R(I) AS(I)=S(I) 1 CONTINUE DO 2 I=14,28 AQ(I)=1.0D0 2 CONTINUE AQ(29)=2.0D0 AQ(30)=3.0D0 AQ(31)=3.0D0 RETURN END SUBROUTINE S27C09(TAU,ALPID) ***** SUBROUTINE TO CALCULATE ALPHA FOR THE IDEAL GAS STATE ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT : TAU=TCR/T **** OUTPUT : ALPID=ALPHA(ID) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION F(1:8) CALL S30C09(F) ALPID=(F(1)/TAU+F(2))/TAU/TAU/TAU + (F(3)+F(4)*TAU)*TAU - +F(5)*LOG(TAU)+F(6)*LOG(EXP(F(7)*TAU)-1.0D0)+F(8) RETURN END SUBROUTINE S28C09(TAU,DATID) ***** SUBROUTINE TO CALCULATE D(ALPHA)/D(TAU) ***** FOR THE IDEAL GAS STATE ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT : TAU=TCR/T **** OUTPUT : DATID=D(ALPHA)/D(TAU)(ID) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION F(1:8) CALL S30C09(F) DATID=-(4.0D0*F(1)/TAU+(3.0D0*F(2)))/TAU/TAU/TAU/TAU - +F(3)+2.0D0*F(4)*TAU+F(5)/TAU - +F(6)*F(7)*EXP(F(7)*TAU)/(EXP(F(7)*TAU)-1.0D0) RETURN END SUBROUTINE S29C09(TAU,DAT2ID) ***** SUBROUTINE TO CALCULATE D2(ALPHA)/D(TAU)2 ***** ***** FOR THE IDEAL GAS STATE ***** ***** LISTED IN TABLE 4.D ON PAGE 39. ***** **** INPUT : TAU=TCR/T **** OUTPUT : DAT2ID=D2(ALPHA)/D(TAU)2(ID) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION F(1:8) CALL S30C09(F) DAT2ID=(20.0D0*F(1)/TAU+(12.0D0*F(2)))/TAU/TAU/TAU/TAU/TAU - +2.0D0*F(4) - F(5)/TAU/TAU - -F(6)*F(7)*F(7)*EXP(F(7)*TAU) - /(EXP(F(7)*TAU)-1.0D0)**2 RETURN END SUBROUTINE S30C09(AF) ***** SUBROUTINE TO SUPPLY THE COEFFICIENTS F(I) OF EQ.(4.3) ***** ***** LISTED IN TABLE 4.C ON PAGE 38. ***** **** OUTPUT : AF(I) ( I=1...8 ) FOR F IMPLICIT DOUBLE PRECISION(A,F) DIMENSION AF(1:8),F(1:8) DATA F(1)/ 3.0717001D-06/, F(2)/-5.2985762D-05/, - F(3)/-1.6372517D+01/, F(4)/ 3.6884682D-05/, - F(5)/ 2.5011231D+00/, F(6)/ 1.0127670D+00/, - F(7)/ 8.9057501D+00/, F(8)/ 4.3887271D+00/ DO 1 I=1,8 AF(I)=F(I) 1 CONTINUE RETURN END C======================================================================= SUBROUTINE S40C09(IHS,PI,Z,TH1,TH2,Z1,Z2,TH) *** IHS : IHS=1 FOR F64C09(FP,FH) *** IHS=2 FOR F65C09(FP,FS) AND *** F71CO9(FP,FS), F79C09(FP,FS), F80C09(FP,FS) *** TH : OUTPUT=T/TCR=1.0D0/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=TH1+DTH TAU2=1.0D0/TH2 CALL S5C09(PI,TAU2,OMEGA2) IF (OMEGA2.LT.-1.0D+10) GO TO 1000 IF (IHS.EQ.1) THEN CALL S10C09(OMEGA2,TAU2,Z2) ELSE IF (IHS.EQ.2) THEN CALL S8C09(OMEGA2,TAU2,Z2) END IF 1 CONTINUE 1000 TH=-1.0D+10 RETURN END SUBROUTINE S90C09(FP,FT,ILL) ***** SUBROUTINE TO CHECK IF THE RANGE OF ARGUMENTS:(P,T) IS PROPER***** ***** FOR F25C09:HPT, F35C09:SPT, F44C09:UPT, F51C09:VPT ****** ***** FOR F18C09:CPPT, F77C09:CVPT, F82C09:AKPT, F83C09:WPT ****** ***** FOR F84C09:AJTPT ****** ***** FOR F90C09:BSPT, F91C09:BTPT, F92C09:BPPT, F93C09:BVPT ****** IMPLICIT REAL(F) ILL=10000 IF ((FP.LT.2.3899E-03).OR.(FP.GT.200.01E+0)) RETURN FTMIN=F69C09(FP)-0.01 FTMAX=26.85+0.04 IF ((FT.GT.FTMIN).AND.(FT.LT.FTMAX)) THEN ILL=0 END IF RETURN END C================================================================= SUBROUTINE S97C09(FUN) *** LEVEL 1 ERROR MESSAGE *** CHARACTER FUN*6, MSG*125 INTEGER KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN MSG='**** NO CONVERGENCE AT '//FUN//' FOR FLUORINE ****' WRITE(6,1000) MSG 1000 FORMAT(1H ,5X,A) END IF RETURN END SUBROUTINE S98C09(IARG,ARG1,ARG2,NARG1,NARG2,NFUN) *** IARG=1 FOR ONE ARGUMENT (SECOND ARUMENT IS DUMMY) *** IARG=2 FOR TWO 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 FLUORINE', - ' WHEN ',A1,' =', 1PE14.7,' ****') 2020 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR FLUORINE', - ' WHEN ',A1,' =',1PE14.7,' AND ',A1,' =',1PE14.7,' ****') RETURN END SUBROUTINE S99C09(FUN) *** LEVEL 3 ERROR MESSAGE *** CHARACTER FUN*6, MSG*125 INTEGER KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN MSG='**** FUNCTION '//FUN//' UNAVAILABLE FOR FLUORINE ****' WRITE(6,3000) MSG 3000 FORMAT(1H ,5X,A) END IF RETURN END C----------------------------------------------------------------------- REAL FUNCTION T90(T) * Conversion of temperature scale 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