c  and the composition of one phase
c
c  inputs:
c        t--temperature [K]
c        x--composition [array of mol frac] (phase specified by kph)
c      kph--phase flag: 1 = input x is liquid composition (bubble point)
c                       2 = input x is vapor composition (dew point)
c                       3 = input x is liquid comp (freezing point)
c                       4 = input x is vapor comp (sublimation point)
c  outputs:
c        p--pressure [kPa]
c     rhol--molar density [mol/L] of saturated liquid
c     rhov--molar density [mol/L] of saturated vapor
c     xliq--liquid phase composition [array of mol frac]
c     xvap--vapor phase composition [array of mol frac]
c     ierr--error flag:   0 = successful
c                         1 = T < Tmin
c                         8 = x out of range
c                         9 = T and x out of range
c                       120 = CRITP did not converge
c                       121 = T > Tcrit
c                       122 = TPRHO-liquid did not converge (pure fluid)
c                       123 = TPRHO-vapor did not converge (pure fluid)
c                       124 = pure fluid iteration did not converge
c           following 3 error codes are advisory--iteration will either
c           converge on later guess or error out (ierr = 128)
c                      -125 = TPRHO did not converge for parent ph (mix)
c                      -126 = TPRHO did not converge for incipient (mix)
c                      -127 = composition iteration did not converge
c                       128 = mixture iteration did not converge
c     herr--error string if ierr<>0 (character*255)



      subroutine SATP (p,x,kph,t,rhol,rhov,xliq,xvap,ierr,herr)
c
c  iterate for saturated liquid and vapor states given temperature
c  and the composition of one phase
c
c  inputs:
c        p--pressure [kPa]
c        x--composition [array of mol frac] (phase specified by kph)
c      kph--phase flag:  1 = input x is liquid composition
c                        2 = input x is vapor composition
c
c  outputs:
c        t--temperature [K]
c     rhol--molar density [mol/L] of saturated liquid
c     rhov--molar density [mol/L] of saturated vapor
c     xliq--liquid phase composition [array of mol frac]
c     xvap--vapor phase composition [array of mol frac]
c     ierr--error flag:  0 = successful
c                        4 = P < 0
c                        8 = x out of range
c                       12 = P and x out of range
c                      140 = CRITP did not converge
c                      141 = P > Pcrit
c                      142 = TPRHO-liquid did not converge (pure fluid)
c                      143 = TPRHO-vapor did not converge (pure fluid)
c                      144 = pure fluid iteration did not converge
c           following 3 error codes are advisory--iteration will either
c           converge on later guess or error out (ierr = 148)
c                     -144 = Raoult's law (mixture intial guess) did
c                            not converge
c                     -145 = TPRHO did not converge for parent ph (mix)
c                     -146 = TPRHO did not converge for incipient (mix)
c                     -147 = composition iteration did not converge
c                      148 = mixture iteration did not converge
c     herr--error string if ierr<>0 (character*255)



      subroutine SATD (rho,x,kph,kr,t,p,rhol,rhov,xliq,xvap,ierr,herr)
c
c  iterate for temperature and pressure given a density along the
c  saturation boundary and the composition
c
c  inputs:
c      rho--molar density [mol/L]
c        x--composition [array of mol frac]
c      kph--flag specifying desired root for multi-valued inputs
c           has meaning only for water at temps close to its triple point
c          -1 = return middle root (between 0 and 4 C)
c           1 = return highest temperature root (above 4 C)
c           3 = return lowest temperature root (along freezing line)
c  outputs:
c        t--temperature [K]
c        p--pressure [kPa]
c     rhol--molar density [mol/L] of saturated liquid
c     rhov--molar density [mol/L] of saturated vapor
c     xliq--liquid phase composition [array of mol frac]
c     xvap--vapor phase composition [array of mol frac]
c       kr--phase flag: 1 = input state is liquid
c                       2 = input state is vapor in equilibrium with liq
c                       3 = input state is liq in equilibrium with solid
c                       4 = input state is vap in equilibrium with solid
c     ierr--error flag:   0 = successful
c                         2 = D > Dmax
c                         8 = x out of range
c                        10 = D and x out of range
c                       160 = CRITP did not converge
c                       161 = SATD did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  N.B. kr = 3,4 presently working only for pure components
c
c  either (rhol,xliq) or (rhov,xvap) will correspond to the input state
c  with the other pair corresponding to the other phase in equilibrium
c  with the input state

      subroutine SATH (h,x,kph,nroot,k1,t1,p1,d1,k2,t2,p2,d2,ierr,herr)
