c  begin file ftn_pas.f
c
c  This file contains the Fortran routine which serves as the link to
c  the Pascal-based graphical interface, along with subsidiary routines.
c
c  contained here are:
c     subroutine REFP32 (cstr,mode32,nin32,xin,nout32,xout)
c     subroutine REFPRP6 (cstr,mode,nin,xin,nout,xout)
c     subroutine FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c     subroutine SATOUT (tl,tv,pl,pv,Dl,Dv,xl,xv,xout)
c     subroutine VLOUT (tl,tv,pl,pv,Dl,Dv,xl,xv,xout)
c     subroutine SPOUT (t,D,x,kph,xout)
c  following routines are analogous to above, but with arguments
c  rearranged for compatibility with version 6.1
c     subroutine REFP61 (mode,xin,intin,cstrin,xblk,xliq,xvap,ierri,
c    &                   nwarn,iwarn,cstr)
c     subroutine FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,
c    &                   xblk,xliq,xvap)
c     subroutine STOUT2 (tl,tv,pl,pv,Dl,Dv,xl,xv,xblk,xliq,xvap)
c     subroutine VLOUT2 (tl,tv,pl,pv,Dl,Dv,xl,xv,xliq,xvap)
c     subroutine SPOUT2 (t,D,x,kph,xblk,xliq,xvap)
c
c ======================================================================
c ======================================================================
c
      subroutine REFP32 (cstr,mode32,nin32,xin,nout32,xout)
c
c  This subroutine is the main link between the Fortran property routines
c  and the outside world via a DLL.  This routine is set up for a 32-bit
c  DLL.  It calls the original "interface" routine REFPRP6, which was set
c  up for a 16-bit DLL.
c
c  inputs:
c    see comments for REFPRP6, immediately following
c  outputs:
c    see comments for REFPRP6, immediately following
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  06-10-97  MM, original version; based on VLOUT
c  06-16-97  MM, corrections to MS$ATTRIBUTES for use with DEC Fortran
c
c  compiler switches for use with Digital Visual Fortran; these should be 
c  treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: REFP32
!MS$ATTRIBUTES ALIAS :'REFP32' :: REFP32
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      integer*2 mode,nin,nout
      integer*4 mode32,nin32,nout32
      character*255 cstr
      dimension xin(100),xout(100)
c
      mode=mode32
      nin=nin32
      call REFPRP6 (cstr,mode,nin,xin,nout,xout)
      nout32=nout
c
      RETURN
      end                                             !subroutine REFP32
c
c ======================================================================
c
      subroutine REFPRP6 (cstr,mode,nin,xin,nout,xout)
c
c  This subroutine is the main link between the Fortran property
c  routines and the Pascal-based graphical user interface.  It calls the
c  appropriate property routine(s) and returns values to the interface.
c
c  inputs:
c     cstr--used to pass text info (e.g. file names) to the Fortran
c           property routines, (applicable only for certain modes)
c           [character*255 with last character the ASCII null (char(0))]
c        if mode = -4:  cstr(1:3) = mixing rule code (e.g. 'VDW')
c                       cstr(4:83) = file containing mixture data
c           mode = -3:  cstr(1:3) = code for overall/mix EOS (e.g. 'HMX')
c                       cstr(4:6) = code for component 1 EOS (e.g. 'BWR')
c                       cstr(7:9) = etc. for each of the nc components
c           mode = -1:  cstr consists of various data separated by '|'
c                       --path to the fluids directory (e.g. C:\REFPROP\fluids\)
c                       --file containing mixture data (e.g. HMX.bnc)
c                       --file(s) containing fluid data for each of the nc
c                         components (e.g. R134a.fld)
c                       for example, to specify a mixture of R32 and
c                       R134a, with mixture data in file HMX.bnc, this
c                       string might be:
c                cstr = 'C:\REFPROP\fluids\|HMX.bnc|R32.fld|R134a.fld|'//char(0)
c     mode--specifies the type of calculation to be done
c           mode = -5:  get mixture interaction coef      (call GETKTV)
c                  -4:  set mixture interaction coef      (call SETKTV)
c                  -3:  specify non-standard model        (call SETMOD)
c                  -2:  specify non-standard ref state    (call SETREF)
c                  -1:  specify fluid(s), ref state, etc. (call SETUP)
c                   0:  calculate properties as specified by xin(1)
c      nin--number of meaningful elements in the 'xin' array
c           this is the number of independent variables (e.g. T,P,rho,h)
c           plus the number of components in the mixture
c           (never used in practice, but retained for compatibility)
c      xin--input variables [array of dimension 100]
c           xin(1): code specifying calculation (with routine noted)
c                   1:  thermo props as f(T,P,x)          (call TPFLSH)
c                   2:  thermo props as f(T,rho,x)        (call TDFLSH)
c                   3:  thermo props as f(P,h,x)          (call PHFLSH)
c                   4:  thermo props as f(P,s,x)          (call PSFLSH)
c                   5:  thermo props as f(T,rho,x)        (call THERM)
c                       N.B. mode = 5 does no error checking and should be used
c                            only when conditions are known to be single-phase
c                   6:  satn thermo props as f(T,xl=xv)   (call SATT)
c                   7:  satn thermo props as f(P,xl=xv)   (call SATP)
c                   8:  critical parameters as f(x)       (call CRITP)
c                   9:  transport propertes as f(t,rho,x) (call TRNPRP)
c                  10:  fugacities as f(t,rho,x)          (call FGCTY)
c                  11:  surface tension as f(t,rho,x)     (call SURFT)
c                  12:  bubble point props as f(t,xl)     (call SATT)
c                  13:  dew point props as f(t,xv)        (call SATT)
c                  14:  bubble point props as f(p,xl)     (call SATP)
c                  15:  dew point props as f(p,xv)        (call SATP)
c                  16:  thermo props as f(T,q,x)          (call TQFLSH)
c                  17:  thermo props as f(p,q,x)          (call PQFLSH)
c                  18:  molecular weight as f(x)          (wm = XMOL)
c                  19:  convert mole to mass fraction     (call XMASS)
c                  20:  convert mass to mole fraction     (call XMOLE)
c                  21:  CAS # and name for component i    (call NAME)
c                  22:  fixed properties for component i  (call INFO)
c                  23:  limits of specified model as f(x) (call LIMITX)
c                  24:  limits for a single component     (call LIMITK)
c                  25:  single-phase T-P flash            (call TPRHO)
c                  26:  single-phase P-h flash            (call PHFL1)
c                  27:  single-phase P-s flash            (call PSFL1)
c                  31:  dP/d(rho) as f(T,rho,x)           (call DPDD)
c                  32:  d2P/d(rho)2 as f(T,rho,x)         (call DPDD2)
c                  33:  dP/dT as f(T,rho,x)               (call DPDT)
c                  34:  d(rho)/dP as f(T,rho,x)           (call DDDP)
c                  35:  d(rho)/dT as f(T,rho,x)           (call DDDT)
c         xin(2..nin):  independent variables in order as noted above
c                       e.g. T,P,x(1)..x(nc)   for xin(1) = 1
c                       e.g. T,rho,x(1)..x(nc) for xin(1) = 2,9,10,11,31-35
c                       e.g. T,x(1)..x(nc)     for xin(1) = 6
c                       e.g. x(1)..x(nc)       for xin(1) = 8,18,19,20
c
c      but, in the case of the setup calls (mode < 0) xin is different:
c        if mode = -5:  xin(1) = component i
c                       xin(2) = component j (for the pair i,j)
c        if mode = -4:  xin(1) = component i
c                       xin(2) = component j (for the pair i,j)
c                       xin(3..6) = mixing coefficients 1..4
c        if mode = -3:  xin(1) = number of component(s) in mixture
c                       xin(2) = code for type of model to be set
c                                1:  equation of state ('EOS')
c                                2:  viscosity ('ETA')
c                                3:  thermal conductivity ('TCX')
c                                4:  surface tension ('STN')
c        if mode = -2:  xin(1) = code for reference state
c                                1:  normal boiling point ('NBP')
c                                2:  ASHRAE convention ('ASH')
c                                3:  IIR convention ('IIR')
c                                4:  other ('OTH')
c                       xin(2) = composition flag
c                                1:  ref state applied to pure components
c                                2:  ref state applied to mixture of
c                                    composition specified by xin(8..7+nc)
c                       xin(3) = Tref
c                       xin(4) = Pref
c                       xin(5) = href
c                       xin(6) = sref
c                       xin(7) = number of component(s)
c                       xin(8..7+nc) = composition(s)
c        if mode = -1:  xin(1) = number of component(s) in mixture
c                       xin(2) = code for reference state
c                                1:  normal boiling point ('NBP')
c                                2:  ASHRAE convention ('ASH')
c                                3:  IIR convention ('IIR')
c                                4:  other ('OTH') [must use mode = -2]
c                       xin(3) = code for equation of state [obsolete]
c      in the case of the NAME and INFO calls (mode = 21, 22) xin is different:
c        if xin(1) = 21:  xin(2) = component i
c        if xin(1) = 22:  xin(2) = component i
c      in the case of the call to LIMITX xin is:
c        if xin(1) = 23:  xin(2) = code for type of model to be set
c                                1:  equation of state ('EOS')
c                                2:  viscosity ('ETA')
c                                3:  thermal conductivity ('TCX')
c                                4:  surface tension ('STN')
c                         xin(3) = temperature of mixture
c                         xin(4) = density of mixture
c                         xin(5) = pressure of mixture
c                         xin(6..5+nc) = composition(s)
c      in the case of the call to LIMITK xin is:
c        if xin(1) = 24:  xin(2) = code for type of model to be set
c                                1:  equation of state ('EOS')
c                                2:  viscosity ('ETA')
c                                3:  thermal conductivity ('TCX')
c                                4:  surface tension ('STN')
c                         xin(3) = component number
c                         xin(4) = temperature of mixture
c                         xin(5) = density of mixture
c                         xin(6) = pressure of mixture
c      in the case of the call to TPRHO xin is:
c        if xin(1) = 25:  xin(2) = temperature
c                         xin(3) = pressure
c                         xin(4..3+nc) = composition(s)
c                         xin(4+nc) = phase flag:  1 = find liquid solution
c                                                  2 = find vapor solution
c                         xin(5+nc) = guess flag:  1 = initial guess for density provided
c                         xin(6+nc) = initial guess for density (ignored if xin(5+nc)<>1)
c
c  outputs:
c     cstr--used to pass error/warning messages to GUI,
c           except for:
c             mode = -5, cstr(1:3) = mixing rule code (e.g. 'VDW')
c                        cstr(4..) = documentation for the binary parameters
c           xin(1) = 21, cstr(1:12)   = short name
c                        cstr(13:92)  = full name
c                        cstr(93:104) = CAS number
c     mode--error/warning message code (0 indicates successful call)
c     nout--number of meaningful elements in the 'xout' array
c     xout--calculated (output) values
c           there are 15 possible properties + the compositions, for most
c           cases, they are returned in the xout array in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c           the first 15 + nc elements are the bulk properties
c           the next 15 + nc elements are liquid-phase properties
c           the next 15 + nc elements are vapor-phase properties
c
c           but there are exceptions to the above order:
c               mode = -1,  SETUP:  xout(1) = ncomp
c               mode = -2, SETREF:  nothing in xout
c               mode = -3, SETMOD:  nothing in xout
c               mode = -4, SETKTV:  nothing in xout
c               mode = -5, GETKTV:  xout(1..4) = mixing coefficients 1..4
c               mode = 21,   NAME:  nothing in xout
c             xin(1) =  8,  CRITP:  xout(1) = Tcrit
c                                   xout(2) = Pcrit
c                                   xout(3) = Dcrit
c             xin(1) =  9, TRNPRO:  xout(1) = eta (viscosity)
c                                   xout(2) = tcx (thermal cond)
c             xin(1) = 10,  FGCTY:  xout(1..nc) = the nc fugacities
c             xin(1) = 11,  SURFT:  xout(1) = the surface tension [N/m]
c             xin(1) = 18,   WMOL:  xout(1) = the molecular weight [g/mol]
c             xin(1) = 19,  XMASS:  xout(1..nc) = the composition on basis opposite of input
c                      20,  XMOLE:  xout(nc+1)  = the molecular weight [g/mol]
c             xin(1) = 21,   NAME:  nothing in xout (all info in string cstr)
c             xin(1) = 22,   INFO:  xout(1) = molecular weight [g/mol]
c                                   xout(2) = triple point temperature [K]
c                                   xout(3) = normal boiling point temperature [K]
c                                   xout(4) = critical temperature [K]
c                                   xout(5) = critical pressure [kPa]
c                                   xout(6) = critical density [mol/L]
c                                   xout(7) = compressibility at critical point [pc/(R*Tc*Dc)]
c                                   xout(8) = accentric factor [-]
c                                   xout(9) = dipole moment [debye]
c                                  xout(10) = gas constant [J/mol-K]
c             xin(1) = 23, LIMITX:  xout(1) = lower temperature limit [K]
c                                   xout(2) = upper temperature limit [K]
c                                   xout(3) = upper density limit [mol/L]
c                                   xout(4) = upper pressure limit [kPa]
c             xin(1) = 24, LIMITK:  same as for xin(1) = 23, LIMITX, except
c                                   limits valid only for specified component
c             xin(1) = 31, DPDD:    xout(1) = the derivative dP/d(rho)
c                      32, DPDD2:   xout(1) = the derivative d2P/d(rho)2
c                      33, DPDT:    xout(1) = the derivative dP/dT
c                      34, DDDP:    xout(1) = the derivative d(rho)/dP
c                      35, DDDT:    xout(1) = the derivative d(rho)/dT
c
c           but note that viscosity, thermal cond, and surface tension
c           are not calculated by the thermodynamic routines and these
c           properties are returned in different locations by the
c           TRNPRP and SURFT routines, so that elements 12, 13, 14 in
c           the 'normal' output array are always set to -9.999d6 (flag
c           for an uncalculated property)
c
c  written by S.A. Klein, NIST Thermophysics Division, Boulder, Colorado
c  (8/95 through 7/96)
c  07-31-96  MM, document inputs, outputs
c  08-02-96  MM, major restructuring of code
c  08-06-96  MM, if error on SATT/SATP set mode=ierr before RETURN
c  11-14-96  MM, change nmxpar from 4 to 6
c  02-19-97  MM, add calls to NAME, INFO, direct call to THERM
c  03-25-97  MM, add calls to LIMITX, LIMITK, XMASS, XMOLE
c  03-26-97  MM, add DEF (default) option for reference states
c  04-18-97  MM; fix bug in mass<->mole conversion (modes 19,20)
c  05-16-97  MM, get special-case flags from /FLAGS/
c  06-06-97  MM, add calls to single-phase T-p, p-h, p-s flashes
c  08-20-97  MM, if mode out of range, also write value of nin--this avoids
c                compiler warning that nin not used
c  08-29-97  MM, add modes 31 - 35 for various derivatives
c  12-02-97  MM, initialize ierr=0, herr=hnull for modes which do not return errors
c                move definition of hrefs and htypes to BDDLL
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (nmxpar=6)       !number of binary mixture parameters
      parameter (nrefst=6)       !number of reference state choices
      integer*2 mode,nin,nout
      character*1 htab,hnull
      character*3 hrf,hcomp(ncmax),hmodij
      character*3 hrefs,htypes,htype,hmix
      character*8  hfij(nmxpar)
      character*12 hcas,hname
      character*80 hfmix,hfiles(ncmax),hdir,hmxrul,hn80
      character*255 cstr,hbinp
      dimension xin(100),xout(100)
      dimension x(ncmax),xl(ncmax),xv(ncmax),xjunk(ncmax),
     &          fug(ncmax),fij(nmxpar),xkg(ncmax),xmol(ncmax)
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
c  pointers to the xout array (to aid in adding possible future props)
      common /IPOINT/
     &     it,ip,irho,ivol,ie,ih,is,icv,icp,iw,iq,ieta,itcx,ist,iz,ix
