*================== PAIRV81 ==1990.9.28===============================* *================== PAAIR ====1997.3.10===============================* C==================================================================== C May 2, 2001: modified by Ryo Akasaka for ver.12.1 C The following 18 functions were added but not implemented. C AKPD(8A), AKPDD(8B), AKTD(8C), AKTDD(8D), CVPD(7A), CVTD(7B), C EPSPD(2A), EPSPDD(2B), EPSTD(2C), EPSTDD(2D), GAMPD(9A), C GAMTD(9B), TPH2(6H), TPS2(6S), WPD(8E), WPDD(8F), WTD(8G), C WTDD(8H) C C------------------------------------------------- F8A = AKPD REAL FUNCTION AKPD(P) CHARACTER FUN*6 COMMON/UNIT/KPA,MESS FUN='AKPD' IF (MESS.NE.0) CALL S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(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 S99AIR(FUN) WTDD=-1.0E+30 RETURN END C C==================================================================== *------------------------------------------------- F1AIR = AIPPT REAL FUNCTION AIPPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99AIR('AIPPT') AIPPT=-1.0E+30 RETURN END *------------------------------------------------- F2AIR = ALAPP REAL FUNCTION ALAPP(P) REAL P,PI PI=P CALL S99AIR('ALAPP') ALAPP=-1.0E+30 RETURN END *------------------------------------------------- F3AIR = ALAPT REAL FUNCTION ALAPT(T) REAL T,TI TI=T CALL S99AIR('ALAPT') ALAPT=-1.0E+30 RETURN END *------------------------------------------------- F4AIR = ALHP REAL FUNCTION ALHP(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) ALHP=F4AIR(PI) IF(ALHP.EQ.-1.0E+10) THEN CALL S97AIR('ALHP') ELSE IF(ALHP.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','ALHP') END IF RETURN END *------------------------------------------------- F5AIR = ALHT REAL FUNCTION ALHT(T) REAL T,TI TI=T CALL S99AIR('ALHT') ALHT=-1.0E+30 RETURN END *------------------------------------------------- F6AIR = ALMPD REAL FUNCTION ALMPD(P) REAL P,PI PI=P CALL S99AIR('ALMPD') ALMPD=-1.0E+30 RETURN END *------------------------------------------------- F7AIR = ALMPDD REAL FUNCTION ALMPDD(P) REAL P,PI PI=P CALL S99AIR('ALMPDD') ALMPDD=-1.0E+30 RETURN END *------------------------------------------------- F8AIR = ALMPT REAL FUNCTION ALMPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) ALMPT=F8AIR(PI,TI) IF(ALMPT.EQ.-1.0E+10) THEN CALL S97AIR('ALMPT') ELSE IF(ALMPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','ALMPT') END IF RETURN END *------------------------------------------------- F9AIR = ALMTD REAL FUNCTION ALMTD(T) REAL T,TI TI=T CALL S99AIR('ALMTD') ALMTD=-1.0E+30 RETURN END *------------------------------------------------- F10AIR = ALMTDD REAL FUNCTION ALMTDD(T) REAL T,TI TI=T CALL S99AIR('ALMTDD') ALMTDD=-1.0E+30 RETURN END *------------------------------------------------- F11AIR = AMUPD REAL FUNCTION AMUPD(P) REAL P,PI PI=P CALL S99AIR('AMUPD') AMUPD=-1.0E+30 RETURN END *------------------------------------------------- F12AIR = AMUPDD REAL FUNCTION AMUPDD(P) REAL P,PI PI=P CALL S99AIR('AMUPDD') AMUPDD=-1.0E+30 RETURN END *------------------------------------------------- F13AIR = AMUPT REAL FUNCTION AMUPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) AMUPT=F13AIR(PI,TI) IF(AMUPT.EQ.-1.0E+10) THEN CALL S97AIR('AMUPT') ELSE IF(AMUPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','AMUPT') END IF RETURN END *------------------------------------------------- F14AIR = AMUTD REAL FUNCTION AMUTD(T) REAL T,TI TI=T CALL S99AIR('AMUTD') AMUTD=-1.0E+30 RETURN END *------------------------------------------------- F15AIR = AMUTDD REAL FUNCTION AMUTDD(T) REAL T,TI TI=T CALL S99AIR('AMUTDD') AMUTDD=-1.0E+30 RETURN END *------------------------------------------------- F16AIR = CPPD REAL FUNCTION CPPD(P) REAL P,PI PI=P CALL S99AIR('CPPD') CPPD=-1.0E+30 RETURN END *------------------------------------------------- F17AIR = CPPDD REAL FUNCTION CPPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) CPPDD=F17AIR(PI) IF(CPPDD.EQ.-1.0E+10) THEN CALL S97AIR('CPPDD') ELSE IF(CPPDD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','CPPDD') END IF RETURN END *------------------------------------------------- F18AIR = CPPT REAL FUNCTION CPPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) CPPT=F18AIR(PI,TI) IF(CPPT.EQ.-1.0E+10) THEN CALL S97AIR('CPPT') ELSE IF(CPPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','CPPT') END IF RETURN END *------------------------------------------------- F19AIR = CPTD REAL FUNCTION CPTD(T) REAL T,TI TI=T CALL S99AIR('CPTD') CPTD=-1.0E+30 RETURN END *------------------------------------------------- F20AIR = CPTDD REAL FUNCTION CPTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) CPTDD=F20AIR(TI) IF(CPTDD.EQ.-1.0E+10) THEN CALL S97AIR('CPTDD') ELSE IF(CPTDD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','CPTDD') END IF RETURN END *------------------------------------------------- F21AIR = 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=F21AIR(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 AIR', - ' 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 *------------------------------------------------- F22AIR = EPSPT REAL FUNCTION EPSPT(P,T) REAL P,PI,T,TI PI=P TI=T CALL S99AIR('EPSPT') EPSPT=-1.0E+30 RETURN END *------------------------------------------------- F23AIR = HPD C************************************************ C FUNCTION FOR FUNDDAMENTAL CONSTANTS C PROPATH VER.7.1, MAY 8, 1990 C USAGE: B=FC(A) C A, B : CHARACTER TYPE VALIABLES C B='28.96' WHEN A='M' C B='287.22' WHEN A='R' C************************************************ REAL FUNCTION FC(A) CHARACTER A*1,MSG*120 COMMON/UNIT/KPA,MESS IF (A.EQ.'M') THEN FC=28.96 ELSE IF (A.EQ.'R') THEN FC=287.22 ELSE FC=-1.E+20 IF (MESS.NE.0) THEN MSG='**** OUT OF RANGE AT FC FOR AIR WHEN A=''' & //A//''' ****' WRITE(6,'(1H ,A)') MSG END IF END IF RETURN END C------------------------------------------------- F89 = FC REAL FUNCTION HPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) HPD=F23AIR(PI) IF(HPD.EQ.-1.0E+10) THEN CALL S97AIR('HPD') ELSE IF(HPD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','HPD') END IF RETURN END *------------------------------------------------- F24AIR = HPDD REAL FUNCTION HPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) HPDD=F24AIR(PI) IF(HPDD.EQ.-1.0E+10) THEN CALL S97AIR('HPDD') ELSE IF(HPDD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','HPDD') END IF RETURN END *------------------------------------------------- F25AIR = HPT REAL FUNCTION HPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) HPT=F25AIR(PI,TI) IF(HPT.EQ.-1.0E+10) THEN CALL S97AIR('HPT') ELSE IF(HPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','HPT') END IF RETURN END *------------------------------------------------- F26AIR = HPX REAL FUNCTION HPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) HPX=F26AIR(PI,X) IF(HPX.EQ.-1.0E+10) THEN CALL S97AIR('HPX') ELSE IF(HPX.EQ.-1.0E+20) THEN CALL S98AIR(2,P,X,'P','X','HPX') END IF RETURN END *------------------------------------------------- F27AIR = HTD REAL FUNCTION HTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) HTD=F27AIR(TI) IF(HTD.EQ.-1.0E+10) THEN CALL S97AIR('HTD') ELSE IF(HTD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','HTD') END IF RETURN END *------------------------------------------------- F28AIR = HTDD REAL FUNCTION HTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) HTDD=F28AIR(TI) IF(HTDD.EQ.-1.0E+10) THEN CALL S97AIR('HTDD') ELSE IF(HTDD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','HTDD') END IF RETURN END *------------------------------------------------- F29AIR = HTX REAL FUNCTION HTX(T,X) REAL T,TI,X,XI TI=T XI=X CALL S99AIR('HTX') HTX=-1.0E+30 RETURN END *------------------------------------------------- F84AIR = IDENTF C************************************************ C FUNCTION FOR IDENTIFICATION OF SUBSTANCE C PROPATH VER.12.1, MAY 2, 2001 C USAGE: B=IDENTF(A) C A, B : CHARACTER TYPE VALIABLES C B='AIR' WHEN A='S' C B='AIR' 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='AIR' ELSE IF (A.EQ.'C') THEN IDENTF='AIR' ELSE IF (A.EQ.'V') THEN IDENTF='12.1' ELSE IDENTF='????????????????????' IF (MESS.NE.0) THEN MSG='**** OUT OF RANGE AT IDENTF FOR AIR WHEN A=''' & //A//''' ****' WRITE(6,'(1H ,A)') MSG END IF END IF RETURN END *------------------------------------------------- F30AIR = PST REAL FUNCTION PST(T) REAL TI,T TI=T CALL S99AIR('PST') PST=-1.0E+30 RETURN END *------------------------------------------------- F31AIR = SIGP REAL FUNCTION SIGP(P) REAL P,PI PI=P CALL S99AIR('SIGP') SIGP=-1.0E+30 RETURN END *------------------------------------------------- F32AIR = SIGT REAL FUNCTION SIGT(T) REAL T,TI TI=T CALL S99AIR('SIGT') SIGT=-1.0E+30 RETURN END *------------------------------------------------- F33AIR = SPD REAL FUNCTION SPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) SPD=F33AIR(PI) IF(SPD.EQ.-1.0E+10) THEN CALL S97AIR('SPD') ELSE IF(SPD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','SPD') END IF RETURN END *------------------------------------------------- F34AIR = SPDD REAL FUNCTION SPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) SPDD=F34AIR(PI) IF(SPDD.EQ.-1.0E+10) THEN CALL S97AIR('SPDD') ELSE IF(SPDD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','SPDD') END IF RETURN END *------------------------------------------------- F35AIR = SPT REAL FUNCTION SPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) SPT=F35AIR(PI,TI) IF(SPT.EQ.-1.0E+10) THEN CALL S97AIR('SPT') ELSE IF(SPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','SPT') END IF RETURN END *------------------------------------------------- F36AIR = SPX REAL FUNCTION SPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) SPX=F36AIR(PI,X) IF(SPX.EQ.-1.0E+10) THEN CALL S97AIR('SPX') ELSE IF(SPX.EQ.-1.0E+20) THEN CALL S98AIR(2,P,X,'P','X','SPX') END IF RETURN END *------------------------------------------------- F37AIR = STD REAL FUNCTION STD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) STD=F37AIR(TI) IF(STD.EQ.-1.0E+10) THEN CALL S97AIR('STD') ELSE IF(STD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','STD') END IF RETURN END *------------------------------------------------- F38AIR = STDD REAL FUNCTION STDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) STDD=F38AIR(TI) IF(STDD.EQ.-1.0E+10) THEN CALL S97AIR('STDD') ELSE IF(STDD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','STDD') END IF RETURN END *------------------------------------------------- F39AIR = STX REAL FUNCTION STX(T,X) REAL T,TI,X,XI TI=T XI=X CALL S99AIR('STX') STX=-1.0E+30 RETURN END *------------------------------------------------- F40AIR = TSP REAL FUNCTION TSP(P) REAL P,PI PI=P CALL S99AIR('TSP') TSP=-1.0E+30 RETURN END *------------------------------------------------- F41AIR = TRPL REAL FUNCTION TRPL(A) CHARACTER A*1,AI*1 AI=A CALL S99AIR('TRPL') TRPL=-1.0E+30 RETURN END *------------------------------------------------- F42AIR = UPD REAL FUNCTION UPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) UPD=F42AIR(PI) IF(UPD.EQ.-1.0E+10) THEN CALL S97AIR('UPD') ELSE IF(UPD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','UPD') END IF RETURN END *------------------------------------------------- F43AIR = UPDD REAL FUNCTION UPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) UPDD=F43AIR(PI) IF(UPDD.EQ.-1.0E+10) THEN CALL S97AIR('UPDD') ELSE IF(UPDD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','UPDD') END IF RETURN END *------------------------------------------------- F44AIR = UPT REAL FUNCTION UPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) UPT=F44AIR(PI,TI) IF(UPT.EQ.-1.0E+10) THEN CALL S97AIR('UPT') ELSE IF(UPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','UPT') END IF RETURN END *------------------------------------------------- F45AIR = UPX REAL FUNCTION UPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) UPX=F45AIR(PI,X) IF(UPX.EQ.-1.0E+10) THEN CALL S97AIR('UPX') ELSE IF(UPX.EQ.-1.0E+20) THEN CALL S98AIR(2,P,X,'P','X','UPX') END IF RETURN END *------------------------------------------------- F46AIR = UTD REAL FUNCTION UTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) UTD=F46AIR(TI) IF(UTD.EQ.-1.0E+10) THEN CALL S97AIR('UTD') ELSE IF(UTD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','UTD') END IF RETURN END *------------------------------------------------- F47AIR = UTDD REAL FUNCTION UTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) UTDD=F47AIR(TI) IF(UTDD.EQ.-1.0E+10) THEN CALL S97AIR('UTDD') ELSE IF(UTDD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','UTDD') END IF RETURN END *------------------------------------------------- F48AIR = UTX REAL FUNCTION UTX(T,X) REAL T,TI,X,XI TI=T XI=X CALL S99AIR('UTX') UTX=-1.0E+30 RETURN END *------------------------------------------------- F49AIR = VPD REAL FUNCTION VPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) VPD=F49AIR(PI) IF(VPD.EQ.-1.0E+10) THEN CALL S97AIR('VPD') ELSE IF(VPD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','VPD') END IF RETURN END *------------------------------------------------- F50AIR = VPDD REAL FUNCTION VPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) VPDD=F50AIR(PI) IF(VPDD.EQ.-1.0E+10) THEN CALL S97AIR('VPDD') ELSE IF(VPDD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','VPDD') END IF RETURN END *------------------------------------------------- F51AIR = VPT REAL FUNCTION VPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) VPT=F51AIR(PI,TI) IF(VPT.EQ.-1.0E+10) THEN CALL S97AIR('VPT') ELSE IF(VPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','VPT') END IF RETURN END *------------------------------------------------- F52AIR = VPX REAL FUNCTION VPX(P,X) REAL P,PI,X INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) VPX=F52AIR(PI,X) IF(VPX.EQ.-1.0E+10) THEN CALL S97AIR('VPX') ELSE IF(VPX.EQ.-1.0E+20) THEN CALL S98AIR(2,P,X,'P','X','VPX') END IF RETURN END *------------------------------------------------- F53AIR = VTD REAL FUNCTION VTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) VTD=F53AIR(TI) IF(VTD.EQ.-1.0E+10) THEN CALL S97AIR('VTD') ELSE IF(VTD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','VTD') END IF RETURN END *------------------------------------------------- F54AIR = VTDD REAL FUNCTION VTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) VTDD=F54AIR(TI) IF(VTDD.EQ.-1.0E+10) THEN CALL S97AIR('VTDD') ELSE IF(VTDD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','VTDD') END IF RETURN END *------------------------------------------------- F55AIR = VTX REAL FUNCTION VTX(T,X) REAL T,TI,X,XI TI=T XI=X CALL S99AIR('VTX') VTX=-1.0E+30 RETURN END *------------------------------------------------- F56AIR = XPH REAL FUNCTION XPH(P,H) REAL P,PI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) XPH=F56AIR(PI,H) IF(XPH.EQ.-1.0E+10) THEN CALL S97AIR('XPH') ELSE IF(XPH.EQ.-1.0E+20) THEN CALL S98AIR(2,P,H,'P','H','XPH') END IF RETURN END *------------------------------------------------- F57AIR = XPS REAL FUNCTION XPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) XPS=F57AIR(PI,S) IF(XPS.EQ.-1.0E+10) THEN CALL S97AIR('XPS') ELSE IF(XPS.EQ.-1.0E+20) THEN CALL S98AIR(2,P,S,'P','S','XPS') END IF RETURN END *------------------------------------------------- F58AIR = XPU REAL FUNCTION XPU(P,U) REAL P,PI,U INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) XPU=F58AIR(PI,U) IF(XPU.EQ.-1.0E+10) THEN CALL S97AIR('XPU') ELSE IF(XPU.EQ.-1.0E+20) THEN CALL S98AIR(2,P,U,'P','U','XPU') END IF RETURN END *------------------------------------------------- F59AIR = XPV REAL FUNCTION XPV(P,V) REAL P,PI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) XPV=F59AIR(PI,V) IF(XPV.EQ.-1.0E+10) THEN CALL S97AIR('XPV') ELSE IF(XPV.EQ.-1.0E+20) THEN CALL S98AIR(2,P,V,'P','V','XPV') END IF RETURN END *------------------------------------------------- F60AIR = XTH REAL FUNCTION XTH(T,H) REAL T,TI,H,HI TI=T HI=H CALL S99AIR('XTH') XTH=-1.0E+30 RETURN END *------------------------------------------------- F61AIR = XTS REAL FUNCTION XTS(T,S) REAL T,TI,S,SI TI=T SI=S CALL S99AIR('XTS') XTS=-1.0E+30 RETURN END *------------------------------------------------- F62AIR = XTU REAL FUNCTION XTU(T,U) REAL T,TI,U,UI TI=T UI=U CALL S99AIR('XTU') XTU=-1.0E+30 RETURN END *------------------------------------------------- F63AIR = XTV REAL FUNCTION XTV(T,V) REAL T,TI,V,VI TI=T VI=V CALL S99AIR('XTV') XTV=-1.0E+30 RETURN END *------------------------------------------------- F64AIR = TPH REAL FUNCTION TPH(P,H) REAL P,PI,H INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TPH=F64AIR(PI,H) IF(TPH.EQ.-1.0E+10) THEN CALL S97AIR('TPH') RETURN ELSE IF(TPH.EQ.-1.0E+20) THEN CALL S98AIR(2,P,H,'P','H','TPH') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPH=TPH+273.15 RETURN END *------------------------------------------------- F65AIR = TPS REAL FUNCTION TPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TPS=F65AIR(PI,S) IF(TPS.EQ.-1.0E+10) THEN CALL S97AIR('TPS') RETURN ELSE IF(TPS.EQ.-1.0E+20) THEN CALL S98AIR(2,P,S,'P','S','TPS') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPS=TPS+273.15 RETURN END *------------------------------------------------- F66AIR = PLDT REAL FUNCTION PLDT(T) REAL T,TI TI=T CALL S99AIR('PLDT') PLDT=-1.0E+30 RETURN END *------------------------------------------------- F67AIR = TLDP REAL FUNCTION TLDP(P) REAL P,PI PI=P CALL S99AIR('TLDP') TLDP=-1.0E+30 RETURN END *------------------------------------------------- F68AIR = PMLT REAL FUNCTION PMLT(T) REAL T,TI TI=T CALL S99AIR('PMLT') PMLT=-1.0E+30 RETURN END *------------------------------------------------- F69AIR = TMLP REAL FUNCTION TMLP(P) REAL P,PI PI=P CALL S99AIR('TMLP') TMLP=-1.0E+30 RETURN END *------------------------------------------------- F70AIR = TPV REAL FUNCTION TPV(P,V) REAL P,PI,V INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TPV=F70AIR(PI,V) IF(TPV.EQ.-1.0E+10) THEN CALL S97AIR('TPV') RETURN ELSE IF(TPV.EQ.-1.0E+20) THEN CALL S98AIR(2,P,V,'P','V','TPV') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TPV=TPV+273.15 RETURN END *------------------------------------------------- F71AIR = HPS REAL FUNCTION HPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) HPS=F71AIR(PI,S) IF(HPS.EQ.-1.0E+10) THEN CALL S97AIR('HPS') ELSE IF(HPS.EQ.-1.0E+20) THEN CALL S98AIR(2,P,S,'P','S','HPS') END IF RETURN END *------------------------------------------------- F72AIR = PSTD REAL FUNCTION PSTD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) PSTD=F72AIR(TI) IF(PSTD.EQ.-1.0E+10) THEN CALL S97AIR('PSTD') RETURN ELSE IF(PSTD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','PSTD') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.2)) RETURN PSTD=PSTD*1.0E+05 RETURN END *------------------------------------------------- F73AIR = PSTDD REAL FUNCTION PSTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) PSTDD=F73AIR(TI) IF(PSTDD.EQ.-1.0E+10) THEN CALL S97AIR('PSTDD') RETURN ELSE IF(PSTDD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','PSTDD') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.2)) RETURN PSTDD=PSTDD*1.0E+05 RETURN END *------------------------------------------------- F74AIR = TSPD REAL FUNCTION TSPD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TSPD=F74AIR(PI) IF(TSPD.EQ.-1.0E+10) THEN CALL S97AIR('TSPD') RETURN ELSE IF(TSPD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','TSPD') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TSPD=TSPD+273.15 RETURN END *------------------------------------------------- F75AIR = TSPDD REAL FUNCTION TSPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TSPDD=F75AIR(PI) IF(TSPDD.EQ.-1.0E+10) THEN CALL S97AIR('SPDD') RETURN ELSE IF(TSPDD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','TSPDD') RETURN END IF IF((KPA.EQ.1).OR.(KPA.EQ.3)) RETURN TSPDD=TSPDD+273.15 RETURN END *------------------------------------------------- F76AIR = CVPDD REAL FUNCTION CVPDD(P) REAL P,PI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) CVPDD=F76AIR(PI) IF(CVPDD.EQ.-1.0E+10) THEN CALL S97AIR('CVPDD') ELSE IF(CVPDD.EQ.-1.0E+20) THEN CALL S98AIR(1,P,P,'P','P','CVPDD') END IF RETURN END *------------------------------------------------- F77AIR = CVPT REAL FUNCTION CVPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) CVPT=F77AIR(PI,TI) IF(CVPT.EQ.-1.0E+10) THEN CALL S97AIR('CVPT') ELSE IF(CVPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','CVPT') END IF RETURN END *------------------------------------------------- F78AIR = CVTDD REAL FUNCTION CVTDD(T) REAL T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS TI=G99AIR(KPA,T) CVTDD=F78AIR(TI) IF(CVTDD.EQ.-1.0E+10) THEN CALL S97AIR('CVTDD') ELSE IF(CVTDD.EQ.-1.0E+20) THEN CALL S98AIR(1,T,T,'T','T','CVTDD') END IF RETURN END *------------------------------------------------- F79AIR = UPS REAL FUNCTION UPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) UPS=F79AIR(PI,S) IF(UPS.EQ.-1.0E+10) THEN CALL S97AIR('UPS') ELSE IF(UPS.EQ.-1.0E+20) THEN CALL S98AIR(2,P,S,'P','S','UPS') END IF RETURN END *------------------------------------------------- F80AIR = VPS REAL FUNCTION VPS(P,S) REAL P,PI,S INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) VPS=F80AIR(PI,S) IF(VPS.EQ.-1.0E+10) THEN CALL S97AIR('VPS') ELSE IF(VPS.EQ.-1.0E+20) THEN CALL S98AIR(2,P,S,'P','S','VPS') END IF RETURN END *------------------------------------------------- F81AIR = PRPT REAL FUNCTION PRPT(P,T) REAL P,PI,T,TI INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) PRPT=F81AIR(PI,TI) IF(PRPT.EQ.-1.0E+10) THEN CALL S97AIR('PRPT') ELSE IF(PRPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','PRPT') END IF RETURN END *------------------------------------------------- F82AIR = AKPT REAL FUNCTION AKPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) AKPT=F82AIR(PI,TI) IF(AKPT.EQ.-1.0E+10) THEN CALL S97AIR('AKPT') ELSE IF(AKPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','AKPT') END IF RETURN END *------------------------------------------------- F83AIR = WPT REAL FUNCTION WPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) WPT=F83AIR(PI,TI) IF(WPT.EQ.-1.0E+10) THEN CALL S97AIR('WPT') ELSE IF(WPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','WPT') END IF RETURN END *------------------------------------------------- F85AIR = PRPD REAL FUNCTION PRPD(P) REAL P,PI PI=P CALL S99AIR('PRPD') PRPD=-1.0E+30 RETURN END *------------------------------------------------- F86AIR = PRPDD REAL FUNCTION PRPDD(P) REAL P,PI PI=P CALL S99AIR('PRPDD') PRPDD=-1.0E+30 RETURN END *------------------------------------------------- F87AIR = PRTD REAL FUNCTION PRTD(T) REAL T,TI TI=T CALL S99AIR('PRTD') PRTD=-1.0E+30 RETURN END *------------------------------------------------- F88AIR = PRTDD REAL FUNCTION PRTDD(T) REAL T,TI TI=T CALL S99AIR('PRTDD') PRTDD=-1.0E+30 RETURN END *------------------------------------------------- F90AIR = BSPT REAL FUNCTION BSPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) BSPT=F90AIR(PI,TI) IF(BSPT.EQ.-1.0E+10) THEN CALL S97AIR('BSPT') RETURN ELSE IF(BSPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','BSPT') RETURN END IF BSPT=BSPT RETURN END *------------------------------------------------- F91AIR = BTPT REAL FUNCTION BTPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) BTPT=F91AIR(PI,TI) IF(BTPT.EQ.-1.0E+10) THEN CALL S97AIR('BTPT') RETURN ELSE IF(BTPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','BTPT') RETURN END IF BTPT=BTPT RETURN END *------------------------------------------------- F92AIR = BPPT REAL FUNCTION BPPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) BPPT=F92AIR(PI,TI) IF(BPPT.EQ.-1.0E+10) THEN CALL S97AIR('BPPT') ELSE IF(BPPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','BPPT') END IF RETURN END *------------------------------------------------- F93AIR = BVPT REAL FUNCTION BVPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) BVPT=F93AIR(PI,TI) IF(BVPT.EQ.-1.0E+10) THEN CALL S97AIR('BVPT') ELSE IF(BVPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','BVPT') END IF RETURN END *------------------------------------------------- F94AIR = AJTPT REAL FUNCTION AJTPT(P,T) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS PI=G98AIR(KPA,P) TI=G99AIR(KPA,T) AJTPT=F94AIR(PI,TI) IF(AJTPT.EQ.-1.0E+10) THEN CALL S97AIR('AJTPT') RETURN ELSE IF(AJTPT.EQ.-1.0E+20) THEN CALL S98AIR(2,P,T,'P','T','AJTPT') RETURN END IF AJTPT=AJTPT RETURN END *------------------------------------------------- F95AIR = GAMPT REAL FUNCTION GAMPT(P,T) REAL T,TI,P,PI TI=T PI=P CALL S99AIR('GAMPT') GAMPT=-1.0E+30 RETURN END *------------------------------------------------- F96AIR = GAMPDD REAL FUNCTION GAMPDD(P) REAL P,PI PI=P CALL S99AIR('GAMPDD') GAMPDD=-1.0E+30 RETURN END *------------------------------------------------- F97AIR = GAMTDD REAL FUNCTION GAMTDD(T) REAL T,TI TI=T CALL S99AIR('GAMTDD') GAMTDD=-1.0E+30 RETURN END *------------------------------------------------- F98AIR = TPSEUP REAL FUNCTION TPSEUP(P) REAL P,PI PI=P CALL S99AIR('TPSEUP') TPSEUP=-1.0E+30 RETURN END *------------------------------------------------- F99AIR = PSBT REAL FUNCTION PSBT(T) REAL T,TI TI=T CALL S99AIR('PSBT') PSBT=-1.0E+30 RETURN END *------------------------------------------------- F100AIR = TSBP REAL FUNCTION TSBP(P) REAL P,PI PI=P CALL S99AIR('TSBP') TSBP=-1.0E+30 RETURN END *------------------------------------------------- G98AIR REAL FUNCTION G98AIR(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 G98AIR=P*PBAR RETURN END *------------------------------------------------- G99AIR REAL FUNCTION G99AIR(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 G99AIR=T-T0K RETURN END ***** FUNCTIONS : FP IN (BAR) AND FCT IN (C) *************************** REAL FUNCTION F4AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F4AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01)) RETURN F4AIR=-1.0E+10 PI=DBLE(FP)/PK CALL S28AIR(1,PI,0.0D0,HPDD) CALL S28AIR(1,PI,1.0D0,HPD) IF((HPDD.LE.-1.0D+10).OR.(HPD.LE.-1.0D+10)) RETURN ALHP=HPDD-HPD F4AIR=REAL(ALHP) RETURN END REAL FUNCTION F8AIR(FP,FCT) **** BY NAGASHIMA, ET AL. (EQ.(5)) **** IMPLICIT DOUBLE PRECISION(A-E,G-H,L,O-Z) PARAMETER(TSTAR=132.5D0,RSTAR=314.3D0,LAMBDA=25.9778D-03) DIMENSION C(-4:1),D(1:5) DATA C(1)/ 0.239503D+00/, C05/ 0.649768D-02/,C(0)/ 1.0D+00/, - C(-1)/-0.192615D+01/,C(-2)/ 0.200383D+01/, - C(-3)/-0.107553D+01/,C(-4)/ 0.229414D+00/ DATA D(1)/ 0.402287D+00/, D(2)/ 0.356603D+00/,D(3)/-0.163159D+00/, - D(4)/ 0.138059D+00/, D(5)/-0.201725D-01/ F8AIR=-1.0E+20 IF((FP.GT.1000.01).OR. - (FCT.LT.-188.16).OR.(FCT.GT.1726.86)) RETURN FT=FCT+273.15 IF((FCT.GE.-114.16).AND.(FCT.LT.-3.36)) THEN FPMIN=0.43E-01*FT**2-0.103E+02*FT+0.649E+03 IF(FP.GT.FPMIN) RETURN ELSE IF((FCT.GE.-143.33).AND.(FCT.LT.-114.16)) THEN IF(FP.GT.33.01) RETURN ELSE IF((FCT.GE.-188.16).AND.(FCT.LT.-143.33)) THEN **** DEW-POINT CURVE BY VASSERMAN, ET AL. (EQ.(6)) **** FPMIN=(10.0**(9.838-527.88/FT-8.4974E-02*FT - +4.5759E-04*FT**2-0.88391E-06*FT**3))*10.0 IF(FP.GT.FPMIN) RETURN END IF F8AIR=-1.0E+10 P=DBLE(FP)*1.0D+05 T=DBLE(FCT)+273.15D0 IF(FP.EQ.0.0) THEN R=0.0D0 ELSE IF(FCT.LT.0.0) THEN CALL S50AIR(P,T,R) ELSE IF((FCT.GE.0.0).AND.(FCT.LT.1000.0)) THEN CALL S51AIR(P,T,R) ELSE IF(FCT.GE.1000.0) THEN CALL S52AIR(P,T,R) END IF IF(R.LT.0.0D0) RETURN END IF RR=R/RSTAR TR=T/TSTAR CT=C(0)+(C(-1)+(C(-2)+(C(-3)+C(-4)/TR)/TR)/TR)/TR DR=(D(1)+(D(2)+(D(3)+(D(4)+D(5)*RR)*RR)*RR)*RR)*RR ALMPT=LAMBDA*(C(1)*TR+C05*SQRT(TR)+CT+DR) F8AIR=REAL(ALMPT) RETURN END REAL FUNCTION F13AIR(FP,FCT) **** BY NAGASHIMA, ET AL. (EQ.(3)) **** IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TSTAR=132.5D0,RSTAR=314.3D0,H=6.16090D-06) DIMENSION A(-4:1),B(1:4) DATA A(1)/ 0.128517D+00/, A05/ 0.260661D+01/,A(0)/-1.0D+00/, - A(-1)/-0.709661D+00/,A(-2)/ 0.662534D+00/, - A(-3)/-0.197846D+00/,A(-4)/ 0.770147D-02/ DATA B(1)/ 0.465601D+00/, B(2)/ 0.126469D+01/,B(3)/-0.511425D+00/, - B(4)/ 0.274600D+00/ F13AIR=-1.0E+20 IF((FP.GT.1000.01).OR. - (FCT.LT.-188.16).OR.(FCT.GT.1726.86))RETURN FT=FCT+273.15 IF((FCT.GE.-140.65).AND.(FCT.LT.-3.36)) THEN FPMIN=0.43E-01*FT**2-0.103E+02*FT+0.649E+03 IF(FP.GT.FPMIN) RETURN ELSE IF((FCT.GE.-188.16).AND.(FCT.LT.-140.65)) THEN **** DEW-POINT CURVE BY VASSERMAN, ET AL. (EQ.(6)) **** FPMIN=(10.0**(9.838-527.88/FT-8.4974E-02*FT - +4.5759E-04*FT**2-0.88391E-06*FT**3))*10.0 IF(FP.GT.FPMIN) RETURN END IF F13AIR=-1.0E+10 P=DBLE(FP)*1.0D+05 T=DBLE(FCT)+273.15D0 IF(FP.EQ.0.0) THEN R=0.0D0 ELSE IF(FCT.LT.0.0) THEN CALL S50AIR(P,T,R) ELSE IF((FCT.GE.0.0).AND.(FCT.LT.1000.0)) THEN CALL S51AIR(P,T,R) ELSE IF(FCT.GE.1000.0) THEN CALL S52AIR(P,T,R) END IF IF(R.LT.0.0D0) RETURN END IF RR=R/RSTAR TR=T/TSTAR AT=A(0)+(A(-1)+(A(-2)+(A(-3)+A(-4)/TR)/TR)/TR)/TR BR=(B(1)+(B(2)+(B(3)+B(4)*RR)*RR)*RR)*RR AMUPT=H*(A(1)*TR+A05*SQRT(TR)+AT+BR) F13AIR=REAL(AMUPT) RETURN END REAL FUNCTION F17AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F17AIR=-1.0E+20 IF((FP.LT.0.01).OR.(FP.GT.37.6626)) RETURN F17AIR=-1.0E+10 PI=DBLE(FP)/PK IF(ABS(PI-1.0D0).LT.1.0D-05) THEN CALL S35AIR(2,1.0D0,1.0D0,CPPDD) ELSE CALL S22AIR(PI,THDD) CALL S31AIR(PI,THDD,DELTDD) CALL S35AIR(2,THDD,DELTDD,CPPDD) END IF IF(CPPDD.LT.0.0D0) RETURN F17AIR=REAL(CPPDD) RETURN END REAL FUNCTION F18AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F18AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK IF(FP.EQ.0.0) THEN IF(THETA.LT.1.6838968D0) THEN CALL S35AIR(2,THETA,0.0D0,CPPT0) ELSE CALL S15AIR(2,THETA,0.0D0,CPPT0) END IF F18AIR=REAL(CPPT0) RETURN END IF CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F18AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S31AIR(PI,THETA,DELTA) CALL S35AIR(2,THETA,DELTA,CPPT) ELSE CALL S11AIR(PI,THETA,DELTA) CALL S15AIR(2,THETA,DELTA,CPPT) END IF IF(CPPT.LT.0.0D0) RETURN F18AIR=REAL(CPPT) RETURN END REAL FUNCTION F20AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F20AIR=-1.0E+20 IF((FCT.LT.-213.16).OR.(FCT.GT.-140.62)) RETURN F20AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN CALL S35AIR(2,1.0D0,1.0D0,CPTDD) ELSE CALL S21AIR(THETA,PIDD) CALL S31AIR(PIDD,THETA,DELTDD) CALL S35AIR(2,THETA,DELTDD,CPTDD) END IF IF(CPTDD.LT.0.0D0) RETURN F20AIR=REAL(CPTDD) RETURN END REAL FUNCTION F21AIR(A) CHARACTER*1 A,B(1:5) DOUBLE PRECISION CRP(1:5),HK,SK DATA B(1)/'H'/,B(2)/'P'/,B(3)/'S'/,B(4)/'T'/,B(5)/'V'/ CALL S34AIR(1,1.0D0,1.0D0,HK) CALL S34AIR(2,1.0D0,1.0D0,SK) CRP(1)=HK CRP(2)=37.6625D0 CRP(3)=SK CRP(4)=132.52D0-273.15D0 CRP(5)=3.19489D-03 DO 10 I=1,5 F21AIR=REAL(CRP(I)) IF(A.EQ.B(I)) RETURN 10 CONTINUE F21AIR=-1.0E+20 RETURN END REAL FUNCTION F23AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F23AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01)) RETURN F23AIR=-1.0E+10 PI=DBLE(FP)/PK CALL S28AIR(1,PI,1.0D0,HPD) IF(HPD.LE.-1.0D+10) RETURN F23AIR=REAL(HPD) RETURN END REAL FUNCTION F24AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F24AIR=-1.0E+20 IF((FP.LT.0.01).OR.(FP.GT.37.6626)) RETURN F24AIR=-1.0E+10 PI=DBLE(FP)/PK IF(ABS(PI-1.0D0).LT.1.0D-05) THEN CALL S34AIR(1,1.0D0,1.0D0,HPDD) ELSE CALL S28AIR(1,PI,0.0D0,HPDD) END IF IF(HPDD.LE.-1.0D+10) RETURN F24AIR=REAL(HPDD) RETURN END REAL FUNCTION F25AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F25AIR=-1.0E+20 CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F25AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S31AIR(PI,THETA,DELTA) CALL S34AIR(1,THETA,DELTA,HPT) ELSE CALL S11AIR(PI,THETA,DELTA) CALL S14AIR(1,THETA,DELTA,HPT) END IF IF(HPT.LE.-1.0D+10) RETURN F25AIR=REAL(HPT) RETURN END REAL FUNCTION F26AIR(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F26AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F26AIR=-1.0E+10 PI=DBLE(FP)/PK X=DBLE(FX) CALL S28AIR(1,PI,1.0D0-X,HPX) IF(HPX.LE.-1.0D+10) RETURN F26AIR=REAL(HPX) RETURN END REAL FUNCTION F27AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F27AIR=-1.0E+20 IF((FCT.LT.-205.16).OR.(FCT.GT.-141.14)) RETURN F27AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK CALL S24AIR(THETA,PID) CALL S28AIR(1,PID,1.0D0,HTD) IF(HTD.LE.-1.0D+10) RETURN F27AIR=REAL(HTD) RETURN END REAL FUNCTION F28AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F28AIR=-1.0E+20 IF((FCT.LT.-213.16).OR.(FCT.GT.-140.62)) RETURN F28AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN CALL S34AIR(1,1.0D0,1.0D0,HTDD) ELSE CALL S21AIR(THETA,PIDD) CALL S31AIR(PIDD,THETA,DELTDD) CALL S34AIR(1,THETA,DELTDD,HTDD) END IF IF(HTDD.LE.-1.0D+10) RETURN F28AIR=REAL(HTDD) RETURN END REAL FUNCTION F33AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F33AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01)) RETURN F33AIR=-1.0E+10 PI=DBLE(FP)/PK CALL S28AIR(2,PI,1.0D0,SPD) IF(SPD.LT.0.0D0) RETURN F33AIR=REAL(SPD) RETURN END REAL FUNCTION F34AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F34AIR=-1.0E+20 IF((FP.LT.0.01).OR.(FP.GT.37.6626)) RETURN F34AIR=-1.0E+10 PI=DBLE(FP)/PK IF(ABS(PI-1.0D0).LT.1.0D-05) THEN CALL S34AIR(2,1.0D0,1.0D0,SPDD) ELSE CALL S28AIR(2,PI,0.0D0,SPDD) END IF IF(SPDD.LT.0.0D0) RETURN F34AIR=REAL(SPDD) RETURN END REAL FUNCTION F35AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F35AIR=-1.0E+20 CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F35AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S31AIR(PI,THETA,DELTA) CALL S34AIR(2,THETA,DELTA,SPT) ELSE CALL S11AIR(PI,THETA,DELTA) CALL S14AIR(2,THETA,DELTA,SPT) END IF IF(SPT.LT.0.0D0) RETURN F35AIR=REAL(SPT) RETURN END REAL FUNCTION F36AIR(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F36AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F36AIR=-1.0E+10 PI=DBLE(FP)/PK X=DBLE(FX) CALL S28AIR(2,PI,1.0D0-X,SPX) IF(SPX.LT.0.0D0) RETURN F36AIR=REAL(SPX) RETURN END REAL FUNCTION F37AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F37AIR=-1.0E+20 IF((FCT.LT.-205.16).OR.(FCT.GT.-141.14)) RETURN F37AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK CALL S24AIR(THETA,PID) CALL S28AIR(2,PID,1.0D0,STD) IF(STD.LT.0.0D0) RETURN F37AIR=REAL(STD) RETURN END REAL FUNCTION F38AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F38AIR=-1.0E+20 IF((FCT.LT.-213.16).OR.(FCT.GT.-140.62)) RETURN F38AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN CALL S34AIR(2,1.0D0,1.0D0,STDD) ELSE CALL S21AIR(THETA,PIDD) CALL S31AIR(PIDD,THETA,DELTDD) CALL S34AIR(2,THETA,DELTDD,STDD) END IF IF(STDD.LT.0.0D0) RETURN F38AIR=REAL(STDD) RETURN END REAL FUNCTION F42AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F42AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01)) RETURN F42AIR=-1.0E+10 PI=DBLE(FP)/PK CALL S28AIR(3,PI,1.0D0,UPD) IF(UPD.LE.-1.0D+10) RETURN F42AIR=REAL(UPD) RETURN END REAL FUNCTION F43AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F43AIR=-1.0E+20 IF((FP.LT.0.01).OR.(FP.GT.37.6626)) RETURN F43AIR=-1.0E+10 PI=DBLE(FP)/PK IF(ABS(PI-1.0D0).LT.1.0D-05) THEN CALL S34AIR(3,1.0D0,1.0D0,UPDD) ELSE CALL S28AIR(3,PI,0.0D0,UPDD) END IF IF(UPDD.LE.-1.0D+10) RETURN F43AIR=REAL(UPDD) RETURN END REAL FUNCTION F44AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F44AIR=-1.0E+20 CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F44AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S31AIR(PI,THETA,DELTA) CALL S34AIR(3,THETA,DELTA,UPT) ELSE CALL S11AIR(PI,THETA,DELTA) CALL S14AIR(3,THETA,DELTA,UPT) END IF IF(UPT.LE.-1.0D+10) RETURN F44AIR=REAL(UPT) RETURN END REAL FUNCTION F45AIR(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F45AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F45AIR=-1.0E+10 PI=DBLE(FP)/PK X=DBLE(FX) CALL S28AIR(3,PI,1.0D0-X,UPX) IF(UPX.LE.-1.0D+10) RETURN F45AIR=REAL(UPX) RETURN END REAL FUNCTION F46AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F46AIR=-1.0E+20 IF((FCT.LT.-205.16).OR.(FCT.GT.-141.14)) RETURN F46AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK CALL S24AIR(THETA,PID) CALL S28AIR(3,PID,1.0D0,UTD) IF(UTD.LE.-1.0D+10) RETURN F46AIR=REAL(UTD) RETURN END REAL FUNCTION F47AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F47AIR=-1.0E+20 IF((FCT.LT.-213.16).OR.(FCT.GT.-140.62)) RETURN F47AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN CALL S34AIR(3,1.0D0,1.0D0,UTDD) ELSE CALL S21AIR(THETA,PIDD) CALL S31AIR(PIDD,THETA,DELTDD) CALL S34AIR(3,THETA,DELTDD,UTDD) END IF IF(UTDD.LE.-1.0D+10) RETURN F47AIR=REAL(UTDD) RETURN END REAL FUNCTION F49AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,VK=3.19489D-03) F49AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01)) RETURN F49AIR=-1.0E+10 PI=DBLE(FP)/PK CALL S25AIR(PI,THD) CALL S27AIR(THD,DELTD) IF(DELTD.LE.0.0D0) RETURN F49AIR=REAL(VK/DELTD) RETURN END REAL FUNCTION F50AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,VK=3.19489D-03) F50AIR=-1.0E+20 IF((FP.LT.0.01).OR.(FP.GT.37.6626)) RETURN F50AIR=-1.0E+10 PI=DBLE(FP)/PK IF(ABS(PI-1.0D0).LT.1.0D-05) THEN DELTDD=1.0D0 ELSE CALL S22AIR(PI,THDD) CALL S31AIR(PI,THDD,DELTDD) IF(DELTDD.LE.0.0D0) RETURN END IF F50AIR=REAL(VK/DELTDD) RETURN END REAL FUNCTION F51AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(VK=3.19489D-03) F51AIR=-1.0E+20 CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F51AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S31AIR(PI,THETA,DELTA) ELSE CALL S11AIR(PI,THETA,DELTA) END IF IF(DELTA.LE.0.0D0) RETURN F51AIR=REAL(VK/DELTA) RETURN END REAL FUNCTION F52AIR(FP,FX) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F52AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FX.LT.0.0).OR.(FX.GT.1.0)) RETURN F52AIR=-1.0E+10 PI=DBLE(FP)/PK X=DBLE(FX) CALL S28AIR(4,PI,1.0D0-X,VPX) IF(VPX.LT.0.0D0) RETURN F52AIR=REAL(VPX) RETURN END REAL FUNCTION F53AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,VK=3.19489D-03) F53AIR=-1.0E+20 IF((FCT.LT.-205.16).OR.(FCT.GT.-140.62)) RETURN F53AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN DELTD=1.0D0 ELSE CALL S27AIR(THETA,DELTD) IF(DELTD.LE.0.0D0) RETURN END IF F53AIR=REAL(VK/DELTD) RETURN END REAL FUNCTION F54AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,VK=3.19489D-03) F54AIR=-1.0E+20 IF((FCT.LT.-213.16).OR.(FCT.GT.-140.62)) RETURN F54AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN DELTDD=1.0D0 ELSE CALL S21AIR(THETA,PIDD) CALL S31AIR(PIDD,THETA,DELTDD) IF(DELTDD.LE.0.0D0) RETURN END IF F54AIR=REAL(VK/DELTDD) RETURN END REAL FUNCTION F56AIR(FP,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,L,O-Z) PARAMETER(PK=37.6625D0) DATA FMIN/0.999999/,FMAX/1.00001/ F56AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FH.LT.F23AIR(FP)*FMIN).OR.(FH.GT.F24AIR(FP)*FMAX)) RETURN F56AIR=-1.0E+10 PI=DBLE(FP)/PK H=DBLE(FH) CALL S29AIR(1,PI,H,LAMBDA) IF((LAMBDA.LT.0.0D0).OR.(LAMBDA.GT.1.0D0)) RETURN XPH=1.0D0-LAMBDA F56AIR=REAL(XPH) RETURN END REAL FUNCTION F57AIR(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,L,O-Z) PARAMETER(PK=37.6625D0) DATA FMIN/0.999999/,FMAX/1.00001/ F57AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FS.LT.F33AIR(FP)*FMIN).OR.(FS.GT.F34AIR(FP)*FMAX)) RETURN F57AIR=-1.0E+10 PI=DBLE(FP)/PK S=DBLE(FS) CALL S29AIR(2,PI,S,LAMBDA) IF((LAMBDA.LT.0.0D0).OR.(LAMBDA.GT.1.0D0)) RETURN XPS=1.0D0-LAMBDA F57AIR=REAL(XPS) RETURN END REAL FUNCTION F58AIR(FP,FU) IMPLICIT DOUBLE PRECISION(A-E,G-H,L,O-Z) PARAMETER(PK=37.6625D0) DATA FMIN/0.999999/,FMAX/1.00001/ F58AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FU.LT.F42AIR(FP)*FMIN).OR.(FU.GT.F43AIR(FP)*FMAX)) RETURN F58AIR=-1.0E+10 PI=DBLE(FP)/PK U=DBLE(FU) CALL S29AIR(3,PI,U,LAMBDA) IF((LAMBDA.LT.0.0D0).OR.(LAMBDA.GT.1.0D0)) RETURN XPU=1.0D0-LAMBDA F58AIR=REAL(XPU) RETURN END REAL FUNCTION F59AIR(FP,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,L,O-Z) PARAMETER(PK=37.6625D0) DATA FMIN/0.999999/,FMAX/1.00001/ F59AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.36.01).OR. - (FV.LT.F49AIR(FP)*FMIN).OR.(FV.GT.F50AIR(FP)*FMAX)) RETURN F59AIR=-1.0E+10 PI=DBLE(FP)/PK V=DBLE(FV) CALL S29AIR(4,PI,V,LAMBDA) IF((LAMBDA.LT.0.0D0).OR.(LAMBDA.GT.1.0D0)) RETURN XPV=1.0D0-LAMBDA F59AIR=REAL(XPV) RETURN END REAL FUNCTION F64AIR(FP,FH) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F64AIR=-1.0E+20 CALL S91AIR(1,64,FP,FH,T,ILL91) IF(ILL91.EQ.10000) RETURN F64AIR=-1.0E+10 IF(ILL91.EQ.1000) RETURN F64AIR=REAL(T-273.15D0) RETURN END REAL FUNCTION F65AIR(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F65AIR=-1.0E+20 CALL S91AIR(2,65,FP,FS,T,ILL91) IF(ILL91.EQ.10000) RETURN F65AIR=-1.0E+10 IF(ILL91.EQ.1000) RETURN F65AIR=REAL(T-273.15D0) RETURN END REAL FUNCTION F70AIR(FP,FV) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0,VK=3.19489D-03) F70AIR=-1.0E+20 IF((FP.LT.0.0099).OR.(FP.GT.4500.02)) RETURN PI=DBLE(FP)/PK VR=DBLE(FV)/VK IF((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(VR-1.0D0).LT.1.0D-05)) - THEN F70AIR=REAL(TK-273.15D0) RETURN END IF THMAX=(1250.2D0+273.15D0)/TK CALL S11AIR(PI,THMAX,DELMIN) IF(DELMIN.LE.0.0D0) RETURN VRMAX=1.0D0/DELMIN IF(FP.LT.1691.0) THEN TH1=(-50.01D0+273.15D0)/TK CALL S11AIR(PI,TH1,DEL1) IF(DEL1.LE.0.0D0) RETURN VR1=1.0D0/DEL1 END IF IF(FP.GE.4000.0) THEN THMIN1=(-2.0D-04*(PI*PK)**2+1.98D0*(PI*PK)-4284.0D0+273.14D0) - /TK CALL S11AIR(PI,THMIN1,DEL1) IF(DEL1.LE.0.0D0) RETURN VRMIN1=1.0D0/DEL1 VRMIN2=VRMIN1 VRMIN3=VRMIN1 ELSE IF((FP.GT.3500.0).AND.(FP.LT.4000.0)) THEN THMIN1=(-6.65D-04*(PI*PK)**2+5.45D0*(PI*PK)-10724.0D0+273.14D0) - /TK CALL S11AIR(PI,THMIN1,DEL1) IF(DEL1.LE.0.0D0) RETURN VRMIN1=1.0D0/DEL1 VRMIN2=VRMIN1 VRMIN3=VRMIN1 ELSE IF((FP.GE.1691.0).AND.(FP.LE.3500.0)) THEN VRMIN1=0.00120398D0/VK VRMIN2=VRMIN1 VRMIN3=VRMIN1 ELSE IF((FP.GE.927.2).AND.(FP.LT.1691.0)) THEN VRMIN1=VR1 VRMIN2=VRMIN1 VRMIN3=VRMIN1 ELSE IF((FP.GT.206.1).AND.(FP.LT.927.21)) THEN VRMIN1=VR1 VRMIN2=0.0014198D0/VK VRMIN3=VRMIN2 ELSE IF((FP.GE.37.6625).AND.(FP.LE.206.1)) THEN VRMIN1=VR1 CALL S31AIR(PI,1.0D0,DELK) IF(DELK.LE.0.0D0) RETURN VRMIN2=1.0D0/DELK VRMIN3=VRMIN2 ELSE IF((FP.GT.36.0).AND.(FP.LT.37.6625)) THEN VRMIN1=VR1 CALL S22AIR(PI,THDD) CALL S31AIR(PI,THDD,DELTDD) IF(DELTDD.LE.0.0D0) RETURN VRMIN2=1.0D0/DELTDD VRMIN3=VRMIN2 ELSE IF((FP.GE.0.4).AND.(FP.LE.36.0)) THEN VRMIN1=VR1 CALL S22AIR(PI,THDD) CALL S31AIR(PI,THDD,DELTDD) IF(DELTDD.LE.0.0D0) RETURN VRMIN2=1.0D0/DELTDD CALL S25AIR(PI,THD) CALL S27AIR(THD,DELTD) IF(DELTD.LE.0.0D0) RETURN VRMIN3=1.0D0/DELTD ELSE IF(FP.LT.0.4) THEN VRMIN1=VR1 CALL S22AIR(PI,THDD) CALL S31AIR(PI,THDD,DELTDD) IF(DELTDD.LE.0.0D0) RETURN VRMIN2=1.0D0/DELTDD VRMIN3=VRMIN2 END IF F70AIR=-1.0E+10 IF((VR.GE.VRMIN1).AND.(VR.LE.VRMAX)) THEN CALL S12AIR(PI,1.0D0/VR,THETA) IF(THETA.LE.0.0D0) RETURN ELSE IF((VR.GE.VRMIN2).AND.(VR.LT.VRMIN1)) THEN CALL S32AIR(PI,1.0D0/VR,THETA) IF(THETA.LE.0.0D0) RETURN ELSE IF((VR.GE.VRMIN3).AND.(VR.LT.VRMIN2)) THEN CALL S22AIR(PI,THDD) CALL S31AIR(PI,THDD,DELTDD) CALL S25AIR(PI,THD) CALL S27AIR(THD,DELTD) IF((DELTDD.LE.0.0D0).OR.(DELTD.LE.0.0D0)) RETURN THETA=THDD-((1.0D0/DELTDD-VR)/(1.0D0/DELTDD-1.0D0/DELTD))* - (THDD-THD) ELSE F70AIR=-1.0E+20 RETURN END IF F70AIR=REAL(THETA*TK-273.15D0) RETURN END REAL FUNCTION F71AIR(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F71AIR=-1.0E+20 CALL S91AIR(2,71,FP,FS,H,ILL91) IF(ILL91.EQ.10000) RETURN F71AIR=-1.0E+10 IF(ILL91.EQ.1000) RETURN F71AIR=REAL(H) RETURN END REAL FUNCTION F72AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0) F72AIR=-1.0E+20 IF((FCT.LT.-205.16).OR.(FCT.GT.-140.72)) RETURN F72AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK CALL S24AIR(THETA,PID) IF(PID.LT.0.0D0) RETURN PSTD=PID*PK F72AIR=REAL(PSTD) RETURN END REAL FUNCTION F73AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0) F73AIR=-1.0E+20 IF((FCT.LT.-213.16).OR.(FCT.GT.-140.62)) RETURN THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN PIDD=1.0D0 ELSE CALL S21AIR(THETA,PIDD) IF(PIDD.LT.0.0D0) RETURN END IF PSTDD=PIDD*PK F73AIR=REAL(PSTDD) RETURN END REAL FUNCTION F74AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0) F74AIR=-1.0E+20 IF((FP.LT.0.4).OR.(FP.GT.37.7435)) RETURN F74AIR=-1.0E+10 PI=DBLE(FP)/PK CALL S25AIR(PI,THD) IF(THD.LT.0.0D0) RETURN TSPD=(THD*TK-273.15D0) F74AIR=REAL(TSPD) RETURN END REAL FUNCTION F75AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0) F75AIR=-1.0E+20 IF((FP.LT.0.01).OR.(FP.GT.37.6626)) RETURN F75AIR=-1.0E+10 PI=DBLE(FP)/PK IF(ABS(PI-1.0D0).LT.1.0D-05) THEN THDD=1.0D0 ELSE CALL S22AIR(PI,THDD) IF(THDD.LT.0.0D0) RETURN END IF TSPDD=(THDD*TK-273.15D0) F75AIR=REAL(TSPDD) RETURN END REAL FUNCTION F76AIR(FP) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0) F76AIR=-1.0E+20 IF((FP.LT.0.01).OR.(FP.GT.37.6626)) RETURN F76AIR=-1.0E+10 PI=DBLE(FP)/PK IF(ABS(PI-1.0D0).LT.1.0D-05) THEN CALL S35AIR(1,1.0D0,1.0D0,CVPDD) ELSE CALL S22AIR(PI,THDD) CALL S31AIR(PI,THDD,DELTDD) CALL S35AIR(1,THDD,DELTDD,CVPDD) END IF IF(CVPDD.LT.0.0D0) RETURN F76AIR=REAL(CVPDD) RETURN END REAL FUNCTION F77AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F77AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK IF(FP.EQ.0.0) THEN IF(THETA.LT.1.6838968D0) THEN CALL S35AIR(1,THETA,0.0D0,CVPT0) ELSE CALL S15AIR(1,THETA,0.0D0,CVPT0) END IF F77AIR=REAL(CVPT0) RETURN END IF CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F77AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S31AIR(PI,THETA,DELTA) CALL S35AIR(1,THETA,DELTA,CVPT) ELSE CALL S11AIR(PI,THETA,DELTA) CALL S15AIR(1,THETA,DELTA,CVPT) END IF IF(CVPT.LT.0.0D0) RETURN F77AIR=REAL(CVPT) RETURN END REAL FUNCTION F78AIR(FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F78AIR=-1.0E+20 IF((FCT.LT.-213.16).OR.(FCT.GT.-140.62)) RETURN F78AIR=-1.0E+10 THETA=(DBLE(FCT)+273.15D0)/TK IF(ABS(THETA-1.0D0).LT.1.0D-05) THEN CALL S35AIR(1,1.0D0,1.0D0,CVTDD) ELSE CALL S21AIR(THETA,PIDD) CALL S31AIR(PIDD,THETA,DELTDD) CALL S35AIR(1,THETA,DELTDD,CVTDD) END IF IF(CVTDD.LT.0.0D0) RETURN F78AIR=REAL(CVTDD) RETURN END REAL FUNCTION F79AIR(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F79AIR=-1.0E+20 CALL S91AIR(2,79,FP,FS,U,ILL91) IF(ILL91.EQ.10000) RETURN F79AIR=-1.0E+10 IF(ILL91.EQ.1000) RETURN F79AIR=REAL(U) RETURN END REAL FUNCTION F80AIR(FP,FS) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) F80AIR=-1.0E+20 CALL S91AIR(2,80,FP,FS,V,ILL91) IF(ILL91.EQ.10000) RETURN F80AIR=-1.0E+10 IF(ILL91.EQ.1000) RETURN F80AIR=REAL(V) RETURN END REAL FUNCTION F81AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) FCPPT=F18AIR(FP,FCT) FAMUPT=F13AIR(FP,FCT) FALMPT=F8AIR(FP,FCT) * F81AIR=-1.0E+20 IF((FCPPT.EQ.-1.0E+20).OR.(FAMUPT.EQ.-1.0E+20) - .OR.(FALMPT.EQ.-1.0E+20)) RETURN F81AIR=-1.0E+10 IF((FCPPT.EQ.-1.0E+10).OR.(FAMUPT.EQ.-1.0E+10) - .OR.(FALMPT.EQ.-1.0E+10)) RETURN * FPRPT=FCPPT*FAMUPT/FALMPT F81AIR=FPRPT RETURN END REAL FUNCTION F82AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0) F82AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK IF(FP.EQ.0.0) THEN IF(THETA.LT.1.6838968D0) THEN CALL S35AIR(1,THETA,0.0D0,CV0) CALL S35AIR(2,THETA,0.0D0,CP0) ELSE CALL S15AIR(1,THETA,0.0D0,CV0) CALL S15AIR(2,THETA,0.0D0,CP0) END IF AKPT0=CP0/CV0 F82AIR=REAL(AKPT0) RETURN END IF CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F82AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S36AIR(1,PI,THETA,AKPT) ELSE CALL S16AIR(1,PI,THETA,AKPT) END IF IF(AKPT.LT.0.0D0) RETURN F82AIR=REAL(AKPT) RETURN END REAL FUNCTION F83AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,GASC=287.22D0) F83AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK IF(FP.EQ.0.0) THEN IF(THETA.LT.1.6838968D0) THEN CALL S35AIR(1,THETA,0.0D0,CV0) CALL S35AIR(2,THETA,0.0D0,CP0) ELSE CALL S15AIR(1,THETA,0.0D0,CV0) CALL S15AIR(2,THETA,0.0D0,CP0) END IF W0=SQRT(GASC*(THETA*TK)*(CP0/CV0)) F83AIR=REAL(W0) RETURN END IF CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F83AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S36AIR(2,PI,THETA,WPT) ELSE CALL S16AIR(2,PI,THETA,WPT) END IF IF(WPT.LT.0.0D0) RETURN F83AIR=REAL(WPT) RETURN END REAL FUNCTION F90AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,GASC=287.22D0) ***** OUTPUT:F90AIR=BSPT IN (1/BAR) ***** F90AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F90AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S36AIR(3,PI,THETA,BSPT) ELSE CALL S16AIR(3,PI,THETA,BSPT) END IF F90AIR=REAL(BSPT*1.0D+05) RETURN END REAL FUNCTION F91AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,GASC=287.22D0) ***** OUTPUT:F91AIR=BTPT IN (1/BAR) ***** F91AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F91AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S36AIR(4,PI,THETA,BTPT) ELSE CALL S16AIR(4,PI,THETA,BTPT) END IF F91AIR=REAL(BTPT*1.0D+05) RETURN END REAL FUNCTION F92AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,GASC=287.22D0) ***** OUTPUT:F92AIR=BPPT IN (1/K)=(1/C) ***** F92AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK IF(FP.EQ.0.0) THEN BP0=1.0D0/(THETA*TK) F92AIR=REAL(BP0) RETURN END IF CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F92AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S36AIR(5,PI,THETA,BPPT) ELSE CALL S16AIR(5,PI,THETA,BPPT) END IF F92AIR=REAL(BPPT) RETURN END REAL FUNCTION F93AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,GASC=287.22D0) ***** OUTPUT:F93AIR=BVPT IN (1/K)=(1/C) ***** F93AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK IF(FP.EQ.0.0) THEN BV0=1.0D0/(THETA*TK) F93AIR=REAL(BV0) RETURN END IF CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F93AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S36AIR(6,PI,THETA,BVPT) ELSE CALL S16AIR(6,PI,THETA,BVPT) END IF F93AIR=REAL(BVPT) RETURN END REAL FUNCTION F94AIR(FP,FCT) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(TK=132.52D0,GASC=287.22D0) ***** OUTPUT:F94AIR=AJTPT IN (K/BAR)=(C/BAR) ***** F94AIR=-1.0E+20 THETA=(DBLE(FCT)+273.15D0)/TK IF(FP.EQ.0.0) THEN AJT0=0.0D0 F94AIR=REAL(AJT0) RETURN END IF CALL S90AIR(FP,FCT,PI,THETA,ILL90) IF(ILL90.NE.0) RETURN F94AIR=-1.0E+10 IF(THETA.LT.1.6838968D0) THEN CALL S36AIR(7,PI,THETA,AJTPT) ELSE CALL S16AIR(7,PI,THETA,AJTPT) END IF F94AIR=REAL(AJTPT*1.0D+05) RETURN END ***** SUBROUTINES (V81) ************************************************ SUBROUTINE S11AIR(PI,THETA,DELTA) *** EQUATION OF STATE, EQ.(28) *** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) PARAMETER(EPS=1.0D-12,ITMAX=10000) DIMENSION K(1:4,0:8),KT(0:8) DATA (K(1,I),I=0,8)/ 5.002039D0, 7.598056D0, 3.465304D0, - 0.453735D0,-0.698108D0, 0.090003D0, - 0.064788D0,-0.743022D0,-0.435251D0/ DATA (K(2,I),I=0,8)/-2.327576D0,-4.680725D0,-2.174357D0, - 3.125956D0, 5.076446D0, 1.586512D0, - 1.352110D0, 3.719013D0, 1.824133D0/ DATA (K(3,I),I=0,8)/-3.085569D0,-6.197521D0,-1.590371D0, - -0.280159D0,-5.881947D0,-6.645001D0, - -6.329591D0,-5.695303D0,-1.930506D0/ DATA (K(4,I),I=0,8)/ 1.518170D0, 3.360568D0,-0.475690D0, - -3.960499D0, 1.588962D0, 7.888349D0, - 7.082484D0, 2.425508D0, 0.0D0/ DELTA=-1.0D+20 IF((PI.LT.0.0D0).OR.(THETA.LE.0.0D0)) RETURN DO 1 I=0,8 KT(I)=K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2 1 CONTINUE XI1=(PI-KT(0))/KT(1) IF((1.0D0+4.0D0*KT(2)*(PI-KT(0))/KT(1)**2).LE.0.0D0) THEN XI2=XI1 ELSE XI2=(2.0D0*(PI-KT(0))/KT(1))/ - (SQRT(1.0D0+4.0D0*KT(2)*(PI-KT(0))/KT(1)**2)+1.0D0) END IF XI=(XI1+XI2)/2.0D0 IF(XI.GE.1.0D0) THEN XI=1.0D0 END IF DO 2 IT=1,ITMAX F1=KT(0)+(KT(1)+(KT(2)+(KT(3)+(KT(4)+(KT(5)+(KT(6) - +(KT(7)+KT(8)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI- PI DF1=KT(1)+(2.0D0*KT(2)+(3.0D0*KT(3)+(4.0D0*KT(4)+(5.0D0*KT(5) - +(6.0D0*KT(6)+(7.0D0*KT(7)+8.0D0*KT(8)*XI) - *XI)*XI)*XI)*XI)*XI)*XI XI=XI-F1/DF1 IF(ABS(-F1/(DF1*XI)).LT.EPS) THEN DELTA=XI+1.0D0 RETURN END IF 2 CONTINUE DELTA=-1.0D+10 RETURN END SUBROUTINE S12AIR(PI,DELTA,THETA) *** EQUATION OF STATE, EQ.(28) *** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) PARAMETER(SIGMA=3.16323D0,EPS=1.0D-12,ITMAX=10000) DIMENSION K(1:4,0:8),KXI(1:4) DATA (K(1,I),I=0,8)/ 5.002039D0, 7.598056D0, 3.465304D0, - 0.453735D0,-0.698108D0, 0.090003D0, - 0.064788D0,-0.743022D0,-0.435251D0/ DATA (K(2,I),I=0,8)/-2.327576D0,-4.680725D0,-2.174357D0, - 3.125956D0, 5.076446D0, 1.586512D0, - 1.352110D0, 3.719013D0, 1.824133D0/ DATA (K(3,I),I=0,8)/-3.085569D0,-6.197521D0,-1.590371D0, - -0.280159D0,-5.881947D0,-6.645001D0, - -6.329591D0,-5.695303D0,-1.930506D0/ DATA (K(4,I),I=0,8)/ 1.518170D0, 3.360568D0,-0.475690D0, - -3.960499D0, 1.588962D0, 7.888349D0, - 7.082484D0, 2.425508D0, 0.0D0/ THETA=-1.0D+20 IF((PI.LE.0.0D0).OR.(DELTA.LE.0.0D0)) RETURN XI=DELTA-1.0D0 DO 1 J=1,4 KXI(J)=0.0D0 DO 1 I=8,0,-1 KXI(J)=K(J,I)+KXI(J)*XI 1 CONTINUE IF(PI.LT.40.0D0) THEN THETA1=PI/(SIGMA*DELTA) ELSE THETA1=12.0D0 END IF THETA2=1.53D0 GO TO 20 15 THETA2=THETA2*1.1D0 GO TO 25 20 THETAW=THETA1 GO TO 26 25 THETAW=THETA2 26 CONTINUE DO 2 IT=1,ITMAX F2=KXI(4)+(KXI(3)+((KXI(2)-PI)+KXI(1)*THETAW)*THETAW)*THETAW DF2=KXI(3)+(2.0D0*(KXI(2)-PI)+3.0D0*KXI(1)*THETAW)*THETAW THETAW=THETAW-F2/DF2 IF(ABS(-F2/(DF2*THETAW)).LT.EPS) THEN IF(THETAW.LT.1.68D0) THEN GO TO 15 END IF CALL S11AIR(PI,THETAW,DELTAW) IF(ABS((DELTA-DELTAW)/DELTA).LT.0.05D0) THEN THETA=THETAW RETURN ELSE GO TO 15 END IF END IF 2 CONTINUE THETA=-1.0D+10 RETURN END SUBROUTINE S13AIR(THETA,DELTA,PI) *** EQUATION OF STATE, EQ.(28) *** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) DIMENSION K(1:4,0:8) DATA (K(1,I),I=0,8)/ 5.002039D0, 7.598056D0, 3.465304D0, - 0.453735D0,-0.698108D0, 0.090003D0, - 0.064788D0,-0.743022D0,-0.435251D0/ DATA (K(2,I),I=0,8)/-2.327576D0,-4.680725D0,-2.174357D0, - 3.125956D0, 5.076446D0, 1.586512D0, - 1.352110D0, 3.719013D0, 1.824133D0/ DATA (K(3,I),I=0,8)/-3.085569D0,-6.197521D0,-1.590371D0, - -0.280159D0,-5.881947D0,-6.645001D0, - -6.329591D0,-5.695303D0,-1.930506D0/ DATA (K(4,I),I=0,8)/ 1.518170D0, 3.360568D0,-0.475690D0, - -3.960499D0, 1.588962D0, 7.888349D0, - 7.082484D0, 2.425508D0, 0.0D0/ PI=-1.0D+20 IF((THETA.LE.0.0D0).OR.(DELTA.LT.0.0D0)) RETURN XI=DELTA-1.0D0 PIW=0.0D0 DO 3 I=8,0,-1 PIW=K(1,I)*THETA+K(2,I)+(K(3,I)+K(4,I)/THETA)/THETA+PIW*XI 3 CONTINUE PI=PIW RETURN END SUBROUTINE S14AIR(IHSU,THETA,DELTA,HSU) *** SPECIFIC ENTHALPY, EQ.(31) *** *** SPECIFIC ENTROPY, EQ.(32) *** IMPLICIT DOUBLE PRECISION(A-H,L-M,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0,VK=3.19489D-03, - GASC=287.22D0,C=2.32888242D0) DIMENSION Q(0:9),M(1:3,0:7),MT(0:7),R(0:8),L(1:3,0:7),LT(0:7) DATA (Q(I),I=0,9)/ 0.05544316D+00, 2.32888242D+00, 0.16858981D+00, - -0.90699946D-01, 0.26894910D-01,-0.44073844D-02, - 0.43749271D-03,-0.26458042D-04, 0.90200000D-06, - -0.13333333D-07/ DATA (M(1,I),I=0,7)/ -0.887155D0, -0.735823D0, -0.004042D0, - 0.021536D0, 0.216772D0, -0.038791D0, - 0.003728D0, 0.082381D0/ DATA (M(2,I),I=0,7)/ -1.160292D0, -1.950900D0, -0.008341D0, - 0.326242D0, -0.529476D0, -0.092374D0, - -0.193293D0, -0.174370D0/ DATA (M(3,I),I=0,7)/ 0.253554D0, 1.439831D0, 0.153749D0, - -0.835323D0, 0.237077D0, 0.423264D0, - 0.383392D0, 0.0D0/ DATA (R(I),I=0,8)/ 1.60629354D+01, 0.33717962D+00,-0.13604992D+00, - 0.35859880D-01,-0.55092305D-02, 0.52499126D-03, - -0.30867716D-04, 0.10308571D-05,-0.15D-07/ DATA (L(1,I),I=0,7)/ -0.504254D0, -0.581309D0, -0.119689D0, - -0.011811D0, 0.041701D0, -0.015496D0, - -0.006717D0, 0.019657D0/ DATA (L(2,I),I=0,7)/ -0.580145D0, -0.975450D0, -0.004171D0, - 0.163121D0, -0.264738D0, -0.046187D0, - -0.096647D0, -0.087185D0/ DATA (L(3,I),I=0,7)/ 0.169036D0, 0.959887D0, 0.102499D0, - -0.556882D0, 0.158052D0, 0.282176D0, - 0.255594D0, 0.0D0/ HSU=-1.0D+20 IF((THETA.LE.0.0D0).OR.(DELTA.LT.0.0D0)) RETURN XI=DELTA-1.0D0 IF((IHSU.EQ.1).OR.(IHSU.EQ.3)) THEN QTHETA=Q(0)+(Q(1)+(Q(2)+(Q(3)+(Q(4)+(Q(5)+(Q(6)+(Q(7) - +(Q(8)+Q(9)*THETA)*THETA)*THETA)*THETA)*THETA)*THETA) - *THETA)*THETA)*THETA DO 1 I=0,7 MT(I)=M(1,I)+M(2,I)/THETA+M(3,I)/THETA**2 1 CONTINUE MTXI=MT(0)+(MT(1)+(MT(2)+(MT(3)+(MT(4)+(MT(5) - +(MT(6)+MT(7)*XI)*XI)*XI)*XI)*XI)*XI)*XI U=(GASC*TK)*(QTHETA+MTXI) IF(IHSU.EQ.3) THEN HSU=U ELSE CALL S13AIR(THETA,DELTA,PI) HSU=U+(PI*PK*1.0D+05)*(VK/DELTA) END IF ELSE IF(IHSU.EQ.2) THEN RTHETA=R(0)+(R(1)+(R(2)+(R(3)+(R(4)+(R(5)+(R(6) - +(R(7)+R(8)*THETA)*THETA)*THETA)*THETA)*THETA)*THETA) - *THETA)*THETA DO 2 I=0,7 LT(I)=L(1,I)+L(2,I)/THETA**2+L(3,I)/THETA**3 2 CONTINUE LTXI=LT(0)+(LT(1)+(LT(2)+(LT(3)+(LT(4)+(LT(5) - +(LT(6)+LT(7)*XI)*XI)*XI)*XI)*XI)*XI)*XI HSU=GASC*(RTHETA+C*LOG(THETA)-LOG(XI+1.0D0)+LTXI) ELSE HSU=-1.0D+20 END IF RETURN END SUBROUTINE S15AIR(IVP,THETA,DELTA,CVP) *** ISOCHORIC SPECIFIC HEAT, EQ.(34) *** *** ISOBARIC SPECIFIC HEAT, EQ.(35) *** IMPLICIT DOUBLE PRECISION(A-H,N-Z) PARAMETER(GASC=287.22D0) DIMENSION P0(0:8),N(1:2,0:7),NT(0:7),O(1:3,0:7),OT(0:7), - P(1:4,0:7),PT(0:7) DATA(P0(I),I=0,8)/ 2.32888242D+00, 0.33717962D+00,-0.27209984D+00, - 0.10757964D+00,-0.22036922D-01, 0.26249563D-02, - -0.1852063D-03, 0.7216D-05, -0.12D-06/ DATA (N(1,I),I=0,7)/ 1.160292D0, 1.950900D0, 0.008341D0, - -0.326242D0, 0.529476D0, 0.092374D0, - 0.193293D0, 0.174370D0/ DATA (N(2,I),I=0,7)/ -0.507108D0, -2.879661D0, -0.307498D0, - 1.670646D0, -0.474155D0, -0.846529D0, - -0.766783D0, 0.0D0/ DATA (O(1,I),I=0,7)/ 2.812431D0, 1.459629D0, 0.488763D0, - -0.233647D0, -0.158869D0, 0.209474D0, - -0.173046D0, -0.244723D0/ DATA (O(2,I),I=0,7)/ 1.734883D0, 1.749717D0, -0.855520D0, - 1.013041D0, 2.294125D0, 1.442073D0, - 2.116783D0, 1.085441D0/ DATA (O(3,I),I=0,7)/ -1.707203D0, -2.071802D0, 2.606722D0, - 1.846914D0, -3.633724D0, -5.236835D0, - -2.727518D0, 0.0D0/ DATA (P(1,I),I=0,7)/ 7.598056D0, 6.930608D0, 1.361205D0, - -2.792432D0, 0.450015D0, 0.388728D0, - -5.201154D0, -3.482008D0/ DATA (P(2,I),I=0,7)/ -4.680725D0, -4.348714D0, 9.377868D0, - 20.305784D0, 7.932560D0, 8.112660D0, - 26.033091D0, 14.593064D0/ DATA (P(3,I),I=0,7)/ -6.197521D0, -3.180742D0, -0.840477D0, - -23.527788D0,-33.225005D0,-37.977546D0, - -39.867121D0,-15.444048D0/ DATA (P(4,I),I=0,7)/ 3.360568D0, -0.951380D0,-11.881497D0, - 6.355848D0, 39.441745D0, 42.494904D0, - 16.978556D0, 0.0D0/ CVP=-1.0D+20 IF((THETA.LE.0.0D0).OR.(DELTA.LT.0.0D0)) RETURN XI=DELTA-1.0D0 CV0=GASC*(P0(0)+(P0(1)+(P0(2)+(P0(3)+(P0(4)+(P0(5) - +(P0(6)+(P0(7)+P0(8)*THETA)*THETA)*THETA)*THETA) - *THETA)*THETA)*THETA)*THETA) DO 1 I=0,7 NT(I)=(N(1,I)+N(2,I)/THETA)/THETA**2 1 CONTINUE NTXI=NT(0)+(NT(1)+(NT(2)+(NT(3)+(NT(4)+(NT(5) - +(NT(6)+NT(7)*XI)*XI)*XI)*XI)*XI)*XI)*XI CV=CV0+GASC*NTXI IF(IVP.EQ.1) THEN CVP=CV ELSE IF(IVP.EQ.2) THEN DO 2 I=0,7 OT(I)=O(1,I)+O(2,I)/THETA**2+O(3,I)/THETA**3 PT(I)=P(1,I)+P(2,I)/THETA+P(3,I)/THETA**2+P(4,I)/THETA**3 2 CONTINUE OTXI=OT(0)+(OT(1)+(OT(2)+(OT(3)+(OT(4)+(OT(5) - +(OT(6)+OT(7)*XI)*XI)*XI)*XI)*XI)*XI)*XI PTXI=PT(0)+(PT(1)+(PT(2)+(PT(3)+(PT(4)+(PT(5) - +(PT(6)+PT(7)*XI)*XI)*XI)*XI)*XI)*XI)*XI CVP=CV+GASC*(OTXI**2/PTXI) END IF RETURN END SUBROUTINE S16AIR(IQ,PI,THETA,Q) *** FROM EQUATION OF STATE, EQ.(28) *** *** IQ = 1 : Q = AK (ISENTROPIC EXPONENT, -) *** IQ = 2 : Q = W (VELOCITY OF SOUND, M/S) *** IQ = 3 : Q = BS (ADIABATIC COMPRESSIBILITY, 1/PA) *** IQ = 4 : Q = BT (ISOTHERMAL COMPRESSIBILITY, 1/PA) *** IQ = 5 : Q = BP (VOLUMETRIC EXPANSION COEFF., 1/K) *** IQ = 6 : Q = BV (PRESSURE COEFFICIENT, 1/K) *** IQ = 7 : Q = AJT(JOULE-THOMSON COEFF., K/PA) ***** NOTE THAT PK IS IN PA AND TK IS IN K ***** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) PARAMETER(PK=37.6625D+05,VK=3.19489D-03,TK= 132.52D0) DIMENSION K(1:4,0:8),KT(0:8),KD(0:8) DATA (K(1,I),I=0,8)/ 5.002039D0, 7.598056D0, 3.465304D0, - 0.453735D0,-0.698108D0, 0.090003D0, - 0.064788D0,-0.743022D0,-0.435251D0/ DATA (K(2,I),I=0,8)/-2.327576D0,-4.680725D0,-2.174357D0, - 3.125956D0, 5.076446D0, 1.586512D0, - 1.352110D0, 3.719013D0, 1.824133D0/ DATA (K(3,I),I=0,8)/-3.085569D0,-6.197521D0,-1.590371D0, - -0.280159D0,-5.881947D0,-6.645001D0, - -6.329591D0,-5.695303D0,-1.930506D0/ DATA (K(4,I),I=0,8)/ 1.518170D0, 3.360568D0,-0.475690D0, - -3.960499D0, 1.588962D0, 7.888349D0, - 7.082484D0, 2.425508D0, 0.0D0/ Q=-1.0D+20 IF ((PI.LT.0.0D0).OR.(THETA.LE.0.0D0)) RETURN IF ((IQ.LE.0).OR.(IQ.GE.8)) RETURN * CALL S11AIR(PI,THETA,DELTA) XI=DELTA-1.0D0 * *** CALCULATION OF D(PI)/D(DELTA) = DPIDD *** DO 1 I=0,8 KT(I)=K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2 1 CONTINUE DPIDD=KT(1)+(2.0D0*KT(2)+(3.0D0*KT(3)+(4.0D0*KT(4) - +(5.0D0*KT(5)+(6.0D0*KT(6)+(7.0D0*KT(7)+8.0D0*KT(8)*XI) - *XI)*XI)*XI)*XI)*XI)*XI * **** IQ=4 : Q=BT **** IF (IQ.EQ.4) THEN Q=1.0D0/(PK*DELTA*DPIDD) RETURN END IF * *** CALCULATIONS OF CV AND CP *** CALL S15AIR(1,THETA,DELTA,CV) CALL S15AIR(2,THETA,DELTA,CP) **** IQ=1 : Q=AK, IQ=2 : Q=W, IQ=3 : Q=BS **** IF (IQ.EQ.1) THEN Q=(CP/CV)*(DELTA/PI)*DPIDD RETURN ELSE IF (IQ.EQ.2) THEN Q=SQRT((CP/CV)*(PK*VK)*DPIDD) RETURN ELSE IF (IQ.EQ.3) THEN Q=(CV/CP)/(PK*DELTA*DPIDD) RETURN END IF * *** CALCULATION OF D(PI)/D(THETA) = DPIDTH *** DO 2 I=0,8 KD(I)=K(1,I)-K(3,I)/THETA**2-2.0D0*K(4,I)/THETA**3 2 CONTINUE DPIDTH=KD(0)+(KD(1)+(KD(2)+(KD(3)+(KD(4)+(KD(5)+(KD(6) - +(KD(7)+KD(8)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI * **** IQ=5 : Q=BP, IQ=6 : Q=BV, IQ=7 : Q=AJT **** IF (IQ.EQ.5) THEN Q=DPIDTH/(TK*DELTA*DPIDD) RETURN ELSE IF (IQ.EQ.6) THEN Q=DPIDTH/(TK*PI) RETURN ELSE IF (IQ.EQ.7) THEN Q=(VK/(DELTA*CP))*((THETA/DELTA)*(DPIDTH/DPIDD)-1.0D0) RETURN END IF END SUBROUTINE S21AIR(THETA,PIDD) *** DEW-POINT CURVE, EQ.(36) *** IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(ALPHA=2.53293D0,BETA=-2.53901D0,GAMMA=0.00609D0, - ZETA=271.6D0) PIDD=-1.0D+20 IF(THETA.LE.0.0D0) RETURN PIDD=1.0D0 IF(THETA.GE.1.0D0) RETURN THETA1=1.0D0-SQRT(ABS(1.0D0-THETA)) PIDD=10.0D0**(ALPHA+BETA/THETA - +GAMMA*THETA1*EXP(ZETA*(1.0D0-1.0D0/THETA))) RETURN END SUBROUTINE S22AIR(PI,THDD) *** DEW-POINT CURVE, EQ.(36) *** IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(ALPHA=2.53293D0,BETA=-2.53901D0,GAMMA=0.00609D0, - ZETA=271.6D0) PARAMETER(EPS=1.0D-12,ITMAX=10000) THDD=-1.0D+20 IF(PI.LE.0.0D0) RETURN THDD=1.0D0 IF(PI.GE.1.0D0) RETURN THDDW=BETA/(LOG10(PI)-ALPHA) DO 1 IT=1,ITMAX THDD1=1.0D0-SQRT(ABS(1.0D0-THDDW)) F=ALPHA+BETA/THDDW+GAMMA*THDD1*EXP(ZETA*(1.0D0-1.0D0/THDDW)) - -LOG10(PI) DF=-BETA/THDDW**2 - +GAMMA*EXP(ZETA*(1.0D0-1.0D0/THDDW)) - *(ZETA*THDD1/THDDW**2-0.5D0/SQRT(ABS(1.0D0-THDDW))) THDDW=THDDW-F/DF IF(ABS(-F/(DF*THDDW)).LT.EPS) THEN THDD=THDDW RETURN END IF 1 CONTINUE THDD=-1.0D+10 RETURN END SUBROUTINE S23AIR(THETA,DPIDTH) *** DEW-POINT CURVE, EQ.(36) *** IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(ALPHA=2.53293D0,BETA=-2.53901D0,GAMMA=0.00609D0, - ZETA=271.6D0,E=2.718281828D0) DPIDTH=-1.0D+20 IF((THETA.LE.0.0D0).OR.(THETA.GE.1.0D0)) RETURN THETA1=1.0D0-SQRT(ABS(1.0D0-THETA)) PIDD=10.0D0**(ALPHA+BETA/THETA - +GAMMA*THETA1*EXP(ZETA*(1.0D0-1.0D0/THETA))) DPIDTH=(PIDD/LOG10(E))*(-BETA/THETA**2 - +GAMMA*EXP(ZETA*(1.0D0-1.0D0/THETA)) - *(ZETA*THETA1/THETA**2+0.5D0/SQRT(ABS(1.0D0-THETA)))) RETURN END SUBROUTINE S24AIR(THETA,PID) *** BUBBLE-POINT CURVE, EQ.(37) *** IMPLICIT DOUBLE PRECISION(A-H,K-Z) PARAMETER(ETA=2.2997D0,KAPPA=-2.30116D0,MU=-41.503D0, - NU=-1.000755D0,TAU=75.0893D0,PSI=-1.00053D0) PID=-1.0D+20 IF(THETA.LE.0.0D0) RETURN PID=1.00215D0 IF(THETA.GT.0.99925D0) RETURN THETA1=3.0D0+TAU*SQRT(ABS(PSI+1.0D0/THETA)) PID=10.0D0**(ETA+KAPPA/THETA - +1.0D-03*THETA1*EXP(MU*(NU+1.0D0/THETA))) RETURN END SUBROUTINE S25AIR(PI,THD) *** BUBBLE-POINT CURVE, EQ.(37) *** IMPLICIT DOUBLE PRECISION(A-H,K-Z) PARAMETER(ETA=2.2997D0,KAPPA=-2.30116D0,MU=-41.503D0, - NU=-1.000755D0,TAU=75.0893D0,PSI=-1.00053D0) PARAMETER(EPS=1.0D-12,ITMAX=10000) THD=-1.0D+20 IF(PI.LE.0.0D0) RETURN THD=0.99925D0 IF(PI.GT.1.00215D0) RETURN THDW=KAPPA/(LOG10(PI)-ETA) IF(THDW.GT.0.99947D0) THEN THDW=0.999D0 END IF DO 1 IT=1,ITMAX THD1=3.0D0+TAU*SQRT(ABS(PSI+1.0D0/THDW)) F=ETA+KAPPA/THDW+1.0D-03*THD1*EXP(MU*(NU+1.0D0/THDW)) - -LOG10(PI) DF=-(KAPPA+1.0D-03*EXP(MU*(NU+1.0D0/THDW)) - *(0.5D0*TAU/SQRT(ABS(PSI+1.0D0/THDW))+MU*THD1))/THDW**2 THDW=THDW-F/DF IF(ABS(-F/(DF*THDW)).LT.EPS) THEN THD=THDW RETURN END IF 1 CONTINUE THD=-1.0D+10 RETURN END SUBROUTINE S26AIR(THETA,DPIDTH) *** BUBBLE-POINT CURVE, EQ.(37) *** IMPLICIT DOUBLE PRECISION(A-H,K-Z) PARAMETER(ETA=2.2997D0,KAPPA=-2.30116D0,MU=-41.503D0, - NU=-1.000755D0,TAU=75.0893D0,PSI=-1.00053D0, - E=2.718281828D0) DPIDTH=-1.0D+20 IF((THETA.LE.0.0D0).OR.(THETA.GE.0.99925D0)) RETURN THETA1=3.0D0+TAU*SQRT(ABS(PSI+1.0D0/THETA)) PID=10.0D0**(ETA+KAPPA/THETA+1.0D-03*THETA1 - *EXP(MU*(NU+1.0D0/THETA))) DPIDTH=-(PID/LOG10(E))*(KAPPA+1.0D-03*EXP(MU*(NU+1.0D0/THETA)) - *(0.5D0*TAU/SQRT(ABS(PSI+1.0D0/THETA))+MU*THETA1))/THETA**2 RETURN END SUBROUTINE S27AIR(THETA,DELTAD) *** DENSITY OF SATURATED LIQUID, EQ.(38A) *** IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(EPS=1.0D-12,ITMAX=10000) DIMENSION G(0:7) DATA (G(I),I=0,7)/ 1.0D0, 0.0D0, -0.178421D0, - 0.893612D0, -2.119651D0, 1.985004D0, - -0.844330D0, 0.136084D0/ DELTAD=-1.0D+20 IF(THETA.LE.0.0D0) RETURN DELTAD=1.0D0 IF(THETA.GE.1.0D0) RETURN XIW=SQRT((THETA-1.0D0)/G(2)) DO 1 IT=1,ITMAX F=G(0)+(G(1)+(G(2)+(G(3)+(G(4)+(G(5)+(G(6)+G(7)*XIW) - *XIW)*XIW)*XIW)*XIW)*XIW)*XIW-THETA DF=G(1)+(2.0D0*G(2)+(3.0D0*G(3)+(4.0D0*G(4)+(5.0D0*G(5) - +(6.0D0*G(6)+7.0D0*G(7)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW XIW=XIW-F/DF IF(ABS(-F/(DF*XIW)).LT.EPS) THEN DELTAD=XIW+1.0D0 RETURN END IF 1 CONTINUE DELTAD=-1.0D+10 RETURN END SUBROUTINE S28AIR(IHSUV,PI,LAMBDA,YP) IMPLICIT DOUBLE PRECISION(A-H,L,O-Z) PARAMETER(PKPA=37.6625D+05,TK=132.52D0,VK=3.19489D-03) YP=-1.0D+20 IF((PI.LE.0.0D0).OR.(LAMBDA.LT.0.0D0).OR.(LAMBDA.GT.1.0D0)) - RETURN YP=-1.0D+10 CALL S22AIR(PI,THDD) CALL S23AIR(THDD,DPDTDD) CALL S31AIR(PI,THDD,DELTDD) CALL S25AIR(PI,THD) CALL S26AIR(THD,DPDTD) CALL S27AIR(THD,DELTD) IF((DELTDD.LE.0.0D0).OR.(DELTD.LE.0.0D0)) RETURN VP=(1.0D0/DELTDD-LAMBDA*(1.0D0/DELTDD-1.0D0/DELTD))*VK IF((IHSUV.EQ.1).OR.(IHSUV.EQ.3)) THEN CALL S34AIR(1,THDD,DELTDD,HPDD) IF(HPDD.LE.-1.0D+10) RETURN HP=HPDD-(PKPA*VK)*(THDD+0.5D0*LAMBDA*(THD-THDD)) - *(1.0D0/DELTDD-1.0D0/DELTD) - *(DPDTDD+0.5D0*LAMBDA*(DPDTD-DPDTDD))*LAMBDA IF(IHSUV.EQ.1) THEN YP=HP ELSE YP=HP-(PI*PKPA)*VP END IF ELSE IF(IHSUV.EQ.2) THEN CALL S34AIR(2,THDD,DELTDD,SPDD) IF(SPDD.LT.0.0D0) RETURN SP=SPDD-(PKPA*VK/TK)*(1.0D0/DELTDD-1.0D0/DELTD) - *(DPDTDD+0.5D0*LAMBDA*(DPDTD-DPDTDD))*LAMBDA YP=SP ELSE IF(IHSUV.EQ.4) THEN YP=VP END IF RETURN END SUBROUTINE S29AIR(IHSUV,PI,HSUVP,LAMBDA) IMPLICIT DOUBLE PRECISION(A-H,L,O-Z) PARAMETER(PKPA=37.6625D+05,TK=132.52D0,VK=3.19489D-03, - EPS=1.0D-07,ITMAX=10000) LAMBDA=-1.0D+20 IF(PI.LE.0.0D0) RETURN LAMBDA=-1.0D+10 CALL S22AIR(PI,THDD) CALL S23AIR(THDD,DPDTDD) CALL S31AIR(PI,THDD,DELTDD) CALL S25AIR(PI,THD) CALL S26AIR(THD,DPDTD) CALL S27AIR(THD,DELTD) IF((DELTDD.LE.0.0D0).OR.(DELTD.LE.0.0D0)) RETURN CALL S28AIR(IHSUV,PI,1.0D0,HSUVD) CALL S28AIR(IHSUV,PI,0.0D0,HSUVDD) IF((HSUVD.LE.-1.0D+10).OR.(HSUVDD.LE.-1.0D+10)) RETURN LAM=(HSUVDD-HSUVP)/(HSUVDD-HSUVD) IF(IHSUV.EQ.1) THEN DO 1 IT=1,ITMAX F=HSUVDD-HSUVP-(PKPA*VK)*(THDD+0.5D0*LAM*(THD-THDD)) - *(1.0D0/DELTDD-1.0D0/DELTD) - *(DPDTDD+0.5D0*LAM*(DPDTD-DPDTDD))*LAM DF=-(PKPA*VK)*(1.0D0/DELTDD-1.0D0/DELTD) - *(THDD*DPDTDD+LAM*(THDD*(DPDTD-DPDTDD)+DPDTDD*(THD-THDD)) - +(3.0D0/4.0D0)*LAM**2*(THD-THDD)*(DPDTD-DPDTDD)) LAM=LAM-F/DF IF(ABS(-F/(DF*LAM)).LT.EPS) THEN LAMBDA=LAM RETURN END IF 1 CONTINUE ELSE IF(IHSUV.EQ.2) THEN DO 2 IT=1,ITMAX F=HSUVDD-HSUVP-(PKPA*VK/TK)*(1.0D0/DELTDD-1.0D0/DELTD) - *(DPDTDD+0.5D0*LAM*(DPDTD-DPDTDD))*LAM DF=-(PKPA*VK/TK)*(1.0D0/DELTDD-1.0D0/DELTD) - *(DPDTDD+LAM*(DPDTD-DPDTDD)) LAM=LAM-F/DF IF(ABS(-F/(DF*LAM)).LT.EPS) THEN LAMBDA=LAM RETURN END IF 2 CONTINUE ELSE IF(IHSUV.EQ.3) THEN CALL S28AIR(4,PI,0.0D0,VPDD) CALL S28AIR(4,PI,1.0D0,VPD) IF((VPDD.LE.0.0D0).OR.(VPD.LE.0.0D0)) RETURN DO 3 IT=1,ITMAX F=HSUVDD-HSUVP-(PKPA*VK)*(THDD+0.5D0*LAM*(THD-THDD)) - *(1.0D0/DELTDD-1.0D0/DELTD) - *(DPDTDD+0.5D0*LAM*(DPDTD-DPDTDD))*LAM - +(PI*PKPA)*(LAM*(VPDD-VPD)) DF=-(PKPA*VK)*(1.0D0/DELTDD-1.0D0/DELTD) - *(THDD*DPDTDD+LAM*(THDD*(DPDTD-DPDTDD)+DPDTDD*(THD-THDD)) - +(3.0D0/4.0D0)*LAM**2*(THD-THDD)*(DPDTD-DPDTDD)) - +(PI*PKPA)*(VPDD-VPD) LAM=LAM-F/DF IF(ABS(-F/(DF*LAM)).LT.EPS) THEN LAMBDA=LAM RETURN END IF 3 CONTINUE END IF IF(IHSUV.EQ.4) THEN LAMBDA=LAM RETURN END IF RETURN END SUBROUTINE S31AIR(PI,THETA,DELTA) *** EQUATION OF STATE, EQ.(43B) *** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) PARAMETER(SIGMA=3.16323D0,EPS=1.0D-12,ITMAX=10000) DIMENSION K(1:5,0:12),KT(0:12) DATA (K(1,I),I=0,12)/ 5.105452D0, 7.460852D0, 2.519800D0, - 0.013958D0, 1.535990D0, 2.819382D0, - -0.548933D0,-2.273959D0,-0.019206D0, - 0.511083D0,-0.061787D0, 0.0D0, - 0.0D0/ DATA (K(2,I),I=0,12) /-2.915652D0,-3.866052D0, 2.701200D0, - 4.692584D0,-4.818280D0,-9.559664D0, - 1.910566D0, 7.839118D0, 0.482412D0, - -1.622166D0, 0.123574D0, 0.0D0, - 0.0D0/ DATA (K(3,I),I=0,12)/ -1.9298D0, -7.6746D0, -10.5642D0, - -2.8946D0, 9.1810D0, 9.5752D0, - -2.9770D0, -7.3840D0, -0.8880D0, - 1.2000D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(4,I),I=0,12)/ 0.7400D0, 4.1268D0, 5.3432D0, - -1.4952D0, -6.3568D0, -1.6092D0, - 2.7360D0, 1.4400D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(5,I),I=0,12)/ -0.13409D0, 1.90804D0, 4.11445D0, - -11.68106D0,-22.57846D0, 29.73481D0, - 54.26586D0,-38.82202D0,-67.30346D0, - 24.84753D0, 41.52300D0, -6.10000D0, - -10.00000D0/ DELTA=-1.0D+20 IF((PI.LT.0.0D0).OR.(THETA.LE.0.0D0)) RETURN DELTA=1.0D0 IF((PI.EQ.1.0D0).AND.(THETA.EQ.1.0D0)) RETURN IF(THETA.LT.1.0D0) THEN DO 1 I=0,12 KT(I)= K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2 1 CONTINUE ELSE DO 2 I=0,12 KT(I)= K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2 - +K(5,I)*(THETA-1.0D0)*(THETA-2.0D0)/THETA**11 2 CONTINUE END IF IF(PI.LT.1.0D0) THEN XIW=PI/(SIGMA*THETA)-1.0D0 ELSE XIW1=(PI-KT(0))/KT(1) IF((1.0D0+4.0D0*KT(2)*(PI-KT(0))/KT(1)**2).LE.0.0D0) THEN XIW2=XIW1 ELSE XIW2= (2.0D0*(PI-KT(0))/KT(1)) - /(SQRT(1.0D0+4.0D0*KT(2)*(PI-KT(0))/KT(1)**2)+1.0D0) END IF XIW=(XIW1+XIW2)/2.0D0 IF(XIW.LT.0.0D0) THEN XIW=0.0D0 END IF END IF DO 3 IT=1,ITMAX F1=KT(0)+(KT(1)+(KT(2)+(KT(3)+(KT(4)+(KT(5)+(KT(6)+(KT(7) - +(KT(8)+(KT(9)+(KT(10)+(KT(11)+KT(12)*XIW) - *XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW-PI DF1=KT(1)+(2.0D0*KT(2)+(3.0D0*KT(3)+(4.0D0*KT(4)+(5.0D0*KT(5) - +(6.0D0*KT(6)+(7.0D0*KT(7)+(8.0D0*KT(8)+(9.0D0*KT(9) - +(10.0D0*KT(10)+(11.0D0*KT(11)+12.0D0*KT(12)*XIW) - *XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW)*XIW XIW=XIW-F1/DF1 IF(ABS(-F1/(DF1*XIW)).LT.EPS) THEN DELTA=XIW+1.0D0 RETURN END IF 3 CONTINUE DELTA=-1.0D+10 RETURN END SUBROUTINE S32AIR(PI,DELTA,THETA) *** EQUATION OF STATE, EQ.(43B) *** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) PARAMETER(SIGMA=3.16323D0,EPS=1.0D-12,ITMAX=10000) DIMENSION K(1:5,0:12),KXI(1:5) DATA (K(1,I),I=0,12)/ 5.105452D0, 7.460852D0, 2.519800D0, - 0.013958D0, 1.535990D0, 2.819382D0, - -0.548933D0,-2.273959D0,-0.019206D0, - 0.511083D0,-0.061787D0, 0.0D0, - 0.0D0/ DATA (K(2,I),I=0,12) /-2.915652D0,-3.866052D0, 2.701200D0, - 4.692584D0,-4.818280D0,-9.559664D0, - 1.910566D0, 7.839118D0, 0.482412D0, - -1.622166D0, 0.123574D0, 0.0D0, - 0.0D0/ DATA (K(3,I),I=0,12)/ -1.9298D0, -7.6746D0, -10.5642D0, - -2.8946D0, 9.1810D0, 9.5752D0, - -2.9770D0, -7.3840D0, -0.8880D0, - 1.2000D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(4,I),I=0,12)/ 0.7400D0, 4.1268D0, 5.3432D0, - -1.4952D0, -6.3568D0, -1.6092D0, - 2.7360D0, 1.4400D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(5,I),I=0,12)/ -0.13409D0, 1.90804D0, 4.11445D0, - -11.68106D0,-22.57846D0, 29.73481D0, - 54.26586D0,-38.82202D0,-67.30346D0, - 24.84753D0, 41.52300D0, -6.10000D0, - -10.00000D0/ THETA=-1.0D+20 IF((PI.LE.0.0D0).OR.(DELTA.LE.0.0D0)) RETURN THETA=1.0D0 IF((PI.EQ.1.0D0).AND.(DELTA.EQ.1.0D0)) RETURN XI=DELTA-1.0D0 DO 1 J=1,5 KXI(J)=0.0D0 DO 1 I=12,0,-1 KXI(J)=K(J,I)+KXI(J)*XI 1 CONTINUE IF(PI.LT.1.0D0) THEN CALL S31AIR(PI,1.0D0,DELTAK) XIPIK=DELTAK-1.0D0 THETAW=PI/(SIGMA*DELTA) DO 2 IT=1,ITMAX IF(XI.GT.XIPIK) THEN F2= KXI(1)*THETAW+KXI(2)+KXI(3)/THETAW+KXI(4)/THETAW**2-PI DF2= KXI(1)-KXI(3)/THETAW**2-2.0D0*KXI(4)/THETAW**3 ELSE F2= KXI(1)*THETAW+KXI(2)+KXI(3)/THETAW+KXI(4)/THETAW**2 - +KXI(5)*(THETAW-1.0D0)*(THETAW-2.0D0)/THETAW**11-PI DF2= KXI(1)-KXI(3)/THETAW**2-2.0D0*KXI(4)/THETAW**3 - -KXI(5)*(9.0D0*THETAW**2-30.0D0*THETAW+22.0)/THETAW**12 END IF THETAW=THETAW-F2/DF2 IF(ABS(-F2/(DF2*THETAW)).LT.EPS) THEN THETA=THETAW RETURN END IF 2 CONTINUE ELSE THETA1=3.4D0 THETA2=0.9091D0 GO TO 20 15 THETA2=THETA2*1.1D0 IF(THETA2.GT.3.4D0) GO TO 1000 GO TO 25 20 THETAW=THETA1 GO TO 26 25 THETAW=THETA2 26 CONTINUE DO 3 IT=1,ITMAX F2= KXI(1)*THETAW+KXI(2)+KXI(3)/THETAW+KXI(4)/THETAW**2 - +KXI(5)*(THETAW-1.0D0)*(THETAW-2.0D0)/THETAW**11-PI DF2= KXI(1)-KXI(3)/THETAW**2-2.0D0*KXI(4)/THETAW**3 - -KXI(5)*(9.0D0*THETAW**2-30.0D0*THETAW+22.0)/THETAW**12 THETAW=THETAW-F2/DF2 IF(ABS(-F2/(DF2*THETAW)).LT.EPS) THEN IF(THETAW.LE.1.0D0) THEN GO TO 15 END IF CALL S31AIR(PI,THETAW,DELTAW) IF(ABS((DELTA-DELTAW)/DELTA).LT.0.05D0) THEN THETA=THETAW RETURN ELSE GO TO 15 END IF END IF 3 CONTINUE END IF 1000 THETA=-1.0D+10 RETURN END SUBROUTINE S33AIR(THETA,DELTA,PI) *** EQUATION OF STATE, EQ.(43B) *** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) DIMENSION K(1:5,0:12) DATA (K(1,I),I=0,12)/ 5.105452D0, 7.460852D0, 2.519800D0, - 0.013958D0, 1.535990D0, 2.819382D0, - -0.548933D0,-2.273959D0,-0.019206D0, - 0.511083D0,-0.061787D0, 0.0D0, - 0.0D0/ DATA (K(2,I),I=0,12) /-2.915652D0,-3.866052D0, 2.701200D0, - 4.692584D0,-4.818280D0,-9.559664D0, - 1.910566D0, 7.839118D0, 0.482412D0, - -1.622166D0, 0.123574D0, 0.0D0, - 0.0D0/ DATA (K(3,I),I=0,12)/ -1.9298D0, -7.6746D0, -10.5642D0, - -2.8946D0, 9.1810D0, 9.5752D0, - -2.9770D0, -7.3840D0, -0.8880D0, - 1.2000D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(4,I),I=0,12)/ 0.7400D0, 4.1268D0, 5.3432D0, - -1.4952D0, -6.3568D0, -1.6092D0, - 2.7360D0, 1.4400D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(5,I),I=0,12)/ -0.13409D0, 1.90804D0, 4.11445D0, - -11.68106D0,-22.57846D0, 29.73481D0, - 54.26586D0,-38.82202D0,-67.30346D0, - 24.84753D0, 41.52300D0, -6.10000D0, - -10.00000D0/ PI=-1.0D+20 IF((THETA.LE.0.0D0).OR.(DELTA.LE.0.0D0)) RETURN PI=1.0D0 IF((THETA.EQ.1.0D0).AND.(DELTA.EQ.1.0D0)) RETURN XI=DELTA-1.0D0 PIW=0.0D0 IF(THETA.LT.1.0D0) THEN DO 4 I=12,0,-1 PIW=K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2+PIW*XI 4 CONTINUE ELSE DO 5 I=12,0,-1 PIW= K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2 - +K(5,I)*(THETA-1.0D0)*(THETA-2.0D0)/THETA**11+PIW*XI 5 CONTINUE END IF PI=PIW RETURN END SUBROUTINE S34AIR(IHSU,THETA,DELTA,HSU) *** SPECIFIC ENTHALPY, EQ.(45) *** *** SPECIFIC ENTROPY, EQ.(46) *** IMPLICIT DOUBLE PRECISION(A-H,L-M,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0,VK=3.19489D-03, - GASC=287.22D0,C=2.49493733D0) DIMENSION Q(0:7),M(1:4,0:11),MT(0:11),R(0:6),L(1:4,0:11),LT(0:11) DATA (Q(I),I=0,7)/-0.00835080D0, 2.49493733D0,-0.00874345D0, - 0.01079545D0,-0.00703999D0, 0.00232117D0, - -0.00032104D0, 0.00001628D0/ DATA (M(1,I),I=0,11)/ -1.149809D0, -0.921734D0, 0.310640D0, - 0.177703D0, -0.051005D0, -0.329657D0, - 0.079745D0, 0.185049D0, -0.073869D0, - 0.004341D0, 0.0D0, 0.0D0/ DATA (M(2,I),I=0,11)/ -0.321854D0, -1.220147D0, -1.206047D0, - -0.211683D0, 0.463008D0, 0.547163D0, - -0.211598D0, -0.296984D0, 0.094840D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (M(3,I),I=0,11)/ -0.371804D0, 0.701815D0, 1.255111D0, - -0.218258D0, -0.654680D0, -0.027314D0, - 0.227616D0, 0.0D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (M(4,I),I=0,11)/ -0.056495D0, -0.042390D0, 0.343987D0, - -0.010949D0, -1.078763D0, 0.305030D0, - 1.777483D0, -0.814247D0, -1.442296D0, - 0.833291D0, 0.439425D0, -0.287394D0/ DATA (R(I),I=0,6)/ 16.29445500D0,-0.01748691D0, 0.01619317D0, - -0.00938666D0, 0.00290147D0,-0.00038525D0, - 0.00001900D0/ DATA (L(1,I),I=0,11)/ -0.556447D0, -0.614002D0, -0.065309D0, - 0.026215D0, -0.007771D0, -0.100410D0, - 0.023981D0, 0.055402D0, -0.025079D0, - 0.002170D0, 0.0D0, 0.0D0/ DATA (L(2,I),I=0,11)/ -0.160927D0, -0.610073D0, -0.603024D0, - -0.105841D0, 0.231504D0, 0.273581D0, - -0.105799D0, -0.148492D0, 0.047420D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (L(3,I),I=0,11)/ -0.247869D0, 0.467877D0, 0.836741D0, - -0.145506D0, -0.436453D0, -0.018209D0, - 0.151744D0, 0.0D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (L(4,I),I=0,11)/ -0.056494D0, -0.042390D0, 0.343987D0, - -0.010949D0, -1.078763D0, 0.305030D0, - 1.777483D0, -0.814247D0, -1.442296D0, - 0.833291D0, 0.439425D0,-0.287394D0/ HSU=-1.0D+20 IF((THETA.LE.0.0D0).OR.(DELTA.LE.0.0D0)) RETURN XI=DELTA-1.0D0 IF((IHSU.EQ.1).OR.(IHSU.EQ.3)) THEN QTHETA=Q(0)+(Q(1)+(Q(2)+(Q(3)+(Q(4)+(Q(5) - +(Q(6)+Q(7)*THETA)*THETA)*THETA)*THETA)*THETA)*THETA) - *THETA IF(THETA.LT.1.0D0) THEN DO 10 I=0,11 MT(I)= M(1,I)+M(2,I)/THETA+M(3,I)/THETA**2 10 CONTINUE ELSE DO 11 I=0,11 MT(I)= M(1,I)+M(2,I)/THETA+M(3,I)/THETA**2 - +M(4,I)*(10.0D0-33.0D0/THETA+24.0D0/THETA**2)/THETA**9 11 CONTINUE END IF MTXI=MT(0)+(MT(1)+(MT(2)+(MT(3)+(MT(4)+(MT(5)+(MT(6)+(MT(7) - +(MT(8)+(MT(9)+(MT(10)+MT(11)*XI) - *XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI U=(GASC*TK)*(QTHETA+MTXI) IF(IHSU.EQ.3) THEN HSU=U ELSE CALL S33AIR(THETA,DELTA,PI) HSU=U+(PI*PK*1.0D+05)*(VK/DELTA) END IF ELSE IF(IHSU.EQ.2) THEN RTHETA=R(0)+(R(1)+(R(2)+(R(3)+(R(4)+(R(5)+R(6)*THETA) - *THETA)*THETA)*THETA)*THETA)*THETA IF(THETA.LT.1.0D0) THEN DO 20 I=0,11 LT(I)= L(1,I)+L(2,I)/THETA**2+L(3,I)/THETA**3 20 CONTINUE ELSE DO 21 I=0,11 LT(I)= L(1,I)+L(2,I)/THETA**2+L(3,I)/THETA**3 - +L(4,I)*(9.0D0-30.0D0/THETA+22.0D0/THETA**2)/THETA**10 21 CONTINUE END IF LTXI=LT(0)+(LT(1)+(LT(2)+(LT(3)+(LT(4)+(LT(5)+(LT(6)+(LT(7) - +(LT(8)+(LT(9)+(LT(10)+LT(11)*XI) - *XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI HSU=GASC*(RTHETA+C*LOG(THETA)-LOG(XI+1.0D0)+LTXI) ELSE HSU=-1.0D+20 END IF RETURN END SUBROUTINE S35AIR(IVP,THETA,DELTA,CVP) *** ISOCHORIC SPECIFIC HEAT, EQ.(48) *** *** ISOBARIC SPECIFIC HEAT, EQ.(49) *** IMPLICIT DOUBLE PRECISION(A-H,N-Z) PARAMETER(GASC=287.22D0) DIMENSION P0(0:6),N(1:3,0:11),NT(0:11),O(1:4,0:11),OT(0:11), - P(1:5,0:11),PT(0:11) DATA(P0(I),I=0,6)/ 2.49493733D0, -0.01748691D0, 0.03238635D0, - -0.02815998D0, 0.01160588D0, -0.00192628D0, - 0.00011400D0/ DATA (N(1,I),I=0,11)/ 0.321854D0, 1.220147D0, 1.206047D0, - 0.211683D0, -0.463008D0, -0.547163D0, - 0.211598D0, 0.296984D0, -0.094840D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (N(2,I),I=0,11)/ 0.743608D0, -1.403630D0, -2.510222D0, - 0.436516D0, 1.309359D0, 0.054628D0, - -0.455231D0, 0.0D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (N(3,I),I=0,11)/ 0.338965D0, 0.254342D0, -2.063925D0, - 0.065692D0, 6.472576D0, -1.830182D0, - -10.664901D0, 4.885479D0, 8.653777D0, - -4.999748D0, -2.636549D0, 1.724361D0/ DATA (O(1,I),I=0,11)/ 2.870576D0, 1.324340D0, 0.092435D0, - -0.084587D0, 0.948208D0, 0.637009D0, - -0.945650D0, -0.332899D0, 0.322100D0, - -0.034740D0, 0.0D0, 0.0D0/ DATA (O(2,I),I=0,11)/ 1.085044D0, 3.230054D0, 2.709741D0, - -1.082232D0, -4.079849D0, -1.303874D0, - 2.977713D0, 1.173993D0, -0.674708D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (O(3,I),I=0,11)/ -0.832140D0, -3.808504D0, -2.199999D0, - 3.881373D0, 3.266938D0, -1.457370D0, - -1.619300D0, 0.0D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (O(4,I),I=0,11)/ 0.075393D0, -1.148202D0, -1.165176D0, - 7.732934D0, 4.961963D0,-21.680568D0, - -8.830792D0, 30.658744D0, 7.183097D0, - -21.153795D0, -2.192802D0, 5.622570D0/ DATA (P(1,I),I=0,11)/ 7.460852D0, 5.039600D0, 0.041874D0, - 6.143960D0, 14.096910D0, -3.293598D0, - -15.917713D0, -0.153648D0, 4.599747D0, - -0.617870D0, 0.0D0, 0.0D0/ DATA (P(2,I),I=0,11)/ -3.866052D0, 5.402400D0, 14.077752D0, - -19.273120D0,-47.798320D0, 11.463396D0, - 54.873826D0, 3.859296D0,-14.599494D0, - 1.235740D0, 0.0D0, 0.0D0/ DATA (P(3,I),I=0,11)/ -7.6746D0, -21.1284D0, -8.6838D0, - 36.7240D0, 47.8760D0, -17.8620D0, - -51.6880D0, -7.1040D0, 10.8000D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (P(4,I),I=0,11)/ 4.1268D0, 10.6864D0, -4.4856D0, - -25.4272D0, -8.0460D0, 16.4160D0, - 10.0800D0, 0.0D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0/ DATA (P(5,I),I=0,11)/ 1.90804D0, 8.22890D0, -35.04318D0, - -90.31384D0, 148.67405D0, 325.59516D0, - -271.75414D0,-538.42768D0, 223.62777D0, - 415.23000D0, -67.10000D0,-120.00000D0/ CVP=-1.0D+20 IF((THETA.LE.0.0D0).OR.(DELTA.LT.0.0D0)) RETURN XI=DELTA-1.0D0 CV0=GASC*(P0(0)+(P0(1)+(P0(2)+(P0(3)+(P0(4)+(P0(5)+P0(6)*THETA) - *THETA)*THETA)*THETA)*THETA)*THETA) IF(THETA.LT.1.0D0) THEN DO 10 I=0,11 NT(I)=(N(1,I)+N(2,I)/THETA)/THETA**2 10 CONTINUE ELSE DO 11 I=0,11 NT(I)=(N(1,I)+N(2,I)/THETA)/THETA**2 - +N(3,I)*(15.0D0-55.0D0/THETA+44.0D0/THETA**2)/THETA**10 11 CONTINUE END IF NTXI=NT(0)+(NT(1)+(NT(2)+(NT(3)+(NT(4)+(NT(5)+(NT(6)+(NT(7) - +(NT(8)+(NT(9)+(NT(10)+NT(11)*XI) - *XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI CV=CV0+GASC*NTXI IF(IVP.EQ.1) THEN CVP=CV ELSE IF(IVP.EQ.2) THEN IF(THETA.LT.1.0D0) THEN DO 20 I=0,11 OT(I)=O(1,I)+O(2,I)/THETA**2+O(3,I)/THETA**3 PT(I)=P(1,I)+P(2,I)/THETA+P(3,I)/THETA**2+P(4,I)/THETA**3 20 CONTINUE ELSE DO 21 I=0,11 OT(I)=O(1,I)+O(2,I)/THETA**2+O(3,I)/THETA**3 - +O(4,I)*(9.0D0-30.0D0/THETA+22.0D0/THETA**2)/THETA**10 PT(I)=P(1,I)+P(2,I)/THETA+P(3,I)/THETA**2+P(4,I)/THETA**3 - +P(5,I)*(THETA-1.0D0)*(THETA-2.0D0)/THETA**12 21 CONTINUE END IF OTXI=OT(0)+(OT(1)+(OT(2)+(OT(3)+(OT(4)+(OT(5)+(OT(6)+(OT(7) - +(OT(8)+(OT(9)+(OT(10)+OT(11)*XI) - *XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI PTXI=PT(0)+(PT(1)+(PT(2)+(PT(3)+(PT(4)+(PT(5)+(PT(6)+(PT(7) - +(PT(8)+(PT(9)+(PT(10)+PT(11)*XI) - *XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI CVP=CV+GASC*(OTXI**2/PTXI) END IF RETURN END SUBROUTINE S36AIR(IQ,PI,THETA,Q) *** FROM EQUATION OF STATE, EQ.(43B) *** *** IQ = 1 : Q = AK (ISENTROPIC EXPONENT, -) *** IQ = 2 : Q = W (VELOCITY OF SOUND, M/S) *** IQ = 3 : Q = BS (ADIABATIC COMPRESSIBILITY, 1/PA) *** IQ = 4 : Q = BT (ISOTHERMAL COMPRESSIBILITY, 1/PA) *** IQ = 5 : Q = BP (VOLUMETRIC EXPANSION COEFF., 1/K) *** IQ = 6 : Q = BV (PRESSURE COEFFICIENT, 1/K) *** IQ = 7 : Q = AJT(JOULE-THOMSON COEFF., K/PA) ***** NOTE THAT PK IS IN PA AND TK IS IN K ***** IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) PARAMETER(PK=37.6625D+05,VK=3.19489D-03,TK= 132.52D0) DIMENSION K(1:5,0:12),KT(0:12),KD(0:12) DATA (K(1,I),I=0,12)/ 5.105452D0, 7.460852D0, 2.519800D0, - 0.013958D0, 1.535990D0, 2.819382D0, - -0.548933D0,-2.273959D0,-0.019206D0, - 0.511083D0,-0.061787D0, 0.0D0, - 0.0D0/ DATA (K(2,I),I=0,12) /-2.915652D0,-3.866052D0, 2.701200D0, - 4.692584D0,-4.818280D0,-9.559664D0, - 1.910566D0, 7.839118D0, 0.482412D0, - -1.622166D0, 0.123574D0, 0.0D0, - 0.0D0/ DATA (K(3,I),I=0,12)/ -1.9298D0, -7.6746D0, -10.5642D0, - -2.8946D0, 9.1810D0, 9.5752D0, - -2.9770D0, -7.3840D0, -0.8880D0, - 1.2000D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(4,I),I=0,12)/ 0.7400D0, 4.1268D0, 5.3432D0, - -1.4952D0, -6.3568D0, -1.6092D0, - 2.7360D0, 1.4400D0, 0.0D0, - 0.0D0, 0.0D0, 0.0D0, - 0.0D0/ DATA (K(5,I),I=0,12)/ -0.13409D0, 1.90804D0, 4.11445D0, - -11.68106D0,-22.57846D0, 29.73481D0, - 54.26586D0,-38.82202D0,-67.30346D0, - 24.84753D0, 41.52300D0, -6.10000D0, - -10.00000D0/ Q=-1.0D+20 IF ((PI.LT.0.0D0).OR.(THETA.LE.0.0D0)) RETURN IF ((IQ.LE.0).OR.(IQ.GE.8)) RETURN * CALL S31AIR(PI,THETA,DELTA) XI=DELTA-1.0D0 * *** CALCULATION OF D(PI)/D(DELTA) = DPIDD *** IF (THETA.LT.1.0D0) THEN DO 1 I=0,12 KT(I)=K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2 1 CONTINUE ELSE DO 2 I=0,12 KT(I)=K(1,I)*THETA+K(2,I)+K(3,I)/THETA+K(4,I)/THETA**2 - +K(5,I)*(THETA-1.0D0)*(THETA-2.0D0)/THETA**11 2 CONTINUE END IF DPIDD=KT(1)+(2.0D0*KT(2)+(3.0D0*KT(3)+(4.0D0*KT(4)+(5.0D0*KT(5) - +(6.0D0*KT(6)+(7.0D0*KT(7)+(8.0D0*KT(8)+(9.0D0*KT(9) - +(10.0D0*KT(10)+(11.0D0*KT(11)+12.0D0*KT(12)*XI) - *XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI * **** IQ=4 : Q=BT **** IF (IQ.EQ.4) THEN Q=1.0D0/(PK*DELTA*DPIDD) RETURN END IF * *** CALCULATIONS OF CV AND CP *** CALL S35AIR(1,THETA,DELTA,CV) CALL S35AIR(2,THETA,DELTA,CP) * IF (IQ.EQ.1) THEN Q=(CP/CV)*(DELTA/PI)*DPIDD RETURN ELSE IF (IQ.EQ.2) THEN Q=SQRT((CP/CV)*(PK*VK)*DPIDD) RETURN ELSE IF (IQ.EQ.3) THEN Q=(CV/CP)/(PK*DELTA*DPIDD) RETURN END IF * *** CALCULATION OF D(PI)/D(THETA) = DPIDTH *** IF (THETA.LT.1.0D0) THEN DO 3 I=0,12 KD(I)=K(1,I)-K(3,I)/THETA**2-2.0D0*K(4,I)/THETA**3 3 CONTINUE ELSE DO 4 I=0,12 KD(I)=K(1,I)-K(3,I)/THETA**2-2.0D0*K(4,I)/THETA**3 - -K(5,I)*(9.0D0*THETA**2-30.0D0*THETA+22.0D0)/THETA**12 4 CONTINUE END IF DPIDTH=KD(0)+(KD(1)+(KD(2)+(KD(3)+(KD(4)+(KD(5) - +(KD(6)+(KD(7)+(KD(8)+(KD(9)+(KD(10)+(KD(11)+KD(12)*XI) - *XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI)*XI **** IQ=5 : Q=BP, IQ=6 : Q=BV, IQ=7 : Q=AJT **** IF (IQ.EQ.5) THEN Q=DPIDTH/(TK*DELTA*DPIDD) RETURN ELSE IF (IQ.EQ.6) THEN Q=DPIDTH/(TK*PI) RETURN ELSE IF (IQ.EQ.7) THEN Q=(VK/(DELTA*CP))*((THETA/DELTA)*(DPIDTH/DPIDD)-1.0D0) RETURN END IF END SUBROUTINE S40AIR(IHS,IREG,PI,Y,TH1,TH2,Y1,Y2,THETA) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(EPS=1.0D-07,ITMAX=10000) DO 10 IT=1,ITMAX THDELT=(Y-Y2)*(TH2-TH1)/(Y2-Y1) IF(ABS(THDELT/TH1).LT.EPS) THEN THETA=TH2 RETURN END IF TH1=TH2 Y1=Y2 TH2=TH2+THDELT IF(IREG.EQ.1) THEN CALL S11AIR(PI,TH2,DELTA2) CALL S14AIR(IHS,TH2,DELTA2,Y2) ELSE IF(IREG.EQ.2) THEN CALL S31AIR(PI,TH2,DELTA2) CALL S34AIR(IHS,TH2,DELTA2,Y2) END IF IF(Y2.LE.-1.0D+10) GO TO 1000 10 CONTINUE 1000 THETA=-1.0D+10 RETURN END SUBROUTINE S50AIR(P,T,R) **** BY VASSERMAN, ET AL. **** IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASC=287.1D0,TCR=132.5D0,VCR=3.16D-03, - EPS=1.0D-12,ITMAX=10000) DIMENSION B1(0:7),B2(0:6),B3(0:6),B4(0:5),B5(0:5),B6(0:4),B7(0:4), - B8(0:3),BT(1:8) DATA (B1(I),I=0,7) - / 0.428634670D+00, -0.847966366D+00, -0.861404401D+00, - 0.111483206D+01, -0.233976796D+01, 0.206954038D+01, - -0.675486780D+00, 0.312761031D-01 / DATA (B2(I),I=0,6) - / -0.143799318D-01, 0.994555861D+00, -0.330041831D+01, - 0.417520713D+01, -0.284529597D+00, -0.663967033D+00, - -0.375382118D+00 / DATA (B3(I),I=0,6) - / 0.254988224D+00, -0.818506542D-01, 0.101659488D+01, - -0.153393797D+01, -0.297568316D+01, 0.220848019D+01, - 0.343492590D+00 / DATA (B4(I),I=0,5) - / -0.622279028D+00, -0.736240568D-01, 0.305277765D+00, - 0.222879604D+01, 0.368788683D+00, -0.624927407D+00 / DATA (B5(I),I=0,5) - / 0.950714847D+00, -0.543639606D+00, -0.104265746D+01, - -0.555836584D+00, -0.138965988D+00, -0.184478015D+00 / DATA (B6(I),I=0,4) - / -0.683183608D+00, 0.849474492D+00, 0.118301745D+00, - 0.259511944D+00, 0.160186152D+00/ DATA (B7(I),I=0,4) - / 0.206032187D+00, -0.273845411D+00, -0.683725069D-02, - -0.121714815D+00, 0.440167769D-01 / DATA (B8(I),I=0,3) - / -0.191877127D-01, 0.663981061D-02, 0.500572273D-01, - -0.250093944D-01 / R=-1.0D+20 IF((P.LE.0.0D0).OR.(T.LE.0.0D0)) RETURN THETA=TCR/T BT(1)=B1(0)+(B1(1)+(B1(2)+(B1(3)+(B1(4)+(B1(5)+(B1(6)+B1(7)*THETA) - *THETA)*THETA)*THETA)*THETA)*THETA)*THETA BT(2)=B2(0)+(B2(1)+(B2(2)+(B2(3)+(B2(4)+(B2(5)+B2(6)*THETA) - *THETA)*THETA)*THETA)*THETA)*THETA BT(3)=B3(0)+(B3(1)+(B3(2)+(B3(3)+(B3(4)+(B3(5)+B3(6)*THETA) - *THETA)*THETA)*THETA)*THETA)*THETA BT(4)=B4(0)+(B4(1)+(B4(2)+(B4(3)+(B4(4)+B4(5)*THETA) - *THETA)*THETA)*THETA)*THETA BT(5)=B5(0)+(B5(1)+(B5(2)+(B5(3)+(B5(4)+B5(5)*THETA) - *THETA)*THETA)*THETA)*THETA BT(6)=B6(0)+(B6(1)+(B6(2)+(B6(3)+B6(4)*THETA) - *THETA)*THETA)*THETA BT(7)=B7(0)+(B7(1)+(B7(2)+(B7(3)+B7(4)*THETA) - *THETA)*THETA)*THETA BT(8)=B8(0)+(B8(1)+(B8(2)+B8(3)*THETA) - *THETA)*THETA IF(P.GT.37.7D0) THEN OMEGA=(2.0D0*P*VCR/(GASC*T))/ - (SQRT(ABS(1.0D0+4.0D0*BT(1)*P*VCR/(GASC*T)))+1.0D0) ELSE OMEGA=(P/(GASC*T))*VCR END IF DO 1 IT=1,ITMAX F=(P/(GASC*T))*(VCR/OMEGA)-1.0D0 - -(BT(1)+(BT(2)+(BT(3)+(BT(4)+(BT(5)+(BT(6)+(BT(7)+BT(8)*OMEGA) - *OMEGA)*OMEGA)*OMEGA)*OMEGA)*OMEGA)*OMEGA)*OMEGA DF=-(P/(GASC*T))*(VCR/OMEGA**2) - -(BT(1)+(2.0D0*BT(2)+(3.0D0*BT(3)+(4.0D0*BT(4)+(5.0D0*BT(5) - +(6.0D0*BT(6)+(7.0D0*BT(7)+8.0D0*BT(8)*OMEGA) - *OMEGA)*OMEGA)*OMEGA)*OMEGA)*OMEGA)*OMEGA) OMEGA=OMEGA-F/DF IF(ABS(-F/(DF*OMEGA)).LT.EPS) THEN R=OMEGA/VCR ITER=IT RETURN END IF 1 CONTINUE R=-1.0D+10 ITER=IT RETURN END SUBROUTINE S51AIR(P,T,R) **** BY VUKALOVICH, ET AL. **** IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER(GASC=287.097D0,EPS=1.0D-12,ITMAX=10000) DIMENSION B1(0:4),B2(0:4),B3(0:3),B4(0:3),BT(1:4) DATA (B1(I),I=0,4)/ 1.261360D0, -1.102233D0, 0.448216D0, - -1.620055D0, 0.762740D0/ DATA (B2(I),I=0,4)/ 0.461681D0, 3.123778D0,-10.751605D0, - 16.262414D0, -7.167238D0/ DATA (B3(I),I=0,3)/ 4.091575D0, 10.963910D0,-32.618341D0, - 16.315865D0/ DATA (B4(I),I=0,3)/-28.019642D0, 69.497330D0,-47.463328D0, - 10.439594D0/ R=-1.0D+20 IF((P.LE.0.0D0).OR.(T.LE.0.0D0)) RETURN TAU=T/304.2D0 BT(1)=B1(0)+(B1(1)+(B1(2)+(B1(3)+B1(4)/TAU) - /TAU)/TAU)/TAU BT(2)=B2(0)+(B2(1)+(B2(2)+(B2(3)+B2(4)/TAU) - /TAU)/TAU)/TAU BT(3)=B3(0)+(B3(1)+(B3(2)+B3(3)/TAU) - /TAU)/TAU BT(4)=B4(0)+(B4(1)+(B4(2)+B4(3)/TAU) - /TAU)/TAU RHO=(P/(GASC*T))*1.0D-03 DO 1 IT=1,ITMAX F=(P/(GASC*T))*(1.0D-03/RHO)-1.0D0 - -(BT(1)+(BT(2)+(BT(3)+BT(4)*RHO)*RHO)*RHO)*RHO DF=-(P/(GASC*T))*(1.0D-03/RHO**2) - -(BT(1)+(2.0D0*BT(2)+(3.0D0*BT(3)+4.0D0*BT(4)*RHO)*RHO)*RHO) RHO=RHO-F/DF IF(ABS(-F/(DF*RHO)).LT.EPS) THEN R=RHO*1.0D+03 ITER=IT RETURN END IF 1 CONTINUE R=-1.0D+10 ITER=IT RETURN END SUBROUTINE S52AIR(P,T,R) **** BY NAGASHIMA, ET AL. **** IMPLICIT DOUBLE PRECISION(A-H,M,O-Z) PARAMETER(M=28.9644D0,UGASC=8.31433D+03,BSTAR=20.5D-06, - EPS=1.0D-12,ITMAX=10000) DIMENSION B(-4:1),C(-3:0) DATA (B(I),I=-4,1,1)/-0.168785D-02,-0.223299D-01,-0.170400D+00, - -0.194783D+01, 0.216059D+01,-0.222430D-01/ DATA (C(I),I=-3,0,1)/ 0.250287D+01,-0.552109D+01, 0.429631D+01, - 0.119665D+01/ R=-1.0D+20 IF((P.LE.0.0D0).OR.(T.LE.0.0D0)) RETURN TAU=T/340.0D0 BT=B(1)*TAU+B(0)+(B(-1)+(B(-2)+(B(-3)+B(-4)/TAU)/TAU)/TAU)/TAU CT=C(0)+(C(-1)+(C(-2)+C(-3)/TAU)/TAU)/TAU RHOB=(P/(UGASC*T))/(1.0D-03) DO 1 IT=1,ITMAX F=(P/(UGASC*T))/(1.0D-03)/RHOB - -1.0D0-BT*(BSTAR*RHOB)-CT*(BSTAR*RHOB)**2 DF=-(P/(UGASC*T))/(1.0D-03)/RHOB**2 - -BT*BSTAR-2.0D0*CT*BSTAR**2*RHOB RHOB=RHOB-F/DF IF(ABS(-F/(DF*RHOB)).LT.EPS) THEN R=RHOB*1.0D-03*M ITER=IT RETURN END IF 1 CONTINUE R=-1.0D+10 ITER=IT RETURN END SUBROUTINE S90AIR(FP,FCT,PI,THETA,ILL) IMPLICIT DOUBLE PRECISION(A-E,G-H,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0,VK=3.19489D-03) PI=DBLE(FP)/PK THETA=(DBLE(FCT)+273.15D0)/TK IF((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(THETA-1.0D0).LT.1.0D-05)) - THEN PI=1.0D0 THETA=1.0D0 ILL=0 RETURN END IF ILL=10000 IF((FCT.LT.-213.16).OR.(FCT.GT.1250.02)) RETURN IF((FP.GE.0.01).AND.(FP.LT.37.6625)) THEN FCTMIN=F75AIR(FP) CALL S22AIR(PI,THDD) IF(THDD.LT.0.0D0) GO TO 1000 IF(ABS(THETA-THDD).LT.1.0D-05) THETA=THDD ELSE IF((FP.GE.37.6625).AND.(FP.LE.206.1)) THEN FCTMIN=-140.631 ELSE IF((FP.GT.206.1).AND.(FP.LT.927.2)) THEN DELMAX=VK/0.001420D0 CALL S32AIR(PI,DELMAX,THMIN) IF(THMIN.LT.0.0D0) GO TO 1000 FCTMIN=REAL(THMIN*TK)-273.16 ELSE IF((FP.GE.927.2).AND.(FP.LE.1691.0)) THEN FCTMIN=-50.01 ELSE IF((FP.GT.1691.0).AND.(FP.LE.3500.0)) THEN DELMAX=VK/0.001204D0 CALL S12AIR(PI,DELMAX,THMIN) IF(THMIN.LT.0.0D0) GO TO 1000 FCTMIN=REAL(THMIN*TK)-273.16 ELSE IF((FP.GT.3500.0).AND.(FP.LT.4000.0)) THEN FCTMIN=-6.65E-04*FP**2+5.45*FP-10724.0 ELSE IF((FP.GE.4000.0).AND.(FP.LE.4500.01)) THEN FCTMIN=-2.0E-04*FP**2+1.98*FP-4284.0 ELSE RETURN END IF IF(FCT.LE.FCTMIN) RETURN ILL=0 RETURN 1000 ILL=1000 RETURN END SUBROUTINE S91AIR(IHS,IFUN,FP,FHS,THUV,ILL) ***** IFUN=64(TPH),65(TPS),71(HPS),79(UPS) OR 80(VPS) ***** IMPLICIT DOUBLE PRECISION(A-E,G-H,L,O-Z) PARAMETER(PK=37.6625D0,TK=132.52D0,VK=3.19489D-03) THUV=-1.0D+20 ILL=10000 IF((FP.LT.0.0099).OR.(FP.GT.4500.02)) RETURN PI=DBLE(FP)/PK HS=DBLE(FHS) CALL S34AIR(IHS,1.0D0,1.0D0,HSK) IF((ABS(PI-1.0D0).LT.1.0D-05).AND.(ABS(HS/HSK-1.0D0).LT.1.0D-05)) - THEN IF((IFUN.EQ.64).OR.(IFUN.EQ.65)) THEN THUV=TK ELSE IF(IFUN.EQ.71) THEN CALL S34AIR(1,1.0D0,1.0D0,HK) THUV=HK ELSE IF(IFUN.EQ.79) THEN CALL S34AIR(3,1.0D0,1.0D0,UK) THUV=UK ELSE IF(IFUN.EQ.80) THEN CALL S31AIR(1.0D0,1.0D0,DELTAK) THUV=VK/DELTAK END IF ILL=0 RETURN END IF THMAX=(1250.02D0+273.15D0)/TK CALL S11AIR(PI,THMAX,DELMIN) CALL S14AIR(IHS,THMAX,DELMIN,HSMAX) IF(HSMAX.LE.-1.0D+10) RETURN IF(FP.LE.1691.0) THEN TH1=(-50.01D0+273.15D0)/TK CALL S11AIR(PI,TH1,DEL1) CALL S14AIR(IHS,TH1,DEL1,HS1) IF(HS1.LE.-1.0D+10) RETURN END IF IF(FP.GE.4000.0) THEN THMIN1=(-2.0D-04*(PI*PK)**2+1.98D0*(PI*PK)-4284.0D0+273.15D0) - /TK CALL S11AIR(PI,THMIN1,DEL1) CALL S14AIR(IHS,THMIN1,DEL1,HS1) IF(HS1.LE.-1.0D+10) RETURN HSMIN1=HS1 HSMIN2=HSMIN1 HSMIN3=HSMIN1 THI1=THMIN1 HSI1=HS1 THI2=THI1 HSI2=HSI1 ELSE IF((FP.GT.3500.0).AND.(FP.LT.4000)) THEN THMIN1=(-6.65D-04*(PI*PK)**2+5.45D0*(PI*PK)-10724.0D0+273.15D0) - /TK CALL S11AIR(PI,THMIN1,DEL1) CALL S14AIR(IHS,THMIN1,DEL1,HS1) IF(HS1.LE.-1.0D+10) RETURN HSMIN1=HS1 HSMIN2=HSMIN1 HSMIN3=HSMIN1 THI1=THMIN1 HSI1=HS1 THI2=THI1 HSI2=HSI1 ELSE IF((FP.GT.1691.0).AND.(FP.LE.3500.0)) THEN DMIN1=VK/0.00120398D0 CALL S12AIR(PI,DMIN1,TH1) CALL S14AIR(IHS,TH1,DMIN1,HS1) IF(HS1.LE.-1.0D+10) RETURN HSMIN1=HS1 HSMIN2=HSMIN1 HSMIN3=HSMIN1 THI1=TH1 HSI1=HS1 THI2=THI1 HSI2=HSI1 ELSE IF((FP.GE.927.2).AND.(FP.LE.1691.0)) THEN HSMIN1=HS1 HSMIN2=HSMIN1 HSMIN3=HSMIN1 THI1=TH1 HSI1=HS1 THI2=THI1 HSI2=HSI1 ELSE IF((FP.GT.206.1).AND.(FP.LT.927.2)) THEN HSMIN1=HS1 DMIN2=VK/0.0014198D0 CALL S32AIR(PI,DMIN2,TH2) CALL S34AIR(IHS,TH2,DMIN2,HS2) IF(HS2.LE.-1.0D+10) RETURN HSMIN2=HS2 HSMIN3=HSMIN2 THI1=TH1 HSI1=HS1 THI2=TH2 HSI2=HS2 ELSE IF((FP.GE.37.6625).AND.(FP.LE.206.1)) THEN HSMIN1=HS1 CALL S31AIR(PI,1.0D0,DELK) CALL S34AIR(IHS,1.0D0,DELK,HS2) IF(HS2.LE.-1.0D+10) RETURN HSMIN2=HS2 HSMIN3=HSMIN2 THI1=TH1 HSI1=HS1 THI2=1.0D0 HSI2=HS2 ELSE IF((FP.GT.36.0).AND.(FP.LT.37.6625)) THEN HSMIN1=HS1 CALL S28AIR(IHS,PI,0.0D0,HSPDD) IF(HSPDD.LE.-1.0D+10) RETURN HSMIN2=HSPDD HSMIN3=HSMIN2 THI1=TH1 HSI1=HS1 CALL S22AIR(PI,THDD) IF(THDD.LT.0.0D0) RETURN THI2=THDD HSI2=HSPDD ELSE IF((FP.GE.0.4).AND.(FP.LE.36.0)) THEN HSMIN1=HS1 CALL S28AIR(IHS,PI,0.0D0,HSPDD) CALL S28AIR(IHS,PI,1.0D0,HSPD) IF((HSPDD.LE.-1.0D+10).OR.(HSPD.LE.-1.0D+10)) RETURN HSMIN2=HSPDD HSMIN3=HSPD THI1=TH1 HSI1=HS1 CALL S22AIR(PI,THDD) IF(THDD.LT.0.0D0) RETURN THI2=THDD HSI2=HSPDD ELSE IF(FP.LT.0.4) THEN HSMIN1=HS1 CALL S28AIR(IHS,PI,0.0D0,HSPDD) IF(HSPDD.LE.-1.0D+10) RETURN HSMIN2=HSPDD HSMIN3=HSMIN2 THI1=TH1 HSI1=HS1 CALL S22AIR(PI,THDD) IF(THDD.LT.0.0D0) RETURN THI2=THDD HSI2=HSPDD END IF THUV=-1.0D+10 ILL=1000 IF((HS.GE.HSMIN1).AND.(HS.LE.HSMAX)) THEN THETA0=THI1 HS0=HSI1 CALL S11AIR(PI,THETA0,DELTA0) CALL S15AIR(2,THETA0,DELTA0,CP0) IF(CP0.LT.0.0D0) RETURN IF(IFUN.EQ.64) THEN THETA1=THETA0+(HS-HS0)/(CP0*TK) ELSE IF((IFUN.EQ.65).OR.(IFUN.EQ.71).OR.(IFUN.EQ.79).OR. - (IFUN.EQ.80)) THEN THETA1=THETA0*EXP((HS-HS0)/CP0) END IF CALL S11AIR(PI,THETA1,DELTA1) CALL S14AIR(IHS,THETA1,DELTA1,HS1) CALL S40AIR(IHS,1,PI,HS,THETA0,THETA1,HS0,HS1,THETAW) IF(THETAW.LT.0.0D0) RETURN THUV=THETAW*TK IF(IFUN.EQ.71) THEN CALL S11AIR(PI,THETAW,DELTAW) CALL S14AIR(1,THETAW,DELTAW,HW) IF(HW.LE.-1.0D+10) RETURN THUV=HW ELSE IF(IFUN.EQ.79) THEN CALL S11AIR(PI,THETAW,DELTAW) CALL S14AIR(3,THETAW,DELTAW,UW) IF(UW.LE.-1.0D+10) RETURN THUV=UW ELSE IF(IFUN.EQ.80) THEN CALL S11AIR(PI,THETAW,DELTAW) IF(DELTAW.LE.-1.0D+10) RETURN THUV=VK/DELTAW END IF ELSE IF((HS.GE.HSMIN2).AND.(HS.LT.HSMIN1)) THEN THETA0=THI2 HS0=HSI2 CALL S31AIR(PI,THETA0,DELTA0) CALL S35AIR(2,THETA0,DELTA0,CP0) IF(CP0.LT.0.0D0) RETURN IF(IFUN.EQ.64) THEN THETA1=THETA0+(HS-HS0)/(CP0*TK) ELSE IF((IFUN.EQ.65).OR.(IFUN.EQ.71).OR.(IFUN.EQ.79).OR. - (IFUN.EQ.80)) THEN THETA1=THETA0*EXP((HS-HS0)/CP0) END IF CALL S31AIR(PI,THETA1,DELTA1) CALL S34AIR(IHS,THETA1,DELTA1,HS1) CALL S40AIR(IHS,2,PI,HS,THETA0,THETA1,HS0,HS1,THETAW) IF(THETAW.LT.0.0D0) RETURN THUV=THETAW*TK IF(IFUN.EQ.71) THEN CALL S31AIR(PI,THETAW,DELTAW) CALL S34AIR(1,THETAW,DELTAW,HW) IF(HW.LE.-1.0D+10) RETURN THUV=HW ELSE IF(IFUN.EQ.79) THEN CALL S31AIR(PI,THETAW,DELTAW) CALL S34AIR(3,THETAW,DELTAW,UW) IF(UW.LE.-1.0D+10) RETURN THUV=UW ELSE IF(IFUN.EQ.80) THEN CALL S31AIR(PI,THETAW,DELTAW) IF(DELTAW.LE.-1.0D+10) RETURN THUV=VK/DELTAW END IF ELSE IF((HS.GE.HSMIN3).AND.(HS.LT.HSMIN2)) THEN CALL S29AIR(IHS,PI,HS,LAMBDA) IF(LAMBDA.LT.0.0D0) RETURN IF((IFUN.EQ.64).OR.(IFUN.EQ.65)) THEN CALL S22AIR(PI,THDD) CALL S25AIR(PI,THD) IF((THDD.LT.0.0D0).OR.(THD.LT.0.0D0)) RETURN THUV=(THDD-LAMBDA*(THDD-THD))*TK ELSE IF(IFUN.EQ.71) THEN CALL S28AIR(1,PI,LAMBDA,HW) IF(HW.LE.-1.0D+10) RETURN THUV=HW ELSE IF(IFUN.EQ.79) THEN CALL S28AIR(3,PI,LAMBDA,UW) IF(UW.LE.-1.0D+10) RETURN THUV=UW ELSE IF(IFUN.EQ.80) THEN CALL S28AIR(4,PI,LAMBDA,VW) IF(VW.LE.-1.0D+10) RETURN THUV=VW END IF ELSE THUV=-1.0D+20 ILL=10000 RETURN END IF ILL=0 RETURN END *-------------------------ERROR MESSAGES FOR LEVEL 1, 2 AND 3------ SUBROUTINE S97AIR(NFUN) *** LEVEL 1 ERROR MESSAGE *** CHARACTER NFUN*6, MSG*125 INTEGER KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN MSG='**** NO CONVERGENCE AT '//NFUN//' FOR AIR ****' WRITE(6,1000) MSG 1000 FORMAT(1H ,5X,A) END IF RETURN END SUBROUTINE S98AIR(IARG,ARG1,ARG2,NARG1,NARG2,NFUN) *** LEVEL 2 ERROR MESSAGE *** *** IARG=1 FOR ONE ARGUMENT (SECOND ARGUMENT IS DUMMY) *** IARG=2 FOR TWO ARGUMENTS CHARACTER NFUN*6, NARG1*1,NARG2*1 INTEGER IARG,KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN IF (IARG.EQ.1) THEN WRITE(6,2000) NFUN,NARG1,ARG1 ELSE IF (IARG.EQ.2) THEN WRITE(6,2010) NFUN,NARG1,ARG1,NARG2,ARG2 END IF END IF 2000 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR AIR', - ' WHEN ',A1,' =', 1PE14.7,' ****') 2010 FORMAT(1H ,5X,'**** OUT OF RANGE AT ',A6,' FOR AIR', - ' WHEN ',A1,' =',1PE14.7,' AND ',A1,' =',1PE14.7,' ****') RETURN END SUBROUTINE S99AIR(NFUN) *** LEVEL 3 ERROR MESSAGE *** CHARACTER NFUN*6, MSG*125 INTEGER KPA,MESS COMMON/UNIT/KPA,MESS IF (MESS.NE.0) THEN MSG='**** FUNCTION '//NFUN//' UNAVAILABLE FOR AIR ****' WRITE(6,3000) MSG 3000 FORMAT(1H ,5X,A) END IF RETURN END REAL FUNCTION T90(T) * Conversion of temperature scal from IPTS-68 to ITS-90 * Range -200C<=T68<=630C anf T68>=1064 REAL T68,T0K REAL*8 CA(8) INTEGER KPA,MESS COMMON/UNIT/KPA,MESS DATA CA/-0.14542, -0.26722, 1.06471, 1.131286, -4.25835, + -1.59924, 7.29176, -3.52573/ IF((KPA.EQ.1).OR.(KPA.EQ.3)) THEN T0K=0.0 ELSE T0K=273.15 END IF T68=T-T0K * Check of argument range IF((-200.0.LE.T68).AND.(T68.LE.630.0)) THEN RT68=T68/630.0 TEMP=T68+RT68*(CA(1)+RT68*(CA(2)+RT68*(CA(3)+RT68*(CA(4) + +RT68*(CA(5)+RT68*(CA(6)+RT68*(CA(7)+RT68*CA(8)))))))) ELSE IF((1064.0.LE.T68).AND.(T68.LE.3000.0)) THEN TEMP=T68-0.25*(T68+273.15)*(T68+273.15)/1337.33/1337.33 ELSE IF(MESS.NE.0) WRITE(6,6020) T 6020 FORMAT(1H ,5X,'**** OUT OF RANGE AT T90', & ' WHEN T68 =',1PE14.7,' ****') T90=-1.0E+20 RETURN END IF T90=TEMP+T0K RETURN END REAL FUNCTION T68(T) * Conversion of temperature scal from ITS-90 to IPTS-68 * Range -200C<=T68<=630C anf T68>=1064 T1=T90(T) T2=T ITER=0 DT=T-T1 1000 CONTINUE ITER=ITER+1 IF(ITER.GE.10000) THEN WRITE(6,*) ' ITERATION EXCEEDED MAXIMUM ' T68=-1.0E+10 RETURN END IF T2=T2+DT T1=T90(T2) DT=T-T1 IF(ABS(DT).GE.1.0E-4) GOTO 1000 T68=T2 RETURN END * Subroutine Program Specifying KPA and MESS * MS-FORTRAN and MS-C Mixed Langage Programing * SUBROUTINE KPAMES(KPAC,MESSC) COMMON/UNIT/KPA,MESS INTEGER KPAC INTEGER MESSC KPA=KPAC MESS=MESSC RETURN END