c
c  iterate for temperature, pressure, and density given an enthalpy along
c  the saturation boundary and the composition
c
c  inputs:
c        h--molar enthalpy [J/mol]
c        x--composition [array of mol frac]
c      kph--flag specifying desired root
c           0 = return all roots along the liquid-vapor line
c           1 = return only liquid VLE root
c           2 = return only vapor VLE roots
c           3 = return liquid SLE root (melting line)
c           4 = return vapor SVE root (sublimation line)
c  outputs:
c    nroot--number of roots.  Set to one for kph=1,3,4 if ierr=0
c       k1--phase of first root (1-liquid, 2-vapor, 3-melt, 4-subl)
c       t1--temperature of first root [K]
c       p1--pressure of first root [kPa]
c       dl--molar density of first root [mol/L]
c       k2--phase of second root (1-liquid, 2-vapor, 3-melt, 4-subl)
c       t2--temperature of second root [K]
c       p2--pressure of second root [kPa]
c       d2--molar density of second root [mol/L]
c     ierr--error flag:   0 = successful
c                         2 = h < hmin
c                         4 = h > hmax
c                         8 = h > htrp (for subl input)
c                       160 = CRITP did not converge
c                       161 = SATH did not converge for one root
c                       162 = SATH did not converge for both roots
c     herr--error string (character*255 variable if ierr<>0)
c
c  The second root is always set as the root in the vapor at temperatures
c  below the maximum enthalpy on the vapor saturation line.  If kph = 2,
c  and only one root is found in the vapor (this occurs when h<hcrit)
c  the state point will be placed in k2,t2,p2,d2.  If kph = 0 and this
c  situation occured, the first root (k1,t1,p1,d1) would be in the liquid
c  (k1 =1, k2 = 2).
c
c  N.B. kph = 3,4 presently working only for pure components

      subroutine SATS (s,x,kph,nroot,k1,t1,p1,d1,k2,t2,p2,d2,
     &                 k3,t3,p3,d3,ierr,herr)
c
c  iterate for temperature, pressure, and density given an entropy along
c  the saturation boundary and the composition
c
c  inputs:
c        s--molar entropy [J/mol-K]
c        x--composition [array of mol frac]
c      kph--flag specifying desired root
c           0 = return all roots along the liquid-vapor line
c           1 = return only liquid VLE root
c           2 = return only vapor VLE roots
c           3 = return liquid SLE root (melting line)
c           4 = return vapor SVE root (sublimation line)
c  outputs:
c    nroot--number of roots.  Set to one for kph=1,3,4 if ierr=0
c       k1--phase of first root (1-liquid, 2-vapor, 3-melt, 4-subl)
c       t1--temperature of first root [K]
c       p1--pressure of first root [kPa]
c       dl--molar density of first root [mol/L]
c       k2--phase of second root (1-liquid, 2-vapor, 3-melt, 4-subl)
c       t2--temperature of second root [K]
c       p2--pressure of second root [kPa]
c       d2--molar density of second root [mol/L]
c       k3--phase of thrid root (1-liquid, 2-vapor, 3-melt, 4-subl)
c       t3--temperature of third root [K]
c       p3--pressure of third root [kPa]
c       d3--molar density of third root [mol/L]
c     ierr--error flag:   0 = successful
c                         2 = s < smin
c                         4 = s > smax
c                         8 = s > strp (for subl input)
c                       160 = CRITP did not converge
c                       161 = SATS did not converge for one root
c                       162 = SATS did not converge for two roots
c                       163 = SATS did not converge for all roots
c     herr--error string (character*255 variable if ierr<>0)
c
c  The second root is always set as the root in the vapor at temperatures
c  below the maximum entropy on the vapor saturation line.  If kph = 2,
c  and only one root is found in the vapor (this occurs when s<scrit)
c  the state point will be placed in k2,t2,p2,d2.  If kph = 0 and this
c  situation occured, the first root (k1,t1,p1,d1) would be in the liquid
c  (k1=1, k2=2).
c
c  The third root is the root with the lowest temperature.
c  For fluids with multiple roots:  When only one root is found in the
c  vapor phase (this happens only at very low temperatures past the
c  region where three roots are located), the value of the root is still
c  placed in k3,t3,p3,d3.  For fluids that never have more than one root
c  (when there is no maximum entropy along the saturated vapor line), the
c  value of the root is always placed in k1,t1,p1,d1.
c
c  N.B. kph = 3,4 presently working only for pure components




      subroutine CSATK (icomp,t,kph,p,rho,csat,ierr,herr)