c  maximum number of properties returned to the GUI on a single call
      common /NPROP/ nprops
c  common block containing flags to GUI (initialized in BDSET in setup.f)
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c  reference state choices and model types
      common /HREFST/ hrefs(nrefst)
      common /HMODTY/ htypes(4)
c
d     write (*,1080) (xin(i),i=1,nin)
d1080 format (1x,'REFPRP--input array:  ',20f14.6)
d     write (*,1255) cstr
d1255 format (1x,'REFPRP--input C_str:  (',a255,')')
c
c  load the output array with flag signifying 'not calculated'; the
c  GUI will crash if any elements in the output array are not defined;
c  this is insurance in case any of the property routines crash
      nout=3*(nprops+nc)
      do 80 i=1,nout
  80  xout(i)=xnotc
c
      if (mode.ge.0) then
        imode=NINT(xin(1))       !mode is passed in real xin array
c
        if (imode.eq.1) then
c
c  temperature-pressure flash (call TPFLSH)
c
          t=xin(2)
          p=xin(3)
          do 110 i=1,nc
  110     x(i)=xin(3+i)
          call TPFLSH (t,p,x,D,Dl,Dv,xl,xv,q,e,h,s,cv,cp,w,ierr,cstr)
          nout=3*(nprops+nc)
          call FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c
        else if (imode.eq.2) then
c
c  temperature-density flash (call TDFLSH)
c
          t=xin(2)
          D=xin(3)
          do 120 i=1,nc
  120     x(i)=xin(3+i)
          call TDFLSH (t,D,x,p,Dl,Dv,xl,xv,q,e,h,s,cv,cp,w,ierr,cstr)
          nout=3*(nprops+nc)
          call FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c
        else if (imode.eq.3) then
c
c  pressure-enthalpy flash (call PHFLSH)
c
          p=xin(2)
          h=xin(3)
          do 130 i=1,nc
  130     x(i)=xin(3+i)
          call PHFLSH (p,h,x,t,D,Dl,Dv,xl,xv,q,e,s,cv,cp,w,ierr,cstr)
          nout=3*(nprops+nc)
          call FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c
        else if (imode.eq.4) then
c
c  pressure-entropy flash (call PSFLSH)
c
          p=xin(2)
          s=xin(3)
          do 140 i=1,nc
  140     x(i)=xin(3+i)
          call PSFLSH (p,s,x,t,D,Dl,Dv,xl,xv,q,e,h,cv,cp,w,ierr,cstr)
          nout=3*(nprops+nc)
          call FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c
        else if (imode.eq.5) then
c
c  temperature-density calculation (call THERM)
c  N.B.  this option does not check for two-phase states
c        it will return nonsense values for meta-stable states
c
          t=xin(2)
          D=xin(3)
          do 150 i=1,nc
          x(i)=xin(3+i)
          xl(i)=x(i)
          xv(i)=x(i)
  150     continue
          call THERM (t,D,x,p,e,h,s,cv,cp,w,hjt)
c  quality of 999 flags FLOUT to set liq & vap properties to bulk
          q=999.0d0
          Dl=D
          Dv=D
          call FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
          nout=3*(nprops+nc)
c  this mode does not return any errors
          ierr=0
          cstr=hnull
c
        else if (imode.eq.6) then
c
c  saturation state at specified t, xl=xv (call SATT)
c
          t=xin(2)
          do 160 i=1,nc
  160     x(i)=xin(2+i)
c  call SATT at both dew point and bubble point
          call SATT (t,x,1,pl,Dl,Djunk,xl,xjunk,ierr,cstr)
          if (ierr.gt.0) then
            mode=ierr
            RETURN
          end if
          call SATT (t,x,2,pv,Djunk,Dv,xjunk,xv,ierr,cstr)
          if (ierr.gt.0) then
            mode=ierr
            RETURN
          end if
          nout=3*(nprops+nc)
          tv=t
          call SATOUT (t,tv,pl,pv,Dl,Dv,xl,xv,xout)
          xout(it)=t
c
        else if (imode.eq.7) then
c
c  saturation state at specified p, xl=xv (call SATP)
c
          p=xin(2)
          do 170 i=1,nc
  170     x(i)=xin(2+i)
c  call SATP at both dew point and bubble point
          call SATP (p,x,1,tl,Dl,Djunk,xl,xjunk,ierr,cstr)
          if (ierr.gt.0) then
            mode=ierr
            RETURN
          end if
          call SATP (p,x,2,tv,Djunk,Dv,xjunk,xv,ierr,cstr)
          if (ierr.gt.0) then
            mode=ierr
            RETURN
          end if
          nout=3*(nprops+nc)
          pv=p
          call SATOUT (tl,tv,p,pv,Dl,Dv,xl,xv,xout)
          xout(ip)=p
c
        else if (imode.eq.8) then
c
c  calculate critical parameters (call CRITP)
c
          do 180 i=1,nc
  180     x(i)=xin(1+i)
          call CRITP (x,tc,pc,Dc,ierr,cstr)
          nout=3
          xout(1)=tc
          xout(2)=pc
          xout(3)=Dc
c
        else if (imode.eq.9) then
c
c  transport properties (viscosity, thermal conductivity) (call TRNPRP)
c
          t=xin(2)
          rho=xin(3)
          do 190 i=1,nc
  190     x(i)=xin(3+i)
          call TRNPRP (t,rho,x,eta,tcx,ierr,cstr)
          nout=2
          xout(1)=eta
          xout(2)=tcx
c
        else if (imode.eq.10) then
c
c  calculate fugacity for each of the nc components (call FGCTY)
c
          t=xin(2)
          rho=xin(3)
          do 200 i=1,nc
  200     x(i)=xin(3+i)
          call FGCTY (t,rho,x,fug)
          nout=nc
          ierr=0
          cstr=hnull
          do 201 i=1,nc
          xout(i)=fug(i)
  201     continue
c
        else if (imode.eq.11) then
c
c  surface tension (call SURFT)
c
          t=xin(2)
          rho=xin(3)
          do 210 i=1,nc
  210     x(i)=xin(3+i)
          call SURFT (t,rho,x,sigma,ierr,cstr)
          nout=1
          xout(1)=sigma
c
        else if (imode.ge.12 .and. imode.le.15) then
c
c  saturation state with xl and xv in equilibrium
c
          do 220 i=1,nc
  220     x(i)=xin(2+i)
c
          if (imode.eq.12) then
c
c  saturation state at bubble point
c  specified t,xl; xv is equilibium composition  (call SATT)
c
            t=xin(2)
            call SATT (t,x,1,p,Dl,Dv,xl,xv,ierr,cstr)
            if (ierr.gt.0) then
              mode=ierr
              RETURN
            end if
c
          else if (imode.eq.13) then
c
c  saturation state at dew point
c  specified t,xv; xl is equilibium composition  (call SATT)
c
            t=xin(2)
            call SATT (t,x,2,p,Dl,Dv,xl,xv,ierr,cstr)
            if (ierr.gt.0) then
              mode=ierr
              RETURN
            end if
c
          else if (imode.eq.14) then
c
c  saturation state at bubble point
c  specified p,xl; xv is equilibium composition  (call SATP)
c
            p=xin(2)
            call SATP (p,x,1,t,Dl,Dv,xl,xv,ierr,cstr)
            if (ierr.gt.0) then
              mode=ierr
              RETURN
            end if
c
          else if (imode.eq.15) then
c
c  saturation state at dew point
c  specified p,xv; xl is equilibium composition  (call SATP)
c
            p=xin(2)
            call SATP (p,x,2,t,Dl,Dv,xl,xv,ierr,cstr)
            if (ierr.gt.0) then
              mode=ierr
              RETURN
            end if
c
          end if
          nout=3*(nprops+nc)
          tl=t
          tv=t
          pl=p
          pv=p
          call SATOUT (tl,tv,pl,pv,Dl,Dv,xl,xv,xout)
          xout(it)=t
          xout(ip)=p
c
        else if (imode.eq.16) then
c
c  temperature-quality flash (call TQFLSH)
c
          t=xin(2)
          q=xin(3)
          do 260 i=1,nc
  260     x(i)=xin(3+i)
          kq=2    !input quality is on mass basis
          call TQFLSH (t,q,x,kq,p,D,Dl,Dv,xl,xv,e,h,s,cv,cp,w,ierr,cstr)
          nout=3*(nprops+nc)
          call FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c
        else if (imode.eq.17) then
c
c  pressure-quality flash (call PQFLSH)
c
          p=xin(2)
          q=xin(3)
          do 270 i=1,nc
  270     x(i)=xin(3+i)
          kq=2    !input quality is on mass basis
          call PQFLSH (p,q,x,kq,t,D,Dl,Dv,xl,xv,e,h,s,cv,cp,w,ierr,cstr)
          nout=3*(nprops+nc)
          call FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c
        else if (imode.eq.18) then
c
c  calculate molecular weight (wm = WMOL)
c
          do 280 i=1,nc
  280     x(i)=xin(1+i)
          nout=1
          ierr=0
          cstr=hnull
          xout(1)=WMOL(x)
          nout=1
