c  These routines accept both single-phase and two-phase states as the
c  input; if the phase is known, the specialized routines are faster
c
c  inputs--two of the following as indicated by the first two letters of
c          the subroutine name:
c        t--temperature [K]
c        p--pressure [kPa]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [[J/mol-K]
c        q--vapor quality [basis specified by kq]
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q < 0 or q > 1 are not allowed and will result in warning

c  additional input--required for all routines
c        z--overall (bulk) composition [array of mol frac]

c  additional input--only for TQFLSH and PQFLSH
c       kq--flag specifying units for input quality
c           kq = 1 quality on MOLAR basis [moles vapor/total moles]
c           kq = 2 quality on MASS basis [mass vapor/total mass]
c

c  outputs--one, two, or all of the following, depending on the inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c        D--overall (bulk) molar density [mol/L]
c
c  additional outputs--common to all routines
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c           if only one phase is present, Dl = Dv = D
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c           if only one phase is present, x = y = z
c
c  additional output--common to all routines except TQFLSH and PQFLSH
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c           q < 0 indicates subcooled (compressed) liquid
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q > 1 indicates superheated vapor
c           q = 998 superheated vapor, but quality not defined (t > Tc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c
c  additional outputs--common to all routines, except that input
c                      quantities are not repeated
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful
c                        1 = Tin < Tmin
c                        4 = Pin < 0
c                        5 = T and P out of range
c                        8 = x out of range (component and/or sum < 0
c                            or > 1)
c                        9 = x and T out of range
c                       12 = x out of range and P < 0
c                       13 = x and T and P out of range
c     herr--error string (character*255 variable if ierr<>0)



SINGLE-PHASE FLASH SUBROUTINES

These routines accept only single-phase states as inputs.  They will be
faster than the corresponding general routines, but will fail if called
with an incorrect phase specification.  The phase-specific subroutines
also do not check limits, so may fail if called outside the range of the
equation of state.  The following single-phase routines are available.

      subroutine TEFL1 (t,e,z,Dmin,Dmax,D,ierr,herr)
      subroutine THFL1 (t,h,z,Dmin,Dmax,D,ierr,herr)
      subroutine TSFL1 (t,s,z,Dmin,Dmax,D,ierr,herr)
      subroutine HSFL1 (h,s,z,Dmin,Dmax,t,D,ierr,herr)
c
c  inputs--two of the following as indicated by the first two letters of c          the subroutine name:
c        t--temperature [K]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [[J/mol-K]

c  additional inputs--required for all routines
c        z--overall (bulk) composition [array of mol frac]
c     Dmin--lower bound on density [mol/L]
c     Dmax--upper bound on density [mol/L]
c
c  outputs:
c        t--temperature [K] (present only for HSFL1)
c        D--molar density [mol/L]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)


Single-phase flash routines (continued)

      subroutine PEFL1 (p,e,z,kph,t,D,ierr,herr)
      subroutine PHFL1 (p,h,z,kph,t,D,ierr,herr)
      subroutine PSFL1 (p,s,z,kph,t,D,ierr,herr)
c
c  inputs--two of the following as indicated by the first two letters of
c          the subroutine name:
c        p--pressure [kPa]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [[J/mol-K]

c  additional inputs--required for all routines
c        z--overall (bulk) composition [array of mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c
c  outputs:
c        t--temperature [K]
c        D--molar density [mol/L]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)


      subroutine DEFL1 (D,e,z,t,ierr,herr)
      subroutine DHFL1 (D,h,z,t,ierr,herr)
      subroutine DSFL1 (D,s,z,t,ierr,herr)
      subroutine PDFL1 (p,D,z,t,ierr,herr)
c
c  inputs--two of the following as indicated by the first two letters of
c          the subroutine name:
c        p--pressure [kPa]
c        D--density [mol/L]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [[J/mol-K]

c  additional inputs--required for all routines
c        z--overall (bulk) composition [array of mol frac]
c
c  outputs:
c        t--temperature [K]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)

The single-phase temperature-pressure flash is called many times by
other routines and has been optimized for speed and requires a specific
calling sequence.

      subroutine TPRHO (t,p,x,kph,kguess,rho,ierr,herr)
c
c  iterate for density as a function of temperature, pressure, and
c  composition for a specified phase
c
c  inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c        x--composition [array of mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c                 N.B.:  0 = stable phase--NOT ALLOWED (use TPFLSH)
c   kguess--input flag:  1 = first guess for rho provided
c                        0 = no first guess provided
c      rho--first guess for molar density [mol/L], only if kguess = 1
c
c  outputs:
c      rho--molar density [mol/L]
c     ierr--error flag:  0 = successful
c                      200 = CRITP did not converge
c                      201 = illegal input (kph <= 0)
c                      202 = liquid-phase iteration did not converge
c                      203 = vapor-phase iteration did not converge
c     herr--error string (character*255 variable if ierr<>0)

TWO-PHASE FLASH SUBROUTINES

These routines accept only two-phase (liquid + vapor) states as inputs.
They will be faster than the corresponding general routines, but will
fail if called with an incorrect phase specification.  The
phase-specific subroutines also do not check limits, so may fail if
called outside the range of the equation of state.  The following
two-phase routines are available.

      subroutine TPFL2 (t,p,z,Dl,Dv,x,y,q,ierr,herr)
      subroutine DHFL2 (D,h,z,t,p,Dl,Dv,x,y,q,ierr,herr)
      subroutine DSFL2 (D,s,z,t,p,Dl,Dv,x,y,q,ierr,herr)
      subroutine DEFL2 (D,e,z,t,p,Dl,Dv,x,y,q,ierr,herr)

c  inputs--two of the following as indicated by the first two letters of c          the subroutine name:
c        t--temperature [K]
c        p--pressure [kPa]
c        D--bulk molar density [mol/L]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [[J/mol-K]
c
c  additional inputs--required for all routines
c        z--overall (bulk) composition [array of mol frac]
c     Dmin--lower bound on density [mol/L]
c     Dmax--upper bound on density [mol/L]
c
c  outputs:
c        t--temperature [K] (not present for TPFL2)
c        p--pressure [kPa] (not present for TPFL2)
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)



In the following two-phase flash routines, there is the option to pass
the dew and bubble point conditions as inputs; if these data are known
(from a previous call to SATT or SATP, for example), these two-phase
routines will be significantly faster than the corresponding general
routines described in Section C.5.1.  Otherwise, the general routines
will be more reliable.

      subroutine TDFL2 (t,D,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  p,Dl,Dv,x,y,q,ierr,herr)
      subroutine TEFL2 (t,e,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  p,Dl,Dv,x,y,q,ierr,herr)
      subroutine THFL2 (t,h,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  p,Dl,Dv,x,y,q,ierr,herr)
      subroutine TSFL2 (t,s,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  p,Dl,Dv,x,y,q,ierr,herr)
      subroutine PDFL2 (p,d,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,Dl,Dv,x,y,q,ierr,herr)
      subroutine PEFL2 (p,e,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
    &                  t,Dl,Dv,x,y,q,ierr,herr)
      subroutine PHFL2 (p,h,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,Dl,Dv,x,y,q,ierr,herr)
      subroutine PSFL2 (p,s,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
    &                  t,Dl,Dv,x,y,q,ierr,herr)
      subroutine TQFL2 (t,q,z,kq,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  p,Dl,Dv,x,y,ierr,herr)
      subroutine PQFL2 (p,q,z,kq,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,Dl,Dv,x,y,ierr,herr)

c  inputs--two of the following as indicated by the first two letters of c          the subroutine name:
c        t--temperature [K]
c        p--pressure [kPa]
c        D--overall (bulk) molar density [mol/L]
c        e--internal energy [J/mol]
c        h--enthalpy [J/mol]
c        s--entropy [[J/mol-K]
c
c  additional inputs
c        z--overall (bulk) composition [array of mol frac]
c     ksat--flag for bubble and dew point limits
c           0 = dew and bubble point limits computed within routine
c           1 = must provide values for following:


c  additional additional input--only for TQFL2 and PQFL2
c       kq--flag specifying units for input quality
c           kq = 1 quality on MOLAR basis [moles vapor/total moles]
c           kq = 2 quality on MASS basis [mass vapor/total mass]
c
c  additional inputs if ksat = 1
c     tbub--bubble point temperature [K] at (p,x=z)
c     tdew--dew point temperature [K] at (p,y=z)
c    --or--
c     pbub--bubble point pressure [kPa] at (t,x=z)
c     pdew--dew point pressure [kPa] at (t,y=z)
c    --and--
c    Dlbub--liquid density [mol/L] at bubble point
c    Dvdew--vapor density [mol/L] at dew point
c     ybub--vapor composition [array of mol frac] at bubble point
c     xdew--liquid composition [array of mol frac] at dew point
c
c  outputs--one of the following, depending on the inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c
c  additional outputs--common to all routines
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c           (not present for TQFL2 and PQFL2)
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)



      subroutine DQFL2 (d,q,z,kq,t,p,Dl,Dv,x,y,ierr,herr)
c
c  flash calculation given bulk density, quality, and composition
c
c  This routine accepts only two-phase states as input.
c
c  inputs:
c        d--overall (bulk) molar density [mol/L]
c        h--overall (bulk) molar enthalpy [J/mol]
c        z--overall (bulk) composition [array of mol frac]
c       kq--flag specifying units for input quality
c           kq = 1 quality on MOLAR basis [moles vapor/total moles]
c           kq = 2 quality on MASS basis [mass vapor/total mass]
c
c  outputs:
c        t--temperature [K]
c        p--pressure [kPa]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)



THERMODYNAMIC PROPERTY SUBROUTINES AS F(T,?,X)

The following routines provide thermodynamic properties as a function of
temperature, density, and composition.  Typically, one or more of these
will be called after finding the temperature and/or density with a call
to one of the saturation or flash routines.  Note that these routines
assume that valid inputs are supplied--no range checking is performed.

      subroutine CRITP (x,tcrit,pcrit,Dcrit,ierr,herr)
c
c  critical parameters as a function of composition
c
c  input:
c        x--composition [array of mol frac]