c
c  compute the heat capacity along the saturation line as a function of
c  temperature for a given component
c
c  csat can be calculated two different ways:
c     Csat = Cp - T(DvDT)(DPDTsat)
c     Csat = Cp - beta/rho*hvap/(vliq - vvap)
c     where beta is the volulme expansivity
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      kph--phase flag: 1 = liquid calculation
c                       2 = vapor calculation
c  outputs:
c        p--saturation pressure [kPa]
c      rho--saturation molar density [mol/L]
c     csat--saturation heat capacity [J/mol-K]


Two similar routines are provided for calculating surface tension.
SURTEN is more efficient if the liquid and vapor density and composition
are known (e.g. from a previous call to SATT).  If these are not known,
then SURFT may be used.


      subroutine SURFT (t,rhol,xl,sigma,ierr,herr)
c
c  compute surface tension
c
c  inputs:
c        t--temperature [K]
c       xl--composition of liquid phase [array of mol frac]
c  outputs:
c     rhol--molar density of liquid phase [mol/L]
c    sigma--surface tension [N/m]
c     ierr--error flag:   0 = successful
c                         1 = T < Tmin
c                         8 = x out of range
c                         9 = T and x out of range
c                       120 = CRITP did not converge
c                       121 = T > Tcrit
c                       122 = TPRHO-liquid did not converge in SATT
c                       123 = TPRHO-vapor did not converge in SATT
c                       124 = SATT pure fluid iteration did not converge
c                       128 = SATT mixture iteration did not converge
c     herr--error string if ierr<>0 (character*255)


      subroutine SURTEN (t,rhol,rhov,xl,xv,sigma,ierr,herr)
c
c  compute surface tension
c
c  inputs:
c        t--temperature [K]
c     rhol--molar density of liquid phase [mol/L]
c     rhov--molar density of vapor phase [mol/L]
c           if either rhol or rhov < 0 call SATT to find densities
c       xl--composition of liquid phase [array of mol frac]
c       xv--composition of liquid phase [array of mol frac]
c           (xv is optional input if rhol < 0 or rhov < 0)
c  outputs:
c    sigma--surface tension [N/m]
c     ierr--error flag:   0 = successful
c                         [all error codes identical to those for SURFT]
c     herr--error string if ierr<>0 (character*255)

FLASH SUBROUTINES

So-called "flash" calculations involve a determination of the
thermodynamic state given two independent variables plus composition.
In addition to the inputs of temperature and pressure which is the most
common flash calculation, REFPROP provides routines for virtually all
combinations of temperature, pressure, or density as the first variable
and density, internal energy, enthalpy, entropy, or quality as the
second variable.  The combination of enthalpy and entropy is also
supported.

Flash subroutines are provided for cases where the state is known to be
single phase (liquid or vapor), known to be two-phase (liquid plus
vapor), and also for the general case where the phase is not known.
Because of the many combinations and their parallel structure, these
routines are described in groups.  The first two letters of the
subroutine name indicate the independent variables, where

      T = temperature [K]
      P = pressure [kPa]
      D = density [mol/L]
      E = internal energy [J/mol]
      H = enthalpy [J/mol]
      S = entropy [J/mol-K]
      Q = vapor quality [moles vapor/total moles]
          or [kg vapor/total kg] depending on the
          value of the input flag kq


GENERAL FLASH SUBROUTINES

For cases where the phase is not known, the following routines are available.

      subroutine TPFLSH (t,p,z,D,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
      subroutine TDFLSH (t,D,z,p,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
      subroutine TEFLSH (t,e,z,kr,p,D,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
      subroutine THFLSH (t,h,z,kr,p,D,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr)
      subroutine TSFLSH (t,s,z,kr,p,D,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr)
      subroutine PDFLSH (p,D,z,t,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
      subroutine PEFLSH (p,e,z,t,D,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
      subroutine PHFLSH (p,h,z,t,D,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr)
      subroutine PSFLSH (p,s,z,t,D,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr)
      subroutine DHFLSH (D,h,z,t,p,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr)
      subroutine DSFLSH (D,s,z,t,p,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr)
      subroutine DEFLSH (D,e,z,t,p,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
      subroutine HSFLSH (h,s,z,t,p,D,Dl,Dv,x,y,q,e,cv,cp,w,ierr,herr)
      subroutine TQFLSH (t,q,z,kq,p,D,Dl,Dv,x,y,e,h,s,cv,cp,w,ierr,herr)
      subroutine PQFLSH (p,q,z,kq,t,D,Dl,Dv,x,y,e,h,s,cv,cp,w,ierr,herr)
c
c  flash calculation given two independent variables and bulk composition
c