C********************************************************************* C PROPATH --- MOIST AIR V10.1 (MARF-SS.EXE) by T.FUJITA (AUG.1996) * C********************************************************************* CHARACTER LSS(7)*78, LME(15)*63, LUP(7)*6, LUT(2)*3, LUR(2)*3, & LUX(2)*9, LUV*11, LUH(3)*12, LUS(3)*14, LL, L(10) DIMENSION UP(7),UT(2),UR(2),UX(2),UHS(3) * DATA LSS(1)/'+--------------------------------------------------- * &------------------------+'/, * & LSS(2)/'| ##### ##### #### ##### ## ####### # # * & #### # |'/, * & LSS(3)/'| # # # # # # # # # # # # # * & # # ## |'/, * & LSS(4)/'| # # # # # # # # # # # # # * & # # # |'/, * & LSS(5)/'| ##### ##### # # ##### # # # ###### * & ##### # |'/, * & LSS(6)/'| # # # # # # ###### # # # * & # # |'/, * & LSS(7)/'| # # # # # # # # # # # * & # ## # |'/, * & LSS(8)/'| # # # #### # # # # # # * & ### ## ### |'/ * DATA LSS(9)/'| * & |'/, * & LSS(10)/'| #### ### # # #### # ###### ### * &# # # #### ####### |'/, * & LSS(11)/'| # # # ## # # # # # # * & # # # # # # |'/, * & LSS(12)/'| # # # # # # # # # # * & # # # # # |'/, * & LSS(13)/'| #### # # # # # # ##### ### * &# ###### # # # |'/, * & LSS(14)/'| # # # # # # ### # # * & # # # # # # |'/, * & LSS(15)/'| # # # # ## # # # # # * & # # # # # # |'/, * & LSS(16)/'| #### ### # # #### ###### ###### ### * &# # # #### # |'/ * DATA LSS(17)/'| * & |'/, * & LSS(18)/'| ##### ##### #### #### ##### * &## # # |'/, * & LSS(19)/'| # # # # # # # # # # # * & # ## ## |'/, * & LSS(20)/'| # # # # # # # # # # # * & # # # # # |'/, * & LSS(21)/'| ##### ##### # # # ##### # * & # # # # # |'/, * & LSS(22)/'| # # # # # # ### # # ## * &#### # # # |'/, * & LSS(23)/'| # # # # # # # # # # * & # # # |'/, * & LSS(24)/'| # # # #### #### # # # * & # # # |'/ * DATA LSS(25)/'+--------------------------------------------------- * &------------------------+'/, * & LSS(26)/'| Application Program -- SS.FO * &R -- |'/, * & LSS(27)/'| Version 12.1 : May, 7 * &99 |'/, * & LSS(28)/'| Copyright PROPATH GROUP * & |'/, * & LSS(29)/'+-------------------------------------------------- * &-------------------------+'/ DATA LSS(1)/'+--------------------------------------------------- &--------------+'/, & LSS(2)/'| PROPATH Version 12.1 & |'/, & LSS(3)/'| A Program Package for Thermophysical Properties o &f Fluids |'/, & LSS(4)/'| & |'/, & LSS(5)/'| Application Program: |'/, & LSS(6)/'| Copyright PROPATH Group, May 7, 2001 & |'/, & LSS(7)/'+--------------------------------------------------- &--------------+'/ DATA LME/'======================================================== &=======', &'| SS: An Application Program for PROPATH Providing |', &'| Thermophysical Properties of MOIST AIR as Real Fluids |', &'| Ver.12.1 |', &'| No. PROPATH Functions |', &'| 1 --(P,T,WB)--> DPA RHA DSA RWA XA VA HA SA |', &'| 2 --(P,T,DP)--> WBB RHB DSB RWB XB VB HB SB |', &'| 3 --(P,T,RH)--> WBC DPC DSC RWC XC VC HC SC |', &'| 4 --(P,T,X)---> WBD DPD RHD DSD RWD VD HD SD |', &'| 5 --(P,T,H)---> WBE DPE RHE DSE RWE XE VE SE |', &'| 6 --(P,X,H)---> TF WBF DPF RHF DSF RWF VF SF |', &'| 7 ------------> PST(T) ENHFAC(P,T) |', &'| 8 ------------> FC |', &'| 9 ------------> Change System of Unit |', &'| 0 ------------> Quit |'/ DATA LUP/'[Pa] ','[kPa] ','[MPa] ','[bar] ','[ata] ','[atm] ', & '[mmHg]'/, LUT/'[K]','[C]'/, LUR/'[-]','[%]'/, & LUX/'[kg/kgDA] ','[g/kgDA] '/, LUV/'[m**3/kgDA] '/, & LUH/'[J/kgDA] ','[kJ/kgDA] ','[kcal/kgfDA]'/, & LUS/'[J/(kgDA*K)] ','[kJ/(kgDA*K)] ','[kcal/kgfDA*K]'/ DATA UP(1)/101325./,UP(2)/101.325/,UP(3)/0.101325/, & UP(4)/1.01325/,UP(5)/1.03323/,UP(6)/1./,UP(7)/760./, & UT(1)/273.15/,UT(2)/0./, UR(1)/1./,UR(2)/100./, & UX(1)/1./,UX(2)/1000./, UHS(1)/4186.8/,UHS(2)/4.1868/,UHS(3)/1./ DATA JP/1/,JT/1/,JR/1/,JX/1/,JHS/1/ WRITE(*,30) (LSS(J),J=1,7) 30 FORMAT(1H ,A78) WRITE(*,'(A)')' --- Hit RETURN Key ---' READ(*,'(A)') LL 1000 WRITE(*,31) (LME(J),J=1,4),LME(1),LME(5) 1,LME(1),(LME(J),J=6,15),LME(1) 31 FORMAT(1H ,A63) WRITE(*,'(A)')' Input No. ====> ' DO 1001 JK=1,10 1001 L(JK)=' ' READ(*,*) QNM NM=IFIX(QNM) IF(NM*(NM-9).LT.0) GOTO(100,200,300,400,500,600,700,800),NM IF(NM.EQ.0) THEN WRITE(*,'(A)')' ***** See You Again ! *****' GOTO 999 ENDIF IF(NM.NE.9) THEN WRITE(*,'(A)')' ***** Invalid No. *****' GOTO 1000 ENDIF 900 CALL SUNIT(JP,JT,JR,JX,JHS) 901 WRITE(*,'(A)')' Input No. ====> ' READ(*,*) QNU NU=IFIX(QNU) IF(NU.EQ.0) GOTO 1000 IF(NU*(NU-6).GT.0) GOTO 901 GOTO(910,920,930,940,950,960),NU 910 WRITE(*,'(A)')' ===================' WRITE(*,'(A)')' Unit for Pressure' WRITE(*,'(A)')' ===================' WRITE(*,'(A)')' 1 ---> [Pa]' WRITE(*,'(A)')' 2 ---> [kPa]' WRITE(*,'(A)')' 3 ---> [MPa]' WRITE(*,'(A)')' 4 ---> [bar]' WRITE(*,'(A)')' 5 ---> [ata]' WRITE(*,'(A)')' 6 ---> [atm]' WRITE(*,'(A)')' 7 ---> [mmHg]' WRITE(*,'(A)')' ===================' 911 WRITE(*,'(A)')' Input No. ====> ' READ(*,*) QJP JP=IFIX(QJP) IF(JP*(JP-8).GE.0) GOTO 911 GOTO 900 920 WRITE(*,'(A)')' ==============================' WRITE(*,'(A)')' Unit for Temperature' WRITE(*,'(A)')' ==============================' WRITE(*,'(A)')' 1 ---> [K]' WRITE(*,'(A)')' 2 ---> [C]' WRITE(*,'(A)')' ==============================' 921 WRITE(*,'(A)')' Input No. ====> ' READ(*,*) QJT JT=IFIX(QJT) IF(JT*(JT-3).GE.0) GOTO 921 GOTO 900 930 WRITE(*,'(A)')' ============================' WRITE(*,'(A)')' Unit for Relative Humidity' WRITE(*,'(A)')' and Degree of Saturation' WRITE(*,'(A)')' ============================' WRITE(*,'(A)')' 1 ---> [-]' WRITE(*,'(A)')' 2 ---> [%]' WRITE(*,'(A)')' ============================' 931 WRITE(*,'(A)')' Input No. ====> ' READ(*,*) QJR JR=IFIX(QJR) IF(JR*(JR-3).GE.0) GOTO 931 GOTO 900 940 WRITE(*,'(A)')' =========================' WRITE(*,'(A)')' Unit for Humidity Ratio' WRITE(*,'(A)')' =========================' WRITE(*,'(A)')' 1 ---> [kg/kgDA]' WRITE(*,'(A)')' 2 ---> [g/kgDA]' WRITE(*,'(A)')' =========================' 941 WRITE(*,'(A)')' Input No. ====> ' READ(*,*) QJX JX=IFIX(QJX) IF(JX*(JX-3).GE.0) GOTO 941 GOTO 900 950 WRITE(*,'(A)')' ==========================' WRITE(*,'(A)')' Unit for Specific Volume' WRITE(*,'(A)')' ==========================' WRITE(*,'(A)')' [m**3/kgDA]' WRITE(*,'(A)')' ==========================' WRITE(*,'(A)')' --- Hit RETURN Key ---' READ(*,'(A)') LL GOTO 900 960 WRITE(*,'(A)')' ========================================' WRITE(*,'(A)')' Unit for Specific Enthalpy and Entropy' WRITE(*,'(A)')' ========================================' WRITE(*,'(A)')' 1 ---> [J/kgDA] [J/(kgDA*K)]' WRITE(*,'(A)')' 2 ---> [kJ/kgDA] [kJ/(kgDA*K)]' WRITE(*,'(A)')' 3 ---> [kcal/kgfDA] [kcal/kgfDA*K]' WRITE(*,'(A)')' ========================================' 961 WRITE(*,'(A)')' Input No. ====> ' READ(*,*) QJHS JHS=IFIX(QJHS) IF(JHS*(JHS-4).GE.0) GOTO 961 GOTO 900 100 DO 101 J=1,2 101 L(J)=' ' DO 102 J=3,10 102 L(J)='A' WRITE(*,'(A)')' =======================================' WRITE(*,'(A)')' Calculation of DPA, RHA, DSA' WRITE(*,'(A)')' RWA, XA, VA, HA, SA' WRITE(*,'(A)')' =======================================' 110 WRITE(*,32) 'Input P',LUP(JP),' (P=0 : Quit) ====> ' 32 FORMAT(1H ,A,A,A) READ(*,*) XP IF(XP.LE.0.) GOTO 1000 WRITE(*,32) 'Input T',LUT(JT),' ====> ' READ(*,*) XT WRITE(*,32) 'Input WB',LUT(JT),' ====> ' READ(*,*) XWB YP=XP*UP(4)/UP(JP) YT=XT+UT(2)-UT(JT) YWB=XWB+UT(2)-UT(JT) CALL SCAL(NM,YP,YT,YWB,YDP,YRH,YDS,YRW,YX,YV,YH,YS,YFS) GOTO 666 200 L(1)=' ' DO 201 J=2,10 201 L(J)='B' L(3)=' ' WRITE(*,'(A)')' =======================================' WRITE(*,'(A)')' Calculation of WBB, RHB, DSB' WRITE(*,'(A)')' RWB, XB, VB, HB, SB' WRITE(*,'(A)')' =======================================' 210 WRITE(*,32) 'Input P',LUP(JP),' (P=0 : Quit) ====> ' READ(*,*) XP IF(XP.LE.0.) GOTO 1000 WRITE(*,32) 'Input T',LUT(JT),' ====> ' READ(*,*) XT WRITE(*,32) 'Input DP',LUT(JT),' ====> ' READ(*,*) XDP YP=XP*UP(4)/UP(JP) YT=XT+UT(2)-UT(JT) YDP=XDP+UT(2)-UT(JT) CALL SCAL(NM,YP,YT,YWB,YDP,YRH,YDS,YRW,YX,YV,YH,YS,YFS) GOTO 666 300 L(1)=' ' DO 301 J=2,10 301 L(J)='C' L(4)=' ' WRITE(*,'(A)')' =======================================' WRITE(*,'(A)')' Calculation of WBC, DPC, DSC' WRITE(*,'(A)')' RWC, XC, VC, HC, SC' WRITE(*,'(A)')' =======================================' 310 WRITE(*,32) 'Input P',LUP(JP),' (P=0 : Quit) ====> ' READ(*,*) XP IF(XP.LE.0.) GOTO 1000 WRITE(*,32) 'Input T',LUT(JT),' ====> ' READ(*,*) XT WRITE(*,32) 'Input RH',LUR(JR),' ====> ' READ(*,*) XRH YP=XP*UP(4)/UP(JP) YT=XT+UT(2)-UT(JT) YRH=XRH*UR(1)/UR(JR) CALL SCAL(NM,YP,YT,YWB,YDP,YRH,YDS,YRW,YX,YV,YH,YS,YFS) GOTO 666 400 L(1)=' ' DO 401 J=2,10 401 L(J)='D' L(7)=' ' WRITE(*,'(A)')' ========================================' WRITE(*,'(A)')' Calculation of WBD, DPD, RHD, DSD' WRITE(*,'(A)')' RWD, VD, HD, SD' WRITE(*,'(A)')' ========================================' 410 WRITE(*,32) 'Input P',LUP(JP),' (P=0 : Quit) ====> ' READ(*,*) XP IF(XP.LE.0.) GOTO 1000 WRITE(*,32) 'Input T',LUT(JT),' ====> ' READ(*,*) XT WRITE(*,32) 'Input X',LUX(JX),' ====> ' READ(*,*) XX YP=XP*UP(4)/UP(JP) YT=XT+UT(2)-UT(JT) YX=XX*UX(1)/UX(JX) CALL SCAL(NM,YP,YT,YWB,YDP,YRH,YDS,YRW,YX,YV,YH,YS,YFS) GOTO 666 500 L(1)=' ' DO 501 J=2,10 501 L(J)='E' L(9)=' ' WRITE(*,'(A)')' ========================================' WRITE(*,'(A)')' Calculation of WBE, DPE, RHE, DSE' WRITE(*,'(A)')' RWE, XE, VE, SE' WRITE(*,'(A)')' ========================================' 510 WRITE(*,32) 'Input P',LUP(JP),' (P=0 : Quit) ====> ' READ(*,*) XP IF(XP.LE.0.) GOTO 1000 WRITE(*,32) 'Input T',LUT(JT),' ====> ' READ(*,*) XT WRITE(*,32) 'Input H',LUH(JHS),' ====> ' READ(*,*) XH YP=XP*UP(4)/UP(JP) YT=XT+UT(2)-UT(JT) YH=XH*UHS(1)/UHS(JHS) CALL SCAL(NM,YP,YT,YWB,YDP,YRH,YDS,YRW,YX,YV,YH,YS,YFS) GOTO 666 600 DO 601 J=1,10 601 L(J)='F' L(7)=' ' L(9)=' ' WRITE(*,'(A)')' ========================================' WRITE(*,'(A)')' Calculation of TF, WBF, DPF, RHF, DSF' WRITE(*,'(A)')' RWF, VF, SF' WRITE(*,'(A)')' ========================================' 610 WRITE(*,32) 'Input P',LUP(JP),'(P=0 : Quit) ====> ' READ(*,*) XP IF(XP.LE.0.) GOTO 1000 WRITE(*,32) 'Input X',LUX(JX),' ====> ' READ(*,*) XX WRITE(*,32) 'Input H',LUH(JHS),' ====> ' READ(*,*) XH YX=XX*UX(1)/UX(JX) YH=XH*UHS(1)/UHS(JHS) YP=XP*UP(4)/UP(JP) CALL SCAL(NM,YP,YT,YWB,YDP,YRH,YDS,YRW,YX,YV,YH,YS,YFS) 666 YPS=F49L01(YT) YP=XP IF(YT.NE.-1.E20) YT=YT-UT(2)+UT(JT) IF(YWB.NE.-1.E20) YWB=YWB-UT(2)+UT(JT) IF(YDP.NE.-1.E20) YDP=YDP-UT(2)+UT(JT) IF(YRH.NE.-1.E20) YRH=YRH*UR(JR)/UR(1) IF(YDS.NE.-1.E20) YDS=YDS*UR(JR)/UR(1) IF(YX.NE.-1.E20) YX=YX*UX(JX)/UX(1) IF(YH.NE.-1.E20) YH=YH*UHS(JHS)/UHS(1) IF(YS.NE.-1.E20) YS=YS*UHS(JHS)/UHS(1) IF(YPS.NE.-1.E20) YPS=YPS*UP(JP)/UP(4) WRITE(*,33) YP, LUP(JP) WRITE(*,34) L(1),YT, LUT(JT), L(6),YRW,LUR(1) WRITE(*,35) L(2),YWB,LUT(JT), L(7), YX,LUX(JX) WRITE(*,36) L(3),YDP,LUT(JT), L(8), YV,LUV WRITE(*,37) L(4),YRH,LUR(JR), L(9), YH,LUH(JHS) WRITE(*,38) L(5),YDS,LUR(JR), L(10),YS,LUS(JHS) IF(NM.NE.6) THEN WRITE(*,39) YFS,LUR(1),YPS,LUP(JP) ELSE WRITE(*,40) YFS,LUR(1),YPS,LUP(JP) ENDIF WRITE(*,'(A)')' ' 33 FORMAT(' P =',1PE13.5,1H ,A) 34 FORMAT(' T',A1,'=',1PE13.5,1H ,A,' RW',A1,'=',E13.5,1H ,A) 35 FORMAT(' WB',A1,'=',1PE13.5,1H ,A,' X',A1,'=',E13.5,1H ,A) 36 FORMAT(' DP',A1,'=',1PE13.5,1H ,A,' V',A1,'=',E13.5,1H ,A) 37 FORMAT(' RH',A1,'=',1PE13.5,1H ,A,' H',A1,'=',E13.5,1H ,A) 38 FORMAT(' DS',A1,'=',1PE13.5,1H ,A,' S',A1,'=',E13.5,1H ,A) 39 FORMAT(' Note: ENHFAC(P,T)=',1PE13.5,1H ,A &,4X,'PS(T)=',1PE13.5,1H ,A) 40 FORMAT(' Note: ENHFAC(P,TF)=',1PE13.5,1H ,A &,4X,'PS(TF)=',1PE13.5,1H ,A) GOTO(110,210,310,410,510,610),NM 700 WRITE(*,'(A)')' ====================================' WRITE(*,'(A)')' Calculation of PST(T), ENHFAC(P,T)' WRITE(*,'(A)')' ====================================' 710 WRITE(*,32) 'Input P',LUP(JP),' (P=0 : Quit) ====> ' READ(*,*) XP IF(XP.EQ.0.) GOTO 1000 WRITE(*,32) 'Input T',LUT(JT),' ====> ' READ(*,*) XT YP=XP*UP(4)/UP(JP) YT=XT+UT(2)-UT(JT) YPS=F49L01(YT) YFS=F50L01(YP,YT) YP=XP YT=XT IF(YPS.NE.-1.E20) YPS=YPS*UP(JP)/UP(4) WRITE(*,41) YP,LUP(JP),YT,LUT(JT) WRITE(*,42) YPS,LUP(JP),YFS,LUR(1) WRITE(*,'(A)')' ' 41 FORMAT(' P=',1PE13.5,1H ,A,' T=',1PE13.5,1H ,A) 42 FORMAT(' PST=',1PE13.5,1H ,A,' ENHFAC=',1PE13.5,1H ,A) GOTO 710 800 CALL SFC(JHS) WRITE(*,'(A)')' --- Hit RETURN Key ---' READ(*,'(A)') LL GOTO 1000 999 STOP END C--SUNIT SUBROUTINE SUNIT(JP,JT,JR,JX,JHS) CHARACTER*6 LUP(7) CHARACTER*3 LUT(2),LUR(2) CHARACTER*9 LUX(2) CHARACTER*11 LUV CHARACTER*12 LUH(3) CHARACTER*14 LUS(3) DATA LUP/'[Pa] ','[kPa] ','[MPa] ','[bar] ','[ata] ','[atm] ', & '[mmHg]'/, LUT/'[K]','[C]'/, LUR/'[-]','[%]'/, & LUX/'[kg/kgDA] ','[g/kgDA] '/, LUV/'[m**3/kgDA] '/, & LUH/'[J/kgDA] ','[kJ/kgDA] ','[kcal/kgfDA]'/, & LUS/'[J/(kgDA*K)] ','[kJ/(kgDA*K)] ','[kcal/kgfDA*K]'/ WRITE(*,30) WRITE(*,31) WRITE(*,30) WRITE(*,32) LUP(JP) WRITE(*,33) LUT(JT) WRITE(*,34) LUR(JR) WRITE(*,35) LUX(JX) WRITE(*,36) LUV WRITE(*,37) LUH(JHS),LUS(JHS) WRITE(*,38) WRITE(*,30) WRITE(*,39) 30 FORMAT(' ============================================') 31 FORMAT(' No. Unit (Current) ') 32 FORMAT(' 1 ---> Pressure ',A) 33 FORMAT(' 2 ---> Temperature ',A) 34 FORMAT(' 3 ---> Relative Humidity ',A,/ & ' Degree of Saturation ') 35 FORMAT(' 4 ---> Humidity Ratio ',A) 36 FORMAT(' 5 ---> Specific Volume ',A) 37 FORMAT(' 6 ---> Enthalpy ',A,/ & ' Entropy ',A) 38 FORMAT(' 0 ---> Return to the Function Menu') 39 FORMAT(' (DA : of Dry Air)') RETURN END C--SFC SUBROUTINE SFC(JE) CHARACTER LUE(3)*13 DIMENSION UE(3) DATA R/8.31441/,AM/28.9645/,WM/18.01528/ DATA LUE/' [J/(kg*K)] ',' [kJ/(kg*K)] ',' [kcal/kgf*K]'/ DATA UE(1)/4186.8/,UE(2)/4.1868/,UE(3)/1./ YAM=AM YWM=WM RA=R/AM RW=R/WM YRA=RA*UE(JE)/UE(2) YRW=RW*UE(JE)/UE(2) WRITE(*,'(A)')' ================================================== &====' WRITE(*,'(A)')' Thermodynamic Properties of MOIST AIR as Real Flu &ids' WRITE(*,'(A)')' PROPATH VER.10.1, August 28, 1996' WRITE(*,'(A)')' ================================================== &====' WRITE(*,'(A)')' Fundamental Constants' WRITE(*,30) YAM,YRA,LUE(JE) WRITE(*,31) YWM,YRW,LUE(JE) 30 FORMAT(' Dry Air : M=',1PE13.5,' [-]',6X,'R=',E13.5,A) 31 FORMAT(' Water Vapor : M=',1PE13.5,' [-]',6X,'R=',E13.5,A) WRITE(*,'(A)')' Temperature Scale' WRITE(*,'(A)')' based on the International Practical Temperat &ure Scale of 1968' WRITE(*,'(A)')' References' WRITE(*,'(A)')' [1] R.W.Hyland & A.Wexler : ASHRAE Transactions, &89(2A), 1983, pp.520-535.' WRITE(*,'(A)')' "Formulations for the Thermodynamic Propertie &s of Dry Air from 173.15 K' WRITE(*,'(A)')' to 473.15 K, and of Saturated Moist Air from & 173.15 K to 372.15 K, at ' WRITE(*,'(A)')' Pressures to 5 MPa" ' WRITE(*,'(A)')' [2] R.W.Hyland & A.Wexler : ASHRAE Transactions, &89(2A), 1983, pp.500-519.' WRITE(*,'(A)')' "Formulations for the Thermodynamic Propertie &s of the Saturated Phases ' WRITE(*,'(A)')' of H2O from 173.15 K to 473.15 K" ' RETURN END C--SCAL SUBROUTINE SCAL(M,P,T,WB,DP,RH,DS,RW,X,V,H,S,FS) C M=1(T,WB), 2(T,DP), 3(T,RH), 4(T,X), 5(T,H), 6(X,H) C 1: F1 =DPA F22=RHA F6 =DSA F16=RWA F45=XA F34=VA F12=HA F27=SA C 2: F40=WBB F23=RHB F7 =DSB F17=RWB F46=XB F35=VB F13=HB F28=SB C 3: F41=WBC F2 =DPC F8 =DSC F18=RWC F47=XC F36=VC F14=HC F29=SC C 4: F42=WBD F3 =DPD F24=RHD F9 =DSD F19=RWD F37=VD F15=HD F30=SD C 5: F43=WBE F4 =DPE F25=RHE F10=DSE F20=RWE F48=XE F38=VE F31=SE C 6: F33=TF F44=WBF F5 =DPF F26=RHF F11=DSF F21=RWF F39=VF F32=SF DATA EM/0.621978/ IF((P.LE.0.).OR.(P.GT.50.)) GOTO 99 PK=P*1.E5 IF(M.NE.6) THEN TK=G1L01(T) IF(TK.EQ.-1.E20) GOTO 99 ENDIF IF(M.GE.5) HK=H*1.E-3 GOTO(1,2,3,4,5,6),M 1 WBK=G1L01(WB) IF(WBK.EQ.-1.E20) GOTO 99 CALL S2L01(PK,TK,WBK,RW,EX,EH) GOTO 9 2 DPK=G1L01(DP) IF(DPK.EQ.-1.E20) GOTO 99 RW=G5L01(PK,TK,DPK) GOTO 9 3 RW=G6L01(PK,TK,RH) GOTO 9 4 IF(X.LT.0.) GOTO 99 RW=X/(X+EM) CALL S1L01(3,PK,TK,E1,E2,RW,EX,E5,E6,E7,E8) IF(EX.EQ.-1.E20) RW=-1.E20 GOTO 9 5 RW=G9L01(PK,TK,HK) GOTO 9 6 IF(X.LT.0.) GOTO 99 TK=G10L01(PK,X,HK) T=G2L01(TK) IF(T.EQ.-1.E20) GOTO 99 GOTO 4 9 IF(RW.LT.0.) GOTO 99 GOTO(10,10,10,40,50,60),M 10 CALL S1L01(7,PK,TK,RWS,XS,RW,X,V,HK,SK,FS) GOTO 70 40 CALL S1L01(7,PK,TK,RWS,XS,RW,EX,V,HK,SK,FS) GOTO 70 50 CALL S1L01(7,PK,TK,RWS,XS,RW,X,V,EH,SK,FS) GOTO 80 60 CALL S1L01(7,PK,TK,RWS,XS,RW,EX,V,EH,SK,FS) GOTO 80 70 H=HK IF(HK.NE.-1.E20) H=HK*1.E3 80 S=SK IF(SK.NE.-1.E20) S=SK*1.E3 IF(M.NE.1) THEN WBK=G8L01(PK,TK,RW) WB=G2L01(WBK) ENDIF IF(M.NE.2) THEN DPK=G4L01(PK,TK,RW) DP=G2L01(DPK) ENDIF IF(M.NE.3) THEN RH=-1.E20 IF(RWS.GT.0.) RH=RW/RWS ENDIF DS=-1.E20 IF(XS.GT.0.) DS=X/XS RETURN 99 IF(M.NE.1) WB=-1.E20 IF(M.NE.2) DP=-1.E20 IF(M.NE.3) RH=-1.E20 IF((M-4)*(M-6).NE.0) X=-1.E20 IF(M.LE.4) H=-1.E20 RW=-1.E20 DS=-1.E20 V=-1.E20 S=-1.E20 FS=-1.E20 RETURN END C*********************************************************FUNCTION*** C--F49 -------------> PST [bar] FUNCTION F49L01(ZT) T=G1L01(ZT) IF(T.EQ.-1.E20) GOTO 99 F49L01=G3L01(T)*1.E-5 RETURN 99 F49L01=-1.E20 RETURN END C--F50 -------------> ENHFAC [-] FUNCTION F50L01(ZP,ZT) IF((ZP.LE.0.).OR.(ZP.GT.50.)) GOTO 99 P=ZP*1.E5 T=G1L01(ZT) IF(T.EQ.-1.E20) GOTO 99 CALL S1L01(1,P,T,E1,E2,E3,E4,E5,E6,E7,F50L01) RETURN 99 F50L01=-1.E20 RETURN END C--G1 FUNCTION G1L01(TC68) C Temperature scale, from IPTS-68 [C] to TTS [K] IF((TC68+100.)*(TC68-200.).GT.0.) GOTO 99 T=TC68+273.15 IF(TC68.GE.0.1) T=T-0.4931358+0.46094296E-2*T & -0.13746454E-4*T**2+0.12743214E-7*T**3 G1L01=T RETURN 99 G1L01=-1.E20 RETURN END C--G2 FUNCTION G2L01(TTS) C Temperature scale, from TTS [K] to IPTS-68 [C] IF(TTS.EQ.-1.E20) GOTO 99 G2L01=TTS-273.15 IF(G2L01.GE.0.1) G2L01=G2L01+0.4931358-0.46094296E-2*TTS & +0.13746454E-4*TTS**2-0.12743214E-7*TTS**3 IF((G2L01+100.)*(G2L01-200.).LE.0.) RETURN 99 G2L01=-1.E20 RETURN END C--G3 FUNCTION G3L01(T) C Saturation vapor pressure of pure liquid water and ice [Pa] C Temperature(TTS) T=173.15 to 473.15[K] DIMENSION CG(6),CM(7) DATA CG/-0.58002206E4, 0.13914993E1, -0.48640239E-1, & 0.41764768E-4, -0.14452093E-7, 0.65459673E1/ DATA CM/-0.56745359E4, 0.63925247E1, -0.96778430E-2, & 0.62215701E-6, 0.20747825E-8, -0.94840240E-12, & 0.41635019E1/ IF(T.GE.273.16) THEN G3L01=CG(6)*ALOG(T) DO 10 I=1,5 10 G3L01=G3L01+CG(I)*T**(I-2) ELSE G3L01=CM(7)*ALOG(T) DO 20 I=1,6 20 G3L01=G3L01+CM(I)*T**(I-2) ENDIF G3L01=EXP(G3L01) RETURN END C--G4 FUNCTION G4L01(P,T,RW) C Dew point temperature DP[K] DP=T CALL S1L01(3,P,DP,RWS,E2,RW,X,E5,E6,E7,E8) IF(X.LE.0.) GOTO 99 IF((RWS.GT.0.).AND.(ABS(1.-RW/RWS).LT.1.E-6)) GOTO 50 DP=173.15 CALL S1L01(1,P,DP,RWS,E2,E3,E4,E5,E6,E7,E8) IF(RW*(RWS-RW).GT.0.) GOTO 99 C approximate value PWK=P*RW*1.E-3 A=ALOG(PWK) IF(T.LE.273.15) GOTO 1 DPC=6.54+A*(14.526+A*(0.7389+0.09486*A)) &+0.4569*PWK**0.1984 IF(DPC.GE.0.) GOTO 2 1 DPC=6.09+A*(12.608+0.4959*A) 2 DP=DPC+273.15 IF(DP.GT.T) DP=T IF(DP.LT.173.15) DP=173.15 C exact value K=0 10 K=K+1 IF(K.GT.40) GOTO 99 CALL S1L01(1,P,DP,RWA,E2,E3,E4,E5,E6,E7,E8) IF(RWA.LE.0.) THEN DP=DP-0.01 GOTO 10 ENDIF IF(RWA.LT.RW) THEN DP1=DP 20 DP2=DP1+0.01 CALL S1L01(1,P,DP2,RW2,E2,E3,E4,E5,E6,E7,E8) IF(RW2.LT.0.) GOTO 40 IF(RW2.LT.RW) THEN DP1=DP2 GOTO 20 ENDIF ELSE DP2=DP 30 DP1=DP2-0.01 CALL S1L01(1,P,DP1,RW1,E2,E3,E4,E5,E6,E7,E8) IF(RW1.LT.0.) GOTO 99 IF(RW1.GT.RW) THEN DP2=DP1 GOTO 30 ENDIF ENDIF 40 K=K+1 IF(K.GT.80) GOTO 99 DP=(DP1+DP2)/2. CALL S1L01(1,P,DP,RWA,E2,E3,E4,E5,E6,E7,E8) IF((DP2-DP1.LT.1.E-4).OR.(ABS(1.-RW/RWA).LT.1.E-6)) GOTO 50 IF((RWA.LT.0.).OR.(RWA.GT.RW)) THEN DP2=DP ELSE DP1=DP ENDIF GOTO 40 50 G4L01=DP RETURN 99 G4L01=-1.E20 RETURN END C--G5 FUNCTION G5L01(P,T,DP) C Mole fraction of water vapor RW[-] IF(T.LT.DP) GOTO 99 CALL S1L01(1,P,DP,G5L01,E2,E3,E4,E5,E6,E7,E8) RETURN 99 G5L01=-1.E20 RETURN END C--G6 FUNCTION G6L01(P,T,RH) C Mole fraction of water vapor RW[-] IF(RH*(RH-1.).GT.0.) GOTO 99 CALL S1L01(1,P,T,RWS,E2,E3,E4,E5,E6,E7,E8) IF(RWS.EQ.-1.E20) GOTO 99 G6L01=RWS*RH RETURN 99 G6L01=-1.E20 RETURN END C--G7 FUNCTION G7L01(P,T) C Specific enthalpy of liquid water and ice HC[kJ/kg] DIMENSION CF(8),CL(7),CD(5) DATA CF/-0.2403360201E4, -0.140758895E1, 0.1068287657, & -0.2914492351E-3, 0.373497936E-6, -0.21203787E-9, & -0.3424442728E1, 0.1619785E-1/ DATA CL/-0.11411380E4, 0.41930463E1, -0.8134865E-4, & 0.1451133E-6, -0.1005230E-9, -0.563473, -0.036/ DATA CD/-0.647595E3,0.274292,0.2910583E-2,0.1083437E-5,0.107E-5/ PS=G3L01(T) T2=T**2 T3=T2*T T4=T3*T T5=T4*T C For liquid water and ice C VC[m3/kg]: Specific volume, DVCDT[m3/(kg K)]: dVC/dT IF(T.GT.273.15) THEN FA=CF(1)+CF(2)*T+CF(3)*T2+CF(4)*T3+CF(5)*T4+CF(6)*T5 FB=CF(2)+2.*CF(3)*T+3.*CF(4)*T2+4.*CF(5)*T3+5.*CF(6)*T4 VC=(CF(7)+CF(8)*T)/FA DVCDT=(CF(8)-VC*FB)/FA DPSDT=PS*(0.58002206E4/T2+0.65459673E1/T-0.48640239E-1 & +0.83529536E-4*T-0.43356279E-7*T2) HC=CL(1)+CL(2)*T+CL(3)*T2+CL(4)*T3+CL(5)*T4 & +CL(6)*10.**(CL(7)*(T-273.15))-0.01214+T*VC*DPSDT*1.E-3 HC=HC+(VC-T*DVCDT)*(P-PS)*1.E-3 ELSE HC=CD(1)+CD(2)*T+CD(3)*T2+CD(4)*T3+CD(5)*PS HC=HC+(CD(5)-0.371611E-12*T2)*(P-PS) ENDIF G7L01=HC RETURN END C--G8 FUNCTION G8L01(P,T,RW) C Wet-bulb temperature WB[K] WB=T CALL S1L01(5,P,T,RWS,E2,RW,X,E5,H,E7,E8) IF(X.LT.0.) GOTO 99 IF((RWS.GT.0.).AND.(ABS(1.-RW/RWS).LT.1.E-6)) GOTO 50 WB2=WB WB=173.15 CALL S2L01(P,T,WB,RWC,E5,E6) IF(RWC.GT.RW) GOTO 99 WB1=WB C approximate value 1 WB=(WB1+WB2)/2. IF(WB2-WB1.LT.0.1) GOTO 5 CALL S1L01(1,P,WB,RWS,E2,E3,E4,E5,E6,E7,E8) IF(RWS.LT.0.) THEN WB2=WB GOTO 1 ENDIF RWC=G11L01(P,T,WB) IF(ABS(RW-RWC).LT.RW*1.E-3) GOTO 5 IF(RWC.GT.RW) THEN WB2=WB ELSE WB1=WB ENDIF GOTO 1 C exact value 5 CALL S1L01(2,P,WB,RWS,XS,E3,E4,E5,HS,E7,E8) IF(RWS.LT.0.) GOTO 99 HC=G7L01(P,WB) HDEF=HS+(X-XS)*HC IF(HDEF.LT.H) THEN WB1=WB 10 WB2=WB1+0.1 CALL S1L01(2,P,WB2,RWS,XS,E3,E4,E5,HS,E7,E8) IF(RWS.LT.0.) GOTO 99 HC=G7L01(P,WB) HDEF=HS+(X-XS)*HC IF(HDEF.LT.H) THEN WB1=WB2 GOTO 10 ENDIF ELSE WB2=WB 20 WB1=WB2-0.1 CALL S1L01(2,P,WB1,RWS,XS,E3,E4,E5,HS,E7,E8) IF(RWS.LT.0.) GOTO 99 HC=G7L01(P,WB) HDEF=HS+(X-XS)*HC IF(HDEF.GT.H) THEN WB2=WB1 GOTO 20 ENDIF ENDIF K=0 30 K=K+1 IF(K.GT.40) GOTO 99 WB=(WB1+WB2)/2. CALL S1L01(2,P,WB,RWS,XS,E3,E4,E5,HS,E7,E8) HC=G7L01(P,WB) HDEF=HS+(X-XS)*HC IF((WB2-WB1.LT.1.E-4).OR.(ABS(HDEF-H).LT.1.E-4)) GOTO 50 IF((ABS(H).GE.10.).AND.(ABS(1.-HDEF/H).LT.1.E-6)) GOTO 50 IF(HDEF.LT.H) THEN WB1=WB ELSE WB2=WB ENDIF GOTO 30 50 G8L01=WB RETURN 99 G8L01=-1.E20 RETURN END C--G9 FUNCTION G9L01(P,T,H) C Mole fraction of water vapor RW[-] RW=0. CALL S1L01(5,P,T,RWS,XS,RW,E4,E5,H0,E7,E8) IF(H.LT.H0) GOTO 99 TC=T-273.15 X=(H-TC)/(2501.+1.805*TC) IF(X.LT.0.) GOTO 50 RWA=X/(0.622+X) IF(RWS.LT.0.) RWS=0.99 RW1=AMAX1(0.,0.9*RWA) RW2=AMIN1(RWS,1.1*RWA) K=0 10 K=K+1 IF(K.GT.40) GOTO 99 RW=(RW1+RW2)/2. CALL S1L01(5,P,T,RWS,E2,RW,E4,E5,H1,E7,E8) IF(RWS.LT.0.) RWS=0.99 IF((RW2-RW1.LT.RWS*1.E-6).OR.(ABS(H1-H).LT.1.E-4)) GOTO 50 IF((ABS(H).GE.10.).AND.(ABS(1.-H1/H).LT.1.E-6)) GOTO 50 IF(H1.LT.H) THEN RW1=RW ELSE RW2=RW ENDIF GOTO 10 50 G9L01=RW RETURN 99 G9L01=-1.E20 RETURN END C--G10 FUNCTION G10L01(P,X,H) C Dry-bulb temperature of moist air T[K] DATA EM/0.621978/ C approximate value TC=(H-2501.*X)/(1.+1.805*X) IF((TC.LT.-105.).OR.(TC.GT.205.)) GOTO 99 IF(TC.GT.200.) TC=200. IF(TC.LT.-100.) TC=-100. T=G1L01(TC) C exact value RW=X/(X+EM) CALL S1L01(5,P,T,E1,E2,RW,E4,E5,H1,E7,E8) IF(ABS(H1-H).LT.1.E-4) GOTO 50 IF((ABS(H).GE.10.).AND.(ABS(1.-H1/H).LT.1.E-6)) GOTO 50 K=0 IF(H1.LT.H) THEN 10 K=K+1 T=T+0.2 CALL S1L01(5,P,T,E1,E2,RW,E4,E5,H1,E7,E8) IF((K.GT.25).AND.(H1.EQ.-1.E20)) GOTO 99 IF(H1.LT.H) GOTO 10 T2=T T1=T-0.2 ELSE 20 K=K+1 IF(K.GT.25) GOTO 99 T=T-0.2 CALL S1L01(5,P,T,E1,E2,RW,E4,E5,H1,E7,E8) IF(H1.GT.H) GOTO 20 T1=T T2=T+0.2 ENDIF 30 K=K+1 IF(K.GT.50) GOTO 99 T=(T1+T2)/2. CALL S1L01(5,P,T,E1,E2,RW,E4,E5,H1,E7,E8) IF((T2-T1.LT.1.E-4).OR.(ABS(H1-H).LT.1.E-4)) GOTO 50 IF((ABS(H).GE.10.).AND.(ABS(1.-H1/H).LT.1.E-6)) GOTO 50 IF(H1.LT.H) THEN T1=T ELSE T2=T ENDIF GOTO 30 50 G10L01=T RETURN 99 G10L01=-1.E20 RETURN END C--G11 FUNCTION G11L01(P,T,WB) C Mole fraction of water vapor (approximate value) RW[-] TC=T-273.15 WBC=WB-273.15 PS=G3L01(WB) XS=0.622*PS/(P-PS) HA=(1.+1.805*XS)*(TC-WBC) HB=2501.+1.805*TC-4.186*WBC IF(WBC.LE.0.) HB=HB+334.+2.1*WBC X=XS-HA/HB IF(X.LT.0.) X=0. G11L01=X/(0.622+X) RETURN END C********************************************************SUBROUTINE*** C--S1 SUBROUTINE S1L01(N,P,T,RWS,XS,RW,X,V,H,S,FS) C Input: N,P,T Input: RW if N=3-7 C Output: N=1/FS,RWS,XS; 2/1+HS; 3/1+X; 4/3+V; 5/4+H; 6/4+S; 7/all DIMENSION CF(8),CJ(7),CI(7),CA(6),CD(6),CG(6),CK(7) DATA R/8.31441/,AM/28.9645/,WM/18.01528/,EM/0.621978/ DATA CF/-0.2403360201E4, -0.140758895E1, 0.1068287657, & -0.2914492351E-3, 0.373497936E-6, -0.21203787E-9, & -0.3424442728E1, 0.1619785E-1/ DATA CJ/ 0.5088496E2, 0.6163813, 0.1459187E-2, & 0.2008438E-4, -0.5847727E-7, 0.4104110E-9, & 0.1967348E-1/ DATA CI/ 0.50884917E2, 0.62590623, 0.13848668E-2, & 0.21603427E-4, -0.72087667E-7, 0.46545054E-9, & 0.19859983E-1/ DATA HMAD/-0.79141982E4/, HMWD/0.3599417E5/ DATA CA/ 0.63290874E1, 0.28709015E2, 0.26431805E-2, & -0.10405863E-4, 0.18660410E-7, -0.97843331E-11/ DATA CD/-0.5008E-2, 0.32491829E2, 0.65576345E-2, & -0.26442147E-4, 0.51751789E-7, -0.31541624E-10/ DATA SMAD/-0.196125465E3/, SMWD/-0.6331449E2/ DATA CG/ 0.34373874E2, 0.52863609E-2, -0.15608795E-4, & 0.24880547E-7, -0.12230416E-10, 0.28709015E2/ DATA CK/ 0.2196603E1, 0.19743819E-1, -0.70128225E-4, & 0.14866252E-6, -0.14524437E-9, 0.55663583E-13, & 0.32284652E2/ IF(N.GE.3) THEN IF(RW.LT.0.) GOTO 999 RW0=RW ENDIF C Volume-series virial and cross-virial coefficients C WV[m3/mol]: Molar volume of saturated liquid water C WK[1/Pa]: Isothermal compressibility of liquid water C WE[1/Pa]: Henry's law constant for air dissolved in liquid water IF(T.GE.273.16) THEN ROU=CF(1) DO 10 K=2,6 10 ROU=ROU+CF(K)*T**(K-1) ROU=ROU/(CF(7)+CF(8)*T) WV=1.E-3*WM/ROU TC=T-273.15 WK=CJ(1) IF(T.LE.373.15) THEN DO 20 K=2,6 20 WK=WK+CJ(K)*TC**(K-1) WK=WK/(1.+CJ(7)*TC)*1.E-11 ELSE DO 30 K=2,6 30 WK=WK+CI(K)*TC**(K-1) WK=WK/(1.+CI(7)*TC)*1.E-11 ENDIF TAU=1.E3/T TAU2=TAU**2 BO=-(0.0512*TAU+0.1076)/2./0.0005943 CO=(-0.147*TAU2+0.8447*TAU-1.)/0.0005943 EO=10.**(BO+SQRT(BO*BO+CO)) BN=-(0.019*TAU+0.03741)/2./0.1021 CN=(-0.1482*TAU2+0.851*TAU-1.)/0.1021 EN=10.**(BN+SQRT(BN*BN+CN)) WE=(0.22/EO+0.78/EN)*1.E-4/101325. ELSE C WV, WK and WE for ice WV=WM*(0.1070003E-5-0.249936E-10*T+0.371611E-12*T**2) WK=(8.875+0.0165*T)*1.E-11 WE=0. ENDIF C Virial and cross-virial coefficients RT=R*T*1.E6 T2=T**2 T3=T2*T T4=T3*T BAA=0.349568E2-0.668772E4/T-0.210141E7/T2+0.924746E8/T3 CAAA=0.125975E4-0.190905E6/T+0.632467E8/T2 BWW1=0.70E-8-0.147184E-8*EXP(1734.29/T) BWW=RT*BWW1 CWWW1=0.104E-14-0.335297E-17*EXP(3645.09/T) CWWW=RT**2*(CWWW1+BWW1**2) BAW=0.32366097E2-0.141138E5/T-0.1244535E7/T2-0.2348789E10/T4 CAAW=0.482737E3+0.105678E6/T-0.656394E8/T2+0.294442E11/T3 &-0.319317E13/T4 CAWW=-0.10728876E2+0.347802E4/T-0.383383E6/T2+0.33406E8/T3 CAWW=-1.E6*EXP(CAWW) C Mole fraction RWS[-] and Humidity ratio XS[kg/kgDA] at saturation C FS[-]: Enhancement factor PS=G3L01(T) IF(PS.GE.P) GOTO 110 PRT=P/RT PRT2=PRT**2 PSRT=PS/RT PSRT2=PSRT**2 FS=1. DO 100 K=1,2 RWS=FS*PS/P RAS=1.-RWS RAS2=RAS**2 RAS3=RAS2*RAS RAS4=RAS3*RAS RWS2=RWS**2 RWS3=RWS2*RWS FS=((1.+WK*PS)*(PRT-PSRT)-WK*(P*PRT-PS*PSRT)/2.)*WV*1.0E6 FS=FS+ALOG(1.-WE*RAS*P)+RAS2*PRT*BAA-2.*RAS2*PRT*BAW FS=FS-((1.-RAS2)*PRT-PSRT)*BWW+RAS3*PRT2*CAAA FS=FS+3.*RAS2*(RWS-RAS)*PRT2/2.*CAAW-3.*RAS2*RWS*PRT2*CAWW FS=FS-((1.+2.*RAS)*RWS2*PRT2-PSRT2)/2.*CWWW FS=FS-RAS2*(1.-3.*RAS)*RWS*PRT2*BAA*BWW FS=FS-2.*RAS3*(2.-3.*RAS)*PRT2*BAA*BAW FS=FS+6.*RAS2*RWS2*PRT2*BWW*BAW-3.*RAS4*PRT2/2.*BAA*BAA FS=FS-2.*RAS2*RWS*(1.-3.*RAS)*PRT2*BAW*BAW FS=FS-(PSRT2-(1.+3.*RAS)*RWS3*PRT2)/2.*BWW*BWW 100 FS=EXP(FS) IF((FS.LT.1.).OR.(FS*PS.GE.P)) GOTO 110 IF(RWS.GT.0.99) GOTO 110 XS=EM*RWS/RAS IF(N.EQ.1) RETURN GOTO 200 110 FS=-1.E20 RWS=-1.E20 XS=-1.E20 IF(N.LE.2) GOTO 999 C Humidity ratio X[kg/kgDA] 200 IF(N.GE.3) THEN URWS=ABS(RWS) IF(RW.GT.AMIN1(URWS,0.99)) GOTO 999 RA=1.-RW X=EM*RW/RA ENDIF IF(N.EQ.3) RETURN C Specific volume V[m3/kgDA] C VM[cm3/mol]: Molar volume C BM, CM: Virial and cross-virial coefficients for mixture IF(N.EQ.2) THEN RW=RWS RA=RAS ENDIF RTP=RT/P RA2=RA**2 RA3=RA2*RA RW2=RW**2 RW3=RW2*RW RAW=RA*RW BM=RA2*BAA+2.*RAW*BAW+RW2*BWW CM=RA3*CAAA+3.*RAW*(RA*CAAW+RW*CAWW)+RW3*CWWW VM=RTP DO 500 K=1,5 500 VM=RTP*(1.+BM/VM+CM/VM**2) V=VM/(RA*AM)*1.E-3 IF(N.EQ.4) RETURN C Temperature derivatives of virial and cross-virial coefficients T5=T4*T DBAA=0.668772E4/T2+0.420282E7/T3-0.277427E9/T4 DCAAA=0.190905E6/T2-0.126493E9/T3 DBWW1=0.25526E-5/T2*EXP(1734.29/T) DBWW=RT*(DBWW1+BWW1/T) DCWWW1=0.122219E-13/T2*EXP(3645.09/T) DCWWW=RT**2*(DCWWW1+2.*(BWW1*DBWW1+(CWWW1+BWW1**2)/T)) DBAW=0.141138E5/T2+0.248907E7/T3+0.93951568E10/T5 DCAAW=-0.105678E6/T2+0.131279E9/T3-0.883326E11/T4+0.127727E14/T5 DCAWW=CAWW*(-0.347802E4/T2+0.766765E6/T3-0.100218E9/T4) DBM=(RA2*DBAA+2.*RAW*DBAW+RW2*DBWW) DCM=(RA3*DCAAA+3.*RAW*(RA*DCAAW+RW*DCAWW)+RW3*DCWWW) C Specific enthalpy H[kJ/kgDA] C HMA0, HMW0: Ideal-gas (zero pressure) enthalpy IF(N.NE.6) THEN HMA0=CA(1) HMW0=CD(1) DO 600 K=2,6 HMA0=HMA0+CA(K)*T**(K-1) 600 HMW0=HMW0+CD(K)*T**(K-1) H=RA*(HMA0+HMAD)+RW*(HMW0+HMWD)+R*T*((BM-T*DBM)/VM & +(CM-T*DCM/2.)/VM**2) H=H/(AM*RA) IF((N.EQ.2).OR.(N.EQ.5)) RETURN ENDIF C Specific entropy of moist air S[kJ/(kgDA K)] C SMA0, SMW0: Ideal-gas (zero pressure) entropy RP=R*ALOG(0.101325E6) SMA0=RP+CG(1)+CG(6)*ALOG(T) DO 700 K=2,5 700 SMA0=SMA0+CG(K)*T**(K-1) SMW0=RP+CK(1)+CK(7)*ALOG(T) DO 710 K=2,6 710 SMW0=SMW0+CK(K)*T**(K-1) S=RA*(SMA0+SMAD)+RW*(SMW0+SMWD)-R*ALOG(P)+RA*R*ALOG(VM/RTP/RA) &-R*((BM+T*DBM)/VM+(CM+T*DCM)/2./VM**2) IF(RW.GT.0.) S=S+RW*R*ALOG(VM/RTP/RW) S=S/(AM*RA) RETURN 999 FS=-1.E20 RWS=-1.E20 XS=-1.E20 IF(N.GE.3) RW=RW0 X=-1.E20 V=-1.E20 H=-1.E20 S=-1.E20 RETURN END C--S2 SUBROUTINE S2L01(P,T,WB,RW,X,H) C Mole fraction of water vapor RW[-], Humidity ratio X[kg/kgDA], C and Specific enthalpy H[kJ/kgDA] CALL S1L01(2,P,WB,RWS,XS,E3,E4,E5,HS,E7,E8) IF(RWS.LE.0.) GOTO 99 IF(T-WB) 99,100,101 101 RW=0. CALL S1L01(5,P,T,E1,E2,RW,X,E5,H,E7,E8) HC=G7L01(P,WB) HDEF=HS+(X-XS)*HC IF(ABS(HDEF-H).LT.1.E-4) RETURN IF((ABS(H).GE.10.).AND.(ABS(1.-HDEF/H).LT.1.E-6)) RETURN IF(HDEF.LT.H) GOTO 99 C approximate value RW=G11L01(P,T,WB) IF(RW.GT.RWS) RW=RWS IF(RW.LT.0.) RW=0. C exact value CALL S1L01(5,P,T,E1,E2,RW,X,E5,H,E7,E8) HDEF=HS+(X-XS)*HC IF(ABS(HDEF-H).LT.1.E-4) RETURN IF((ABS(H).GE.10.).AND.(ABS(1.-HDEF/H).LT.1.E-6)) RETURN DRW=RWS/1000. K=0 IF(HDEF.GT.H) THEN RW1=RW 10 RW2=RW1+DRW IF(RW2.GE.RWS) THEN RW2=RWS GOTO 30 ENDIF RW=RW2 CALL S1L01(5,P,T,E1,E2,RW,X,E5,H,E7,E8) HDEF=HS+(X-XS)*HC IF(HDEF.GT.H) THEN RW1=RW2 GOTO 10 ENDIF ELSE RW2=RW 20 RW1=RW2-DRW IF(RW1.LE.0.) THEN RW1=0. GOTO 30 ENDIF RW=RW1 CALL S1L01(5,P,T,E1,E2,RW,X,E5,H,E7,E8) HDEF=HS+(X-XS)*HC IF(HDEF.LT.H) THEN RW2=RW1 GOTO 20 ENDIF ENDIF 30 K=K+1 IF(K.GT.40) GOTO 99 RW=(RW1+RW2)/2. CALL S1L01(5,P,T,E1,E2,RW,X,E5,H,E7,E8) HDEF=HS+(X-XS)*HC IF((RW2-RW1.LT.RWS*1.E-6).OR.(ABS(HDEF-H).LT.1.E-4)) RETURN IF((ABS(H).GE.10.).AND.(ABS(1.-HDEF/H).LT.1.E-6)) RETURN IF(HDEF.GT.H) THEN RW1=RW ELSE RW2=RW ENDIF GOTO 30 100 RW=RWS X=XS H=HS RETURN 99 RW=-1.E20 X=-1.E20 H=-1.E20 RETURN END