c
        else if (imode.eq.19) then
c
c  convert composition in mole fraction to mass fraction (call XMASS)
c
          do 290 i=1,nc
  290     xmol(i)=xin(1+i)
          nout=nc+1
          ierr=0
          cstr=hnull
          call XMASS (xmol,xkg,wmix)
          do 292 i=1,nc
  292     xout(i)=xkg(i)
          xout(nc+1)=wmix
c
        else if (imode.eq.20) then
c
c  convert composition in mass fraction to mole fraction (call XMOLE)
c
          do 300 i=1,nc
  300     xkg(i)=xin(1+i)
          nout=nc+1
          ierr=0
          cstr=hnull
          call XMOLE (xkg,xmol,wmix)
          do 302 i=1,nc
  302     xout(i)=xmol(i)
          xout(nc+1)=wmix
c
        else if (imode.eq.21) then
c
c  return CAS and names for a component (call NAME)
c
          icomp=NINT(xin(2))
          call NAME (icomp,hname,hn80,hcas)
          cstr(1:12)=hname
          cstr(13:92)=hn80
          cstr(93:104)=hcas
          cstr(105:105)=hnull
          ierr=0
          nout=0
c
        else if (imode.eq.22) then
c
c  return various fixed points for a component (call INFO)
c
          icomp=NINT(xin(2))
          call INFO (icomp,wm,ttp,tnbp,tc,pc,Dc,Zc,acf,dip,Rgas)
          xout(1)=wm
          xout(2)=ttp
          xout(3)=tnbp
          xout(4)=tc
          xout(5)=pc
          xout(6)=Dc
          xout(7)=Zc
          xout(8)=acf
          xout(9)=dip
          xout(10)=Rgas
          ierr=0
          nout=10
          cstr=hnull
c
        else if (imode.eq.23) then
c
c  return limits to specified model (call LIMITX)
c
          itype=NINT(xin(2))
          htype=htypes(itype)    !code for model (see data stmt)
          t=xin(3)
          D=xin(4)
          p=xin(5)
          do 330 i=1,nc
  330     x(i)=xin(5+i)
          call LIMITX (htype,t,D,p,x,tmin,tmax,Dmax,pmax,ierr,cstr)
          nout=4
          xout(1)=tmin
          xout(2)=tmax
          xout(3)=Dmax
          xout(4)=pmax
c
        else if (imode.eq.24) then
c
c  return limits to specified model for a given component (call LIMITK)
c
          itype=NINT(xin(2))
          htype=htypes(itype)    !code for model (see data stmt)
          icomp=NINT(xin(3))
          t=xin(4)
          D=xin(5)
          p=xin(6)
          call LIMITK (htype,icomp,t,D,p,tmin,tmax,Dmax,pmax,ierr,cstr)
          nout=4
          xout(1)=tmin
          xout(2)=tmax
          xout(3)=Dmax
          xout(4)=pmax
c
        else if (imode.eq.25) then
c
c  single-phase temperature-pressure flash (call TPRHO)
c  N.B.  use TPFLSH (mode = 1) when phase is not known
c
          t=xin(2)
          p=xin(3)
          do 350 i=1,nc
  350     x(i)=xin(3+i)
          kph=NINT(xin(4+nc))
          kguess=NINT(xin(5+nc))
          if (kguess.eq.1) D=xin(6+nc)
          call TPRHO (t,p,x,kph,kguess,D,ierr,cstr)
          nout=3*(nprops+nc)
          call SPOUT (t,D,x,kph,xout)
c
        else if (imode.eq.26) then
c
c  single-phase pressure-entropy flash (call PHFL1)
c  N.B.  use PHFLSH (mode = 3) when phase is not known
c
          p=xin(2)
          h=xin(3)
          do 360 i=1,nc
  360     x(i)=xin(3+i)
          kph=NINT(xin(4+nc))
          t=xin(5+nc)             !initial guess for temperature
          D=xin(6+nc)             !initial guess for density
          call PHFL1 (p,h,x,kph,t,D,ierr,cstr)
          nout=3*(nprops+nc)
          call SPOUT (t,D,x,kph,xout)
c
        else if (imode.eq.27) then
c
c  single-phase pressure-entropy flash (call PSFL1)
c  N.B.  use PSFLSH (mode = 4) when phase is not known
c
          p=xin(2)
          s=xin(3)
          do 370 i=1,nc
  370     x(i)=xin(3+i)
          kph=NINT(xin(4+nc))
          t=xin(5+nc)             !initial guess for temperature
          D=xin(6+nc)             !initial guess for density
          call PSFL1 (p,s,x,kph,t,D,ierr,cstr)
          nout=3*(nprops+nc)
          call SPOUT (t,D,x,kph,xout)
c
        else if (imode.eq.31) then
c
c  derivative of pressure w.r.t. density at constant temperature (call DPDD)
c
          t=xin(2)
          rho=xin(3)
          do 410 i=1,nc
  410     x(i)=xin(3+i)
          call DPDD (t,rho,x,dpdrho)
          nout=1
          xout(1)=dpdrho
c  this mode does not return any errors
          ierr=0
          cstr=hnull
c
        else if (imode.eq.32) then
c
c  2nd derivative of pressure w.r.t. density at constant T (call DPDD2)
c
          t=xin(2)
          rho=xin(3)
          do 420 i=1,nc
  420     x(i)=xin(3+i)
          call DPDD2 (t,rho,x,dp2dD2)
          nout=1
          xout(1)=dp2dD2
c  this mode does not return any errors
          ierr=0
          cstr=hnull
c
        else if (imode.eq.33) then
c
c  derivative of pressure w.r.t. temperature at constant density (call DPDT)
c
          t=xin(2)
          rho=xin(3)
          do 430 i=1,nc
  430     x(i)=xin(3+i)
          call DPDT (t,rho,x,dpt)
          nout=1
          xout(1)=dpt
c  this mode does not return any errors
          ierr=0
          cstr=hnull
c
        else if (imode.eq.34) then
c
c  derivative of density w.r.t. pressure at constant temperature (call DDDP)
c
          t=xin(2)
          rho=xin(3)
          do 440 i=1,nc
  440     x(i)=xin(3+i)
          call DDDP (t,rho,x,drhodp)
          nout=1
          xout(1)=drhodp
c  this mode does not return any errors
          ierr=0
          cstr=hnull
c
        else if (imode.eq.35) then
c
c  derivative of density w.r.t. temperature at constant pressure (call DDDT)
c
          t=xin(2)
          rho=xin(3)
          do 450 i=1,nc
  450     x(i)=xin(3+i)
          call DDDT (t,rho,x,drhodt)
          nout=1
          xout(1)=drhodt
c  this mode does not return any errors
          ierr=0
          cstr=hnull
c
        else
          ierr=99
          write (cstr,1099) imode,hnull
 1099     format ('[DLL error 99] mode passed to DLL is out of range:',
     &            i6,a1)
        end if
c
      else if (mode.eq.-1) then
c
c  call the main setup routine to specify mixture components, the file
c  containing the mixing parameters, and one of the standard reference
c  states  (call SETUP)
c
        ncc=NINT(xin(1))   !number of components
        iref=NINT(xin(2))  !code for reference state (see data stmt)
c       ieos=NINT(xin(3))  !code for eqn of state (obsolete)
c
c  parse the input cstr, first pick off path for the fluids directory
c  (inputs are separated by '|')
        do 820 i=1,254
        if (cstr(i:i).eq.'|') then
          lendir=i-1           !length of the directory string
          hdir=cstr(1:lendir)
          nexti=i+1
          goto 824
        endif
  820   continue
c  continue parsing cstr and pick off file name containing mixture info
  824   i0=nexti
        do 840 i=i0,254
        if (cstr(i:i).eq.'|') then
          hfmix=hdir(1:lendir)//cstr(nexti:i-1)
          nexti=i+1
          goto 844
        endif
  840   continue
c  finally, pick off the names for the '.fld' files for each component
  844   i0=nexti
        jcomp=0                !counter for mixture component
        do 860 i=i0,254
        if (cstr(i:i).eq.'|') then
          jcomp=jcomp+1
          if (jcomp.gt.ncc) goto 880
          hfiles(jcomp)=hdir(1:lendir)//cstr(nexti:i-1)
          nexti=i+1
        endif
  860   continue
  880   continue
c  define the reference state and call the main setup routine
        hrf=hrefs(iref)
        call SETUP (ncc,hfiles,hfmix,hrf,ierr,cstr)
        nout=1
        xout(1)=REAL(ncc)
c
      else if (mode.eq.-2) then
c
c  set the reference state for enthalpy and entropy (call SETREF)
c
        iref=NINT(xin(1))
        hrf=hrefs(iref)   !code for reference state (see data stmt)
        ixflag=NINT(xin(2)) !flag:  apply ref state to pures or mix
c  following inputs are meaninful only if hrf='OTH'
        t0=xin(3)
        p0=xin(4)
        h0=xin(5)
        s0=xin(6)
        ncc=NINT(xin(7))   !number of components
        do 920 i=1,ncc
  920   x(i)=xin(7+i)
        call SETREF (hrf,ixflag,x,h0,s0,t0,p0,ierr,cstr)
        nout=0
c  hrf may have been reset by SETREF, pass as output
c       nout=1
c       do 922 i=1,nrefst
c       if (hrf.eq.hrefs(i)) then
c         xout(1)=REAL(i)
c       end if
c 922   continue
c
      else if (mode.eq.-3) then
c
c  set the model (call SETMOD)
c
        ncc=NINT(xin(1))
        itype=NINT(xin(2))
        htype=htypes(itype)    !code for model (see data stmt)
        hmix=cstr(1:3)         !pick off overall/mixture model
        do 940 i=1,ncc
        hcomp(i)=cstr(i*3+1:i*3+3)  !pick off component models
  940   continue
        call SETMOD (ncc,htype,hmix,hcomp,ierr,cstr)
        nout=0
c
      else if (mode.eq.-4) then
c
c  set the mixing parameters (call SETKTV)
c
        hmodij=cstr(1:3)     !3-letter code for mixing rule
        hfmix=cstr(4:83)     !file containing mixture info
        icomp=NINT(xin(1))
        jcomp=NINT(xin(2))
        do 944 i=1,nmxpar
        fij(i)=xin(i+2)
  944   continue
        mode=4               !unknown meaning
        ierr=4
        call SETKTV (icomp,jcomp,hmodij,fij,hfmix,ierr,cstr)
        nout=0
c
      else if (mode.eq.-5) then
c
c  retrieve the mixing parameters (call GETKTV)
c
        icomp=NINT(xin(1))
        jcomp=NINT(xin(2))
        call GETKTV (icomp,jcomp,hmodij,fij,hfmix,hfij,hbinp,hmxrul)
        cstr=hmodij//hbinp
c  the following should reflect any changes in "parameter (nmxpar=6)"
        nout=nmxpar
        do 980 i=1,nout
  980   xout(i)=fij(i)
        ierr=0
c
      else
        ierr=99
        write (cstr,1098) mode,nin,hnull
 1098   format ('[DLL error 99] mode and/or nin passed to DLL is/are ',
     &          'out of range:  mode =',i5,'; nin ='i5,a1)
c
      end if
      mode=ierr
c
      RETURN
      end                                            !subroutine REFPRP6
c
c ======================================================================
c
      subroutine FLOUT (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xout)
c
c  take bulk-phase results of a flash calculation and, depending on
c  value of quality, calculate liquid- and/or vapor-phase properties and
c  write all quantities to the xout array for return to the GUI
c
c  inputs:
c        t--temperature (K)
c        p--pressure (kPa)
c        D--overall (bulk) molar density (mol/L)
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 array (mol frac) for bulk composition
c       xl--composition array (mol frac) for liquid phase
c       xv--composition array (mol frac) for vapor phase
c           if only one phase is present, xl = xv = x
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 = 999 indicates supercritical state (t > Tc)
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, -9.992d is returned
c
c  outputs:
c     xout--calculated (output) values [array of 100 elements]
c           there are 15 possible properties + the compositions,
c           they are returned in the xout array in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c           the first 15 + nc elements are the bulk properties
c           the next 15 + nc elements are liquid-phase properties
c           the next 15 + nc elements are vapor-phase properties
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  08-01-96  MM, original version
c  05-16-97  MM, revise flags for special cases, get from /FLAGS/
c  06-16-97  MM, add infinite flag for V,s for rho < d-20
c  07-14-97  MM, add explicit check for saturated liquid, sat vapor
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      dimension x(ncmax),xl(ncmax),xv(ncmax)
      dimension xout(100)
      common /Gcnst/ R
      common /NCOMP/ nc
