      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
C   2 V(I)=V(I-1)*TZ/T
      V(I)=V(I-1)*TZ/T
    2 CONTINUE
      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
C   4 B2TT=B2TT+BQ(I)*(I-2)**2*V(I-1)/T/T
      B2TT=B2TT+BQ(I)*(I-2)**2*V(I-1)/T/T
    4 CONTINUE
      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
C  11 IF(DD.LE.0.D0) DD=1.D-8
      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
C 19  IF(L.LE.30) GO TO 9
      IF(L.LE.30) GO TO 9
      CONTINUE
      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)
C  8  CPI=CPI+C(I)*(I-6)*(I-5)*TT**(I-6)
      CPI=CPI+C(I)*(I-6)*(I-5)*TT**(I-6)
   8  CONTINUE
      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
C  4  B=B+A(I)*W**((Z+1.D0)/2.D0)
      B=B+A(I)*W**((Z+1.D0)/2.D0)
   4  CONTINUE
      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
C  4  QT(I)=QT(I-1)*V
      QT(I)=QT(I-1)*V
   4  CONTINUE
      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
C 10  CVR=CVR+G(I)*D2F/GASCON
      CVR=CVR+G(I)*D2F/GASCON
  10  CONTINUE
      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)
C  4  B=B+Y
      B=B+Y
   4  CONTINUE
      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
