FUNCTION BASE(D,T) C C THIS FUNCTION CALCULATES Z (=PBASE/(DRT) OF EQ. Q) (CALLED BASE), C AND ALSO ABASE,GBASE,SBASE,UBASE,HBASE,CVBASE, AND 1/(DRT) * DP/DT C FOR THE BASE FCT (CALLED DPDTB). C THE AB,GB,SB,UB,HB AND CVB ARE CALCULATED IN DIMENSIONLESS UNITS C AB/RT, GB/RT, SB/R, ETC. C G1,G2 AND GF ARE THE ALPHA, BETA AND GAMMA OF EQ 2, WHICH ARE C SUPPLIED BY THE BLOCKDATA ROUTINE. B1 AND B2 ARE THE "EXCLUDED C VOLUME" AND "2ND VIRIAL" (EQS 3 AND 4) SUPPLIED BY THE SUBROUTINE C BB(T). WHICH ALSO SUPPLIES THE 1ST AND 2ND DERIVATIVES WITH C RESPECT TO T (B1T,B2T,B1TT,B2TT). C IMPLICIT REAL*8 (A-H,O-Z) COMMON /ELLCON/ G1,G2,GF,B1,B2,B1T,B2T,B1TT,B2TT COMMON /BASEF/ AB,GB,SB,UB,HB,CVB,DPDTB COMMON /ACONST/ WM,GASCON,TZ,A,Z,DZ,Y,UREF,SREF C Y=.25D0*B1*D X=1.D0-Y Z0=(1.D0+G1*Y+G2*Y*Y)/X**3 Z=Z0+4.D0*Y*(B2/B1-GF) DZ0=(G1+2.D0*G2*Y)/X**3 + 3.D0*(1.D0+G1*Y+G2*Y*Y)/X**4 DZ=DZ0+4.D0*(B2/B1-GF) AB = -DLOG(X)-(G2-1.D0)/X+28.16666667D0/X/X+4.D0*Y*(B2/B1-GF) 1 +15.166666667D0 + DLOG(D*T*GASCON/.101325D0) GB = AB + Z BASE = Z BB2TT=T*T*B2TT UB= -T*B1T*(Z-1.D0-D*B2)/B1-D*T*B2T HB=Z+UB CVB=2.D0*UB+(Z0-1.D0)*((T*B1T/B1)**2-T*T*B1TT/B1) 1 - D*(BB2TT - GF*B1TT*T*T) -(T*B1T/B1)**2*Y*DZ0 DPDTB=BASE/T + BASE*D/Z*(DZ*B1T/4.D0+B2T-B2/B1*B1T) SB = UB - AB RETURN END SUBROUTINE BB(T) C C THIS SUBROUTINE CALCULATES THE B'S OF EQS. 3,4 USING COEFFICIENTS C FROM BLOCKDATA, CALCULATING ALSO THE FIRST AND SECOND DERIVATIVE C WITH RESPECT TO TEMP. THE B'S CALCULATED HERE ARE IN CM3/G. C IMPLICIT REAL*8(A-H,O-Z) COMMON /ELLCON/ G1,G2,GF,B1,B2,B1T,B2T,B1TT,B2TT COMMON /ACONST/ WM,GASCON,TZ,AA,Z,DZ,Y,UREF,SREF COMMON /BCONST/ BP(10),BQ(10) DIMENSION V(10) C V(1)=1. DO 2 I=2,10 2 V(I)=V(I-1)*TZ/T B1=BP(1)+BP(2)*DLOG(1./V(2)) B2=BQ(1) B1T=BP(2)*V(2)/TZ B2T=0. B1TT=0. B2TT=0. DO 4 I=3,10 B1=B1+BP(I)*V(I-1) B2=B2+BQ(I)*V(I-1) B1T=B1T-(I-2)*BP(I)*V(I-1)/T B2T=B2T-(I-2)*BQ(I)*V(I-1)/T B1TT=B1TT+BP(I)*(I-2)**2*V(I-1)/T/T 4 B2TT=B2TT+BQ(I)*(I-2)**2*V(I-1)/T/T B1TT=B1TT-B1T/T B2TT=B2TT-B2T/T RETURN END BLOCK DATA C C THE BLOCKDATA SUBROUTINE SUPPLIES MOST OF THE FIXED PARAMETERS C USED IN THE REST OF THE ROUTINES. P IS THE B(I) OF TABLE I, C Q IS THE B(I) OF TABLE I, G1,G2 AND FG ARE THE ALPHA, BETA C AND GAMMA OF EQ 2, AND G,II,JJ ARE THE G(I), K(I) AND C L(I) OF EQ 5. C IMPLICIT REAL*8(A-H,O-Z) COMMON /ACONST/ WM,GASCON,TZ,AA,ZB,DZB,YB,UREF,SREF COMMON /NCONST/ G(40),II(40),JJ(40),NC COMMON /ELLCON/ G1,G2,GF,B1,B2,B1T,B2T,B1TT,B2TT COMMON /BCONST/ BP(10),BQ(10) COMMON /ADDCON/ ATZ(4),ADZ(4),AAT(4),AAD(4) DATA ATZ/2*64.D1,641.6D0,27.D1/,ADZ/3*.319D0,1.55D0/,AAT/2*2.D4 1,4.D4,25.D0/,AAD/34.D0,4.D1,3.D1,1.05D3/ DATA WM/18.0152D0/,GASCON/.461522D0/,TZ/647.073D0/,AA/1.D0/,NC/36/ DATA UREF,SREF/-4328.4549774261D0,7.6180720166752D0/ DATA G1,G2,GF/11.D0,44.333333333333D0,3.5D0/ DATABP/.7478629D0,-.3540782D0,2*0.D0,.7159876D-2,0.D0,-.3528426D-2 1,3*0.D0/,BQ/1.1278334D0,0.D0,-.5944001D0,-5.010996D0,0.D0 2,.63684256D0,4*0.D0/ DATA G/-.53062968529023D3,.22744901424408D4,.78779333020687D3 1,-.69830527374994D2,.17863832875422D5,-.39514731563338D5 2,.33803884280753D5,-.13855050202703D5,-.25637436613260D6 3,.48212575981415D6,-.34183016969660D6, .12223156417448D6 4,.11797433655832D7,-.21734810110373D7,.10829952168620D7 5,-.25441998064049D6,-.31377774947767D7,.52911910757704D7 6,-.13802577177877D7,-.25109914369001D6, .46561826115608D7 7,-.72752773275387D7,.41774246148294D6,.14016358244614D7 8,-.31555231392127D7,.47929666384584D7,.40912664781209D6 9,-.13626369388386D7, .69625220862664D6,-.10834900096447D7 A,-.22722827401688D6,.38365486000660D6,.68833257944332D4 B,.21757245522644D5,-.26627944829770D4,-.70730418082074D5 C,-.225D0,-1.68D0,.055D0,-93.0D0/ DATA II/4*0,4*1,4*2,4*3,4*4,4*5,4*6,4*8,2*2,0,4,3*2,4/ DATA JJ/2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7 1,2,3,5,7,1,3*4,0,2,0,0/ END SUBROUTINE CORR(T,P,DL,DV,DELG) C C SUBROUTINE CORR WILL CALCULATE, FOR AN INPUT T AND P AT OR NEAR C THE VAPOR PRESSURE, THE CORRESPONDING LIQUID AND VAPOR DENSITIES C AND ALSO DELG = (GL-GV)/RT FOR USE IN CALCULATING THE CORRECTION C TO THE VAPOR PRESSURE IN PCORR OR VAPOR TEMPERATURE IN TCORR FOR C DELG = 0. NOTE THAT DELG IS DEFINED AS ZERO FOR ALL TEMPERATURES C ABOVE 646.3K. C IMPLICIT REAL*8(A-H,O-Z) COMMON /QQQQ/ Q00,Q11 COMMON /ACONST/ WM,GASCON,TZ,AA,ZB,DZB,YB,UREF,SREF COMMON /FCTS/ AD,GD,SD,UD,HD,CVD,CPD,DPDT,DVDT,CJTT,CJTH DATA TCRIT / 647.1260000001D0 / DATA PCRIT / 22.055D0 / C RT=GASCON*T CALL BB(T) IF(T.GT.646.3D0) GO TO 101 DLIQ=DL IF(DL.LE.0.) DLIQ=1.11-.0004*T CALL DFIND(DL,P,DLIQ,T,DQ) CALL THERM(DL,T) GL=GD DVAP=DV IF(DV.LE.0.) DVAP=P/RT CALL DFIND(DV,P,DVAP,T,DQ) IF(DV.LT.5.D-7) DV=5.D-7 CALL THERM(DV,T) GV=GD DELG = GL-GV RETURN 101 P = PCRIT DELG=0 IF(T .GT. TCRIT) RETURN TAU=.657128D0*(1.-T/647.126D0)**.325 DL=.322D0+TAU DV=.322D0-TAU ZDUM = BASE(DV,T) CALL QQ(T,DV) P=RT*DV*ZB + Q00 RETURN END SUBROUTINE DFIND(DOUT,P,D,T,DPD) C C ROUTINE TO FIND DENSITY CORRESPONDING TO INPUT PRESSURE P(MPA), AND C TEMPERATURE T(K), USING INITIAL GUESS DENSITY D(G/CM3). THE OUTPUT C DENSITY IS IN G/CM3, ALSO, AS A BYPRODUCT, DP/DRHO IS CALCULATED C ("DPD", MPA CM3/G) C IMPLICIT REAL*8(A-H,O-Z) COMMON /QQQQ/ Q0,Q5 COMMON /ACONST/ WM,GASCON,TZ,AA,Z,DZ,Y,UREF,SREF C DD=D RT=GASCON*T IF(DD.LE.0.D0) DD=1.D-8 IF(DD.GT.1.9D0) DD=1.9D0 L=0 9 L = L + 1 11 IF(DD.LE.0.D0) DD=1.D-8 IF(DD.GT.1.9D0) DD=1.9D0 CALL QQ(T,DD) PP =RT*DD*BASE(DD,T)+Q0 DPD=RT*(Z+Y*DZ)+Q5 C THE FOLLOWING 3 LINES CHECK FOR NEGATIVE DP/DRHO, AND IF SO ASSUME C GUESS TO BE IN 2-PHASE REGION, AND ADJUST GUESS ACCORDINGLY. IF(DPD.GT.0.D0) GO TO 13 IF(D.GE..2967D0) DD=DD*1.02D0 IF(D.LT..2967D0) DD=DD*.98D0 IF(L.LE.10) GO TO 9 13 DPDX=DPD*1.1D0 IF(DPDX.LT..1D0) DPDX=.1D0 DP=DABS(1.D0-PP/P) IF(DP.LT.1.D-8) GO TO 20 IF(D.GT..3 .AND. DP.LT.1.D-7) GO TO 20 IF(D.GT..7 .AND. DP.LT.1.D-6) GO TO 20 X=(P-PP)/DPDX IF(DABS(X).GT..1D0) X=X*.1D0/DABS(X) DD=DD+X IF(DD.LE.0.) DD=1.D-8 19 IF(L.LE.30) GO TO 9 WRITE(*,25) P, T, D 20 CONTINUE DOUT=DD 25 FORMAT('***NONCONVERGENCE IN DFIND'/' P =',G15.7,2X, 1 'T =',G15.7, 2X,'D =',G15.7 ) RETURN END SUBROUTINE IDEAL(T) C C THIS SUBROUTINE CALCULATES THE THERMODYNAMIC PROPERTIES FOR C WATER IN THE IDEAL GAS STATE FROM FUNCTION OF H.W. WOOLLEY C IMPLICIT REAL*8 (A-H,O-Z) COMMON /IDF/ AI,GI,SI,UI,HI,CVI,CPI DIMENSION C(18) DATA C /.19730271018D2,.209662681977D2,-.483429455355D0 1,.605743189245D1,22.56023885D0,-9.87532442D0,-.43135538513D1 2,.458155781D0,-.47754901883D-1,.41238460633D-2,-.27929052852D-3 3,.14481695261D-4,-.56473658748D-6,.16200446D-7,-.3303822796D-9 4,.451916067368D-11,-.370734122708D-13,.137546068238D-15/ C TT=T/1.D2 TL=DLOG(TT) GI=-(C(1)/TT+C(2))*TL HI=(C(2)+C(1)*(1.D0-TL)/TT) CPI=C(2)-C(1)/TT DO 8 I=3,18 GI=GI-C(I)*TT**(I-6) HI=HI+C(I)*(I-6)*TT**(I-6) 8 CPI=CPI+C(I)*(I-6)*(I-5)*TT**(I-6) AI=GI-1.D0 UI=HI-1.D0 CVI=CPI-1.D0 SI=UI-AI RETURN END SUBROUTINE PCORR(T,P,DL,DV) C C SUBROUTINE PCORR WILL CALCULATE THE VAPOR PRESSURE P AND THE LIQUID C AND VAPOR DENSITIES CORRESPONDING TO THE INPUT T, CORRECTED SUCH C THAT GL-GV=0. THE FUNCTION PS IS REQUIRED WHICH WILL GIVE A C REASONABLY GOOD APPROXIMATION TO THE VAPOR PRESSURE TO BE USED AS C THE STARTING POINT FOR THE ITERATION. NOTE - NO LIMIT ON THE C NUMBER OF ITERATIONS. C IMPLICIT REAL*8 (A-H,O-Z) COMMON /ACONST/ WM,GASCON,TZ,AA,ZB,DZB,YB,UREF,SREF C P = PS(T) 2 CALL CORR(T,P,DL,DV,DELG) DP=DELG*GASCON*T/((1./DV-1./DL)+1D-13) P = P+DP IF(DABS(DELG).LT.1.D-8) RETURN GO TO 2 END FUNCTION PS(T) C C THIS FUNCTION CALCULATES AN APPROXIMATION TO THE VAPOR PRESSURE, C PS, AS A FUNCTION OF THE INPUT TEMPERATURE. THE VAPOR PRESSURE C CALCULATED AGREES WITH THE VAPOR PRESSURE PREDICTED BY THE SURFACE C TO WITHIN .02% TO WITHIN A DEGREE OR SO OF THE CRITICAL TEMPERATURE C AND CAN SERVE AS AN INITIAL GUESS FOR FURTHER REFINEMENT BY C IMPOSING THE CONDITION THAT G1=GV. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(8) DATA A/-7.8889166D0,2.5514255D0,-6.716169D0 1,33.239495D0,-105.38479D0,174.35319D0,-148.39348D0 2,48.631602D0/ C IF(T.GT.314.D0) GO TO 2 PL=6.3573118D0-8858.843D0/T+607.56335D0*T**(-.6) PS=.1D0*DEXP(PL) RETURN 2 IF(T .GT. 646.3D0) GO TO 5 V=T/647.25D0 W=DABS(1.D0-V) B=0.D0 DO 4 I=1,8 Z=I 4 B=B+A(I)*W**((Z+1.D0)/2.D0) Q=B/V PS=22.093D0*DEXP(Q) RETURN 5 PS = -147.7854939D0 + 0.262453402D0*T RETURN END SUBROUTINE QQ(T,D) C C THIS ROUTINE CALCULATES, FOR A GIVEN T(K) AND D(G/CM3), THE C RESIDUAL CONTRIBUTIONS TO PRESSURE (Q), HELMHOLTZ FCT (AR), C DP/DRHP (Q5), AND ALSO THE THE GIBBS FUNCTION, ENTROPY,INTERNAL C ENERGY, ENTHALPY, ISOCHORIC HEAT CAPACITY AND DPDT. (EQ 5). C TERMS 37 THRU 39 ARE THE ADDITIONAL TERMS AFFECTING ONLY THE C IMMEDIATE VICINITY OF THE CRITICAL POINT, AND TERM 40 IS THE C ADDITIONAL TERM IMPROVING THE LOW T, HIGH P REGION. C IMPLICIT REAL*8(A-H,O-Z) COMMON /NCONST/ G(40),II(40),JJ(40),N COMMON /ACONST/ WM,GASCON,TZ,AA,Z,DZ,Y,UREF,SREF COMMON /ADDCON/ ATZ(4),ADZ(4),AAT(4),AAD(4) COMMON /RESF/ AR,GR,SR,UR,HR,CVR,DPDTR COMMON /QQQQ/ Q,Q5 DIMENSION QR(11),QT(10),QZR(9),QZT(9) EQUIVALENCE (QR(3),QZR(1)),(QT(2),QZT(1)) C RT=GASCON*T QR(1)=0.D0 Q5=0.D0 Q=0.D0 AR=0.D0 DADT=0.D0 CVR=0.D0 DPDTR=0.D0 EX0=-AA*D EX0=DMAX1(EX0,-225.D0) EX0=DMIN1(EX0, 225.D0) E=DEXP(EX0) Q10=D*D*E Q20=1.D0-E QR(2)=Q10 V=TZ/T QT(1)=T/TZ DO 4 I=2,10 QR(I+1)=QR(I)*Q20 4 QT(I)=QT(I-1)*V DO 10 I=1,N K=II(I)+1 L=JJ(I) ZZ=K QP=G(I)*AA*QZR(K-1)*QZT(L) Q=Q+QP Q5 = Q5 + AA*(2./D-AA*(1.-E*(K-1)/Q20))*QP AR=AR+G(I)*QZR(K)*QZT(L)/Q10/ZZ/RT DFDT=Q20**K*(1-L)*QZT(L+1)/TZ/K D2F=L*DFDT DPT=DFDT*Q10*AA*K/Q20 DADT=DADT+G(I)*DFDT DPDTR=DPDTR+G(I)*DPT 10 CVR=CVR+G(I)*D2F/GASCON QP=0.D0 Q2A=0.D0 DO 20 J=37,40 IF(G(J).EQ.0.D0) GO TO 20 K=II(J) KM=JJ(J) DDZ = ADZ(J-36) DEL = D/DDZ - 1.D0 IF(DABS(DEL).LT.1.D-10) DEL=1.D-10 DD = DEL*DEL EX1 = -AAD(J-36)*DEL**K EX1=DMAX1(EX1,-225.D0) EX1=DMIN1(EX1, 225.D0) DEX = DEXP(EX1)*DEL**KM ATT = AAT(J-36) TX = ATZ(J-36) TAU = T/TX-1.D0 EX2 = -ATT*TAU*TAU EX2=DMAX1(EX2,-225.D0) EX2=DMIN1(EX2, 225.D0) TEX = DEXP(EX2) Q10 =DEX*TEX QM = KM/DEL - K*AAD(J-36)*DEL**(K-1) FCT=QM*D**2*Q10/DDZ Q5T = FCT*(2./D+QM/DDZ)-(D/DDZ)**2*Q10*(KM/DEL/DEL+ 1 K*(K-1)*AAD(J-36)*DEL**(K-2)) Q5 = Q5 + Q5T*G(J) QP = QP + G(J)*FCT DADT = DADT - 2.D0*G(J)*ATT*TAU*Q10/TX DPDTR = DPDTR - 2.D0*G(J)*ATT*TAU*FCT/TX Q2A = Q2A + T*G(J)*(4.D0*ATT*EX2+2.D0*ATT)*Q10/TX/TX AR = AR + Q10*G(J)/RT 20 CONTINUE SR=-DADT/GASCON UR=AR+SR CVR=CVR+Q2A/GASCON Q=Q+QP RETURN END SUBROUTINE SECVIR(T,VIR) C C THIS SUBROUTINE CALCULATES THE SECOND VIRIAL IN CM3/G C AT TEMPERATURE T IN K. C THIS SUBROUTINE IS FOR CONVENIENCE ONLY, AND NOT USED BY PROGRAM. C IMPLICIT REAL*8(A-H,O-Z) COMMON /NCONST/ G(40),II(40),JJ(40),NC COMMON /ELLCON/ G1,G2,GF,BB1,BB2,B1T,B2T,B1TT,B2TT COMMON /ACONST/ WM,GASCON,TC,AA,Z,DZ,Y,UREF,SREF C CALL BB(T) V=TC/T VIR=BB2 DO 3 J=1,NC IF(II(J).NE.0) GO TO 3 L=JJ(J) VIR=VIR+G(J)*V**(L-1)/T/GASCON 3 CONTINUE RETURN END SUBROUTINE TCORR(T,P,DL,DV) C C SUBROUTINE TCORR WILL CALCULATE THE VAPOR TEMPERATURE T AND THE C LIQUID AND VAPOR DENSITIES CORRESPONDING TO THE INPUT P, C CORRECTED SUCH THAT GL - GV = 0. C SUBROUTINE TCORR IS SIMILAR TO "PCORR" EXCEPT THAT THE VAPOR C TEMPERATURE CORRESPONDING TO THE INPUT VAPOR PRESSURE IS FOUND. C FUNCTIONS CALLED ARE TSAT AND TDPSDT WHICH GIVES AN APPROXIMATION C TO T(SAT) AS A FUNCTION OF T*DP(SAT)/DT. NOTE - NO LIMIT ON C NUMBER OF ITERATIONS. C IMPLICIT REAL*8 (A-H,O-Z) COMMON /ACONST/ WM,GASCON,TZ,AA,ZB,DXB,YB,UREF,SREF C T = TSAT(P) IF(T.EQ.0.D0) RETURN 2 CALL CORR(T,P,DL,DV,DELG) DP=DELG*GASCON*T/((1./DV-1./DL)+1D-13) T = T*(1.-DP/TDPSDT(T)) IF(DABS(DELG).LT.1.D-8) RETURN GO TO 2 END FUNCTION TDPSDT(T) C C THIS FUNCTION CALCULATES T*(DPS/DT), AND IS USED BY THE FUNCTION C TSAT AND TCORR. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(8) DATA A / -7.8889166D0,2.5514255D0,-6.716169D0 1,33.239495D0,-105.38479D0,174.35319D0,-148.39348D0 2,48.631602D0/ C V=T/647.25 W=1.-V B=0. C=0. DO 4 I=1,8 Z=I Y=A(I)*W**((Z+1.)/2.) C=C+Y/W*(.5-.5*Z-1./V) 4 B=B+Y Q=B/V TDPSDT=22.093D0*DEXP(Q)*C RETURN END SUBROUTINE THERM(D,T) C C THIS SUBROUTINE CALCULATES THE THERMODYNAMIC FUNCTIONS IN C DIMENSIONLESS UNITS (AD=A-RT, GD=G/RT, SD=S/R, UD=U/RT, C HD=H/RT, CVD=CV/R, AND CPD=CP/R) C IMPLICIT REAL*8(A-H,O-Z) COMMON /ACONST/ WM,GASCON,TZ,AA,ZB,DZB,Y,UREF,SREF COMMON /QQQQ/ QP, QDP COMMON /BASEF/AB,GB,SB,UB,HB,CVB,DPDTB COMMON /RESF/AR,GR,SR,UR,HR,CVR,DPDTR COMMON /IDF/ AI,GI,SI,UI,HI,CVI,CPI COMMON /FCTS/AD,GD,SD,UD,HD,CVD,CPD,DPDT,DVDT,CJTT,CJTH C CALL IDEAL(T) RT = GASCON*T Z = ZB + QP/RT/D DPDD = RT*(ZB+Y*DZB) + QDP AD = AB + AR + AI - UREF/T + SREF GD = AD + Z UD = UB + UR + UI - UREF/T DPDT = RT*D*DPDTB + DPDTR CVD = CVB + CVR + CVI CPD = CVD + T*DPDT**2/(D*D*DPDD*GASCON) HD = UD + Z SD = SB + SR + SI - SREF DVDT=DPDT/DPDD/D/D CJTT=1./D-T*DVDT CJTH=-CJTT/CPD/GASCON RETURN END FUNCTION TSAT(P) C C THIS FUNCTION CALCULATES THE SATURATION TEMPERATURE FOR A GIVEN C PRESSURE, BY AN ITERATIVE PROCESS USING PS(T) AND TDPSDT. C IMPLICIT REAL*8 (A-H,O-Z) DATA PCRIT / 22.055D0 / DATA TCRIT / 647.1260000001D0 / C TSAT=647.126D0 IF(P.GT.PCRIT) RETURN K=0 PL=2.302585D0+DLOG(P) C PL=LOGE(10)+LOGE(P) TO CONVERT EQUATION FROM BARS TO MPA TG=372.83+PL*(27.7589+PL*(2.3819+PL*(.24834+PL*.0193855))) IF(TG .GT. 647.D0) TG = 647.D0 1 IF(TG.LT.273.15D0) TG=273.15D0 IF(K.LT.8) GO TO 2 WRITE(*,25) K,P,PP,TG,DP GO TO 8 2 K=K+1 PP=PS(TG) IF(ABS(1.-PP/P).LT.1.D-5) GO TO 8 DP=TDPSDT(TG) TG = TG*(1.+(P-PP)/DP) IF(TG .GT. 647.126D0) TG = 647.126D0 GO TO 1 8 TSAT=TG RETURN 25 FORMAT( '***NONCONVERGENCE IN TSAT(P)' / ' K =',I4,2X,'P =', 1 F10.4,2X,'PP =',F10.4,2X,'TG =',F10.4,2X,'DP =',F10.4 ) END FUNCTION TTI(T) C C FUNCTION TO CONVERT INTERNAL TEMPERATURES IN DEG K TO EXTERNAL UNITS C REAL*8 T,TTI,FT,FD,FP,FH,NT,ND,NP,NH COMMON /UNITS/ IT,ID,IP,IH,FT,FD,FP,FH COMMON /CUNITS/ NT,ND,NP,NH C GO TO (5,6,7,8),IT 5 TTI=T RETURN 6 TTI=T-273.15D0 RETURN 7 TTI=T*1.8D0 RETURN 8 TTI=T*1.8D0-459.67D0 RETURN END FUNCTION TTT(T) C C FUNCTION TO CONVERT INPUT TEMPERATURES IN EXTERNAL UNITS TO DEG K C REAL*8 T,TTT,FT,FD,FP,FH,NT,ND,NP,NH COMMON /UNITS/ IT,ID,IP,IH,FT,FD,FP,FH COMMON /CUNITS/ NT,ND,NP,NH C GO TO (1,2,3,4),IT 1 TTT=T FT=1. RETURN 2 TTT=T+273.15D0 FT=1. RETURN 3 TTT=T/1.8D0 FT=.5555555555556D0 RETURN 4 TTT=(T+459.67D0)/1.8D0 FT=.5555555555556D0 RETURN END SUBROUTINE UNIT(IFC) C C SUBROUTINE UNIT QUERIES THE USER AS TO HIS CHOICE OF UNITS AND SETS C INTERNAL PARAMETERS APPROPRIATELY. THE INTERNAL UNITS OF THE C PROGRAM ARE TEMPERATURE IN K, DENSITY IN G/CM3, ALL OTHER C QUANTITIES ARE CALCULATED IN DIMENSIONLESS UNITS AND DIMENSIONED C AT OUTPUT TIME. C IMPLICIT REAL*8 (A-H,O-Z) COMMON /UNITS/ IT,ID,IP,IH,FT,FD,FP,FH COMMON /CUNITS/ NT,ND,NP,NH CHARACTER *8 NT,ND,NP,NH,NNT,NND,NNP,NNH CHARACTER *12A1, A2, A3, A4 DIMENSION FFD(4),FFP(5),FFH(6),NNT(4),NND(4),NNP(5),NNH(6) DATA FFD/1.D-3,1.D0,.0180152D0,.0160184634D0/ DATA FFP/1.D0,10.D0,9.869232667D0,145.037738D0,10.1971621D0/ DATA FFH/2*1.D0,18.0152D0,.238845897D0,4.30285666D0,.429922614D0/ DATA NNT / 'K', 'C', 'R', 'F' / DATA NND / 'KG/M3', 'G/CM3', 'MOL/L', 'LB/FT3' / DATA NNP / 'MPA ' , 'BAR ' , 'ATM ' , 'PSI ' , 'KG/CM2' / DATA NNH/ 'KJ/KG', 'J/G ' , 'J/MOL', 'CAL/G', 'CA/MOL', 'BTU/LB'/ DATA A1, A2, A3, A4 / 'TEMPERATURE', 'DENSITY' 1, 'PRESSURE', 'ENERGY' / C IF(IFC .NE. 5) GO TO 301 WRITE(*,11) A1 30 WRITE(*,12) 301 READ(IFC,*,END=99) IT IF(IT.EQ.0) STOP IF(IT .GT. 4 .OR. IT .LT. 1) GO TO 30 NT=NNT(IT) IF(IFC .NE. 5) GO TO 311 WRITE(*,11) A2 31 WRITE(*,13) 311 READ(IFC,*,END=99) ID IF(ID.GT.4 .OR. ID.LT.1) GOTO 31 ND=NND(ID) FD=FFD(ID) IF(IFC .NE. 5) GO TO 321 WRITE(*,11) A3 32 WRITE(*,14) 321 READ(IFC,*,END=99) IP IF(IP.GT.5 .OR. IP.LT.1) GOTO 32 NP=NNP(IP) FP=FFP(IP) IF(IFC .NE. 5) GO TO 331 WRITE(*,11) A4 33 WRITE(*,15) 331 READ(IFC,*,END=99) IH IF(IH.GT.6 .OR. IH.LT.1) GOTO 33 NH=NNH(IH) FH=FFH(IH) IF(IFC .NE. 5) WRITE(*,16) NT, ND, NP, NH RETURN 99 STOP 'END NBS/NRC STEAM PROPERTIES - TO RESTART ENTER Y RETURN' 11 FORMAT(' ENTER UNITS CHOSEN FOR ',A12) 12 FORMAT(' CHOOSE FROM 1=DEG K, 2=DEG C, 3=DEG R, 4=DEG F') 13 FORMAT(' CHOOSE FROM 1=KG/M3, 2=G/CM3, 3=MOL/L, 4=LB/FT3') 14 FORMAT(' CHOOSE FROM 1=MPA, 2=BAR, 3=ATM, 4=PSIA, 5=KG/CM2') 15 FORMAT(' CHOOSE FROM 1=KJ/KG, 2=J/G, 3=J/MOL, 4=CALORIES/G, 5=CALO 1RIES/MOL, 6=BTU/LB') 16 FORMAT(' UNITS: TEMPERATURE= ',A1,' DENSITY= ',A6, 1 ' PRESSURE= ',A6,' ENERGY= ',A7) END