c  pointers to the xout array (to aid in adding possible future props)
      common /IPOINT/
     &     it,ip,irho,ivol,ie,ih,is,icv,icp,iw,iq,ieta,itcx,ist,iz,ix
c  maximum number of properties returned to the GUI on a single call
      common /NPROP/ nprops
c  common block containing flags to GUI
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  set bulk-phase outputs
c
      xout(it)=t
      xout(ip)=p
      xout(irho)=D
      if (D.ge.1.0d-20) then
        xout(ivol)=1.0d0/D
        xout(iz)=p/(D*R*t)
        xout(is)=s
      else
c  density is zero, volume & entropy are infinite
        xout(ivol)=xinf
        xout(iz)=1.0d0
        xout(is)=xinf
      end if
      xout(ie)=e
      xout(ih)=h
      xout(icv)=cv
      xout(icp)=cp
      xout(iw)=w
      xout(iq)=q
c  the transport properties and surface tension are not calc here
      xout(ieta)=xnotc
      xout(itcx)=xnotc
      xout(ist)=xnotc
      do 100 i=1,nc
      xout(ix+i-1)=x(i)
  100 continue
c
      nval=nprops+nc
      if (abs(999.0d0-q).lt.1.0d-6) then
c
c  super-critical state (indicated by q = 999) set liquid and vapor
c  properties equal to bulk-phase properties
c
        do 199 j=1,nval
        xout(nval+j)=xout(j)
        xout(2*nval+j)=xout(j)
  199   continue
c
      else if (ABS(q).le.1.0d-8 .or. ABS(1.0d0-q).le.1.0d-8) then
c
c  saturation state, calculate the liquid and vapor-phase properties
c
        tv=t     !temperature and pressure for both phases always
        pv=p     !the same for a flash calculation
        call VLOUT (t,tv,p,pv,Dl,Dv,xl,xv,xout)
c
      else if (q.lt.0.0d0) then
c
c  sub-cooled liquid state, set liquid properties equal to bulk-phase properties; 
c  vapor properties are set to flag indicating subcooled
c
        do 240 j=1,nval
        xout(nval+j)=xout(j)
        xout(2*nval+j)=xsubc
  240   continue
c
      else if (q.gt.1.0d0) then
c
c  superheated vapor state, set vapor properties equal to bulk-phase
c  properties and liquid properties to flag indicating superheated
c
        do 280 j=1,nval
        xout(nval+j)=xsuph
        xout(2*nval+j)=xout(j)
  280   continue
c
      else
c
c  two-phase state, calculate the liquid- and vapor-phase properties
c
        tv=t     !temperature and pressure for both phases always
        pv=p     !the same for a flash calculation
        call VLOUT (t,tv,p,pv,Dl,Dv,xl,xv,xout)
c  bulk-phase heat capacity and speed of sound not defined for two-phase
        xout(icv)=x2ph
        xout(icp)=x2ph
        xout(iw)=x2ph
      end if
c
      xout(nval+iq)=xnotd    !quality is not defined for liq, vap
      xout(2*nval+iq)=xnotd
c
      RETURN
      end                                              !subroutine FLOUT
c
c ======================================================================
c
      subroutine SATOUT (tl,tv,pl,pv,Dl,Dv,xl,xv,xout)
c
c  take results of a saturation calculation and calculate liquid and
c  vapor-phase properties and write all quantities to the xout array
c  for return to the GUI
c
c  inputs:
c       tl--temperature (K) of the liquid phase
c       tv--temperature (K) of the vapor phase
c       pl--pressure (kPa) of the liquid phase
c       pv--pressure (kPa) of the vapor phase
c       Dl--molar density (mol/L) of the liquid phase
c       Dv--molar density (mol/L) of the vapor phase
c       xl--composition array (mol frac) for liquid phase
c       xv--composition array (mol frac) for vapor phase
c
c  outputs:
c     xout--calculated (output) values [array of 100 elements]
c           there are 15 possible properties + the compositions,
c           they are returned in the xout array in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c           the first 15 + nc elements are the bulk properties (which
c             are not applicable in this case; 8.888d8 is returned)
c           the next 15 + nc elements are liquid-phase properties
c           the next 15 + nc elements are vapor-phase properties
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  08-02-96  MM, original version
c  05-16-97  MM, revise flags for special cases, get from /FLAGS/
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      dimension xl(ncmax),xv(ncmax)
      dimension xout(100)
      common /NCOMP/ nc
c  pointers to the xout array (to aid in adding possible future props)
      common /IPOINT/
     &     it,ip,irho,ivol,ie,ih,is,icv,icp,iw,iq,ieta,itcx,ist,iz,ix
c  maximum number of properties returned to the GUI on a single call
      common /NPROP/ nprops
c  flags to GUI
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  set bulk-phase outputs
c
c     nval=nprops+nc
c  the bulk-phase properties are not defined
      do 100 i=1,nprops
      xout(i)=xnotd
  100 continue
c  set the bulk compositions equal to 'not defined' (this is an
c  arbitrary choice, but should return something to prevent GUI crash)
      do 120 i=1,nc
      xout(ix+i-1)=xl(i)
  120 continue
c
c  calculate the liquid- and vapor-phase properties
c
      call VLOUT (tl,tv,pl,pv,Dl,Dv,xl,xv,xout)
c
      RETURN
      end                                             !subroutine SATOUT
c
c ======================================================================
c
      subroutine VLOUT (tl,tv,pl,pv,Dl,Dv,xl,xv,xout)
c
c  calculate liquid and vapor-phase properties and write all quantities
c  to the xout array for return to the GUI
c
c  inputs:
c       tl--temperature (K) of the liquid phase
c       tv--temperature (K) of the vapor phase
c       pl--pressure (kPa) of the liquid phase
c       pv--pressure (kPa) of the vapor phase
c       Dl--molar density (mol/L) of the liquid phase
c       Dv--molar density (mol/L) of the vapor phase
c       xl--composition array (mol frac) for liquid phase
c       xv--composition array (mol frac) for vapor phase
c
c  outputs:
c     xout--calculated (output) values [array of 100 elements]
c           there are 15 possible properties + the compositions,
c           they are returned in the xout array in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c           the first 15 + nc elements are the bulk properties; they
c             are not set in this routine
c           the next 15 + nc elements are liquid-phase properties
c           the next 15 + nc elements are vapor-phase properties
c           note:  the quality is not set in this routine
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  08-02-96  MM, original version
c  05-16-97  MM, revise flags for special cases, get from /FLAGS/
c  06-16-97  MM, add infinite flag for V,s for rho < d-20
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      dimension xl(ncmax),xv(ncmax)
      dimension xout(100)
      common /Gcnst/ R
      common /NCOMP/ nc
c  pointers to the xout array (to aid in adding possible future props)
      common /IPOINT/
     &     it,ip,irho,ivol,ie,ih,is,icv,icp,iw,iq,ieta,itcx,ist,iz,ix
c  maximum number of properties returned to the GUI on a single call
      common /NPROP/ nprops
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  calculate the liquid-phase properties
c
      nval=nprops+nc
      call THERM (tl,Dl,xl,ptherm,el,hl,sl,cvl,cpl,wl,hjt)
      xout(nval+it)=tl
      xout(nval+ip)=pl
      xout(nval+irho)=Dl
      if (Dl.ge.1.0d-20) then
        xout(nval+ivol)=1.0d0/Dl
        xout(nval+iz)=pl/(Dl*R*tl)
        xout(nval+is)=sl
      else
c  density is zero, volume & entropy are infinite
        xout(nval+ivol)=xinf
        xout(nval+iz)=1.0d0
        xout(nval+is)=xinf
      end if
      xout(nval+ie)=el
      xout(nval+ih)=hl
      xout(nval+icv)=cvl
      xout(nval+icp)=cpl
      xout(nval+iw)=wl
      xout(nval+iq)=xnotd    !quality is not defined
      do 100 i=1,nc
      xout(nval+ix+i-1)=xl(i)
  100 continue
c
c  calculate the vapor-phase properties
c
      call THERM (tv,Dv,xv,ptherm,ev,hv,sv,cvv,cpv,wv,hjt)
      xout(2*nval+it)=tv
      xout(2*nval+ip)=pv
      xout(2*nval+irho)=Dv
      if (Dv.ge.1.0d-20) then
        xout(2*nval+ivol)=1.0d0/Dv
        xout(2*nval+iz)=pv/(Dv*R*tv)
        xout(2*nval+is)=sv
      else
c  density is zero, volume & entropy are infinite
        xout(2*nval+ivol)=xinf
        xout(2*nval+iz)=1.0d0
        xout(2*nval+is)=xinf
      end if
      xout(2*nval+ie)=ev
      xout(2*nval+ih)=hv
      xout(2*nval+icv)=cvv
      xout(2*nval+icp)=cpv
      xout(2*nval+iw)=wv
      xout(2*nval+iq)=xnotd  !quality is not defined
      do 200 i=1,nc
      xout(2*nval+ix+i-1)=xv(i)
  200 continue
c  the transport properties and surface tension are not calc here
      do 240 i=1,2
      iprop=i*nval
      xout(iprop+ieta)=xnotc
      xout(iprop+itcx)=xnotc
      xout(iprop+ist)=xnotc
  240 continue
c
      RETURN
      end                                              !subroutine VLOUT
c
c ======================================================================
c
      subroutine SPOUT (t,D,x,kph,xout)
