Category : Science and Education
Archive   : MIPROPS2.ZIP
Filename : HELIUM.FOR

PROGRAM HELIUM
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
COMMON/LIM/TUL,TLL,PUL,TCC,DCC,PCC
1000 FORMAT(//' THIS PROGRAM PROVIDES THE THERMODYNAMIC PROPERTIES
1OF HELIUM',/,' FROM 2 TO 1500 K (-456 TO 2240 F)'
A,/' WITH PRESSURES TO 100 MPA (14503 PSIA)'/)
WRITE(*,1000)
CALL INFO
IP=3
CALL FDATA
TTP=TLL
PTP=.005D0
EM=4.0026
1010 FORMAT(I1)
1 CONTINUE
120 WRITE(*,1040)
1040 FORMAT(' FOR ENGINEERING UNITS ENTER "0", FOR METRIC ENTER "1"',/)
READ(*,1010)IU
WRITE(*,1050)
1050 FORMAT(' FOR SATURATION PROPERTIES ENTER "0", FOR FLUID ENTER '
A '"1"',/)
READ(*,1010)IC
WRITE(*,1060)
1060 FORMAT(' FOR A SINGLE POINT ENTER "0", FOR A TABLE ENTER "1"',/)
READ(*,1010)IV
160 IF(IC.EQ.0)GO TO 240
IF(IV.EQ.1)GO TO 330
170 IF(IU.EQ.0)GO TO 180
WRITE(*,1080)
READ(*,*)PI,D,T
GO TO 190
180 WRITE(*,1070)
1070 FORMAT(' ENTER PRESSURE(PSIA), DENSITY(LB/CU FT), AND TEMPERATURE'
A '(F)',/)
READ(*,*)P,D,T
PI=(P/14.695949D0)*.101325
D=D*16.01846371D0/EM
T=(T-32.D0)/1.8D0+273.15D0
190 IF(PI.LE.0.0D0.AND.D.LE.0.0D0)GO TO 110
IF(PI.GT.0.0D0.AND.D.GT.0.0D0)GO TO 200
IF(PI.LE.0.0D0.AND.T.NE.0.0D0)GO TO 220
IF(PI.NE.0.0D0.AND.T.NE.0.0D0)GO TO 210
GO TO 1
1080 FORMAT(' ENTER PRESSURE(MPA), DENSITY(MOL/L), AND TEMPERATURE(K)'
A /)
200 P=PI/.101325
T=FIND T(P,D)
CALL LIMITS(PI,T,IL)
IF(IL.LE.0)GO TO 170
GO TO 230
210 P=PI/.101325
CALL LIMITS(PI,T,IL)
IF(IL.LE.0)GO TO 170
D=FIND D(P,T)
GO TO 230
220 P=FIND P(D,T)
PI=P*.101325
CALL LIMITS(PI,T,IL)
IF(IL.LE.0)GO TO 170
230 CALL REPRO(PI,D,T,IU,IV,IC,IP,TF,TDEL)
GO TO 170
240 WRITE(*,1090)
1090 FORMAT(' FOR SATURATED LIQUID ENTER "0", FOR VAPOR ENTER "1"',/)
READ(*,1010)IP
IF(IV.EQ.1)GO TO 330
WRITE(*,1095)
1095 FORMAT(' TO ENTER WITH TEMPERATURE ENTER "0", FOR PRESSURE "1"',/)
READ(*,1010)II
IF(II.EQ.1)GO TO 290
250 IF(IU.EQ.1)GO TO 260
1100 FORMAT(' ENTER A TEMPERATURE IN DEGREES F, ENTER 0 FOR RESTART',/)
WRITE(*,1100)
READ(*,*)TI
IF(TI.EQ.0.0D0)GO TO 110
T=(TI-32.D0)/1.8D0+273.15D0
GO TO 270
260 WRITE(*,1110)
1110 FORMAT(' ENTER A TEMPERATURE(K) ,ENTER 0 FOR RESTART',/)
READ(*,*)T
270 IF(T.LT..000001D0)GO TO 110
IF(T.GT.TCC.OR.T.LT.TTP)GO TO 280
P=VPN(T)
IF(IP.EQ.0)P=P+.00001D0
IF(IP.EQ.1)P=P-.00001D0
D=FIND D(P,T)
PI=P*.101325
CALL RE PRO(PI,D,T,IU,IV,IC,IP,TF,TDEL)
GO TO 250
280 X1=(TTP-273.15D0)*1.8D0+32.D0
X2=(TCC-273.15D0)*1.8D0+32.D0
WRITE(*,1120)TTP,TCC,X1,X2
1120 FORMAT(' FOR SATURATION ',F6.2,' < TEMP < ',F6.2,' K',/,
A ' OR ',F7.2,' < TEMP < ',F7.2,' F',/)
GO TO 250
290 IF(IU.EQ.1)GO TO 300
WRITE(*,1130)
1130 FORMAT(' ENTER A PRESSURE IN LB/SQ IN 0 RESTRATS PROGRAM',/)
READ(*,*)P
PI=(P/14.695949D0)*.101325D0
P=PI/.101325D0
IF(PI.LE.0.0D0)GO TO 110
GO TO 310
300 WRITE(*,1140)
1140 FORMAT(' ENTER A PRESSURE(MPA) 0 RESTARTS PROGRAM',/)
READ(*,*)PI
IF(PI.LE.0.0D0)GO TO 110
P=PI/.101325
310 IF(PI.GT.PCC.OR.PI.LT.PTP)GO TO 320
T=FIND TV(P)
P=VPN(T)
IF(IP.EQ.1)P=P-.00001D0
IF(IP.EQ.0)P=P+.00001D0
D=FIND D(P,T)
PI=P*.101325
CALL RE PRO(PI,D,T,IU,IV,IC,IP,TF,TDEL)
GO TO 290
320 PTPF=PTP*14.695949/.101325
PCCF=PCC*14.695949/.101325
WRITE(*,1150) PTP,PCC,PTPF,PCCF
1150 FORMAT(' YOUR INPUT PRESSURE IS OUTSIDE THE RANGE OF SATURATION'
A ' PRESSURES'/' FOR THIS FLUID THE RANGE IS ',F6.5,' TO ',F6.3,
B' MPA,'/' OR ',F6.5,' TO ',F6.3,' PSIA'/' TRY AGAIN',/)
GO TO 290
330 IF(IC.EQ.1)GO TO 370
IF(IU.EQ.1)GO TO 340
WRITE(*,1160)
1160 FORMAT(' ENTER A STARTING TEMPERATURE, A FINAL TEMPERATURE'
A /' AND A TEMPERATURE INCREMENT, IN DEGREES F AND IN THAT ORDER'/
B ' ENTER 0,0,0 TO RESTART'/)
READ(*,*)TS,TF,TDEL
IF(TDEL.LE.0.0D0)GO TO 110
TS=(TS-32.D0)/1.8D0+273.15D0
TF=(TF-32.D0)/1.8D0+273.15D0
TDEL=TDEL/1.8D0
IF(TS.LT.TTP.OR.TS.GT.TCC)GO TO 360
IF(TF.LT.TTP.OR.TF.GT.TCC)GO TO 360
GO TO 350
340 WRITE(*,1170)
1170 FORMAT(' ENTER A STARTING TEMPERATURE, A FINAL TEMPERATURE'
1/' AND A TEMPERATURE INCREMENT IN KELVINS AND IN THAT ORDER',/
2' ENTER 0,0,0 TO RESTART PROGRAM'/)
READ(*,*)TS,TF,TDEL
IF(TDEL.LE.0.0D0)GO TO 110
IF(TS.LT.TTP.OR.TS.GT.TCC)GO TO 360
IF(TF.LT.TTP.OR.TF.GT.TCC)GO TO 360
350 T=TS
P=VPN(T)
IF(IP.EQ.0)P=P+.00001D0
IF(IP.EQ.1)P=P-.00001D0
D=FIND D(P,T)
PI=P*.101325
CALL RE PRO(PI,D,T,IU,IV,IC,IP,TF,TDEL)
GO TO 330
360 X1=(TTP-273.15D0)*1.8D0+32.D0
X2=(TCC-273.15D0)*1.8D0+32.D0
WRITE(*,1180)TTP,TCC,X1,X2
1180 FORMAT(' FOR SATURATION, ',F6.2,' < TEMP < ',F6.2,' K',/,
A ,13X,'OR, ',F7.1,' < TEMP < ',F7.1,' F. TRY AGAIN.',/)
GO TO 330
370 IF(IU.EQ.1)GO TO 380
WRITE(*,1190)
1190 FORMAT(' ENTER PRESSURE(PSIA), STARTING TEMPERATURE(F), FINAL '
A 'TEMPERATURE(F)'/' AND A TEMPERATURE INCREMENT, IN THAT ORDER',/
B ' ENTER 0,0,0,0 TO RESTART PROGRAM'/)
READ(*,*)P,TS,TF,TDEL
IF(TDEL.LE.0.0D0)GO TO 110
PI=(P/14.695949D0)*.101325D0
T=(TS-32.D0)/1.8D0+273.15D0
TF=(TF-32.D0)/1.8D0+273.15D0
TDEL=TDEL/1.8D0
P=PI/.101325
CALL LIMITS(PI,T,IL)
IF(IL.LE.0)GO TO 370
CALL LIMITS(PI,TF,IL)
IF(IL.LE.0)GO TO 370
GO TO 390
380 WRITE(*,1200)
1200 FORMAT(' ENTER A PRESSURE(MPA), A STARTING TEMPERATURE(K), A '
1'FINAL TEMPERATURE AND A'/' TEMPERATURE INCREMENT'
2' AND IN THAT ORDER, TO RESTART PROGRAM ENTER "0,0,0,0"',/)
READ(*,*)PI,TS,TF,TDEL
IF(TDEL.LE.0.0D0)GO TO 110
T=TS
P=PI/.101325
CALL LIMITS(PI,T,IL)
IF(IL.LE.0)GO TO 370
CALL LIMITS(PI,TF,IL)
IF(IL.LE.0)GO TO 370
390 D=FIND D(P,T)
CALL RE PRO(PI,D,T,IU,IV,IC,IP,TF,TDEL)
GO TO 370
110 CONTINUE
WRITE(*,1210)
1210 FORMAT(//' FOR MORE PROPERTIES ENTER 0, TO TERMINATE ENTER 1')
READ(*,1010)IM
IF(IM.EQ.0)GO TO 1
999 CONTINUE
STOP
END
SUBROUTINE REPRO(PI,D,T,IU,IV,IC,IP,TF,TDEL)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
N=500
EM=4.0026
P=PI/.101325
IF(IV.EQ.0)TF=T-1.D0
IF(IU.EQ.0)GO TO 100
WRITE(*,1000)
WRITE(*,1010)
WRITE(*,1020)
GO TO 110
100 WRITE(*,1000)
WRITE(*,1030)
WRITE(*,1040)
110 CONTINUE
DO 210 I=1,N
IF(I.EQ.1)GO TO 120
D=FIND D(P,T)
120 H=ENTHAL(P,D,T)
E=H-101.325D0*P/D
S=ENTROP(D,T)
W=SOUND(D,T)
CPP=CP(D,T)
CVV=CV(D,T)
150 V=VISC(D,T)
TH=THERM(D,T)*100.D0
EPS=FDIEL(D)
IF(IU.EQ.0)GO TO 160
PO=P*.101325
WRITE(*,2030) PO,T,D,E,H,S,CVV,CPP,W,V,TH,EPS
GO TO 200
160 H=H/(2.324445D0*EM)
E=E/(2.324445D0*EM)
S=S/(4.184001D0*EM)
CPP=CPP/(4.184001D0*EM)
CVV=CVV/(4.184001D0*EM)
W=W*3.280840D0
PO=P*14.695949D0
DO=D*EM/16.01846371D0
TO=T*1.8D0-459.67D0
190 V=VISC(D,T)*.0067196897D0
TH=THERM(D,T)*.000578176D0
EPS=FDIEL(P,D,T)
WRITE(*,3030) PO,TO,DO,E,H,S,CVV,CPP,W,V,TH,EPS
200 T=T+TDEL
IF(T.GT.TF+.01D0)GO TO 220
IF(IC.NE.0)GO TO 210
P=VPN(T)
IF(IP.EQ.0)P=P+.00001D0
IF(IP.EQ.1)P=P-.00001D0
210 CONTINUE
220 CONTINUE
WRITE(*,1000)
1000 FORMAT(' ')
RETURN
1010 FORMAT(3X,'P',5X,'T',6X,'DEN',5X,'E',6X,'H',6X,'S',4X,'CV',4X,
A 'CP',3X,'SOUND',2X,'VISC',2X,'COND',2X,'DIEL')
1020 FORMAT(3X,'MPA',3X,'K',6X,'MOL/L',2X,'=== J/MOL ==',2X,'==== ',
A 'J/MOL-K ===',2X,'M/S',3X,'PA-S',1X,'MW/M-K',2X,'==',/,61X,'E+6')
1030 FORMAT(3X,'P',5X,'T',7X,'DENS',3X,'E',8X,'H',5X,'S',5X,'CV',
A 4X,'CP',2X,'SOUND',1X,'VISC',2X,'COND',2X,'DIEL')
1040 FORMAT(3X,'PSIA',2X,'F',7X,'LB/',3X,'==== BTU/ ===',1X,'======',
A ' BTU/ =====',2X,'F/S',3X,'LB/',2X,'BTU/',2X,'==',/,17X,'CU FT',
B 6X,'LB',14X,'LB-F',13X,'FT-S',1X,'FT-HR-F',/,62X,'E+5')
2030 FORMAT(F8.3,F7.2,F7.3,2F7.0,F6.1,F5.1,F6.1,F6.0,F6.1,F6.1,F8.5)
3030 FORMAT(F8.1,F7.2,F7.3,2F7.1,3F6.3,F6.0,F5.2,F6.4,F8.5)
END
SUBROUTINE INFO
WRITE(*,101)
101 FORMAT(' WRITTEN BY'//
A' ROBERT D. MCCARTY'/
1' THERMOPHYSICS DIVISION'/
2' CENTER FOR CHEMICAL ENGINEERING'/
3' NATIONAL BUREAU OF STANDARDS'/
4' BOULDER, COLORADO'//
5' DISTRIBUTED BY'//
A' THE OFFICE OF STANDARD REFERENCE DATA'/
6' NATIONAL BUREAU OF STANDARDS, WASHINGTON, DC'/)
WRITE(*,100)
100 FORMAT(/ ' WHEN THE PROGRAM ASKS FOR A PRESSURE, DENSITY AND'
*' TEMPERATURE,'/ ' ENTER ANY TWO AND A ZERO FOR THE THIRD.'/
* ' TO TERMINATE THE PROGRAM ENTER ZERO FOR ALL THREE.'/)
RETURN
END
SUBROUTINE FDATA(IF)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION B(27,3)
COMMON/COEF/B
COMMON/LIM/TUL,TLL,PUL,TCC,DCC,PCC
OPEN(5,FILE='HELIUM.COF')
DO 60 J=1,3
DO 60 I=1,27
60 READ(5,101)B(I,J)
101 FORMAT(D20.13)
READ(5,102)TUL,TLL,PUL,TCC,DCC,PCC
CLOSE(5,STATUS='KEEP')
102 FORMAT(F10.8,E14.8,3F8.2,2F8.4)
RETURN
END
SUBROUTINE LIMITS(PI,T,IL)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
COMMON/LIM/TUL,TLL,PUL,TCC,DCC,PCC
P=PI/.101325
IF(PI.GT. PUL)GO TO 10
IF(T .GT. TUL.OR.T .LT. TLL)GO TO 12
PM=PMELT(T)
IF(P .GT. PM) GO TO 20
IL=1
RETURN
10 PULF=PUL*14.6959496/.101325
WRITE(*,11)PUL,PULF
11 FORMAT(' THE INPUT PRESSURE IS OUT OF THE RANGE OF THIS EQUATION '
1/' THE RANGE FOR THIS EQUATION IS FROM 0 TO ',F6.0,' MPA'
2/' OR ',F7.0,' PSIA')
IL=0
RETURN
12 TLLF= (TLL-273.15D0)*1.8D0+32.D0
TULF= (TUL-273.15D0)*1.8D0+32.D0
WRITE(*,13) TLL,TUL,TLLF,TULF
13 FORMAT(' THE INPUT TEMPERATURE IS OUT OF RANGE'
A /' THE RANGE FOR THIS EQUATION IS ',F6.2,' K TO ',F6.0,' K',/,
B 27X,' OR ',F8.2,' F TO ',F6.0,' F')
IL=0
RETURN
20 TM=TMELT(P)
TF=(TM-273.15D0)*1.8D0+32.D0
WRITE(*,22) TM,TF
22 FORMAT(' SOLID PHASE DETECTED.',/,' FOR THIS PRESSURE, TEMP'
A ' SHOULD EXCEED ',F8.3,' K, OR',F9.3,' F')
IL=0
END
DOUBLE PRECISION FUNCTION CV(D,T)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
C CALCULATES CV(J/(MOL*K)). INPUT DENS(MOL/L) AND TEMP(K).
CALL PROPS(PP,D,T,6)
CV=PP
END
DOUBLE PRECISION FUNCTION ENTHAL(P,D,T)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
DD=D
TT=T
CALL PROPS(UD,DD,TT,5)
ENTHAL=(UD)
END
DOUBLE PRECISION FUNCTION ENTROP(D,T)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
C CALCULATES ENTROPY(J/(MOL-K), FROM INPUT OF DENSITY(MOL/L) AND TEMP(K).
DD=D
TT=T
CALL PROPS(SD,DD,TT,4)
ENTROP=(SD)
END
DOUBLE PRECISION FUNCTION FINDP(D,T)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION G(32),VP(9)
COMMON/DATA/G,R,GAMMA,VP,DTP,PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC
DD=D
TT=T
IF(TT.LT.TCC)GO TO 10
1 CALL PROPS(PP,DD,TT,1)
FINDP=PP
RETURN
10 P=VPN(TT)
DV=FINDD(P-.0001D0,TT)
DL=FINDD(P+.0001D0,TT)
IF(DD.LE.DV.OR.DD.GE.DL)GO TO 1
WRITE(*,100)DV,DL,DD
CALL PROPS(PP,DV,TT,1)
FINDP=PP
D=DV
100 FORMAT(' THE STATE POINT YOU HAVE SPECIFIED CORRESPONDS TO A '
1/' DENSITY IN THE LIQUID VAPOR COEXISTENCE REGION'
2/' THE DENSITY OF THE SATURATED VAPOR IS ',F6.4,' MOLES/LITER'
3/' THE DENSITY OF THE SATURATED LIQUID IS ',F8.4,' MOLES/LITER'
4/' AND THE INPUT DENSITY IS ',F8.4,' MOLES/LITER'
5/' SATURATED VAPOR IS ASSUMED')
END
DOUBLE PRECISION FUNCTION FINDT(P,D)
IMPLICIT REAL*8(A-H)
IMPLICIT REAL*8(O-Z)
IMPLICIT INTEGER*4(I-N)
C RETURNS TEMPERATURE(K), FROM THE 32-TERM MBWR EQN OF STATE.
C INPUT IS PRESSURE(MPA) AND DENSITY(MOL/L).
DIMENSION G(32),VP(9)
COMMON/DATA/G,R,GAMMA,VP,DTP,PCC,PTP,TCC,TTP,TUL,TLL,PUL,DCC
PP=P
DD=D
IF(P.GE.PCC)GO TO 1
TSAT=FINDTV(PP)
DV=FINDD(PP-.00001D0,TSAT)
DL=FINDD(PP+.0001D0,TSAT)
IF(DD.GT.DV.AND.DD.LT.DL)GO TO 30
TT=TSAT
GO TO 2
1 TT=TCC
2 DO 10 I=1,10
CALL PROPS(P2,DD,TT,1)
IF(DABS(PP-P2)-1.D-7*PP)20,20,11
11 CALL PROPS(DP,DD,TT,3)
CORR=(P2-PP)/DP
IF(DABS(CORR)-1.D-5)20,20,10
10 TT=TT-CORR
20 FINDT=TT
RETURN
30 FINDT=TSAT
D=DV
WRITE(*,100)DV,DL,DD
100 FORMAT(' THE STATE POINT YOU HAVE SPECIFIED CORRESPONDS TO'
1/' A DENSITY IN THE LIQUID VAPOR COEXISTENCE REGION'
2/' DENSITY OF THE SATURATED VAPOR IS',F8.4,' MOLES/LITER'
3/' DENSITY OF THE SATURATED LIQUID IS',F8.4,' MOLES/LITER'
4/' INPUT DENSITY IS',F8.4,' MOLES/LITER'
5/' SATURATED VAPOR CONDITIONS ARE ASSUMED')
END
DOUBLE PRECISION FUNCTION FIND D(PI,TI)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
T=TI
P=PI
IF(T.LT.5.2D0)GO TO 6
PM=1001.D0
PM=1.0D+5
IF(T.LT.15.D0)PM=PMELT(T)
IF(PM.LT.P)GO TO 30
IF(T.GT.100.D0)GO TO 1
PC=2.2449D0+(T-5.2014D0)*1.76D0
IF(P.LT.PC)GO TO 1
PM=200.D0+(T-5.2D0)*12.31D0
D=17.399D0+((PM-PC)/(PM-PC+1.D0))*2.33D0*17.399D0
GO TO 7
2 D=.0001D0
IF(T.LT.4.2D0)GO TO 7
1 VB=VIRB(T)
RT=0.0820558D0*T
P1=RT/P
D=1./(P1+VB)
GO TO 7
6 IF(P.LT.VPN(T))GO TO 2
DS=46.18D0+(T-2.D0)*4.02D0
DL=SATL(T)*1000.D0/4.0026D0
DEL=DS-DL
PM=PMELT(T)
IF(P.GT.PM)GO TO 30
D=DL+DEL*P/PM
7 DO 10 I=1,50
CALL PROPS(P2,D,T,1)
IF(DABS(P-P2)-1.D-7*P)20,20,8
8 CALL PROPS(DP,D,T,2)
CORR=(P2-P)/DP
IF(DABS(CORR)-1.D-7*D)20,20,10
10 D=D-CORR
FIND D=0.D0
RETURN
20 FIND D=D
RETURN
30 FIND D=0.D0
RETURN
END
SUBROUTINE PROPS(PP,DD,TT,ICON)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION A(26),B(27,3)
COMMON/COEF/B
GO TO(1000,3000,2000,4000,5000,6000)ICON
1000 KP=1
GO TO 10
C ENTRY DPDT
2000 KP=2
GO TO 10
C ENTRY DPDD
3000 KP=3
GO TO 10
C ENTRY ENTROP
4000 KP=4
GO TO 10
C ENTRY ENTHAL
5000 KP=5
GO TO 10
C ENTRY CV
6000 KP=6
10 K1=0
KH=1
K=1
KK=1
K4=1
IF(ID.NE.0)GO TO 20
IF(TT.GE.15.D0)GO TO 20
IF(TT.GT.10.D0)GO TO 30
IF(DD.GT.17.3987D0)GO TO 40
11 I=1
T=TT
D=DD
8 GO TO (9,100,200,300,400,500)KP
9 D2=D*D
D3=D2*D
D4=D3*D
D5=D4*D
GAMMA=B(27,I)
EX=DEXP(D2*GAMMA)
EXD3=EX*D3
EXD5=EX*D5
M=I
N=1
A(N)=D5*D
N=N+1
A(N)=A(N-1)/T
N=N+1
DO 2 I=1,6
FI=I
A(N)=D5*T**(.75D0-FI/4.D0)
2 N=N+1
DO 3 I=1,4
FI=I
A(N)=D4*T**(1.5D0-FI)
3 N=N+1
DO 1 I=1,8
FI=I
A(N)=D3*T**(1.5D0-FI/2.D0)
N=N+1
1 CONTINUE
DO 4 I=1,3
FI=I
A(N)=EXD3*T**(1.D0-FI)
4 N=N+1
DO 5 I=1,3
FI=I
A(N)=EXD5*T**(1.D0-FI)
5 N=N+1
N=N-1
I=M
7 P=0.D0
DO 15 J=1,N
15 P=P+B(J,I)*A(J)
P=P+.0820558D0*D*T*(1.+VIRB(T)*D)
IF(KH.LT.1)GO TO 413
GO TO(50,50,30,40)K
100 D2=D*D
D3=D**3
D4=D3*D
D5=D4*D
D6=D5*D
T2=T*T
T3=T2*T
T4=T**4
M=I
GAMMA=B(27,M)
EX=DEXP(D2*GAMMA)
N=1
R=.0820558D0
A(N)=0.0D0
N=N+1
A(N)=(-1.D0)*D6/T2
N=N+1
DO 102 I=1,6
FI=I
A(N)=D5*T**(.75D0-FI/4.D0-1.D0)*(.75D0-FI/4.D0)
102 N=N+1
DO 103 I=1,4
FI=I
A(N)=D4*T**(1.5D0-FI/1.D0-1.D0)*(1.5D0-FI)
103 N=N+1
DO 101 I=1,8
FI=I
A(N)=D3*T**(1.5D0-FI/2.D0-1.D0)*(1.5D0-FI/2.D0)
101 N=N+1
DO 104 I=1,3
FI=I
A(N)=EX*D3*T**(1.D0-FI-1.D0)*(1.D0-FI)
104 N=N+1
DO 105 I=1,3
FI=I
A(N)=EX*D5*T**(1.D0-FI-1.D0)*(1.D0-FI)
105 N=N+1
N=N-1
P=0
DO 115 J=1,N
115 P=P+A(J)*B(J,M)
P=P+R*D*(1.D0+VIRB(T)*D)+R*D*T*DBDT(T)*D
I=M
GO TO(50,50,30,40)K
200 D2=D*D
D3=D2*D
D4=D3*D
D5=D4*D
M=I
GAMMA=B(27,M)
EX=DEXP(D2*GAMMA)
DEX=GAMMA*2.D0*D*EX
N=1
R=0.0820558D0
A(N)=6.D0*D5
N=N+1
A(N)=A(N-1)/T
N=N+1
DO 202 I=1,6
FI=I
A(N)=5.D0*D4*T**(.75D0-FI/4.D0)
202 N=N+1
DO 203 I=1,4
FI=I
A(N)=D3*4.D0*T**(1.5D0-FI)
203 N=N+1
DO 201 I=1,8
FI=I
A(N)=D2*3.D0*T**(1.5D0-FI/2.D0)
201 N=N+1
DO 204 I=1,3
FI=I
A(N)=(DEX*D3+3.D0*D2*EX)*T**(1.D0-FI)
204 N=N+1
DO 205 I=1,3
FI=I
A(N)=(DEX*D5+5.D0*D4*EX)*T**(1.D0-FI)
205 N=N+1
N=N-1
P=0
DO 215 J=1,N
215 P=P+A(J)*B(J,M)
I=M
P=P+R*T*(1.D0+2.D0*D*VIRB(T))
GO TO(50,50,30,40)K
300 D2=D*D
D3=D2*D
D4=D3*D
N=1
R=.0820558D0
M=I
GAMMA=B(27,M)
EX=DEXP(D2*GAMMA)
A(N)=0.0D0
N=N+1
A(N)=(D4*D/5.D0)*T**(-2.D0)*(-1.D0)
N=N+1
DO 302 I=1,6
FI=I
A(N)=(D4/4.D0)*T**(.75D0-FI/4.D0-1.D0)*(.75D0-FI/4.D0)
302 N=N+1
DO 303 I=1,4
FI=I
A(N)=(D3/3.D0)* T**(1.5D0-FI-1.D0)*(1.5D0-FI)
303 N=N+1
DO 301 I=1,8
FI=I
A(N)=(D2/2.D0)*T**(1.5D0-FI/2.D0-1.D0)*(1.5D0-FI/2.D0)
301 N=N+1
DO 304 I=1,3
FI=I
A(N)=(EX/(2.D0*GAMMA))*T**(1.D0-FI-1.D0)*(1.D0-FI)
304 N=N+1
DO 305 I=1,3
FI=I
A(N)=(D2*EX/(2.D0*GAMMA)-EX/(2.D0*GAMMA**2))*
1T**(1.D0-FI-1.D0)*(1.D0-FI)
305 N=N+1
N=N-1
SINT=D*R*(VIRB(T)+T*DBDT(T))
DO 306 I=1,N
306 SINT=SINT+B(I,M)*A(I)
N=21
EX=1.D0
D2=0.0D0
DO 310 I=1,3
FI=I
A(N)=(EX/(2.D0*GAMMA))*T**(1.D0-FI-1.D0)*(1.D0-FI)
310 N=N+1
DO 311 I=1,3
FI=I
A(N)=(D2*EX/(2.D0*GAMMA)-EX/(2.D0*GAMMA**2))*
1T**(1.D0-FI-1.D0)*(1.D0-FI)
311 N=N+1
N=N-1
DO 312 I=21,N
312 SINT=SINT-B(I,M)*A(I)
P=(9.371658D0+5.193043D0*DLOG(T/4.2144D0)-25.31469D0
1*(SINT+R*DLOG(R*T*D)))
P=P*4.0026D0
I=M
GO TO(50,50,30,40)K
400 KH=0
GO TO 9
413 PP=P
KH=1
D2=D*D
D3=D*D2
D4=D3*D
N=1
R=.0820558D0
M=I
GAMMA=B(27,M)
EX=DEXP(D2*GAMMA)
A(N)=(D4*D)/5.D0
N=N+1
A(N)=(D4*D/5.D0)*(2.D0/T)
N=N+1
DO 402 I=1,6
FI=I
A(N)=(D4/4.D0)*(T**(.75D0-FI/4.D0)-T**(.75D0-FI/4.D0)
1*(.75D0-FI/4.D0))
402 N=N+1
DO 403 I=1,4
FI=I
A(N)=(D3/3.D0)*(T**(1.5D0-FI)-T**(1.5D0-FI)*(1.5D0-FI))
403 N=N+1
DO 401 I=1,8
FI=I
A(N)=(D2/2.D0)*(T**(1.5D0-FI/2.D0)-T**(1.5D0-FI/2.D0)
1*(1.5D0-FI/2.D0))
401 N=N+1
DO 404 I=1,3
FI=I
A(N)=(EX/(2.D0*GAMMA))*(T**(1.D0-FI)-T**(1.D0-FI)*(1.D0-FI))
404 N=N+1
DO 405 I=1,3
FI=I
A(N)=(D2*EX/(2.D0*GAMMA)-EX/(2.D0*GAMMA**2))
1*(T**(1.D0-FI)-T**(1.D0-FI)*(1.D0-FI))
405 N=N+1
N=N-1
HINT=R*T*T*D*DBDT(T) *(-1.D0)
DO 406 I=1,N
406 HINT=HINT+B(I,M)*A(I)
N=21
D2=0.0D0
EX=1.D0
DO 410 I=1,3
FI=I
A(N)=(EX/(2.D0*GAMMA))*(T**(1.D0-FI)-T**(1.D0-FI)*(1.D0-FI))
410 N=N+1
DO 411 I=1,3
FI=I
A(N)=(D2*EX/(2.D0*GAMMA)-EX/(2.D0*GAMMA**2))*(T**(1.D0-FI)
1-T**(1.D0-FI)*(1.-FI))
411 N=N+1
N=N-1
DO 412 I=21,N
412 HINT=HINT-B(I,M)*A(I)
P=21.8228D0+5.193043D0*(T-4.2144D0)+25.31469D0
1*(HINT+PP/D-R*T)
P=P*4.0026D0
I=M
GO TO(50,50,30,40)K
500 D2=D*D
D3=D2*D
D4=D3*D
N=1
R=.0820558D0
M=I
GAMMA=B(27,M)
EX=DEXP(D2*GAMMA)
SINT=T*D*R*(2.*DBDT(T )+D2DBDT2(T)*T)
A(N)=0.0D0
N=N+1
A(N)=(D4*D/5.D0)*T**(-2.D0)*(-1.D0)*(-2.D0)
N=N+1
DO 502 I=1,6
FI=I
A(N)=(D4/4.D0)*T**(.75D0-FI/4.D0-1.D0)
1*(.75D0-FI/4.D0)*(.75D0-FI/4.D0-1.D0)
502 N=N+1
DO 503 I=1,4
FI=I
A(N)=(D3/3.D0)* T**(1.5D0-FI-1.D0)*(1.5D0-FI)*(1.5D0-FI-1.D0)
503 N=N+1
DO 501 I=1,8
FI=I
A(N)=(D2/2.D0)*T**(1.5D0-FI/2.D0-1.D0)*(1.5D0-FI/2.D0)
1*(1.5D0-FI/2.D0-1.D0)
501 N=N+1
DO 504 I=1,3
FI=I
A(N)=(EX/(2.D0*GAMMA))*T**(1.D0-FI-1.D0)*(1.D0-FI)
1*(1.D0-FI-1.D0)
504 N=N+1
DO 505 I=1,3
FI=I
A(N)=(D2*EX/(2.D0*GAMMA)-EX/(2.D0*GAMMA**2))*
1T**(1.D0-FI-1.D0)*(1.D0-FI)*(1.D0-FI-1.D0)
505 N=N+1
N=N-1
DO 506 I=1,N
506 SINT=SINT+B(I,M)*A(I)
P=SINT
EX=1.D0
D2=0
N=21
DO 510 I=1,3
FI=I
A(N)=(EX/(2.D0*GAMMA))*T**(1.D0-FI-1.D0)*(1.D0-FI)
1*(1.D0-FI-1.D0)
510 N=N+1
DO 511 I=1,3
FI=I
A(N)=(D2*EX/(2.D0*GAMMA)-EX/(2.D0*GAMMA**2))*
1T**(1.D0-FI-1.D0)*(1.D0-FI)*(1.D0-FI-1.D0)
511 N=N+1
N=N-1
DO 512 I=21,N
512 P=P-B(I,M)*A(I)
P=5.193043D0*(3.D0/5.D0)*4.0026D0-P*101.3278D0
I=M
GO TO(50,50,30,40)K
20 K=2
I=3
T=TT
D=DD
GO TO 8
30 K=3
GO TO(33,34,35)KK
33 D=DD
T=TT
KK=2
IF(DD.GT.17.3987D0)GO TO 40
I=1
GO TO 8
34 PTI=P
I=3
KK=3
T=TT
D=DD
GO TO 8
35 PTIII=P
38 F=(15.D0-T)/5.D0
P=F*PTI+(1.D0-F)*PTIII
IF(KH.LT.1)GO TO 413
PP=P
RETURN
40 IF(K.EQ.3)K1=3
GO TO(41,42,43,44)K4
41 K=4
I=2
K4=2
D=DD
IF(K1.EQ.0)T=TT
GO TO 8
42 PIID=P
D=17.3987D0
IF(T.LT.5.2014D0)D=SATL(T)*1000.D0/4.0026D0
K4=3
GO TO 8
43 PIIDC=P
I=1
K4=4
GO TO 8
44 PIDC=P
P=PIDC+(PIID-PIIDC)
K4=1
IF(K1.EQ.3)GO TO 30
PP=P
RETURN
50 PP=P
RETURN
END
DOUBLE PRECISION FUNCTION CP(D,T)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
CALL PROPS(DPDD,D,T,2)
CALL PROPS(DPDT,D,T,3)
CP=CV(D,T)+(T*(DPDT**2)/((D**2)*DPDD))*101.3278D0
RETURN
END
DOUBLE PRECISIONFUNCTION SOUND(D,T)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
CALL PROPS(DPDD,D,T,2)
SOUND=((CP(D,T)/CV(D,T))*(DPDD*25311.D0))**.5D0
RETURN
END
DOUBLE PRECISION FUNCTION P MELT(TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION A(5)
DATA A/33.28D0,-44.156D0,31.799D0,-4.8159D0,.30313D0/
T=TT
IF(T.LE.5.2D0)GO TO 7
PMELT=-17.80D0+17.31457D0*T**1.555414D0
PMELT=PMELT*.98066D0/1.01325D0
RETURN
7 P=0.0D0
DO 9 I=1,5
9 P=P+A(I)*T**(I-1)
PMELT=P*9.80665D0/10.1325D0
RETURN
END
DOUBLE PRECISION FUNCTION VIRB(T)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION A(9)
DATA A/-5.0815710041D-7,-1.1168680862D-4,1.1652480354D-2,
1 7.4474587998D-2,-5.3143174768D-1,-9.5759219306D-1,
2 3.9374414843D0,-5.1370239224D0,2.0804456338D0/
C COEFFICIENTS FROM PROGRAM 5/28/70-1630
C THIS SUB PROGRAM CALCULATES THE SECOND VIRIAL COEFFICIENT FOR
C HELIUM. THE RANGE IS FROM 2 TO 1500 DEG K. INPUT IS TEMPERATURE
C IN DEGREES KELVIN
C UNITS ARE ATM, DEG KELVIN, AND MOLES/LITER,4/3/69-1253,R.D.MCCARTY
C REVISED 2/12/70-925
1 B=0.0D0
DO 5 I=1,9
FI=I
5 B=B+T**(1.5D0-FI/2.D0)*A(I)
VIRB=B
RETURN
END
DOUBLE PRECISION FUNCTION DBDT (T)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION A(9),V(45)
C THIS SUB PROGRAM CALCULATES THE SECOND VIRIAL COEFFICIENT FOR
C HELIUM. THE RANGE IS FROM 2 TO 1500 DEG K. INPUT IS TEMPERATURE
C IN DEGREES KELVIN
DATA A/-5.0815710041D-7,-1.1168680862D-4,1.1652480354D-2,
1 7.4474587998D-2,-5.3143174768D-1,-9.5759219306D-1,
2 3.9374414843D0,-5.1370239224D0,2.0804456338D0/
C UNITS ARE ATM, DEG KELVIN, AND MOLES/LITER,5/28/70-1630,R.D.MCCARTY
1 B=0.0D0
DO 5 I=1,9
FI=I
5 B=B+T**(.5D0-FI/2.D0)*A(I)*(1.5D0-FI/2.D0)
DBDT=B
RETURN
END
DOUBLE PRECISION FUNCTION D2DBDT2(T)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION A(9)
DATA A/-5.0815710041D-7,-1.1168680862D-4,1.1652480354D-2,
1 7.4474587998D-2,-5.3143174768D-1,-9.5759219306D-1,
2 3.9374414843D0,-5.1370239224D0,2.0804456338D0/
C THIS SUB PROGRAM CALCULATES THE SECOND VIRIAL COEFFICIENT FOR
C HELIUM. THE RANGE IS FROM 2 TO 1500 DEG K. INPUT IS TEMPERATURE
C IN DEGREES KELVIN
1 B=0.0D0
DO 5 I=1,9
FI=I
5 B=B+T**(.5D0-FI/2.D0-1.D0)*(1.5D0-FI/2.D0)*(.5D0-FI/2.D0)*A(I)
D2DBDT2=B
RETURN
END
DOUBLE PRECISION FUNCTION FINDTV(PP)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
P=PP
IF(P.LT..842105D0)GO TO 12
T=5.D0
PCAL=VPN(T)
GO TO 13
12 PCAL=.049737D0
IF(DABS(P-PCAL)-.0000001D0*PP)11,11,1
1 T=2.1720D0
13 DO 10 I=1,50
DP=DPDTVP(T)
DEL=(PCAL-P)/DP
T=T-DEL
PCAL=VPN(T)
IF(DABS(P-PCAL)-.0000001D0*P)11,11,2
2 IF(DABS(DEL)-.0000001D0*T)11,11,10
10 CONTINUE
WRITE(*,100)
11 FINDTV=T
RETURN
100 FORMAT(' TEMPERATURE ITTERATION FAILED AT T=',E14.7)
END
DOUBLE PRECISION FUNCTION DPDTVP(TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION C(12),D(14)
DATA C/-3.9394635287D0,141.27497598D0,-1640.7741565D0
1,11974.557102D0,-55283.309818D0,
1166219.56504D0,-325212.82840D0,398843.22750D0,
2-277718.06992D0,83395.204183D0,0.D0,0.D0/
DATA D/-49.510540356D0,651.9236417D0,-3707.5430856D0
1,12880.673491D0,
1 -30048.545554D0,49532.267436D0,-59337.558548D0,52311.296025D0,
2-33950.233134D0,16028.674003D0,-5354.1038967D0,1199.0301906D0,
3 -161.46362959D0,9.8811553386D0/
P=0.0D0
T=TT-DELT(TT)
IF(T-2.1720D0)10,10,1
1 DO 5 I=1,10
5 P=P+C(I)*T**(1-I)*(2-I)
DPDTVP=P*VPN(T)
RETURN
10 DO 15 I=1,14
15 P=P+D(I)*T**(1-I)*(2-I)
DPDTVP=P*VPN(T)
RETURN
END
DOUBLE PRECISION FUNCTION VPN(TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION C(12),D(14)
DATA C/-3.9394635287,141.27497598,-1640.7741565,11974.557102,
1-55283.309818,166219.56504,-325212.82840,398843.22750,
2-277718.06992,83395.204183,0.D0,0.D0/
DATA D/-49.510540356,651.9236417,-3707.5430856,12880.673491,
1 -30048.545554,49532.267436,-59337.558548,52311.296025,
2-33950.233134,16028.674003,-5354.1038967,1199.0301906,
3 -161.46362959,9.8811553386/
T=TT
T=T-DELT(T)
P=0.0D0
IF(T-2.1720D0)10,10,1
1 DO 5 I=1,10
5 P=P+C(I)*T**(2-I)
VPN=DEXP(P)/.76D+6
RETURN
10 DO 15 I=1,14
15 P=P+D(I)*T**(2-I)
VPN=DEXP(P)/.76D+6
RETURN
END
DOUBLE PRECISION FUNCTION DELT(TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
T=TT
DELT=.001D0+.002D0*T
RETURN
END
DOUBLE PRECISION FUNCTION SATV(TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION GV(6)
DATA GV/-.069267495322D0,-.1292532553D0,.29347470712D0
1,-.40806658212D0,.35809505624D0,-.11315580397D0/
DATA DC/.06964D0/
DATA TC/5.2014D0/
T=TT
DCAL=DC
R=(1.D0/TC)
DO 1 I=1,6
FI=I
1 DCAL=DCAL+GV(I)*R**(FI/3.D0)
SATV=DCAL
RETURN
END
DOUBLE PRECISION FUNCTION SATL(TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DIMENSION GL(9)
DATA GL/.12874326484D0,-.43128217346D0,1.7851911824D0
1,-3.3509624489D0,3.0344215824D0,-1.0981289602D0,0.D0,0.D0,0.D0/
DATA DC/.06964D0/
DATA TC/5.2014D0/
T=TT
DCAL=DC
R=(1.D0-T/TC)
DO 2 I=1,6
FI=I
2 DCAL=DCAL+GL(I)*R**(FI/3.D0)
SATL=DCAL
RETURN
END
DOUBLE PRECISION FUNCTION SURFTEN(T)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
SURFTEN=.5308D0*(1.-T/5.2014D0)
RETURN
END
DOUBLE PRECISION FUNCTION R INDEX(D,W)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DD=D*4.0026D0/1000.D0
ALG=0.123396D0-0.0014*DD+33701.617944/D0W**2-12325284955.D0/W**4
FAC=ALG*DD/0.95555D0
R=((2.D0*FAC+1.D0)/(FAC-1.D0))
IF(R.LT.0.0D0)R=R*(-1.)
R=R**(.5D0)
RINDEX=R
RETURN
END
DOUBLE PRECISION FUNCTION TMELT(PP)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
P=PP*1.01325D0/.98066D0
IF(P-17.80D0.LT.0.0D0)GO TO 1
TMELT=((P-17.80D0)/17.3145D0)**(1./1.555414D0)
RETURN
1 TMELT=2.D0
RETURN
END
DOUBLE PRECISION FUNCTION FDIEL(D)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
DD=D*4.0026D0/1000.D0
ALPHA=.123396D0-.0014D0*DD
ALPHA=ALPHA*DD *1.04652D0
EPS=(1.D0+2.D0*ALPHA)/(1.D0-ALPHA)
FDIEL=EPS
RETURN
END
DOUBLE PRECISION FUNCTION THERM(DD,TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
C HELIUM
C THIS ROUTINE CALCULTATES THERMAL CONDUCTIVITY AND VISCOSITY
C FOR AN INPUT OF DEGREES KELVIN AND DENSITY IN MOLES PER LITER
C THE RANGE OF TEMPERATURE IS FROM 2 TO 2000 K
C FOR TEMPERATURES BELOW 300 K FORMULAS OFD VINCE ARP AND GE STEWARD
C ARE USED, FOR TEMPERATURES ABOVE 300 THE DILUTE GAS OF A CRITICAL
C COMPILATION FROM ENGLAND IS USED FOR BOTH VISCOSITY AND
C THERMAL CONDUCTIVITY AND THE EXCESS FUNCTIONS FROM THE ROUTINES BY
C ARP AND STEWART) THE EXCESS FUNCTIONS ARE CALCULATED FOR TEMPS
C ABOVE 300 K WITH THE TEMPERATURE DEPENDENCE HELD AT 300 K
C FOR TEMPS BELOW 300 K TO 100 K THE VISCOSITY EXCESS IS CALC
C FROM STEWARTS ROUTINE BUT THE DILUTE GAS VALUES ARE TAKEN FROM
C THE ENGLISH CORRELATION FOR TEMPS BETWEEN 100 AND 110 T
C DILUTE GAS CALCULATION IS AVARAGED
1 D=DD
T=TT
RHO=D*4.0026D-3
IF(T.LT.300.D0)GO TO 5
THO30=VISCX(300.D0)*.00781736D0
THO300=CONZ(300.D0)
DEL300=DELC(300.D0,RHO)
THO=VISCX(T)*.00781736D0+THO300-THO30
THE=THO300*DEL300-THO300
THERM=THO+THE
RETURN
5 THERM=CONZ(T)*DELC(T,RHO)+CRITIC(T,RHO)
C OUTPUT IN MW/CM.K
RETURN
END
DOUBLE PRECISION FUNCTION VISC(DD,TT)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
D=DD
T=TT
IF(T.LT.100.D0)GO TO 10
IF(T.LT.300.D0)GO TO 8
ETAO=VISCX(T)
ETO30=VISCX(300.D0)
ETO300=VISCDT(0.0D0,300.D0)
ETE300=VISCDT(D,300.D0)-ETO300
VISC=ETAO+ETE300
C OUTPUT UNITS ARE MICROPOISE
RETURN
8 IF(T.LT.110.D0)GO TO 9
ETAO=VISCX(T)
ETEB=VISCDT(D,T)-VISCDT(0.0D0,T)
VISC=ETAO+ETEB
RETURN
9 ETA1=VISCDT(0.0D0,100.D0)
ETA2=VISCX(110.D0)
ETAO=ETA1+(ETA2-ETA1)*(T-100.D0)/10.D0
VISC=ETAO+VISCDT(D,T)-VISCDT(0.0D0,T)
RETURN
10 VISC=VISCDT(D,T)
RETURN
END
DOUBLE PRECISION FUNCTION VISCX( T )
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
VISCX=196.D0*T**.71938D0*DEXP( 12.451D0/T-295.67D0/T/T-4.1249D0)
RETURN
END
DOUBLE PRECISION FUNCTION DELC(TEMP, RHO)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
C K=KZERO*EXPF(B(T)*RHO + C(T)*RHO**2)
C THIS PROGRAM RETURNS EXPF(B(T)*RHO + C(T)*RHO**2)
1 BB=DLOG(TEMP)
CC=1.D0/TEMP
BETTY=4.7470660612D0-5.3641468153D0*BB+3.4639703698D0*BB**2
2-1.0702455443D0*BB**3+0.1571349306D0*BB**4-.00892140047D0*BB**5
B=DEXP(BETTY)
C=2.2109006708D0+187.74174808D0*CC-1281.0947055D0*CC*CC
3+3645.2393216D0*CC**3-3986.6937948D0*CC**4
DELC=DEXP(B*RHO+C*RHO*RHO)
RETURN
END
DOUBLE PRECISION FUNCTION CONZ(TEMP)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
C KZERO IN MILLIWATTS/CM-K, T IN KELVIN9 22 JUNE 71.
1 ANNE=DLOG(TEMP)
PAT=-4.3611622157D0+1.9250159286D0*ANNE-0.52544120165D0*ANNE**2
1+.090045763885D0*ANNE**3-.0054773874708D0*ANNE**4
CONZ=DEXP(PAT)
RETURN
END
DOUBLE PRECISION FUNCTION VISCDT(DGC,T)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
C W.G.STEWARD,S DATA 23 JUNE 71
C INPUT UNITS ARE KELVIN AND MOL/LITER
C OUTPUT UNITS ARE MICROPOISE
TL=DLOG(T)
R=DGC*4.0026D0/1000.D0
ANNE=-0.135311743D0/TL+1.00347841D0+1.20654649D0*TL
1-0.149564551D0*TL*TL+0.0125208416D0*TL**3
BETTY=R*(-47.5295259D0/TL+87.6799309D0-42.0741589D0*TL
1+8.33128289D0*TL*TL-0.589252385D0*TL**3)
CAROL=R*R*(547.309267D0/TL-904.870586D0+431.404928D0*TL
1-81.4504854D0*TL*TL+5.37008433D0*TL**3)
DAGMAR=R**3*(-1684.39324D0/TL+3331.08630D0-1632.19172D0*TL
1+308.804413D0*TL*TL-20.2936367D0*TL**3)
VISCDT=DEXP(ANNE+ BETTY+ CAROL+ DAGMAR)
RETURN
END
DOUBLE PRECISION FUNCTION CRITIC(TEMP,RHO)
IMPLICIT REAL*8 (A-H)
IMPLICIT REAL*8 (O-Z)
IMPLICIT INTEGER*4(I-N)
C CRITICAL ANOMALY FOR HE THERM. CON., SCALED FROM H-2
C T IN KELVIN, REQUIRES DENSITY IN GRAMS/CC AND CP IN JOULES/MOLE
C THIS DECK OF 18 SEPT 70, I HAVE USED MCCARTY"S HE DECKS OF 7/18/70
4 T=TEMP
5 DML=RHO/0.0040026D0
6 IF(T .GE. 11.83D0) GO TO 11
IF(RHO.GT.0.12D0) GO TO 11
7 CP1=CP(DML,T)
8 CP2=CP(DML,11.83D0)
9 CRITIC=0.0026D0*(CP1-CP2)/4.0026D0
10 IF(CRITIC)11,12,12
11 CRITIC=0.0D0
12 RETURN
END