c
c  calculate properties for the specified phase and write all quantities
c  to the xout array for return to the GUI; values for the "other" phase
c  are set to "not calculated;" the bulk properties are set to those of
c  the specified phase
c
c  inputs:
c        t--temperature [K]
c        D--molar density [mol/L]
c        x--composition array [mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c
c  outputs:
c     xout--calculated (output) values [array of 100 elements]
c           there are 15 possible properties + the compositions,
c           they are returned in the xout array in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c           the first 15 + nc elements are the bulk properties
c           the next 15 + nc elements are liquid-phase properties
c           the next 15 + nc elements are vapor-phase properties
c           note:  the quality is not set in this routine
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  06-06-97  MM, original version; based on VLOUT
c  06-16-97  MM, add infinite flag for V,s for rho < d-20
c  10-31-97  MM, composition array for bulk phase not set properly (all x = x(1))
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      dimension x(ncmax)
      dimension xout(100)
      common /Gcnst/ R
      common /NCOMP/ nc
c  pointers to the xout array (to aid in adding possible future props)
      common /IPOINT/
     &     it,ip,irho,ivol,ie,ih,is,icv,icp,iw,iq,ieta,itcx,ist,iz,ix
c  maximum number of properties returned to the GUI on a single call
      common /NPROP/ nprops
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  calculate the properties of the specified phase
c
      if (kph.eq.1) then
c  liquid phase
        nval=nprops+nc
      else
c  vapor phase
        nval=2*(nprops+nc)
      end if
      call THERM (t,D,x,ptherm,e,h,s,cv,cp,w,hjt)
      xout(nval+it)=t
      xout(nval+ip)=ptherm
      xout(nval+irho)=D
      if (D.ge.1.0d-20) then
        xout(nval+ivol)=1.0d0/D
        xout(nval+iz)=ptherm/(D*R*t)
        xout(nval+is)=s
      else
c  density is zero, volume & entropy are infinite
        xout(nval+ivol)=xinf
        xout(nval+iz)=1.0d0
        xout(nval+is)=xinf
      end if
      xout(nval+ie)=e
      xout(nval+ih)=h
      xout(nval+icv)=cv
      xout(nval+icp)=cp
      xout(nval+iw)=w
      xout(nval+iq)=xnotc    !quality is not calculated
      do 100 i=1,nc
      xout(nval+ix+i-1)=x(i)
  100 continue
c  the transport properties and surface tension are not calc here
      do 240 i=1,2
      iprop=i*nval
      xout(iprop+ieta)=xnotc
      xout(iprop+itcx)=xnotc
      xout(iprop+ist)=xnotc
  240 continue
c  set the bulk properties to those of the specified phase
      do 300 i=1,nprops
      xout(i)=xout(nval+i)
  300 continue
      do 320 i=1,nc
      xout(ix+i-1)=xout(nval+ix+i-1)
  320 continue
c
      RETURN
      end                                              !subroutine SPOUT
c
c ======================================================================
c
      subroutine REFP61 (mode,xin,intin,cstrin,xblk,xliq,xvap,xextra,
     &                   ierri,nwarn,iwarn,cstr)
c
c  This subroutine is the main link between the Fortran property
c  routines and the Pascal-based graphical user interface.  It calls the
c  appropriate property routine(s) and returns values to the interface.
c  
c  This version is very similar to REFPRP6 (see above), but with inputs and
c  outputs restructured for compatibility with Refprop version 6.1 and later.
c
c  N.B. This routine is still under development, and all arguments are
c       subject to change.
c
c  inputs:
c     mode--specifies the type of calculation to be done
c           mode = -5:  get mixture interaction coef      (call GETKTV)
c                  -4:  set mixture interaction coef      (call SETKTV)
c                  -3:  specify non-standard model        (call SETMOD)
c                  -2:  specify non-standard ref state    (call SETREF)
c                  -1:  specify fluid(s), ref state, etc. (call SETUP)
c                   1:  thermo props as f(T,P,x)          (call TPFLSH)
c                   2:  thermo props as f(T,rho,x)        (call TDFLSH)
c                   3:  thermo props as f(P,h,x)          (call PHFLSH)
c                   4:  thermo props as f(P,s,x)          (call PSFLSH)
c                   5:  thermo props as f(T,rho,x)        (call THERM)
c                       N.B. mode = 5 does no error checking and should be used
c                            only when conditions are known to be single-phase
c                   6:  satn thermo props as f(T,xl=xv)   (call SATT)
c                   7:  satn thermo props as f(P,xl=xv)   (call SATP)
c                   8:  critical parameters as f(x)       (call CRITP)
c                   9:  transport propertes as f(t,rho,x) (call TRNPRP)
c                  10:  fugacities as f(t,rho,x)          (call FGCTY)
c                  11:  surface tension as f(t,rho,x)     (call SURFT)
c                  12:  bubble point props as f(t,xl)     (call SATT)
c                  13:  dew point props as f(t,xv)        (call SATT)
c                  14:  bubble point props as f(p,xl)     (call SATP)
c                  15:  dew point props as f(p,xv)        (call SATP)
c                  16:  thermo props as f(T,q,x)          (call TQFLSH)
c                  17:  thermo props as f(p,q,x)          (call PQFLSH)
c                  18:  molecular weight as f(x)          (wm = XMOL)
c                  19:  convert mole to mass fraction     (call XMASS)
c                  20:  convert mass to mole fraction     (call XMOLE)
c                  21:  CAS # and name for component i    (call NAME)
c                  22:  fixed properties for component i  (call INFO)
c                  23:  limits of specified model as f(x) (call LIMITX)
c                  24:  limits for a single component     (call LIMITK)
c                  25:  single-phase T-P flash            (call TPRHO)
c                  26:  single-phase P-h flash            (call PHFL1)
c                  27:  single-phase P-s flash            (call PSFL1)
c      xin--input variables [array of dimension 25]
c         xin(1..nin):  independent variables in order as noted above
c                       e.g. T,P,x(1)..x(nc)   for mode = 1
c                       e.g. T,rho,x(1)..x(nc) for mode = 2,9,10,11,31-35
c                       e.g. T,x(1)..x(nc)     for mode = 6
c                       e.g. x(1)..x(nc)       for mode = 8,18,19,20
c
c      but, in the case of the setup calls (mode < 0) xin is different:
c        if mode = -5:  xin(1) = component i
c                       xin(2) = component j (for the pair i,j)
c        if mode = -4:  xin(1) = component i
c                       xin(2) = component j (for the pair i,j)
c                       xin(3..6) = mixing coefficients 1..4
c        if mode = -3:  xin(1) = number of component(s) in mixture
c                       xin(2) = code for type of model to be set
c                              1:  equation of state ('EOS')
c                              2:  viscosity ('ETA')
c                              3:  thermal conductivity ('TCX')
c                              4:  surface tension ('STN')
c        if mode = -2:  xin(1) = code for reference state
c                              1:  normal boiling point ('NBP')
c                              2:  ASHRAE convention ('ASH')
c                              3:  IIR convention ('IIR')
c                              4:  other ('OTH')
c                       xin(2) = composition flag
c                              1:  ref state applied to pure components
c                              2:  ref state applied to mixture of
c                                  composition specified by xin(8..7+nc)
c                       xin(3) = Tref
c                       xin(4) = Pref
c                       xin(5) = href
c                       xin(6) = sref
c                       xin(7) = number of component(s)
c                       xin(8..7+nc) = composition(s)
c        if mode = -1:  xin(1) = number of component(s) in mixture
c                       xin(2) = code for reference state
c                              1:  normal boiling point ('NBP')
c                              2:  ASHRAE convention ('ASH')
c                              3:  IIR convention ('IIR')
c                              4:  other ('OTH') [must use mode = -2]
c                       xin(3) = code for equation of state [obsolete]
c      in the case of the NAME and INFO calls (mode = 21, 22) xin is different:
c        if mode = 21:  xin(1) = component i
c        if mode = 22:  xin(1) = component i
c      in the case of the call to LIMITX xin is:
c        if mode = 23:  xin(1) = code for type of model to be set
c                              1:  equation of state ('EOS')
c                              2:  viscosity ('ETA')
c                              3:  thermal conductivity ('TCX')
c                              4:  surface tension ('STN')
c                       xin(2) = temperature of mixture
c                       xin(3) = density of mixture
c                       xin(4) = pressure of mixture
c                       xin(5..4+nc) = composition(s)
c      in the case of the call to LIMITK xin is:
c        if mode = 24:  xin(1) = code for type of model to be set
c                              1:  equation of state ('EOS')
c                              2:  viscosity ('ETA')
c                              3:  thermal conductivity ('TCX')
c                              4:  surface tension ('STN')
c                       xin(2) = component number
c                       xin(3) = temperature of mixture
c                       xin(4) = density of mixture
c                       xin(5) = pressure of mixture
c      in the case of the call to TPRHO xin is:
c        if mode = 25:  xin(1) = temperature
c                       xin(2) = pressure
c                       xin(3..2+nc) = composition(s)
c                       xin(3+nc) = phase flag:  1 = find liquid solution
c                                                2 = find vapor solution
c                       xin(5+nc) = guess flag:  1 = initial guess for density provided
c                       xin(6+nc) = initial guess for density (ignored if xin(5+nc)<>1)
c    intin--integer input array
c   cstrin--used to pass text info (e.g. file names) to the Fortran
c           property routines, (applicable only for certain modes)
c           [character*255 with last character the ASCII null (char(0))]
c        if mode = -4:  cstrin(1:3) = mixing rule code (e.g. 'LJ1')
c                       cstrin(4:83) = file containing mixture data
c           mode = -3:  cstrin(1:3) = code for overall/mix EOS (e.g. 'HMX')
c                       cstrin(4:6) = code for component 1 EOS (e.g. 'BWR')
c                       cstrin(7:9,etc) = ditto for each of the nc components
c           mode = -1:  cstrin consists of various data separated by '|'
c                       --path to the fluids directory (e.g. C:\REFPROP\fluids\)
c                       --file containing mixture data (e.g. HMX.bnc)
c                       --file(s) containing fluid data for each of the nc
c                         components (e.g. R134a.fld)
c                       for example, to specify a mixture of R32 and
c                       R134a, with mixture data in file HMX.bnc, this
c                       string might be:
c                cstrin = 'C:\REFPROP\fluids\|HMX.bnc|R32.fld|R134a.fld|'//char(0)
c
c  outputs:
c     xblk--calculated (output) values for the bulk (overall)
c     xliq--calculated (output) values for the liquid phase
c     xvap--calculated (output) values for the vapor phase
c   xextra--calculated (output) values for an as-yet-undefined third phase
c           there are 15 possible properties + the compositions, for most
c           cases, they are returned in the output array(s) in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c
c           but there are exceptions to the above order:
c               mode = -1,  SETUP:  xblk(1) = ncomp
c               mode = -2, SETREF:  nothing in output arrays
c               mode = -3, SETMOD:  nothing in output arrays
c               mode = -4, SETKTV:  nothing in output arrays
c               mode = -5, GETKTV:  xblk(1..4) = mixing coefficients 1..4
c               mode = 21,   NAME:  nothing in output arrays
c               mode =  8,  CRITP:  xblk(1) = Tcrit
c                                   xblk(2) = Pcrit
c                                   xblk(3) = Dcrit
c               mode =  9, TRNPRO:  viscosity and thermal conductivity returned in
c                                   array specified by cstrin
c               mode = 10,  FGCTY:  the nc fugacities returned in array spec by cstrin
c               mode = 18,   WMOL:  xblk(1) = the molecular weight [g/mol]
c               mode = 19,  XMASS:  xblk(1..nc) = the composition on basis opposite of input
c                      20,  XMOLE:  xblk(nc+1)  = the molecular weight [g/mol]
c               mode = 21,   NAME:  nothing in output arrays (all info in string cstr)
c               mode = 22,   INFO:  xblk(1) = molecular weight [g/mol]
c                                   xblk(2) = triple point temperature [K]
c                                   xblk(3) = normal boiling point temperature [K]
c                                   xblk(4) = critical temperature [K]
c                                   xblk(5) = critical pressure [kPa]
c                                   xblk(6) = critical density [mol/L]
c                                   xblk(7) = compressibility at critical point [pc/(R*Tc*Dc)]
c                                   xblk(8) = accentric factor [-]
c                                   xblk(9) = dipole moment [debye]
c                                  xblk(10) = gas constant [J/mol-K]
c               mode = 23, LIMITX:  xblk(1) = lower temperature limit [K]
c                                   xblk(2) = upper temperature limit [K]
c                                   xblk(3) = upper density limit [mol/L]
c                                   xblk(4) = upper pressure limit [kPa]
c               mode = 24, LIMITK:  same as for mode = 23, LIMITX, except
c                                   limits valid only for specified component
c               mode = 31, DPDD:    xblk(1) = the derivative dP/d(rho)
c                      32, DPDD2:   xblk(1) = the derivative d2P/d(rho)2
c                      33, DPDT:    xblk(1) = the derivative dP/dT
c                      34, DDDP:    xblk(1) = the derivative d(rho)/dP
c                      35, DDDT:    xblk(1) = the derivative d(rho)/dT
c    ierri--error code returned by property subroutine
c    nwarn--number of warning(s) returned by property subroutines
c    iwarn--warning code(s) returned by property subroutine [integer array]
c     cstr--used to pass error messages to GUI,
c           except for:
c             mode = -5, cstr(1:3) = mixing rule code (e.g. 'VDW')
c                        cstr(4..) = documentation for the binary parameters
c             mode = 21, cstr(1:12)   = short name
c                        cstr(13:92)  = full name
c                        cstr(93:104) = CAS number
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  11-06-97  MM, original version; based on REFPRP6
c  11-12-97  MM, add integer input array (for GUI-compatibility with STEAM)
c  11-25-97  MM, add xextra array for future use (e.g. LLVE)
c  12-02-97  MM, initialize ierr=0, herr=hnull for mode=5 (call THERM)
c                move definition of hrefst, htypes to BDDLL2
c
c  compiler switches for use with Digital Visual Fortran; these should be 
c  treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: REFP61
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (nmxpar=6)       !number of binary mixture parameters
      parameter (nrefst=6)       !number of reference state choices
      parameter (narray=80)      !size of output (blk, liq, vap) arrays
      parameter (nwarna=10)      !size of warning array
      integer*4 mode,ierri,nwarn,iwarn(nwarna),intin(narray)
      character*1 htab,hnull
      character*3 hrf,hcomp(ncmax),hmodij
      character*3 hrefs,htypes,htype,hmix
      character*8  hfij(nmxpar)
      character*12 hcas,hname
      character*80 hfmix,hfiles(ncmax),hdir,hmxrul,hn80
      character*255 cstrin,cstr,hbinp
      dimension xin(25),xblk(narray),xliq(narray),xvap(narray),
     &          xextra(narray)
      dimension x(ncmax),xl(ncmax),xv(ncmax),xjunk(ncmax),
     &          fug(ncmax),fij(nmxpar),xkg(ncmax),xmol(ncmax)
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
c  pointers to the output arrays for use with REFP61
      common /POINT2/ it,ip,irho,ivol,ie,ih,is,icv,icp,iw,
     &                iq,iz,ihjt,iA,iG,ikappa,ibeta,idPdD,id2P,idPdT,
     &                idDdT,idDdP,ist,ieta,itcx,idiel,ij1,ij2,ij3,ij4,
     &                ix,ifug
c  maximum number of properties returned to the GUI on a single call
      common /NPROP2/ nprops
c  common block containing flags to GUI (initialized in BDSET in setup.f)
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c  reference state choices and model types
      common /HREFS2/ hrefs(nrefst)
      common /HMODT2/ htypes(4)
c
d     write (*,*) ' REFP61--input mode: ',mode
d     write (*,1080) (xin(i),i=1,10)
d1080 format (1x,'REFP61--input array:  ',10f14.6)
d     write (*,1255) cstr
d1255 format (1x,'REFP61--input C_str:  (',a255,')')
c
c  load the output arrays with flag signifying 'not calculated'; the
c  GUI will crash if any elements in the output array are not defined;
c  this is insurance in case any of the property routines crash
      nout=nprops+2*ncmax     !# properties + compositions + fugacities
      do 80 i=1,nout
      xblk(i)=xnotc
      xliq(i)=xnotc
      xvap(i)=xnotc
  80  continue
c  initialize the error flags
      ierri=0
      nwarn=0
      iwarn(1)=0
      cstr=hnull
c
      if (mode.eq.1) then
c
c  temperature-pressure flash (call TPFLSH)
c
        t=xin(1)
        p=xin(2)
        do 110 i=1,nc
  110   x(i)=xin(2+i)
        call TPFLSH (t,p,x,D,Dl,Dv,xl,xv,q,e,h,s,cv,cp,w,ierr,cstr)
        call FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xblk,xliq,xvap)
c
      else if (mode.eq.2) then
c
c  temperature-density flash (call TDFLSH)
c
        t=xin(1)
        D=xin(2)
        do 120 i=1,nc
  120   x(i)=xin(2+i)
        call TDFLSH (t,D,x,p,Dl,Dv,xl,xv,q,e,h,s,cv,cp,w,ierr,cstr)
        call FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xblk,xliq,xvap)
c
      else if (mode.eq.3) then
c
c  pressure-enthalpy flash (call PHFLSH)
c
        p=xin(1)
        h=xin(2)
        do 130 i=1,nc
  130   x(i)=xin(2+i)
        call PHFLSH (p,h,x,t,D,Dl,Dv,xl,xv,q,e,s,cv,cp,w,ierr,cstr)
        call FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xblk,xliq,xvap)
c
      else if (mode.eq.4) then
c
c  pressure-entropy flash (call PSFLSH)
c
        p=xin(1)
        s=xin(2)
        do 140 i=1,nc
  140   x(i)=xin(2+i)
        call PSFLSH (p,s,x,t,D,Dl,Dv,xl,xv,q,e,h,cv,cp,w,ierr,cstr)
        call FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xblk,xliq,xvap)
c
      else if (mode.eq.5) then
c
c  temperature-density calculation (call THERM)
c  N.B.  this option does not check for two-phase states
c        it will return nonsense values for meta-stable states
c
        t=xin(1)
        D=xin(2)
        do 150 i=1,nc
        x(i)=xin(2+i)
        xl(i)=x(i)
        xv(i)=x(i)
  150   continue
        call THERM (t,D,x,p,e,h,s,cv,cp,w,hjt)
c  quality of 999 flags FLOUT2 to set liq & vap properties to bulk
        q=999.0d0
        Dl=D
        Dv=D
        call FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xblk,xliq,xvap)
c  this mode does not return any errors
        ierr=0
        cstr=hnull
c
      else if (mode.eq.6) then
c
c  saturation state at specified t, xl=xv (call SATT)
c
        t=xin(1)
        do 160 i=1,nc
  160   x(i)=xin(1+i)
c  call SATT at both dew point and bubble point
        call SATT (t,x,1,pl,Dl,Djunk,xl,xjunk,ierr,cstr)
        if (ierr.gt.0) then
          RETURN
        end if
        call SATT (t,x,2,pv,Djunk,Dv,xjunk,xv,ierr,cstr)
        if (ierr.gt.0) then
          RETURN
        end if
        tv=t
        call STOUT2 (t,tv,pl,pv,Dl,Dv,xl,xv,xblk,xliq,xvap)
c  only the temperature of the bulk phase is defined
        xblk(it)=t
c
      else if (mode.eq.7) then
c
c  saturation state at specified p, xl=xv (call SATP)
c
        p=xin(1)
        do 170 i=1,nc
  170   x(i)=xin(1+i)
c  call SATP at both dew point and bubble point
        call SATP (p,x,1,tl,Dl,Djunk,xl,xjunk,ierr,cstr)
        if (ierr.gt.0) then
            RETURN
        end if
        call SATP (p,x,2,tv,Djunk,Dv,xjunk,xv,ierr,cstr)
        if (ierr.gt.0) then
          RETURN
        end if
        pv=p
        call STOUT2 (tl,tv,p,pv,Dl,Dv,xl,xv,xblk,xliq,xvap)
c  only the pressure of the bulk phase is defined
        xblk(ip)=p
c
      else if (mode.eq.8) then
c
c  calculate critical parameters (call CRITP)
c
        do 180 i=1,nc
  180   x(i)=xin(i)
        call CRITP (x,tc,pc,Dc,ierr,cstr)
        xblk(1)=tc
        xblk(2)=pc
        xblk(3)=Dc
c
      else if (mode.eq.9) then
c
c  transport properties (viscosity, thermal conductivity) (call TRNPRP)
c
        t=xin(1)
        rho=xin(2)
        do 190 i=1,nc
  190   x(i)=xin(2+i)
        call TRNPRP (t,rho,x,eta,tcx,ierr,cstr)
c  TRNPRP is unaware of phase--use input string to put results in
c  appropriate output array
        if (cstrin(1:3).eq.'liq') then
          xliq(ieta)=eta
          xliq(itcx)=tcx
        else if (cstrin(1:3).eq.'vap') then
          xvap(ieta)=eta
          xvap(itcx)=tcx
        else
          xblk(ieta)=eta
          xblk(itcx)=tcx
        end if
c
        else if (mode.eq.10) then
c
c  calculate fugacity for each of the nc components (call FGCTY)
c
        t=xin(1)
        rho=xin(2)
        do 200 i=1,nc
  200   x(i)=xin(2+i)
        call FGCTY (t,rho,x,fug)
        nout=nc
        ierr=0
        cstr=hnull
        do 201 i=1,nc
        xblk(ifug+i-1)=fug(i)
  201   continue
c
      else if (mode.eq.11) then
c
c  surface tension (call SURFT)
c
        t=xin(1)
        rho=xin(2)
        do 210 i=1,nc
  210   x(i)=xin(2+i)
        call SURFT (t,rho,x,sigma,ierr,cstr)
        xblk(ist)=sigma
c
      else if (mode.ge.12 .and. mode.le.15) then
c
c  saturation state with xl and xv in equilibrium
c
        do 220 i=1,nc
  220   x(i)=xin(1+i)
c
        if (mode.eq.12) then
c
c  saturation state at bubble point
c  specified t,xl; xv is equilibium composition  (call SATT)
c
          t=xin(1)
          call SATT (t,x,1,p,Dl,Dv,xl,xv,ierr,cstr)
          if (ierr.gt.0) then
            RETURN
          end if
c
        else if (mode.eq.13) then
c
c  saturation state at dew point
c  specified t,xv; xl is equilibium composition  (call SATT)
c
          t=xin(1)
          call SATT (t,x,2,p,Dl,Dv,xl,xv,ierr,cstr)
          if (ierr.gt.0) then
            RETURN
          end if
c
        else if (mode.eq.14) then
c
c  saturation state at bubble point
c  specified p,xl; xv is equilibium composition  (call SATP)
c
          p=xin(1)
          call SATP (p,x,1,t,Dl,Dv,xl,xv,ierr,cstr)
          if (ierr.gt.0) then
            RETURN
          end if
c
        else if (mode.eq.15) then
c
c  saturation state at dew point
c  specified p,xv; xl is equilibium composition  (call SATP)
c
          p=xin(1)
          call SATP (p,x,2,t,Dl,Dv,xl,xv,ierr,cstr)
          if (ierr.gt.0) then
            RETURN
          end if
c
        end if
        tl=t
        tv=t
        pl=p
        pv=p
        call STOUT2 (tl,tv,pl,pv,Dl,Dv,xl,xv,xblk,xliq,xvap)
c  only the temperature and pressure are defined for the bulk phase
        xblk(it)=t
        xblk(ip)=p
c
      else if (mode.eq.16) then
c
c  temperature-quality flash (call TQFLSH)
c
        t=xin(1)
        q=xin(2)
        do 260 i=1,nc
  260   x(i)=xin(2+i)
c  N.B.:  REFPRP6 assumed mass basis (but that option not working then)
c         may need to have GUI supply basis
        kq=1    !input quality is on mole basis
        call TQFLSH (t,q,x,kq,p,D,Dl,Dv,xl,xv,e,h,s,cv,cp,w,ierr,cstr)
        call FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xblk,xliq,xvap)
c
      else if (mode.eq.17) then
c
c  pressure-quality flash (call PQFLSH)
c
        p=xin(1)
        q=xin(2)
        do 270 i=1,nc
  270   x(i)=xin(2+i)
c  N.B.:  REFPRP6 assumed mass basis (but that option not working then)
c         may need to have GUI supply basis
        kq=1    !input quality is on mole basis
        call PQFLSH (p,q,x,kq,t,D,Dl,Dv,xl,xv,e,h,s,cv,cp,w,ierr,cstr)
        call FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,xblk,xliq,xvap)
c
      else if (mode.eq.18) then
c
c  calculate molecular weight (wm = WMOL)
c
        do 280 i=1,nc
  280   x(i)=xin(1+i)
        ierr=0
        cstr=hnull
        xblk(1)=WMOL(x)
        nout=1
c
      else if (mode.eq.19) then
c
c  convert composition in mole fraction to mass fraction (call XMASS)
c
        do 290 i=1,nc
  290   xmol(i)=xin(i)
        ierr=0
        cstr=hnull
        call XMASS (xmol,xkg,wmix)
        do 292 i=1,nc
  292   xblk(i)=xkg(i)
        xblk(nc+1)=wmix
c
      else if (mode.eq.20) then
c
c  convert composition in mass fraction to mole fraction (call XMOLE)
c
        do 300 i=1,nc
  300   xkg(i)=xin(i)
        ierr=0
        cstr=hnull
        call XMOLE (xkg,xmol,wmix)
        do 302 i=1,nc
  302   xblk(i)=xmol(i)
        xblk(nc+1)=wmix
c
      else if (mode.eq.21) then
c
c  return CAS and names for a component (call NAME)
c
        icomp=NINT(xin(1))
        call NAME (icomp,hname,hn80,hcas)
        cstr(1:12)=hname
        cstr(13:92)=hn80
        cstr(93:104)=hcas
        cstr(105:105)=hnull
        ierr=0
c
      else if (mode.eq.22) then
c
c  return various fixed points for a component (call INFO)
c
        icomp=NINT(xin(1))
        call INFO (icomp,wm,ttp,tnbp,tc,pc,Dc,Zc,acf,dip,Rgas)
        xblk(1)=wm
        xblk(2)=ttp
        xblk(3)=tnbp
        xblk(4)=tc
        xblk(5)=pc
        xblk(6)=Dc
        xblk(7)=Zc
        xblk(8)=acf
        xblk(9)=dip
        xblk(10)=Rgas
        ierr=0
        cstr=hnull
c
        else if (mode.eq.23) then
c
c  return limits to specified model (call LIMITX)
c
        itype=NINT(xin(1))
        htype=htypes(itype)    !code for model (see data stmt)
        t=xin(2)
        D=xin(3)
        p=xin(4)
        do 330 i=1,nc
  330   x(i)=xin(4+i)
        call LIMITX (htype,t,D,p,x,tmin,tmax,Dmax,pmax,ierr,cstr)
        xblk(1)=tmin
        xblk(2)=tmax
        xblk(3)=Dmax
        xblk(4)=pmax
c
      else if (mode.eq.24) then
c
c  return limits to specified model for a given compontent (call LIMITK)
c
        itype=NINT(xin(1))
        htype=htypes(itype)    !code for model (see data stmt)
        icomp=NINT(xin(2))
        t=xin(3)
        D=xin(4)
        p=xin(5)
        call LIMITK (htype,icomp,t,D,p,tmin,tmax,Dmax,pmax,ierr,cstr)
        xblk(1)=tmin
        xblk(2)=tmax
        xblk(3)=Dmax
        xblk(4)=pmax
c
      else if (mode.eq.25) then
c
c  single-phase temperature-pressure flash (call TPRHO)
c  N.B.  use TPFLSH (mode = 1) when phase is not known
c
        t=xin(1)
        p=xin(2)
        do 350 i=1,nc
  350   x(i)=xin(2+i)
        kph=NINT(xin(3+nc))
        kguess=NINT(xin(4+nc))
        if (kguess.eq.1) D=xin(5+nc)
        call TPRHO (t,p,x,kph,kguess,D,ierr,cstr)
        call SPOUT2 (t,D,x,kph,xblk,xliq,xvap)
c
      else if (mode.eq.26) then
c
c  single-phase pressure-entropy flash (call PHFL1)
c  N.B.  use PHFLSH (mode = 3) when phase is not known
c
        p=xin(1)
        h=xin(2)
        do 360 i=1,nc
  360   x(i)=xin(2+i)
        kph=NINT(xin(3+nc))
        t=xin(4+nc)             !initial guess for temperature
        D=xin(5+nc)             !initial guess for density
        call PHFL1 (p,h,x,kph,t,D,ierr,cstr)
        call SPOUT2 (t,D,x,kph,xblk,xliq,xvap)
c
      else if (mode.eq.27) then
c
c  single-phase pressure-entropy flash (call PSFL1)
c  N.B.  use PSFLSH (mode = 4) when phase is not known
c
        p=xin(1)
        s=xin(2)
        do 370 i=1,nc
  370   x(i)=xin(2+i)
        kph=NINT(xin(3+nc))
        t=xin(4+nc)             !initial guess for temperature
        D=xin(5+nc)             !initial guess for density
        call PSFL1 (p,s,x,kph,t,D,ierr,cstr)
        call SPOUT2 (t,D,x,kph,xblk,xliq,xvap)
c
      else if (mode.eq.-1) then
c
c  call the main setup routine to specify mixture components, the file
c  containing the mixing parameters, and one of the standard reference
c  states  (call SETUP)
c
        ncc=NINT(xin(1))   !number of components
        iref=NINT(xin(2))  !code for reference state (see data stmt)
c
c  parse the input cstrin, first pick off path for the fluids directory
c  (inputs are separated by '|')
        do 820 i=1,254
        if (cstrin(i:i).eq.'|') then
          lendir=i-1           !length of the directory string
          hdir=cstrin(1:lendir)
          nexti=i+1
          goto 824
        endif
  820   continue
c  continue parsing cstrin and pick off file name containing mixture info
  824   i0=nexti
        do 840 i=i0,254
        if (cstrin(i:i).eq.'|') then
          hfmix=hdir(1:lendir)//cstrin(nexti:i-1)
          nexti=i+1
          goto 844
        endif
  840   continue
c  finally, pick off the names for the '.fld' files for each component
  844   i0=nexti
        jcomp=0                !counter for mixture component
        do 860 i=i0,254
        if (cstrin(i:i).eq.'|') then
          jcomp=jcomp+1
          if (jcomp.gt.ncc) goto 880
          hfiles(jcomp)=hdir(1:lendir)//cstrin(nexti:i-1)
          nexti=i+1
        endif
  860   continue
  880   continue
c  define the reference state and call the main setup routine
        hrf=hrefs(iref)
        call SETUP (ncc,hfiles,hfmix,hrf,ierr,cstr)
        xblk(1)=REAL(ncc)
c
      else if (mode.eq.-2) then
c
c  set the reference state for enthalpy and entropy (call SETREF)
c
        iref=NINT(xin(1))
        hrf=hrefs(iref)   !code for reference state (see data stmt)
        ixflag=NINT(xin(2)) !flag:  apply ref state to pures or mix
c  following inputs are meaninful only if hrf='OTH'
        t0=xin(3)
        p0=xin(4)
        h0=xin(5)
        s0=xin(6)
        do 920 i=1,nc
  920   x(i)=xin(6+i)
        call SETREF (hrf,ixflag,x,h0,s0,t0,p0,ierr,cstr)
c
      else if (mode.eq.-3) then
c
c  set the model (call SETMOD)
c
c  need to specify number of component as SETMOD called before SETUP
        ncc=NINT(xin(1))
        itype=NINT(xin(2))
        htype=htypes(itype)      !code for model (see data stmt)
        hmix=cstrin(1:3)         !pick off overall/mixture model
        do 940 i=1,ncc
        hcomp(i)=cstrin(i*3+1:i*3+3)  !pick off component models
  940   continue
        call SETMOD (ncc,htype,hmix,hcomp,ierr,cstr)
c
      else if (mode.eq.-4) then
c
c  set the mixing parameters (call SETKTV)
c
        hmodij=cstrin(1:3)     !3-letter code for mixing rule
        hfmix=cstrin(4:83)     !file containing mixture info
        icomp=NINT(xin(1))
        jcomp=NINT(xin(2))
        do 944 i=1,nmxpar
        fij(i)=xin(i+2)
  944   continue
        call SETKTV (icomp,jcomp,hmodij,fij,hfmix,ierr,cstr)
c
      else if (mode.eq.-5) then
c
c  retrieve the mixing parameters (call GETKTV)
c
        icomp=NINT(xin(1))
        jcomp=NINT(xin(2))
        call GETKTV (icomp,jcomp,hmodij,fij,hfmix,hfij,hbinp,hmxrul)
        cstr=hmodij//hbinp
        do 980 i=1,nmxpar
  980   xblk(i)=fij(i)
        ierr=0
c
      else
        ierr=99
        write (cstr,1099) mode,hnull
 1099   format ('[DLL error 99] mode passed to DLL is out of range:',
     &          i6,a1)
      end if
c  check for warnings
      ierri=ierr
      if (ierr.lt.0) then
        nwarn=1
        iwarn(1)=ierr
      end if
c
      RETURN
      end                                             !subroutine REFP61
c
c ======================================================================
c
      subroutine FLOUT2 (t,p,D,Dl,Dv,x,xl,xv,q,e,h,s,cv,cp,w,
     &                   xblk,xliq,xvap)
c
c  Take bulk-phase results of a flash calculation and, depending on
c  value of quality, calculate liquid- and/or vapor-phase properties and
c  write all quantities to the output arrays for return to the GUI.
c
c  This routine is similar to FLOUT, but outputs are in separate bulk,
c  liquid, and vapor arrays for compatibility with REFP61.
c
c  inputs:
c        t--temperature (K)
c        p--pressure (kPa)
c        D--overall (bulk) molar density (mol/L)
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 array (mol frac) for bulk composition
c       xl--composition array (mol frac) for liquid phase
c       xv--composition array (mol frac) for vapor phase
c           if only one phase is present, xl = xv = x
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 = 999 indicates supercritical state (t > Tc)
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, -9.992d is returned
c
c  outputs:
c     xblk--calculated (output) values for the bulk (overall)
c     xliq--calculated (output) values for the liquid phase
c     xvap--calculated (output) values for the vapor phase
c           there are 15 possible properties + the compositions, for most
c           cases, they are returned in the output array(s) in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  11-06-97  MM, original version; based on FLOUT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (narray=80)      !size of output (blk, liq, vap) arrays
      dimension x(ncmax),xl(ncmax),xv(ncmax)
      dimension xblk(narray),xliq(narray),xvap(narray)
      common /Gcnst/ R
      common /NCOMP/ nc
c  pointers to the output arrays for use with REFP61
      common /POINT2/ it,ip,irho,ivol,ie,ih,is,icv,icp,iw,
     &                iq,iz,ihjt,iA,iG,ikappa,ibeta,idPdD,id2P,idPdT,
     &                idDdT,idDdP,ist,ieta,itcx,idiel,ij1,ij2,ij3,ij4,
     &                ix,ifug
c  maximum number of properties returned to the GUI on a single call
      common /NPROP2/ nprops
c  common block containing flags to GUI
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  set bulk-phase outputs
c
      xblk(it)=t
      xblk(ip)=p
      xblk(irho)=D
      if (D.ge.1.0d-20) then
        xblk(ivol)=1.0d0/D
        xblk(iz)=p/(D*R*t)
        xblk(is)=s
      else
c  density is zero, volume & entropy are infinite
        xblk(ivol)=xinf
        xblk(iz)=1.0d0
        xblk(is)=xinf
      end if
      xblk(ie)=e
      xblk(ih)=h
      xblk(icv)=cv
      xblk(icp)=cp
      xblk(iw)=w
      xblk(iq)=q
c  the transport properties and surface tension are not calc here
      xblk(ieta)=xnotc
      xblk(itcx)=xnotc
      xblk(ist)=xnotc
      do 100 i=1,nc
      xblk(ix+i-1)=x(i)
  100 continue
c
      nval=nprops+nc
      if (abs(999.0d0-q).lt.1.0d-6) then
c
c  super-critical state (indicated by q = 999) set liquid and vapor
c  properties equal to bulk-phase properties
c
        do 199 j=1,nval
        xliq(j)=xblk(j)
        xvap(j)=xblk(j)
  199   continue
c
      else if (ABS(q).le.1.0d-8 .or. ABS(1.0d0-q).le.1.0d-8) then
c
c  saturation state, calculate the liquid and vapor-phase properties
c
        tv=t     !temperature and pressure for both phases always
        pv=p     !the same for a flash calculation
        call VLOUT2 (t,tv,p,pv,Dl,Dv,xl,xv,xliq,xvap)
c
      else if (q.lt.0.0d0) then
c
c  sub-cooled liquid state, set liquid properties equal to bulk-phase properties; 
c  vapor properties are set to flag indicating subcooled
c
        do 240 j=1,nval
        xliq(j)=xblk(j)
        xvap(j)=xsubc
  240   continue
c
      else if (q.gt.1.0d0) then
c
c  superheated vapor state, set vapor properties equal to bulk-phase
c  properties and liquid properties to flag indicating superheated
c
        do 280 j=1,nval
        xliq(j)=xsuph
        xvap(j)=xblk(j)
  280   continue
c
      else
c
c  two-phase state, calculate the liquid- and vapor-phase properties
c
        tv=t     !temperature and pressure for both phases always
        pv=p     !the same for a flash calculation
        call VLOUT2 (t,tv,p,pv,Dl,Dv,xl,xv,xliq,xvap)
c  bulk-phase heat capacity and speed of sound not defined for two-phase
        xblk(icv)=x2ph
        xblk(icp)=x2ph
        xblk(iw)=x2ph
      end if
c
      xliq(iq)=xnotd    !quality is not defined for liq, vap
      xvap(iq)=xnotd
c
      RETURN
      end                                             !subroutine FLOUT2
c
c ======================================================================
c
      subroutine STOUT2 (tl,tv,pl,pv,Dl,Dv,xl,xv,xblk,xliq,xvap)
c
c  Take results of a saturation calculation and calculate liquid and
c  vapor-phase properties and write all quantities to the output arrays
c  for return to the GUI.
c
c  This routine is similar to SATOUT, but outputs are in separate bulk,
c  liquid, and vapor arrays for compatibility with REFP61.
c
c  inputs:
c       tl--temperature (K) of the liquid phase
c       tv--temperature (K) of the vapor phase
c       pl--pressure (kPa) of the liquid phase
c       pv--pressure (kPa) of the vapor phase
c       Dl--molar density (mol/L) of the liquid phase
c       Dv--molar density (mol/L) of the vapor phase
c       xl--composition array (mol frac) for liquid phase
c       xv--composition array (mol frac) for vapor phase
c
c  outputs:
c     xblk--calculated (output) values for the bulk (overall)
c     xliq--calculated (output) values for the liquid phase
c     xvap--calculated (output) values for the vapor phase
c           there are 15 possible properties + the compositions, for most
c           cases, they are returned in the output array(s) in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  11-06-97  MM, original version; based on SATOUT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (narray=80)      !size of output (blk, liq, vap) arrays
      dimension xl(ncmax),xv(ncmax)
      dimension xblk(narray),xliq(narray),xvap(narray)
      common /NCOMP/ nc
c  pointers to the output arrays for use with REFP61
      common /POINT2/ it,ip,irho,ivol,ie,ih,is,icv,icp,iw,
     &                iq,iz,ihjt,iA,iG,ikappa,ibeta,idPdD,id2P,idPdT,
     &                idDdT,idDdP,ist,ieta,itcx,idiel,ij1,ij2,ij3,ij4,
     &                ix,ifug
c  maximum number of properties returned to the GUI on a single call
      common /NPROP2/ nprops
c  flags to GUI
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  the bulk-phase properties are not defined
      do 100 i=1,nprops
      xblk(i)=xnotd
  100 continue
c  set the bulk compositions equal to 'not defined' (this is an
c  arbitrary choice, but should return something to prevent GUI crash)
      do 120 i=1,nc
      xblk(ix+i-1)=xnotd
  120 continue
c
c  calculate the liquid- and vapor-phase properties
c
      call VLOUT2 (tl,tv,pl,pv,Dl,Dv,xl,xv,xliq,xvap)
c
      RETURN
      end                                             !subroutine STOUT2
c
c ======================================================================
c
      subroutine VLOUT2 (tl,tv,pl,pv,Dl,Dv,xl,xv,xliq,xvap)
c
c  Calculate liquid and vapor-phase properties and write all quantities
c  to the output arrays for return to the GUI.
c
c  This routine is similar to VLOUT, but outputs are in separate bulk,
c  liquid, and vapor arrays for compatibility with REFP61.
c
c  inputs:
c       tl--temperature (K) of the liquid phase
c       tv--temperature (K) of the vapor phase
c       pl--pressure (kPa) of the liquid phase
c       pv--pressure (kPa) of the vapor phase
c       Dl--molar density (mol/L) of the liquid phase
c       Dv--molar density (mol/L) of the vapor phase
c       xl--composition array (mol frac) for liquid phase
c       xv--composition array (mol frac) for vapor phase
c
c  outputs:
c     xliq--calculated (output) values for the liquid phase
c     xvap--calculated (output) values for the vapor phase
c           there are 15 possible properties + the compositions, for most
c           cases, they are returned in the output array(s) in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  11-06-97  MM, original version; based on VLOUT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (narray=80)      !size of output (blk, liq, vap) arrays
      dimension xl(ncmax),xv(ncmax)
      dimension xliq(narray),xvap(narray)
      common /Gcnst/ R
      common /NCOMP/ nc
c  pointers to the output arrays for use with REFP61
      common /POINT2/ it,ip,irho,ivol,ie,ih,is,icv,icp,iw,
     &                iq,iz,ihjt,iA,iG,ikappa,ibeta,idPdD,id2P,idPdT,
     &                idDdT,idDdP,ist,ieta,itcx,idiel,ij1,ij2,ij3,ij4,
     &                ix,ifug
c  maximum number of properties returned to the GUI on a single call
      common /NPROP2/ nprops
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  calculate the liquid-phase properties
c
      call THERM (tl,Dl,xl,ptherm,el,hl,sl,cvl,cpl,wl,hjt)
      xliq(it)=tl
      xliq(ip)=pl
      xliq(irho)=Dl
      if (Dl.ge.1.0d-20) then
        xliq(ivol)=1.0d0/Dl
        xliq(iz)=pl/(Dl*R*tl)
        xliq(is)=sl
      else
c  density is zero, volume & entropy are infinite
        xliq(ivol)=xinf
        xliq(iz)=1.0d0
        xliq(is)=xinf
      end if
      xliq(ie)=el
      xliq(ih)=hl
      xliq(icv)=cvl
      xliq(icp)=cpl
      xliq(iw)=wl
      xliq(iq)=xnotd    !quality is not defined
c  the transport properties and surface tension are not calc here
      xliq(ieta)=xnotc
      xliq(itcx)=xnotc
      xliq(ist)=xnotc
      do 100 i=1,nc
      xliq(ix+i-1)=xl(i)
  100 continue
c
c  calculate the vapor-phase properties
c
      call THERM (tv,Dv,xv,ptherm,ev,hv,sv,cvv,cpv,wv,hjt)
      xvap(it)=tv
      xvap(ip)=pv
      xvap(irho)=Dv
      if (Dv.ge.1.0d-20) then
        xvap(ivol)=1.0d0/Dv
        xvap(iz)=pv/(Dv*R*tv)
        xvap(is)=sv
      else
c  density is zero, volume & entropy are infinite
        xvap(ivol)=xinf
        xvap(iz)=1.0d0
        xvap(is)=xinf
      end if
      xvap(ie)=ev
      xvap(ih)=hv
      xvap(icv)=cvv
      xvap(icp)=cpv
      xvap(iw)=wv
      xvap(iq)=xnotd  !quality is not defined
c  the transport properties and surface tension are not calc here
      xvap(ieta)=xnotc
      xvap(itcx)=xnotc
      xvap(ist)=xnotc
      do 200 i=1,nc
      xvap(ix+i-1)=xv(i)
  200 continue
c
      RETURN
      end                                             !subroutine VLOUT2
c
c ======================================================================
c
      subroutine SPOUT2 (t,D,x,kph,xblk,xliq,xvap)
c
c  Calculate properties for the specified phase and write all quantities
c  to the output arrays for return to the GUI; values for the "other" phase
c  are set to "not calculated;" the bulk properties are set to those of
c  the specified phase.
c
c  inputs:
c        t--temperature [K]
c        D--molar density [mol/L]
c        x--composition array [mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c
c  outputs:
c     xliq--calculated (output) values for the liquid phase
c     xvap--calculated (output) values for the vapor phase
c           there are 15 possible properties + the compositions, for most
c           cases, they are returned in the output array(s) in the order:
c             T, P, rho, V (= 1/rho), e, h, s, Cv, Cp, w, quality,
c             viscosity, th cond, sur ten, Z (= PV/RT), composition(s)
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  11-06-97  MM, original version; based on SPOUT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (narray=80)      !size of output (blk, liq, vap) arrays
      dimension x(ncmax)
      dimension xblk(narray),xliq(narray),xvap(narray)
      common /Gcnst/ R
      common /NCOMP/ nc
c  pointers to the output arrays for use with REFP61
      common /POINT2/ it,ip,irho,ivol,ie,ih,is,icv,icp,iw,
     &                iq,iz,ihjt,iA,iG,ikappa,ibeta,idPdD,id2P,idPdT,
     &                idDdT,idDdP,ist,ieta,itcx,idiel,ij1,ij2,ij3,ij4,
     &                ix,ifug
c  maximum number of properties returned to the GUI on a single call
      common /NPROP2/ nprops
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
c  calculate the properties of the specified phase
c
      if (kph.eq.1) then
c  liquid phase
        call THERM (t,D,x,ptherm,e,h,s,cv,cp,w,hjt)
        xliq(it)=t
        xliq(ip)=ptherm
        xliq(irho)=D
        if (D.ge.1.0d-20) then
          xliq(ivol)=1.0d0/D
          xliq(iz)=ptherm/(D*R*t)
          xliq(is)=s
        else
c  density is zero, volume & entropy are infinite
          xliq(ivol)=xinf
          xliq(iz)=1.0d0
          xliq(is)=xinf
        end if
        xliq(ie)=e
        xliq(ih)=h
        xliq(icv)=cv
        xliq(icp)=cp
        xliq(iw)=w
        xliq(iq)=xnotc    !quality is not calculated
        do 100 i=1,nc
        xliq(ix+i-1)=x(i)
  100   continue
c  the transport properties and surface tension are not calc here
        xliq(ieta)=xnotc
        xliq(itcx)=xnotc
        xliq(ist)=xnotc
      else
c  vapor phase
        call THERM (t,D,x,ptherm,e,h,s,cv,cp,w,hjt)
        xvap(it)=t
        xvap(ip)=ptherm
        xvap(irho)=D
        if (D.ge.1.0d-20) then
          xvap(ivol)=1.0d0/D
          xvap(iz)=ptherm/(D*R*t)
          xvap(is)=s
        else
c  density is zero, volume & entropy are infinite
          xvap(ivol)=xinf
          xvap(iz)=1.0d0
          xvap(is)=xinf
        end if
        xvap(ie)=e
        xvap(ih)=h
        xvap(icv)=cv
        xvap(icp)=cp
        xvap(iw)=w
        xvap(iq)=xnotc    !quality is not calculated
        do 120 i=1,nc
        xvap(ix+i-1)=x(i)
  120   continue
c  the transport properties and surface tension are not calc here
        xvap(ieta)=xnotc
        xvap(itcx)=xnotc
        xvap(ist)=xnotc
      end if
c  set the bulk properties to those of the specified phase
      if (kph.eq.1) then
        do 300 i=1,nprops
        xblk(i)=xliq(i)
  300   continue
        do 320 i=1,nc
        xblk(ix+i-1)=xliq(ix+i-1)
  320   continue
      else
        do 360 i=1,nprops
        xblk(i)=xvap(i)
  360   continue
        do 380 i=1,nc
        xblk(ix+i-1)=xvap(ix+i-1)
  380   continue
      end if
c
      RETURN
      end                                             !subroutine SPOUT2
c
c ======================================================================
c
      block data BDDLL
c
c  This block data initializes the common blocks used in the DLL routines.
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  08-02-96  MM, original version
c  05-16-97  MM, revise flags for special cases, put in /FLAGS/
c  07-15-97  MM, move /FLAGS/ from here to setup.f
c  12-02-97  MM, move reference state choices and model types here
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (nrefst=6)       !number of reference state choices
      character*3 hrefs,htypes
c
c  pointers to the xout array (to aid in adding possible future props)
      common /IPOINT/
     &     it,ip,irho,ivol,ie,ih,is,icv,icp,iw,iq,ieta,itcx,ist,iz,ix
c  maximum number of properties returned to the GUI on a single call
      common /NPROP/ nprops
c  flags to GUI indicating 'not applicable', '2-phase', etc.
c     common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c  reference state choices and model types
      common /HREFST/ hrefs(nrefst)
      common /HMODTY/ htypes(4)
c
      data it,ip,irho,ivol,ie,ih,is,icv,icp,iw,iq,ieta,itcx,ist,iz,ix
     &    / 1, 2,  3,   4,  5, 6, 7,  8,  9,10,11,  12,  13, 14,15,16/
      data nprops /15/
c  the following declarations have been moved to BDSET in setup.f
c     data xnota /-9.991d6/  !flag indicating 'not applicable'
c     data x2ph  /-9.992d6/  !flag indicating '2-phase'
c     data xsubc /-9.993d6/  !flag indicating 'subcooled'
c     data xsuph /-9.994d6/  !flag indicating 'superheated'
c     data xsupc /-9.995d6/  !flag indicating 'supercritical'
c     data xinf  /-9.996d6/  !flag indicating 'infinite'
c     data x7    /-9.997d6/  !not used
c     data xnotd /-9.998d6/  !flag indicating property 'not defined'
c     data xnotc /-9.999d6/  !flag indicating property 'not calculated'
c
      data hrefs /'NBP','ASH','IIR','OTH','DEF','MIX'/
      data htypes/'EOS','ETA','TCX','STN'/
c
      end                                              !block data BDDLL
c
c ======================================================================
c
      block data BDDLL2
c
c  This block data initializes the common blocks used in the DLL routines.
c
c  written by M. McLinden, NIST Physical & Chem Properties Div, Boulder, CO
c  11-07-97  MM, original version; based on BDDLL
c  12-02-97  MM, move reference state choices and model types here
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)        !max number of components in mixture
      parameter (nrefst=6)       !number of reference state choices
      character*3 hrefs,htypes
c
c  pointers to the output arrays for use with REFP61
      common /POINT2/ it,ip,irho,ivol,ie,ih,is,icv,icp,iw,
     &                iq,iz,ihjt,iA,iG,ikappa,ibeta,idPdD,id2P,idPdT,
     &                idDdT,idDdP,ist,ieta,itcx,idiel,ij1,ij2,ij3,ij4,
     &                ix,ifug
c  maximum number of properties returned to the GUI on a single call
      common /NPROP2/ nprops
c  reference state choices and model types
      common /HREFS2/ hrefs(nrefst)
      common /HMODT2/ htypes(4)
c
      data it,ip,irho,ivol,ie,ih,is,icv,icp,iw
     &    / 1, 2,  3,   4,  5, 6, 7,  8,  9,10/
      data iq,iz,ihjt,iA,iG,ikappa,ibeta,idPdD,id2P,idPdT
     &    /11,12, 13, 14,15,  16,   17,    18,  19,  20/
      data idDdT,idDdP,ist,ieta,itcx,idiel,ij1,ij2,ij3,ij4
     &    /21,     22,  23, 24,  25,  26,  27, 28, 29, 30/
      data ix,ifug
     &    /31, 36/    !ifug = 31+ncmax
c
      data hrefs /'NBP','ASH','IIR','OTH','DEF','MIX'/
      data htypes/'EOS','ETA','TCX','STN'/
c
      data nprops /30/     !30 total properties, including 4 spares
c
      end                                             !block data BDDLL2
c
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                     end file ftn_pas.f
c ======================================================================