c  begin file flsh_sub.f
c
c  This file contains iterative flash routines which call the
c  intermediate-level routines
c
c  contained here are:
c     subroutine TPRHO (t,p,x,kph,kguess,rho,ierr,herr)
c     subroutine TPFLSH (t,p,z,D,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
c     subroutine TPFL2 (t,p,z,Dl,Dv,x,y,q,ierr,herr)
c     subroutine TDFLSH (t,D,z,p,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
c     subroutine TDFL2 (t,D,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
c    &                  p,Dl,Dv,x,y,q,ierr,herr)
c     subroutine PDFLSH (p,D,z,t,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
c     subroutine PDFL1 (p,rho,x,t,ierr,herr)
c     subroutine PDFL2 (p,d,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
c    &                  t,Dl,Dv,x,y,q,ierr,herr)
c     subroutine PHFLSH (p,h,z,t,D,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr)
c     subroutine PHFL1 (p,h,x,kph,t,D,ierr,herr)
c     subroutine PHFL2 (p,h,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
c    &                  t,Dl,Dv,x,y,q,ierr,herr)
c     subroutine PSFLSH (p,s,z,t,D,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr)
c     subroutine PSFL1 (p,s,x,kph,t,D,ierr,herr)
c     subroutine PSFL2 (p,s,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
c    &                  t,Dl,Dv,x,y,q,ierr,herr)
c     subroutine PEFLSH (p,e,z,t,D,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
c     subroutine PEFL1 (p,e,x,kph,t,D,ierr,herr)
c     subroutine PEFL2 (p,e,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
c    &                  t,Dl,Dv,x,y,q,ierr,herr)
c     subroutine PBFLSH (p,b,z,ab,t,D,Dl,Dv,x,y,q,e,h,s,cv,cp,w,
c    &                   ierr,herr)
c
c  these routines use the following common blocks from other files
c     common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c     common /NCOMP/ nc,ic
c     common /HCHAR/ htab,hnull
c
c  various arrays are dimensioned with parameter statements
c     parameter (ncmax=20)        !max number of components in mixture
c     parameter (nrefmx=10)       !max number of fluids for transport E
c     parameter (n0=-ncmax-nrefmx,nx=ncmax)
c
c ======================================================================
c ======================================================================
c
      subroutine TPRHO (t,p,x,kph,kguess,rho,ierr,herr)
c
c  iterate for density as a function of temperature, pressure, and
c  composition for a specified phase
c
c***********************************************************************
c  WARNING:
c  Invalid densities will be returned for T & P outside range of validity,
c  i.e., pressure > melting pressure, pressure less than saturation
c  pressure for kph=1, etc.
c
c***********************************************************************
c  inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c        x--composition [array of mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c                 N.B.:  0 = stable phase--NOT ALLOWED (use TPFLSH)
c                            (unless an initial guess is supplied for rho)
c                       -1 = force the search in the liquid phase
c                       -2 = force the search in the vapor phase
c   kguess--input flag:  1 = first guess for rho provided
c                        0 = no first guess provided
c      rho--first guess for molar density [mol/L], only if kguess = 1
c
c  outputs:
c      rho--molar density [mol/L]
c     ierr--error flag:  0 = successful
c                      200 = CRITP did not converge
c                      201 = illegal input (kph <= 0)
c                      202 = liquid-phase iteration did not converge
c                      203 = vapor-phase iteration did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  04-03-95  MM, original version
c  09-11-95  MM, add error string to argument list
c  10-11-95  MM, RETURN if any error detected
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  01-22-95  MM, check that input rho > 0 before using as initial guess
c  01-23-95  MM, use Rackett technique for liquid density initial guess
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  03-07-96  MM, change tolr from 1.0d-6 to 1.0d-8
c  10-16-96  MM, change call from DPRHO to DPDD
c  12-02-96  MM, add check that vapor density guess does not give p < 0
c  01-07-97  MM, error message bug (do not concatenate herr with itself)
c  02-11-97 EWL, reduced tolerance after 10 and after 15 iterations
c  04-22-97  MM, add additional iteration for low-p (allow p = 0)
c  05-16-97  MM, do not apply step limitation for t > 1.5*tc
c  06-03-97 EWL, add second order Newton's solution to increase convergence speed
c  06-05-97 EWL, remove checks in liquid solution for large jumps in density
c                change initial guess in liquid phase to 2*Dc
c  07-15-97  MM, renumber and add detail to error messages
c  09-22-97 EWL, modify test to select liquid or vapor iteration
c  09-23-97  MM, increase itmax from 20 to 30 to aid in near-critical convergence
c  10-01-97  MM, add compiler switch to allow access by DLL
c  11-13-97 EWL, revert to 1st order Newton's method if 2nd order gives bad result
c  11-18-97 EWL, ditto for liquid-phase iteration
c  12-31-97  MM, check derivative dP/dD for liquid-phase initial guess
c  02-11-98  MM, check for t > Tc before switching liq to vapor iteration
c  03-31-98 EWL, check for Dguess>Dmax in liquid iteration
c  08-14-98  MM, do not check input rho<0 unless kguess=1 (@stmt 10)
c  09-30-98 EWL, if liquid iteration fails and t>tc, call vapor iteration
c  12-18-98 EWL, increase max value of fvdpl to 2 on check of 1st order Newton's
c  12-22-98 EWL, recalculate R for mixtures based on values for pure fluids
c  02-02-99 EWL, increase max value of fvdpl to 25 on check of 1st order Newton's
c  11-09-00 EWL, allow kph=0 if an initial guess is supplied for rho
c  01-17-02 EWL, allow kph=-1 and -2 to force searches in liquid or vapor
c  02-13-02 EWL, change pressure bound where the ideal gas law kicks in
c  11-02-04 EWL, check for p=0 and kph=1
c  07-20-05 EWL, check for very large rho
c  10-03-05 EWL, change itmax from 30 to 40
c  06-07-06 EWL, slight changes to check for T=Tc or P=Pc
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: TPRHO
c     dll_export TPRHO
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport E
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      character*1 htab,hnull
      character*255 herr,herr1
      dimension x(ncmax)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
      common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
      itmax=40
      tolr=1.0d-8
      ierr=0
      herr=' '
c     write (*,1001) t,p,(x(i),i=1,5)
c1001 format (1x,' TPRHO--input t,p,x:           ',f10.5,f14.6,5f10.6)
c
c  determine which iteration to use:
c  For liquids (or fluids above the critical pressure) the iteration is
c  carried out in transformed coordinates of log(V).  For vapor (or
c  fluids at supercritical temperatures but pressures below the critical
c  value) the iteration is in terms of log(V) and log(p).  For very low
c  pressure vapors a simple iteration in rho and p is used.  The iteration
c  has converged when the pressure calculated from the equation of state
c  agrees with the input pressure within a convergence tolerance.
c
      call Rmix (x)
      if (p.lt.1.0d-8) then
        if (kph.ne.1) then
c  ideal-gas
          if (t.gt.0.d0) rho=p/(R*t)
          RETURN
        elseif (ABS(p).lt.1.d-16) then
          rho=0.d0
          RETURN
        endif
      end if
c
c  call Peng-Robinson routines if in use
c
      i=-1
      call PREOS(i)
      if (i.eq.1) then
        call TPRHOPR (t,p,x,rho1,rho2)
        if (rho1.gt.1.d8) goto 100
        if (rho2.lt.1.d-12) then
          rho=rho1
        else
          if (abs(kph).eq.1) then
            rho=rho1
          elseif (abs(kph).eq.2) then
            rho=rho2
          else
            if (ABS(rho-rho1).lt.ABS(rho-rho2)) then
              rho=rho1
            else
              rho=rho2
            endif
          endif
        endif
c  if the pressure is less than 1 kPa and the density is in the liquid phase,
c  continue on with iterative routine to make sure that the cubic solver got
c  the right answer
        if (p.gt.1.d0 .or. rho.lt.1.d0 .or. kph.eq.2) RETURN
      endif
c
      call CRITP (x,tc,pc,rhoc,ierr,herr1)
      if (ierr.ne.0) then
        ierr=200
        write (herr,1200) herr1(1:236),hnull
 1200   format ('[TPRHO error 200] ',a236,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
      vc=1.0d0/rhoc
      if (vc.le.0.d0) vc=1.d0
      vclog=LOG(vc)
      if (kph.eq.-1) then
        lliq=.true.
      elseif (kph.eq.-2) then
        lliq=.false.
c  EWL modification to aid near-critical convergence;
c  following approximates the pressure along the critical isochore
      elseif (p.lt..99*pc*(1.0d0+7.0d0*(t/tc-1.0d0)) .and. t.gt.tc) then
c  use 'vapor phase' iteration
        lliq=.false.
      else if (p.ge.pc .or. t.ge.tc) then
c  use 'liquid phase' iteration
        lliq=.true.
      else if (kph.le.0) then
        if (kguess.ne.0) then
          lliq=.true.
          if (rho.lt.rhoc) lliq=.false.
        else
c  illegal input to TPRHO (two-phase state)
          ierr=201
          herr='[TPRHO error 201] illegal input to TPRHO (kph <= 0); '
     &           //'use TPFLSH instead'//hnull
          call ERRMSG (ierr,herr)
          RETURN
        endif
      else
c  use input phase specification to find solution
        if (kph.eq.1) then
          lliq=.true.
        else
          lliq=.false.
        end if
      end if
c
c  set initial guess for density
      kgrho=1      !flag is set to 0 if input rho<0
      vlog=0.d0
 10   continue
      if (kguess.eq.1) then
        if (rho.gt.0.0d0) then
c  use input density as initial guess
          vlog=LOG(1.0d0/rho)
        else
c  set flag to trigger initial guess algorithm below
          kgrho=0
        end if
      end if
      if (kguess.eq.0 .or. kgrho.eq.0) then
c  initial guess based on region
        if (p.gt.pc .and. t.gt.tc) then
          vlog=LOG(0.5d0*vc)
        else if (lliq) then
c         vlog=LOG(0.5d0*vc)
c  initial guess for sub-critical liquid density based on modified
c  Rackett technique (Reid, Prausnitz, and Poling (1987), Properties
c  of Gases and Liquids, 4th edition) except that Zra = Zc
          Rtp=R*tc/pc
          Zc=vc/Rtp
c         vlog=LOG(Rtp*Zc**(1.0d0+(1.0d0-t/tc)**(2.0d0/7.0d0)))
          Dguess=1.0d0/(Rtp*Zc**(1.0d0+ABS(1.0d0-t/tc)**(2.0d0/7.0d0)))
          call LIMITS ('EOS',x,tmin,tmax,Dmax,pmax)
          if (Dguess.gt.Dmax) Dguess=Dmax
          iguess=1
 20       continue
          call DPDD (t,Dguess,x,dpdrho)
          if (dpdrho.le.0.0d0 .and. iguess.le.8) then
c  initial guess is in two-phase region, set to higher density
            Dguess=1.1d0*Dguess
            iguess=iguess+1
            goto 20
          end if
          vlog=LOG(1.0d0/Dguess)
          if (vlog.gt.LOG(0.5d0*vc)) vlog=LOG(0.5d0*vc)
        else
          vlog=LOG(R*t/p)
        end if
        if (.not.lliq) vlog=LOG(R*t/p)
      end if
c     write (*,1082) kph,kguess,t,p,1.0/EXP(vlog),(x(i),i=1,nc)
c1082 format (1x,' TPRHO--kph,kguess,t,p,rho,x: ',2i4,3e14.6,5f10.5)
c
c  enter Newton's method iteration, separate loops for liquid and vapor
c
      rho2=0.0d0
      if (lliq) then
c  liquid phase iteration
        do it=1,itmax
          vlog0=vlog
c         write (*,*) ' TPRHO (liquid) vlog: ',vlog
          rho=1.0d0/exp(vlog)
          if (rho.gt.1.d8) goto 100
          call PRESS (t,rho,x,p2)
          call DPDD (t,rho,x,dpd)
          call DPDD2 (t,rho,x,dpd2)
c         write (*,1010) it,t,x(1),p,p2,dpd,rho
c1010     format (1x,'it,t,x(1),p,p2,dpd,rho (liquid): ',i3,2f8.3,4e16.8)
          if (dpd.lt.0.0d0) then
c  unstable portion of two-phase region, make another guess
c           write (*,1012) it,1.0/exp(vlog)
c1012       format (1x,'unstable 2-phase in liquid, it, rho = ',i4,d16.8)
            vlog=vlog-0.10d0
          else
            dpd2=dpd2+dpd/rho              !2nd order Newton's method
            fvdp=(p-p2)/(dpd2*(p-p2)/2.0d0/dpd+dpd)/rho
            if (ABS(fvdp).gt.100.0d0) then
c  revert to 1st order Newton's method if too large step from 2nd order
              dpdlv=-rho*dpd
              fvdp=(p2-p)/dpdlv            !1st order Newton's method
            endif
c           write (*,1014) dpdlv,dpd,fvdp
c1014       format (1x,' TPRHO (liquid) dpdlv,dpd,fvdp: ',3d16.8)
c  if calculation has not converged after 10 or 15 iterations, loosen
c  tolerance (original tolr was too tight near critical)
            if (it.eq.10) tolr=tolr*10.d0
            if (it.eq.15) tolr=tolr*10.d0
c  The rho-rho2 check is only important for very low pressures on the
c  liquid surface (propane or R124).  In this case, the change in density
c  required to get the correct pressure is less than machine precision.
            if (ABS(fvdp/p).lt.0.0001d0*tolr .or.
     &          ABS(rho-rho2).lt.1.0d-11) then
c  iteration has converged
              rho=1.0d0/exp(vlog-fvdp)
c             call PRESS (t,rho,x,p2)        !convergence testing only
c             write (*,1016) it,t,p,p2,rho
c1016         format (1x,8x,'% TPRHO converged (liq); it,t,p,p2,rho  ',
c    &              i3,f8.3,2e18.10,e18.10)
c             write (*,*) ' TPRHO final rho-liq: ',rho
              RETURN
            else
c  next guess
              vlog=vlog0-fvdp
c  do not allow too great of a change in density between iterations
c  for liquid states
              if (ABS(vlog-vlog0).gt.0.1d0 .and. t.lt.1.5d0*tc) then
                vlog=vlog0+SIGN(0.1d0,vlog-vlog0)
              end if
              if (vlog.gt.vclog .and. t.lt.tc) then
                vlog=0.5d0*(vlog0+vclog)
              end if
c  switch to the "vapor" iteration
              if (vlog.lt.-5.0d0 .and. t.ge.tc) then
c               write (*,*) ' TPRHO--switching to the vapor iteration'
                lliq=.false.
                goto 10
              end if
            end if
            rho2=rho
          end if
        enddo
c  iteration has not converged
        rho=1.0d0/exp(vlog)
        if (t.ge.tc) then
          lliq=.false.
          goto 10
        endif
        goto 100
c
c      else if (lowp) then
c   low-pressure vapor iteration
c   disabled--not required
c        do it=1,itmax
c          call PRESS (t,rho,x,p2)
c          call DPDD (t,rho,x,dpd)
c          delp=p2-p
c          if (ABS(delp).lt.0.0001*tolr) then
c   iteration has converged
c            write (*,1024) it,t,p,p2,delp,rho
c 1024       format (9x,'% TPRHO converged (low-p); it,t,p,p2,delp,rho  ',
c     &              i3,f8.3,3e18.10,e18.10)
c            RETURN
c          else
c   next guess
c            rho=rho-delp/dpd
c          end if
c        enddo
c   iteration has not converged
c        rho=p/(R*t)
c        ierr=1
c        write (herr,1002) t,p,rho,hnull
c        call ERRMSG (ierr,herr)
c 1002   format (' ERROR--TPRHO has not converged (low-p),  t,p,rho =',
c     &          f8.3,2e14.6,a1)
c        RETURN
c
      else
c  vapor phase iteration
        plog=LOG(p)
        do it=1,itmax
          vlog0=vlog
          rho=1.0d0/exp(vlog)
          call PRESS (t,rho,x,p2)
          call DPDD (t,rho,x,dpd)
          call DPDD2 (t,rho,x,dpd2)
c         write (*,1030) it,t,x(1),p,p2,dpd,rho
c1030 format (1x,'it,t,x(1),p,p2,dpd,rho (vapor):  ',i3,2f8.3,4e18.10)
          if (dpd.lt.0.0d0 .or. p2.le.0.0d0) then
c  unstable portion of two-phase region or p2<0, make another guess
c         write (*,1032) it,1.0/exp(vlog)
c1032     format (1x,'unstable 2-phase in vapor or p2 < 0, it, rho = ',
c    &            i4,d16.8)
            vlog=vlog+0.10d0
          else

            dpd2=dpd2/p2-(dpd/p2)**2+dpd/rho/p2
            fvdpl=(plog-LOG(p2))/(p2*dpd2*(plog-LOG(p2))
     &            /2.0d0/dpd+dpd/p2)/rho      !2nd order Newton's method
c  if 2nd order Newton's method gives unreasonable result, revert to
c  first order method (large value of fvdpl corresponds to huge change
c  in next guess for volume)
            if (ABS(fvdpl).gt.1.0d0) then
              dpdlv=-rho*dpd
              fvdpl=(LOG(p2)-plog)*p2/dpdlv   !1st order Newton's method
            endif
c  the following used was original set at '.gt.2.d0', but was changed to
c  '.gt.25.d0'.  Vapor densities were failing to converg for a mixture of
c  47% methane, 16% ethane, 7% propane, and 30% i-butane.  This may now cause
c  other areas to fail.
            if (ABS(fvdpl).gt.25.d0) fvdpl=0.1d-5
c  if calculation has not converged after 10 or 15 iterations, loosen
c  tolerance (original tolr was too tight near critical)
            if (it.eq.10) tolr=tolr*10.d0
            if (it.eq.15) tolr=tolr*10.d0
            if (ABS(fvdpl).lt.0.0001d0*tolr) then
c  iteration has converged
              rho=1.0d0/exp(vlog)
              call PRESS (t,rho,x,p2)        !convergence testing only
c              write (*,1034) it,t,p,p2,rho
c1034       format (1x,8x,'% TPRHO converged (vap); it,t,p,p2,rho  ',
c    &              i3,f8.3,2e18.10,e18.10)
c              write (*,*) ' TPRHO final rho-vap: ',rho
              RETURN
            else
c  next guess
              vlog=vlog-fvdpl
c             if (ABS(vlog-vlog0).gt.0.5d0) then
c               vlog=vlog0+SIGN(0.5d0,vlog-vlog0)
c             end if
c             if (vlog.lt.vclog .and. t.lt.tc) then
c               vlog=0.5d0*(vlog0+vclog)
c             end if
            end if
          end if
        enddo
c  iteration has not converged
        rho=1.0d0/exp(vlog)
        ierr=203
        write (herr1,1203) t,p/1000.d0,rho,(x(j),j=1,MIN(nc,5))
 1203   format ('[TPRHO error 203] vapor iteration has not ',
     &          'converged for T =',g11.5,' K, P =',g11.5,
     &          ' MPa, rho (last guess) = ',g11.5,
     &          ' mol/L, x (mol frac) =',5(0pf8.5))
        herr=herr1(1:254)//hnull
        call ERRMSG (ierr,herr)
c       write (*,*) ' TPRHO final rho-vap (not converged): ',rho
        RETURN
c
      end if
      RETURN
c
 100  continue
      ierr=202
      write (herr1,1202) t,p/1000.d0,rho,(x(i),i=1,MIN(nc,5))
 1202 format ('[TPRHO error 202] liquid iteration has not ',
     &        'converged for T =',g11.5,' K, P =',g11.5,
     &        ' MPa, rho (last guess) = ',g11.5,' mol/L,',
     &        ' compositions = ',5(0pf9.5))
      herr=herr1(1:193)//hnull
      call ERRMSG (ierr,herr)
c     write (*,*) ' TPRHO final rho-liq (not converged): ',rho
      RETURN
c
      end                                              !subroutine TPRHO
c
c ======================================================================
c
      subroutine TPFLSH (t,p,z,D,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
c
c  flash calculation given temperature, pressure, and bulk composition
c
c  This routine accepts both single-phase and two-phase states as the
c  input; if the phase is known, the subroutine TPRHO is faster.
c
c  inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c        z--overall (bulk) composition [array of mol frac]
c
c  outputs:
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 of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c           if only one phase is present, x = y = z
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c           q < 0 indicates subcooled (compressed) liquid
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q > 1 indicates superheated vapor
c           q = 998 superheated vapor, but quality not defined (t > Tc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful
c                        1 = Tin < Tmin
c                        4 = Pin < 0
c                        5 = T and P out of range
c                        8 = x out of range (component and/or sum < 0
c                            or > 1)
c                        9 = x and T out of range
c                       12 = x out of range and P < 0
c                       13 = x and T and P out of range
c                      210 = CRITP did not converge
c                      211 = SATT did not converge at bubble point
c                      212 = SATT did not converge at dew point
c                      213 = TPRHO did not converge for liquid state
c                      214 = TPRHO did not converge for vapor state
c                      215 = TPRHO did not convg for supercritical state
c                      216 = TPRHO did not convg for liq in 2-phase it
c                      217 = TPRHO did not convg for vap in 2-phase it
c                      218 = liquid frac for 2-phase it did not converge
c                      219 = composition for 2-phase it did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  09-11-95  MM, original version
c  09-12-95  MM, change argument list (add e, h, s, etc.)
c  09-20-95  MM, define output xl, xv for single-component super-critical
c  09-25-95  MM, rearrange argument list (outputs in order rho, x, q)
c  10-11-95  MM, RETURN if any error detected
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  12-12-95  MM, restructure logic to accommodate mixtures
c                first guess for 2-phase mix from ratio of dew & bubble
c  12-13-95  MM, set undefined output compositions to zero
c                initial rho guess for 1-phase vapor based on p/pdew
c  12-14-95  MM, bug on call to CRITP:  pass z (not x)
c  12-28-95  MM, add full 2-phase mixture iteration using fugacities
c  12-29-95  MM, move 2-phase iteration to separate routine TPFL2
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  05-31-96  MM, check input t,p,z against limits
c  06-03-96  MM, add 'EOS' to calling list for LIMITX
c  01-07-97  MM, error message bug (do not concatenate herr with itself)
c  04-09-97  MM, initialize x,y to z
c  04-22-97  MM, use h, rather than s, to compute q (h converges at p = 0)
c  05-15-97  MM, use V to compute q for comp. liq. (T,h can be double-valued)
c                add q = 998 case for t > Tc, but p < Pc
c  07-15-97  MM, renumber and add detail to error messages
c                get flags for 'not defined' from common /FLAGS/
c  10-01-97  MM, add compiler switch to allow access by DLL
c  02-11-98  MM, check that computed density is within bounds
c  09-15-98 EWL, if T>Tc-5, use critical parameters if call to SATT fails.
c  12-28-98 EWL, if SATT fails (vapor), check for liquid state before failing
c  03-05-99 EWL, add check for region below Ttriple in the vapor phase
c  10-21-99 EWL, do not check satt if pressure is less than .8pc at tc,
c                .1pc at .85tc and the area in between (on a log plot).
c  01-11-00 EWL, remove ierr from line 1008 and several other places in flsh_sub.for
c  01-02-01 EWL, add check for d<dl and d>dv in single phase routine
c  05-21-01 EWL, call PSATT if an ancillary vapor pressure eq. exists (for nc=1)
c  09-11-02 EWL, drop psat by factor of 10 to check for two phase before doing 1 phase
c  12-16-02 EWL, add checks for pseudo-pure fluid
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: TPFLSH
c     dll_export TPFLSH
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport E
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      character*1 htab,hnull
      character*3 hps,hpsk
      character*255 herr,herr1,herr2,herrl
      dimension z(ncmax),x(ncmax),y(ncmax)
      dimension xdew(ncmax),ydew(ncmax),ybub(ncmax),xbub(ncmax)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
      common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
      common /PSMOD/ hps,hpsk(n0:nx)
c  flags to GUI indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      common /FLAGS2/ iamwat,ianc,iwat
      common /prnterr/ iprnterr
c
      icomp=0
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
      endif
      ierr=0
      pbub=0
      pdew=0
      herr=' '
c
c  initialize output liquid and vapor compositions to input values
c  zero output values for undefined components
      if (icomp.eq.0) then
        do i=1,nc
          x(i)=z(i)
          y(i)=z(i)
        enddo
      else
        x(icomp)=1.d0
        y(icomp)=1.d0
      endif
c     if (nc.lt.ncmax) then
c       do i=nc+1,ncmax
c         x(i)=0.0d0
c         y(i)=0.0d0
c       enddo
c     end if
      if (p.lt.1.0d-14) then
c  ideal-gas
        D=p/(R*t)
        Dl=D
        Dv=D
        call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
        RETURN
      endif
c
c  check that input conditions are within limits
c
      rhodum=0.0d0
      call LIMITX ('EOS',t,rhodum,p,z,tmin,tmax,rhomax,pmax,ierrl,herrl)
      if (ierrl.lt.0) then
c  one or inputs are outside limits--if just a warning proceed w/ calc
        write (herr,1006) ierrl,herrl(1:234),hnull
 1006   format ('[TPFLSH warning',i3,'] ',a234,a1)
        call ERRMSG (ierrl,herr)
c  if error (as opposed to warning) set output density to zero and return
      else if (ierrl.gt.0) then
        write (herr,1007) ierrl,herrl(1:236),hnull
 1007   format ('[TPFLSH error',i3,'] ',a236,a1)
        call ERRMSG (ierrl,herr)
        D=0.0d0
        Dl=0.0d0
        Dv=0.0d0
        if (t.gt.0) call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
        q=999.d0     !quality undefined
        ierr=ierrl
        RETURN
      end if
c
      call CRITP (z,tc,pc,rhoc,ierr,herr2)
      if (ierr.ne.0) then
        ierr=210
        write (herr,1008) herr2(1:235),hnull
 1008   format ('[TPFLSH error 210] ',a235,a1)
        call ERRMSG (ierr,herr)
        q=999.d0     !quality undefined
        RETURN
      end if
c
 150  continue
      psat=9999999.d0
      if (t.le.tc) psat=exp(13.85d0*(t/tc)-14.073d0)*pc
      if (icomp.ne.0 .and. t.le.tc) then
        ierr=1
        if (hpsk(icomp).ne.'NBS') call PSATT (t,x,psat2,ierr,herr)
        if (ierr.eq.0) psat=0.95d0*psat2 !Drop 5% to allow for error
      else
c  drop cutoff for single phase vapor check by factor of 10
        psat=psat/10.d0
      endif
      if (icomp.ne.0 .and. t.lt.ttp(icomp)) then
c  temperature less than triple point temperature, but in the gas phase.
        kph=2
        kguess=1
        D=p/(R*t)
        call TPRHO (t,p,z,kph,kguess,D,ierr,herr2)
        Dl=D
        Dv=D
        if (ierr.ne.0) then
          ierr=214
          write (herr,1214) herr2(1:190),hnull
          call ERRMSG (ierr,herr)
          q=999.d0     !quality undefined
          RETURN
        end if
      elseif (t.ge.tc .or. p.lt.psat) then
c  supercritical and vapor states (x = y = z as set above)
c       write (*,*) ' TPFLSH--supercritical'
        kph=2
        kguess=0
        call TPRHO (t,p,z,kph,kguess,D,ierr,herr2)
        call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
        Dl=D
        Dv=D
        if (p.le.pc) then
          q=998.d0     !superheated, but cannot define quality
        else
          q=999.d0     !quality undefined
        end if
        if (ierr.ne.0) then
          ierr=215
          write (herr,1215) herr2(1:175),hnull
 1215     format ('[TPFLSH error 215] supercritical or vapor density ',
     &            'iteration did not converge:  ',a175,a1)
          call ERRMSG (ierr,herr)
c  call two phase states if t is close to tc for a mixture
          if (icomp.eq.0 .and. t.lt.tc*1.1d0) then
            tc=tc*1.1d0
            goto 150
          endif
          RETURN
        end if
c       do i=1,nc
c         x(i)=z(i)
c         y(i)=z(i)
c       enddo
      else
c  subcritical state--call saturation routine to determine liq or vap
        i=iprnterr
        iprnterr=0
        call SATT (t,z,2,pdew,Dldew,Dvdew,xdew,ydew,ierr,herr2)
        iprnterr=i
        iflag=0
        if (ierr.ne.0 .and. t.gt.tc-5) then
          pdew=pc
          Dldew=rhoc*2
          Dvdew=rhoc/2
        elseif (ierr.ne.0 .and. icomp.ne.0) then
          ierr=212
          write (herr,1212) herr2(1:193),hnull
 1212     format ('[TPFLSH error 212] dew point calculation ',
     &            'did not converge:  ',a193,a1)
          call ERRMSG (ierr,herr)
          D=0.0d0
          Dl=0.0d0
          Dv=0.0d0
          q=999.d0     !quality undefined
          RETURN
        else if (ierr.ne.0) then
c  do not exit yet, check if valid liquid state
          iflag=1
          pdew=0
        end if
        if (p.lt.pdew) then
c  single-phase vapor (x = y = z as set above)
c         write (*,*) ' TPFLSH--single-phase vapor'
          kph=2
          kguess=1
          D=Dvdew*p/pdew                     !initial guess for density
          call TPRHO (t,p,z,kph,kguess,D,ierr,herr2)
          Dl=D
          Dv=D
          if (ierr.ne.0) then
            ierr=214
            write (herr,1214) herr2(1:190),hnull
 1214       format ('[TPFLSH error 214] vapor density iteration ',
     &              'did not converge:  ',a190,a1)
            call ERRMSG (ierr,herr)
            q=999.d0     !quality undefined
            RETURN
          end if
          call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
          q=998.d0
          if (ianc.eq.0) then
            call ENTHAL (t,Dldew,xdew,hldew)
            call ENTHAL (t,Dvdew,z,hvdew)
            q=(h-hldew)/(hvdew-hldew)
            if (q.gt.0.0d0 .and. q.lt.1.0d0) q=998.d0
          endif
        else
          if (icomp.ne.0 .and. ianc.eq.0) then
c  special case:  pure component single-phase liquid
c           write (*,*) ' TPFLSH--pure fluid single-phase liquid'
            kph=1
            kguess=1
            D=Dldew
            call TPRHO (t,p,z,kph,kguess,D,ierr,herr2)
            Dl=D
            Dv=D
            if (ierr.ne.0) then
              ierr=213
              write (herr,1213) herr2(1:189),hnull
 1213         format ('[TPFLSH error 213] liquid density iteration ',
     &                'did not converge:  ',a189,a1)
              call ERRMSG (ierr,herr)
              q=999.d0     !quality undefined
              RETURN
            end if
            call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
c  compute quality based on volumes (possible for T,h to be
c  double-valued in compressed liquid)
            q=(1.0d0/D-1.0d0/Dldew)/(1.0d0/Dvdew-1.0d0/Dldew)
          else
c  mixture:  calculate bubble point to determine if liquid or 2-phase
            call SATT (t,z,1,pbub,Dlbub,Dvbub,xbub,ybub,ierr,herr2)
            if (ierr.ne.0 .and. t.gt.tc-5.0d0) then
              pbub=pc
              Dlbub=rhoc*2
              Dvbub=rhoc/2
            elseif (ierr.ne.0) then
              ierr=211
              write (herr,1211) herr2(1:190),hnull
 1211         format ('[TPFLSH error 211] bubble point calculation ',
     &                'did not converge:  ',a190,a1)
              call ERRMSG (ierr,herr)
              D=0.0d0
              Dl=0.0d0
              Dv=0.0d0
              q=999.d0     !quality undefined
              RETURN
            end if
            if (p.ge.pbub) then
c  mixture single-phase liquid
c             write (*,*) ' TPFLSH--mixture single-phase liquid'
              kph=1
              kguess=1
              D=Dlbub
              call TPRHO (t,p,z,kph,kguess,D,ierr,herr2)
              Dl=D
              Dv=D
              if (ierr.ne.0 .or. (d.lt.dlbub .and. d.gt.dvdew)) then
                ierr=213
                write (herr,1213) herr2(1:189),hnull
                call ERRMSG (ierr,herr)
                q=999.d0     !quality undefined
                RETURN
              end if
              call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
c  compute quality based on volumes (possible for T,h to be
c  double-valued in compressed liquid
              q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub)
            else
c
c  two-phase mixture
c             write (*,*) ' TPFLSH--two-phase mixture'
c
c  generate initial guesses for x and y by interpolating input pressure
c  with dew and bubble point pressures
c
              if (iflag.eq.1) then
c  dew point did not converge, location of state point is either two-phase
c  or vapor phase
                ierr=212
                write (herr,1212) herr2(1:193),hnull
                call ERRMSG (ierr,herr)
                D=0.0d0
                Dl=0.0d0
                Dv=0.0d0
                q=999.d0     !quality undefined
                RETURN
              endif
              if (ABS(1.0d0-pbub/pdew).lt.1.0d-6 .and. ianc.eq.0) then
c  possible azeotrope
                q=0.5d0
              else
                q=1.0d0-(p-pdew)/(pbub-pdew)
              end if
              Dl=Dlbub                 !initial guess for liquid density
              Dv=Dvdew                 !initial guess for vapor density
              if (ianc.eq.0) then
                xsum=0.0d0
                ysum=0.0d0       !sums for normalization of compositions
                do i=1,nc
                  x(i)=(1.0d0-q)*z(i)+q*xdew(i)
                  y(i)=q*z(i)+(1.0d0-q)*ybub(i)
                  xsum=xsum+x(i)
                  ysum=ysum+y(i)
                enddo
                do i=1,nc
                  x(i)=x(i)/xsum
                  y(i)=y(i)/ysum
                enddo
                call TPFL2 (t,p,z,Dl,Dv,x,y,q,ierr,herr)
                if (ierr.ne.0) then
c  two-phase iteration did not converge--error message written by TPFL2
                  q=999.d0     !quality undefined
                  RETURN
                end if
c  compute 2-phase properties and load output variables
c  call TPRHO to ensure consistency of t, p, rho
                kguess=1
                call TPRHO (t,p,x,1,kguess,Dl,ierr1,herr1)
                call TPRHO (t,p,y,2,kguess,Dv,ierr1,herr1)
                if (ierr1.ne.0) then
                  herr2='[TPFLSH error] (mix, 2-phase):  '//herr1
                  call ERRMSG (ierr1,herr2)
                  RETURN
                end if
              endif
              call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt)
              call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt)
              cp=xnotd                   !Cp, w not defined for 2-phase
              w=xnotd
              cv=xnotd
c  bulk properties are weighted average of liquid and vapor phases
              alpha=1.0d0-q               !alpha is liq fraction,
              D=1.0d0/(alpha/Dl+q/Dv)     !q is vapor frac
              e=alpha*el+q*ev
              h=alpha*hl+q*hv
              s=alpha*sl+q*sv
            end if
          end if
        end if
      end if
c
c  if limits check resulted in warning (as opposed to error) return
c  that message; do this again in case intermediate iteration did not
c  converge (thereby overwriting any warning message from LIMITX)
c
      if (ierrl.lt.0) then
        ierr=ierrl
        write (herr,2006) ierrl,herrl(1:234),hnull
 2006   format ('[TPFLSH warning',i3,'] ',a234,a1)
        call ERRMSG (ierr,herr)
      end if
c
      if (D.gt.rhomax) then
c  check that computed (output) density is within bounds of EOS
        ierr=2
        write (herr,2008) ierr,t,p/1000.0d0,hnull
 2008   format ('[TPFLSH error',i3,']; input T,p correspond to ',
     &          'a density above the limit of the EOS; T =',g11.5,
     &          ' K, p =',g11.5,' MPa.',a1)
        call ERRMSG (ierr,herr)
      end if
c
      RETURN
      end                                             !subroutine TPFLSH
c
c ======================================================================
c
      subroutine TPFL2 (t,p,z,Dl,Dv,x,y,q,ierr,herr)
c
c  flash calculation given temperature, pressure, and bulk composition
c
c  This routine accepts only two-phase states as input; if the phase is
c  not known use TPFLSH.  Use TPRHO for single-phase states.
c
c  inputs:
c        t--temperature [K]
c        p--pressure [kPa]
c        z--overall (bulk) composition [array of mol frac]
c       Dl--initial guess for molar density [mol/L] of the liquid phase
c       Dv--initial guess for molar density [mol/L] of the vapor phase
c        x--initial guess for composition of liquid phase
c           [array of mol frac]
c        y--initial guess for composition of vapor phase
c           [array of mol frac]
c        q--initial guess for vapor quality on a MOLAR basis
c           [moles vapor/total moles]
c
c  outputs:
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c     ierr--error flag:  (these are also passed up to TPFLSH for output)
c                        0 = successful
c                      216 = TPRHO did not converge for liquid
c                      217 = TPRHO did not converge for vapor
c                      218 = inner loop (liquid frac) did not converge
c                      219 = outer loop (composition) did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  12-29-95  MM, original version, extracted from TPFLSH
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  04-09-97  MM, fix bug:  converged value of q not returned
c  04-11-97  MM, fix divide by zero on saturation boundary
c                adjust tolr --> d-7; allow alpha outside 0,1 by d-5
c  07-15-97  MM, renumber and add detail to error messages
c                get flags for 'not defined' from common /FLAGS/
c  10-01-97  MM, add compiler switch to allow access by DLL
c  11-01-00 EWL, change tol from -7 to -9 so that TPFLSH finds a better root.
c  12-13-00 EWL, changed the calculation of y from x*xk to y*fliq/fvap.  Before
c                this, sum(y) was always equal to one on the first iteration
c                and the program exited prematurely
c  02-20-01 EWL, changed the tolerance after 15 iterations to increase the
c                chance of converging
c  02-26-01 EWL, if dpdrho is too large, liquid root is no good,
c                use SATP to get dl
c  08-09-01 EWL, made the data statement for tolr into a regular executable line
c  09-11-02 EWL, increase itmax to 80
c  09-11-02 EWL, at every 10th iteration, use average of old and new compositions
c                to get next composition.  Helps when x and y are jumping back and forth
c  07-22-03 EWL, modify how alpha(3) is calculated
c  08-23-06 EWL, changed the tolerance after 50 iterations
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: TPFL2
c     dll_export TPFL2
c
      parameter (ncmax=20)        !max number of components in mixture
      character*1 htab,hnull
      character*255 herr,herr2
      dimension z(ncmax),x(ncmax),y(ncmax),xs(ncmax),ys(ncmax)
      dimension xbub(ncmax),xdew(ncmax)
      dimension alpha(3),falpha(2),fliq(ncmax),fvap(ncmax),xk(ncmax)
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,xerr,xnotd,xnotc
      data itmax /80/
c
      tolr=1.0d-9
      ierr=0
      herr=' '
c
c  begin outer iteration loop for the vapor composition
c
      alpha(1)=1.0d0-q
      kguess=1
      j=1
      do 200 iy=1,itmax
      if (iy.eq.12) tolr=tolr*100.d0
      if (iy.eq.24) tolr=tolr*100.d0
      if (iy.eq.50) tolr=tolr*100.d0
c  compute densities and fugacities for each phase
c     write (*,1001) t,p,(x(i),i=1,5)
c1001 format (1x,' TPFL2 call TPRHO; liq t,p,x:  ',f10.5,f14.6,5f10.6)
      call TPRHO (t,p,x,1,kguess,Dl,ierr,herr2)
      call DPDD (t,Dl,x,dpdrho)
      if (ABS(dpdrho).gt.1.d7) then
        call SATP (p,x,1,tt,Dl,Dvdew,xbub,xdew,ierr2,herr2)
      endif
      if (ierr.ne.0) then
c  try modifying the composition to get a workable solution.
        x(1)=x(1)+.01d0
        do i=2,nc
          x(i)=x(i)-.01d0/DBLE(nc-1)
        enddo
        if (x(1).lt.1.d0) goto 200
        ierr=216
        write (herr,1216) herr2(1:170),hnull
 1216   format ('[TPFLSH error 216] liquid density iteration for ',
     &          '2-phase state did not converge:  ',a170,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
      call FGCTY (t,Dl,x,fliq)
c     write (*,1002) t,p,(y(i),i=1,5)
c1002 format (1x,' TPFL2 call TPRHO; vap t,p,x:  ',f10.5,f14.6,5f10.6)
      call TPRHO (t,p,y,2,kguess,Dv,ierr,herr2)
      if (ierr.ne.0) then
        do i=1,nc
          y(i)=0.8d0*y(i)+0.2d0*ys(i)
        enddo
        goto 200
c       ierr=217
c       write (herr,1217) herr2(1:170),hnull
c1217   format ('[TPFLSH error 217] vapor density iteration for ',
c    &          '2-phase state did not converge:  ',a170,a1)
c       call ERRMSG (ierr,herr)
c       RETURN
      end if
      call FGCTY (t,Dv,y,fvap)
c  compute k-factors (assume constant for a given guess of y)
      do i=1,nc
        if (z(i).gt.0.0d0) then
          xk(i)=fliq(i)*y(i)/(x(i)*fvap(i))
        else
c  zero composition for component i
          xk(i)=0.0d0
        end if
      enddo
c     write (*,1004) iy,(xk(i),i=1,nc)
c1004 format (1x,' TPFL2 ity,k-factors: ',i4,5e18.10)
c
c  begin inner iteration loop for the liquid fraction (alpha)
c
      j=1                      !iteration flag
c     lfapos=.false.           !flags for reguli-falsi iteration
c     lfaneg=.false.
      do i=1,nc
        xs(i)=x(i)
      enddo
      do ix=1,itmax
c       write (*,1040) ix,(x(i),i=1,nc)
c1040 format (1x,' TPFL2 itx,old liq comps:  ',i4,10x,5f16.12)
        xsum=0.0d0
        do i=1,nc
          if (z(i).gt.0.0d0) then
c  compute new guess for liquid composition
            x(i)=z(i)/(alpha(j)+xk(i)*(1.0d0-alpha(j)))
          else
c  zero composition for component i
            x(i)=0.0d0
          end if
          xsum=xsum+x(i)
        enddo
c  normalize liquid composition
        do i=1,nc
          x(i)=x(i)/xsum
        enddo
c       write (*,1050) ix,xsum,(x(i),i=1,nc)
c1050 format (1x,' TPFL2 itx,xsum,new comps: ',i4,6f16.12)
        falpha(j)=1.0d0-xsum
c       write (*,1056) ix,alpha(j),xsum,(x(i),i=1,3)
c1056 format (1x,' TPFL2:',i4,35x,2f10.6,3f10.6)
        if (ABS(falpha(j)).lt.1.0d-2*tolr) then
c  inner loop has converged
c  (tolerance on inner loop must be tighter than outer loop)
          goto 640
        else
c  store variables for possible use in reguli-falsi
c  reguli-falsi actually seems to slow convergence--comment out
c       if (falpha(j).lt.0.0d0) then
c         lfaneg=.true.
c         aneg=alpha(j)
c         faneg=falpha(j)
c       else
c         lfapos=.true.
c         apos=alpha(j)
c         fapos=falpha(j)
        end if
c  generate next guess for liquid fraction
        if (j.eq.1) then
c  for second iteration, move 5 % of way to alpha = 1.0
c         alpha(2)=alpha(1)+0.05d0*(1.0d0-alpha(1))
c  above could result in alpha(2) = alpha(1) if alpha(1) = 1
c  try moving towards alpha = 0.5 !MM 04-11-97
          j=2
          alpha(2)=alpha(1)+SIGN(0.05d0,0.5d0-alpha(1))
        else
c  for subsequent iterations, use secant method
          if (ABS(falpha(2)-falpha(1)).gt.1.0d-2*tolr) then
            alpha(3)=alpha(2)-falpha(2)*(alpha(2)-alpha(1))/
     &               (falpha(2)-falpha(1))
          else
            alpha(3)=alpha(1)+SIGN(tolr,0.5d0-alpha(1))
          end if
c  new guess for alpha must lie within bounds of zero and one
c  but allow small tolerance for slop in TPRHO and thus fugacities.
c         write (*,1060) ix,(alpha(i),i=1,3),(falpha(i),i=1,2)
c1060   format (1x,' TPFL2 itx,alpha,falpha:   ',i4,5f16.12)
c  alpha(3) generally goes out of bounds when the slope of falpha vs. x
c  is near flat.  Use next guess of half of alpha(2).
          if (alpha(3).gt.1.0d0+1.0d2*tolr) then
            alpha(3)=alpha(2)/2.d0
          else if (alpha(3).lt.0.0d0-1.0d2*tolr) then
            alpha(3)=alpha(2)/2.d0
          end if
c  check that new guess is not outside previous bounds
c         if (lfaneg .and. lfapos .and. (alpha(3).gt.MAX(fapos,faneg)
c    &       .or. alpha(3).lt.MIN(fapos,faneg))) then
c  if so, use reguli-falsi
c           alpha(3)=apos-fapos*(apos-aneg)/(fapos-faneg)
c         end if
c  discard oldest iteration
          alpha(1)=alpha(2)
          alpha(2)=alpha(3)
          falpha(1)=falpha(2)
        end if
      enddo
c  inner iteration loop has not converged, issue warning and proceed
      ierr=-218
      write (herr,1218) hnull
 1218 format ('[TPFLSH warning 218] inner iteration loop for liquid ',
     &        'fraction in 2-phase state did not converge:  ',a1)
      call ERRMSG (ierr,herr)
c
c  end of inner iteration loop for liquid fraction
c
 640  continue
      ysum=0.0d0
c  compute next guess for vapor composition, using x(i) from inner loop
      do i=1,nc
c       y(i)=x(i)*xk(i)                !Old way, but incorrect
        if (fliq(i).eq.xerr .or. fvap(i).eq.xerr) goto 210
        ys(i)=y(i)
        y(i)=y(i)*fliq(i)/fvap(i)
        ysum=ysum+y(i)
      enddo
c  normalize vapor compositions
      do i=1,nc
        y(i)=y(i)/ysum
      enddo
c     write (*,1076) iy,t,Dl,Dv,ysum,(y(i),i=1,3)
c1076 format (1x,' TPFL2:',i4,f7.2,2e14.6,10x,f10.6,3f10.6)
      if (ABS(1.0d0-ysum).lt.tolr) then
c  outer loop has converged (also check for warning from inner loop)
        if (ierr.ne.0) then
          ierr=218
          write (herr,2218) hnull
 2218     format ('[TPFLSH error 218] inner iteration loop for liquid ',
     &            'fraction in 2-phase state did not converge:  ',a1)
          call ERRMSG (ierr,herr)
        end if
        if (iy.gt.1) goto 840
      end if
      if (int(iy/10)*10.eq.iy) then
        do i=1,nc
          x(i)=(x(i)+xs(i))/2.d0
          y(i)=(y(i)+ys(i))/2.d0
        enddo
      endif
 200  continue
c  outer iteration loop has not converged, issue error and return
 210  ierr=219
      write (herr,1219) hnull
 1219 format ('[TPFLSH error 219] outer iteration loop for ',
     &        'composition in 2-phase state did not converge:  ',a1)
      call ERRMSG (ierr,herr)
c
c  end of outer iteration loop for vapor composition
c
 840  continue
      q=1.0d0-alpha(j)
c
      RETURN
      end                                              !subroutine TPFL2
c
c ======================================================================
c
      subroutine TDFLSH (t,D,z,p,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
c
c  flash calculation given temperature, bulk density, & bulk composition
c
c  This routine accepts both single-phase and two-phase states as input,
c  if the phase is known, then subroutine THERM is (much) faster.
c
c  inputs:
c        t--temperature [K]
c        D--overall (bulk) molar density [mol/L]
c        z--overall (bulk) composition [array of mol frac]
c
c  outputs:
c        p--pressure [kPa]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c           if only one phase is present, Dl = Dv = D
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c           if only one phase is present, 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 = 998 superheated vapor, but quality not defined (t > Tc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful
c                        1 = Tin < Tmin
c                        2 = Din > Dmax or Din < 0
c                        3 = T and D out of range
c                        8 = x out of range (component and/or sum < 0
c                            or > 1)
c                        9 = x and T out of range
c                       10 = x and D out of range
c                       11 = x and T and D out of range
c                      220 = CRITP did not converge
c                      221 = SATT did not converge at bubble point
c                      222 = SATT did not converge at dew point
c                      223 = SATT (bubble pt) did not converge for 2-ph
c                      224 = SATT (dew pt) did not converge for 2-phase
c                      225 = TPFL2 did not converge
c                      226 = 2-phase iteration did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  09-20-95  MM, original version
c  09-25-95  MM, rearrange argument list (outputs in order rho, x, q)
c  10-11-95  MM, RETURN if any error detected
c  11-08-95  MM, patch to allow testing of GUI with mixtures
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  12-13-95  MM, set undefined output compositions to zero
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  01-07-97  MM, error message bug (do not concatenate herr with itself)
c  01-21-97  MM, return error message for mixtures (not yet implemented)
c  04-09-97  MM, add mixture calculations; restructure cases
c                add limits checks
c  04-22-97  MM, use h, rather than s, to compute q (h converges at p = 0)
c  05-15-97  MM, use V to compute q for comp. liq. (T,h can be double-valued),
c                add q = 998 case for t > Tc, but p < Pc
c  07-15-97  MM, renumber and add detail to error messages
c                get flags for 'not defined' from common /FLAGS/
c  10-01-97  MM, add compiler switch to allow access by DLL
c  10-16-97  MM, fix bug--out of range errors not caught properly
c  12-01-97  MM, error code not passed correctly if inputs out of range
c  02-11-98  MM, check that computed pressure is within bounds
c  03-05-99 EWL, remove check on pressure, and call LIMITX to check for
c                p>pmax or p>melt
c  12-16-02 EWL, add checks for pseudo-pure fluid, calculate p using q
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: TDFLSH
c     dll_export TDFLSH
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport E
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      character*1 htab,hnull
      character*255 herr,herr1,herr2,herrl
      dimension x(ncmax),y(ncmax),z(ncmax)
      dimension xbub(ncmax),xdew(ncmax),ybub(ncmax),ydew(ncmax)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
      common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
c  flags to GUI indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      common /FLAGS2/ iamwat,ianc,iwat
c
      icomp=0
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
      endif
      ierr=0
      herr=' '
c
c  set output liquid and vapor compositions to input values
c  zero output values for undefined components
      if (icomp.eq.0) then
        do i=1,nc
          x(i)=z(i)
          y(i)=z(i)
        enddo
      endif
c     if (nc.lt.ncmax) then
c       do i=nc+1,ncmax
c         x(i)=0.0d0
c         y(i)=0.0d0
c       enddo
c     end if
      if (d.lt.1.0d-14) then
c  ideal-gas
        p=D*R*t
        Dl=D
        Dv=D
        call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
        RETURN
      endif
c
c  check that input conditions are within limits
c
      pdum=0.0d0
      call LIMITX ('EOS',t,D,pdum,z,tmin,tmax,rhomax,pmax,ierrl,herrl)
      if (ierrl.lt.0) then
c  one or inputs are outside limits--if just a warning proceed w/ calc
        write (herr,1006) ierrl,herrl(1:234),hnull
 1006   format ('[TDFLSH warning',i3,'] ',a234,a1)
        call ERRMSG (ierrl,herr)
c  if error (as opposed to warning) set output density to zero and return
      else if (ierrl.gt.0) then
        write (herr,1007) ierrl,herrl(1:236),hnull
 1007   format ('[TDFLSH error',i3,'] ',a236,a1)
        call ERRMSG (ierrl,herr)
        q=999.d0     !quality undefined
        ierr=ierrl
        RETURN
      end if
c
      call CRITP (z,tc,pc,rhoc,ierr,herr2)
      if (ierr.ne.0) then
        ierr=220
        write (herr,1008) herr2(1:235),hnull
 1008   format ('[TDFLSH error 220] ',a235,a1)
        call ERRMSG (ierr,herr)
        q=999.d0     !quality undefined
        RETURN
      end if
c
      if (icomp.ne.0 .and. t.lt.ttp(icomp)) then
c  temperature less than triple point temperature, but in the gas phase.
        call THERM (t,D,z,p,e,h,s,cv,cp,w,hjt)
        Dl=D
        Dv=D
      elseif (t.ge.tc .or. D.lt.1.0d-10) then
c  super-critical state or rho = 0 (x = y = z as set above)
        call THERM (t,D,z,p,e,h,s,cv,cp,w,hjt)
        Dl=D
        Dv=D
        if (p.le.pc) then
          q=998.d0     !superheated, but cannot define quality
        else
          q=999.d0     !quality undefined
        end if
      else
c  sub-critical state--call saturation routine to determine phase
        call SATT (t,z,2,pdew,Dldew,Dvdew,xdew,ydew,ierr,herr1)
        if (ierr.ne.0) then
          ierr=222
          write (herr,1222) herr1(1:193),hnull
 1222     format ('[TDFLSH error 222] dew point calculation ',
     &            'did not converge:  ',a193,a1)
          call ERRMSG (ierr,herr)
          D=0.0d0
          Dl=0.0d0
          Dv=0.0d0
          q=999.d0     !quality undefined
          RETURN
        end if
        if (icomp.ne.0 .and. ianc.eq.0) then
c  if pure, bubble point density provided by call to SATT at dew point
          Dlbub=Dldew
          Dvbub=Dvdew
          pbub=pdew
        else
c  if mixture, call saturation routine again at bubble point
          call SATT (t,z,1,pbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1)
          if (ierr.ne.0) then
            ierr=221
            write (herr,1221) herr1(1:190),hnull
 1221       format ('[TPFLSH error 221] bubble point calculation ',
     &              'did not converge:  ',a190,a1)
            call ERRMSG (ierr,herr)
            D=0.0d0
            Dl=0.0d0
            Dv=0.0d0
            q=999.d0     !quality undefined
            RETURN
          end if
        end if
        if (D.le.Dvdew .or. D.ge.Dlbub) then
c  single-phase (liq or vapor) (pure or mixture) (x = y = z as set above)
          call THERM (t,D,z,p,e,h,s,cv,cp,w,hjt)
          Dl=D
          Dv=D
          call ENTHAL (t,Dvdew,z,hvdew)
c         call ENTHAL (t,Dlbub,z,hlbub)
c         q=(h-hlbub)/(hvdew-hlbub)
c  compute quality based on volumes (possible for T,h to be
c  double-valued in compressed liquid)
          if (abs(dvdew-dlbub).gt.1.d-12)
     &        q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub)
        else
          if (icomp.ne.0) then
c  special case:  pure-fluid two-phase (x = y = z as set above)
            if (abs(dvdew-dlbub).gt.1.d-12)
     &          q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub)
            p=q*pdew+(1.0d0-q)*pbub
            Dl=Dlbub
            Dv=Dvdew
          else
c  general case:  mixture 2-phase
            ksat=1   !bubble and dew point data provided to TDFL2
            call TDFL2 (t,D,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  p,Dl,Dv,x,y,q,ierr,herr)
            if (ierr.ne.0) then
c  two-phase iteration did not converge--error message written by TDFL2
              q=999.d0     !quality undefined
              RETURN
            end if
c           write (*,1225) p,q,(x(i),i=1,3),(y(i),i=1,3)
c1225       format (1x,' TDFLSH--TDFL2 return   p,q,x,y = ',8f12.7)
          end if
c  compute remaining properties for 2-phase states
          call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt)
          call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt)
          if (nc.ge.2) p=ptherm
          e=q*ev+(1.0d0-q)*el
          h=q*hv+(1.0d0-q)*hl
          s=q*sv+(1.0d0-q)*sl
          w=xnotd     !Cp,w not defined for 2-phase states
          cp=xnotd
          cv=xnotd
        end if
      end if
c
      call LIMITX ('EOS',t,D,p,z,tmin,tmax,rhomax,pmax,ierr,herrl)
      if (ierr.lt.0) then
c  one or inputs are outside limits--if just a warning proceed w/ calc
        write (herr,1006) ierr,herrl(1:234),hnull
        call ERRMSG (ierr,herr)
c  if error (as opposed to warning) set output density to zero and return
      else if (ierr.gt.0) then
        write (herr,1007) ierr,herrl(1:236),hnull
        call ERRMSG (ierr,herr)
        q=999.d0     !quality undefined
      end if
c
      RETURN
      end                                             !subroutine TDFLSH
c
c ======================================================================
c
      subroutine TDFL2 (t,D,z,ksat,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  p,Dl,Dv,x,y,q,ierr,herr)
c
c  flash calculation given temperature, bulk density, and composition
c
c  This routine accepts only two-phase states as input; it is intended
c  primarily for use by the general temperature-density flash routine
c  TDFLSH.  It may be called independently if the state is known to be
c  two-phase.  But beware--this routine does not check limits, and it
c  will be significantly faster than TDFLSH only if the bubble and dew
c  point limits can be provided (ksat = 1 option).
c
c  inputs:
c        t--temperature [K]
c        D--overall (bulk) molar density [mol/L]
c        z--overall (bulk) composition [array of mol frac]
c     ksat--flag for bubble and dew point limits
c           0 = dew and bubble point limits computed here
c           1 = must provide values for the following:
c     pbub--bubble point pressure [kPa] at (t,x=z)
c     pdew--dew point pressure [kPa] at (t,y=z)
c    Dlbub--liquid density [mol/L] at bubble point
c    Dvdew--vapor density [mol/L] at dew point
c     ybub--vapor composition [array of mol frac] at bubble point
c     xdew--liquid composition [array of mol frac] at dew point
c
c  outputs:
c        p--pressure [kPa]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: TDFL2
c     dll_export TDFL2
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension z(ncmax),x(ncmax),y(ncmax),xdew(ncmax),ybub(ncmax)
c
      call ABFL2 (t,d,z,ksat,0,'TD',
     &                 tbub,tdew,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                 tt,p,Dl,Dv,x,y,q,ierr,herr)
c
      RETURN
      end                                              !subroutine TDFL2
c
c ======================================================================
c
      subroutine PDFLSH (p,D,z,t,Dl,Dv,x,y,q,e,h,s,cv,cp,w,ierr,herr)
c
c  flash calculation given density, pressure, and bulk composition
c
c  This routine accepts both single-phase and two-phase states as the
c  input; for single-phase calculations, the subroutine PDFL1 is faster.
c
c  inputs:
c        D--overall (bulk) molar density [mol/L]
c        p--pressure [kPa]
c        z--overall (bulk) composition [array of mol frac]
c
c  outputs:
c        t--temperature [K]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c           if only one phase is present, Dl = Dv = D
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c           if only one phase is present, x = y = z
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c           q < 0 indicates subcooled (compressed) liquid
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q > 1 indicates superheated vapor
c           q = 998 superheated vapor, but quality not defined (t > Tc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful
c                      210 = CRITP did not converge
c                      211 = SATT did not converge at bubble point
c                      212 = SATT did not converge at dew point
c                      213 = TPRHO did not converge for liquid state
c                      214 = TPRHO did not converge for vapor state
c                      215 = TPRHO did not convg for supercritical state
c                      216 = TPRHO did not convg for liq in 2-phase it
c                      217 = TPRHO did not convg for vap in 2-phase it
c                      218 = liquid frac for 2-phase it did not converge
c                      219 = composition for 2-phase it did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado
c  03-03-99 EWL, original version
c  03-14-00 EWL, change name from DPFLSH to PDFLSH
c  05-25-00 EWL, set t=-1 in initial call to LIMITX
c  12-16-02 EWL, add checks for pseudo-pure fluid, calculate t using q
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PDFLSH
c     dll_export PDFLSH
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport E
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      character*1 htab,hnull
      character*255 herr,herr2,herrl
      dimension z(ncmax),x(ncmax),y(ncmax)
      dimension xdew(ncmax),ydew(ncmax),ybub(ncmax),xbub(ncmax)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
c  flags to GUI indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      common /FLAGS2/ iamwat,ianc,iwat
      common /CCON/ tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx),
     &              ttp(n0:nx),ptp(n0:nx),dtp(n0:nx),dtpv(n0:nx),
     &              tnbp(n0:nx),dnbpl(n0:nx),dnbpv(n0:nx),
     &              wm(n0:nx),accen(n0:nx),dipole(n0:nx),Reos(n0:nx)
      common /prnterr/ iprnterr
c
      icomp=0
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
      endif
      ierr=0
      herr=' '
      t=0.0d0
      Dl=0.0d0
      Dv=0.0d0
      e=0.0d0
      h=0.0d0
      s=0.0d0
      cv=0.0d0
      cp=0.0d0
      w=0.0d0
      q=999.0d0
c
c  initialize output liquid and vapor compositions to input values
c  zero output values for undefined components
      do i=1,nc
        x(i)=z(i)
        y(i)=z(i)
      enddo
c     if (nc.lt.ncmax) then
c       do i=nc+1,ncmax
c         x(i)=0.0d0
c         y(i)=0.0d0
c       enddo
c     end if

      if (p.lt.1.0d-14 .or. d.lt.1.0d-14) then
c  ideal-gas
        if (D.ne.0.d0) t=p/(R*D)
        RETURN
      end if
c
c  check that input conditions are within limits
c
      call CRITP (z,tc,pc,rhoc,ierr,herr2)
      if (ierr.ne.0) then
        ierr=210
        write (herr,1008) herr2(1:235),hnull
 1008   format ('[PDFLSH error 210] ',a235,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
      t=-1
      call LIMITX ('EOS',t,D,p,z,tmin,tmax,rhomax,pmax,ierrl,herrl)
      if (ierrl.lt.0) then
c  one or inputs are outside limits--if just a warning proceed w/ calc
        write (herr,1006) ierrl,herrl(1:234),hnull
 1006   format ('[PDFLSH warning',i3,'] ',a234,a1)
        call ERRMSG (ierrl,herr)
c  if error (as opposed to warning) set output density to zero and return
      else if (ierrl.gt.0) then
        write (herr,1007) ierrl,herrl(1:236),hnull
 1007   format ('[PDFLSH error',i3,'] ',a236,a1)
        call ERRMSG (ierrl,herr)
        ierr=ierrl
        RETURN
      end if
c
      if (p.ge.pc) then
c  supercritical state (x = y = z as set above)
c       write (*,*) ' PDFLSH--supercritical'
        call PDFL1 (p,D,z,t,ierr,herr2)
        if (ierr.ne.0) then
          ierr=215
          write (herr,1215) herr2(1:183),hnull
 1215     format ('[PDFLSH error 215] supercritical density iteration ',
     &            'did not converge:  ',a183,a1)
          call ERRMSG (ierr,herr)
          RETURN
        end if
        Dl=D
        Dv=D
        call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
      else
c  subcritical state--call saturation routine to determine liq or vap
        i=iprnterr
        iprnterr=0
        ierr=0
        Dvdew=D
        Dldew=D
        tdew=0
c  do not call SATP if P<triple point pressure for pure fluids
        if (icomp.eq.0 .or. p.gt.ptp(icomp) .or. ianc.ne.0) then
          call SATP (p,z,2,tdew,Dldew,Dvdew,xdew,ydew,ierr,herr2)
        endif
        iprnterr=i
        iflag=0
        if (ierr.ne.0 .and. icomp.ne.0) then
          ierr=212
          write (herr,1212) herr2(1:193),hnull
 1212     format ('[PDFLSH error 212] dew point calculation ',
     &            'did not converge:  ',a193,a1)
          call ERRMSG (ierr,herr)
          RETURN
        elseif (ierr.ne.0) then
c  do not exit yet, check if valid liquid state
          iflag=1
        end if
        if (D.le.Dvdew+1.0d-12) then
c  single-phase vapor (x = y = z as set above)
c         write (*,*) ' PDFLSH--single-phase vapor'
          call PDFL1 (p,D,z,t,ierr,herr2)
          if (ierr.ne.0) then
            ierr=214
            write (herr,1214) herr2(1:190),hnull
 1214       format ('[PDFLSH error 214] vapor density iteration ',
     &              'did not converge:  ',a190,a1)
            call ERRMSG (ierr,herr)
            t=0.0d0
            RETURN
          end if
          Dl=D
          Dv=D
          call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
          q=998.d0
          if (tdew.gt.0 .and. ianc.eq.0) then
            call ENTHAL (tdew,Dldew,xdew,hldew)
            call ENTHAL (tdew,Dvdew,z,hvdew)
            q=(h-hldew)/(hvdew-hldew)
          endif
        else
c  mixture:  calculate bubble point to determine if liquid or 2-phase
          if (icomp.ne.0 .and. ianc.eq.0) then
            Dlbub=Dldew
            Dvbub=Dvdew
            tbub=tdew
          else
            call SATP (p,z,1,tbub,Dlbub,Dvbub,xbub,ybub,ierr,herr2)
          endif
          if (ierr.ne.0) then
            ierr=211
            write (herr,1211) herr2(1:190),hnull
 1211       format ('[PDFLSH error 211] bubble point calculation ',
     &              'did not converge:  ',a190,a1)
            call ERRMSG (ierr,herr)
            RETURN
          end if
          if (D.gt.Dlbub) then
c  mixture single-phase liquid
c           write (*,*) ' TPFLSH--mixture single-phase liquid'
            call PDFL1 (p,D,z,t,ierr,herr2)
            if (ierr.ne.0) then
              ierr=213
              write (herr,1213) herr2(1:189),hnull
 1213         format ('[PDFLSH error 213] liquid density iteration ',
     &                'did not converge:  ',a189,a1)
              call ERRMSG (ierr,herr)
              t=0.0d0
              RETURN
            end if
            Dl=D
            Dv=D
            call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
c  compute quality based on volumes
            q=(1.0d0/D-1.0d0/Dlbub)/(1.0d0/Dvdew-1.0d0/Dlbub)
c           write (*,*) ' PDFLSH--comp liq q by volumes:  ',q
          else
c
c  two-phase pure fluid
            if (icomp.ne.0) then
              Dl=Dlbub
              Dv=Dvdew
              q=(1.0d0/D-1.0d0/Dl)/(1.0d0/Dv-1.0d0/Dl)
              t=(1.d0-q)*tbub+q*tdew
              if (ianc.eq.1) call PTANC (t,p,q,d,'D',Dl,Dv)
c  two-phase mixture
c             write (*,*) ' PDFLSH--two-phase mixture'
c
c  generate initial guesses for x and y by interpolating input pressure
c  with dew and bubble point pressures
c
            else
              if (iflag.eq.1) then
c  dew point did not converge, location of state point is either two-phase
c  or vapor phase
                ierr=212
                write (herr,1212) herr2(1:193),hnull
                call ERRMSG (ierr,herr)
                t=0.0d0
                RETURN
              endif
              ksat=1   !bubble and dew point data provided to PDFL2
              call PDFL2 (p,D,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                    t,Dl,Dv,x,y,q,ierr,herr)
              if (ierr.ne.0) then
c  two-phase iteration did not converge--error message written by PDFL2
                q=999.d0     !quality undefined
                RETURN
              end if
c             write (*,1225) p,q,(x(i),i=1,3),(y(i),i=1,3)
c1225         format (1x,' TDFLSH--TDFL2 return   p,q,x,y = ',8f12.7)
            endif
            call THERM (t,Dl,x,ptherm,el,hl,sl,cv,cp,w,hjt)
            call THERM (t,Dv,y,ptherm,ev,hv,sv,cv,cp,w,hjt)
            cp=xnotd                   !Cp, w not defined for 2-phase
            w=xnotd
            cv=xnotd
c  bulk properties are weighted average of liquid and vapor phases
            alpha=1.0d0-q               !alpha is liq fraction,
            e=alpha*el+q*ev
            h=alpha*hl+q*hv
            s=alpha*sl+q*sv
          end if
        end if
      end if
c
c  if limits check resulted in warning (as opposed to error) return
c  that message; do this again in case intermediate iteration did not
c  converge (thereby overwriting any warning message from LIMITX)
c
      call LIMITX ('EOS',t,D,p,z,tmin,tmax,rhomax,pmax,ierrl,herrl)
      if (ierrl.lt.0) then
c  one or inputs are outside limits--if just a warning proceed w/ calc
        ierr=ierrl
        write (herr,2006) ierrl,herrl(1:234),hnull
 2006   format ('[PDFLSH warning',i3,'] ',a234,a1)
        call ERRMSG (ierr,herr)
c  if error (as opposed to warning) set output density to zero and return
      else if (ierrl.gt.0) then
        write (herr,1007) ierrl,herrl(1:236),hnull
        ierr=ierrl
        call ERRMSG (ierrl,herr)
        RETURN
      end if
c
      RETURN
      end                                             !subroutine PDFLSH
c
c ======================================================================
c
      subroutine PDFL1 (p,rho,x,t,ierr,herr)
c
c  iterate for single-phase temperature as a function of density, pressure,
c  and composition
c
c  inputs:
c      rho--molar density [mol/L]
c        p--pressure [kPa]
c        x--composition [array of mol frac]
c
c  outputs:
c        t--temperature [K]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PDFL1
c     dll_export PDFL1
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension x(ncmax)
      call ABFL1 (rho,p,x,0,'DP',0.d0,0.d0,t,pp,dd,ierr,herr)
      RETURN
      end                                              !subroutine PDFL1
c
c ======================================================================
c
      subroutine PDFL2 (p,d,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,Dl,Dv,x,y,q,ierr,herr)
c
c  flash calculation given pressure, bulk density, and composition
c
c  This routine accepts only two-phase states as input; it is intended
c  primarily for use by the general pressure-density flash routine
c  PDFLSH.  It may be called independently if the state is known to be
c  two-phase.  But beware--this routine does not check limits, and it
c  will be significantly faster than PDFLSH only if the bubble and dew
c  point limits can be provided (ksat = 1 option).
c
c  inputs:
c        p--pressure [kPa]
c        d--overall (bulk) density [mol/L]
c        z--overall (bulk) composition [array of mol frac]
c     ksat--flag for bubble and dew point limits
c           0 = dew and bubble point limits computed here
c           1 = must provide values for the following:
c     tbub--bubble point temperature [K] at (p,x=z)
c     tdew--dew point temperature [K] at (p,y=z)
c    Dlbub--liquid density [mol/L] at bubble point
c    Dvdew--vapor density [mol/L] at dew point
c     ybub--vapor composition [array of mol frac] at bubble point
c     xdew--liquid composition [array of mol frac] at dew point
c
c  outputs:
c        t--temperature [K]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PDFL2
c     dll_export PDFL2
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension z(ncmax),x(ncmax),y(ncmax),xdew(ncmax),ybub(ncmax)
c
      call ABFL2 (p,d,z,ksat,0,'PD',
     &                  tbub,tdew,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,pp,Dl,Dv,x,y,q,ierr,herr)
c
      RETURN
      end                                              !subroutine PDFL2
c
c ======================================================================
c
      subroutine PHFLSH (p,h,z,t,D,Dl,Dv,x,y,q,e,s,cv,cp,w,ierr,herr)
c
c  flash calculation given pressure, bulk enthalpy, and bulk composition
c
c  inputs:
c        p--pressure [kPa]
c        h--overall (bulk) enthalpy [J/mol]
c        z--composition [array of mol frac]
c
c  outputs:
c        t--temperature [K]
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 of mol frac] for liquid phase
c        y--composition [array of mol frac] for vapor phase
c           if only one phase is present, x = y = z
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c           q < 0 indicates subcooled (compressed) liquid
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q > 1 indicates superheated vapor
c           q = -998 subcooled liquid, but quality not defined (p > Pc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c        e--overall (bulk) internal energy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful (see PBFLSH for others)
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  09-20-95  MM, original version
c  09-25-95  MM, rearrange argument list (outputs in order rho, x, q)
c  10-11-95  MM, correct: t not returned for two-phase pure-component
c  10-11-95  MM, RETURN if any error detected
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  12-13-95  MM, set undefined output compositions to zero
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  01-07-97  MM, error message bug (do not concatenate herr with itself)
c  04-30-97  MM, restructure logic to parallel TDFLSH
c  05-15-97  MM, add q = -998 case for p > Pc, but t < Tc
c  05-21-97  MM, call ENTHAL at rho = 1d-10 to define hmax
c  07-15-97  MM, renumber and add detail to error messages
c                get flags for 'not defined' from common /FLAGS/
c  10-01-97  MM, add compiler switch to allow access by DLL
c  10-16-97  MM, evaluate hmin at p = 100 kPa, rather than p = pmax
c  02-11-98  MM, check that computed density is within bounds
c  05-25-00 EWL, call limitx with tdum=-1 to avoid check on tmelt
c  05-26-00 EWL, check for p<ptrp after calling SATP
c  06-06-00  MM, input h may have been changed by call to THERM
c  08-20-00 EWL, change initial guess from tc to tc+10
c  01-18-01 EWL, remove code and call generic PBFLSH
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PHFLSH
c     dll_export PHFLSH
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension x(ncmax),y(ncmax),z(ncmax)
      call PBFLSH (p,h,z,'PH',t,D,Dl,Dv,x,y,q,e,hh,s,cv,cp,w,ierr,herr)
      RETURN
      end                                             !subroutine PHFLSH
c
c ======================================================================
c
      subroutine PHFL1 (p,h,x,kph,t,D,ierr,herr)
c
c  flash calculation given pressure, enthalpy, and composition.  This routine
c  accepts only single-phase inputs, it is intended primarily for use with
c  the more general flash routine PHFLSH.
c
c  inputs:
c        p--pressure [kPa]
c        h--enthalpy [J/mol]
c        x--composition [array of mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c        t--initial guess for temperature [K]
c        D--initial guess for molar density [mol/L]
c
c  outputs:
c        t--temperature [K]
c        D--molar density [mol/L]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PHFL1
c     dll_export PHFL1
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension x(ncmax)
      call ABFL1 (p,h,x,kph,'PH',0.d0,0.d0,t,pp,D,ierr,herr)
      RETURN
      end                                              !subroutine PHFL1
c
c ======================================================================
c
      subroutine PHFL2 (p,h,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,Dl,Dv,x,y,q,ierr,herr)
c
c  flash calculation given pressure, bulk enthalpy, and composition
c
c  This routine accepts only two-phase states as input; it is intended
c  primarily for use by the general pressure-enthalpy flash routine
c  PHFLSH.  It may be called independently if the state is known to be
c  two-phase.  But beware--this routine does not check limits, and it
c  will be significantly faster than PHFLSH only if the bubble and dew
c  point limits can be provided (ksat = 1 option).
c
c  inputs:
c        p--pressure [kPa]
c        h--overall (bulk) molar enthalpy [J/mol]
c        z--overall (bulk) composition [array of mol frac]
c     ksat--flag for bubble and dew point limits
c           0 = dew and bubble point limits computed here
c           1 = must provide values for the following:
c     tbub--bubble point temperature [K] at (p,x=z)
c     tdew--dew point temperature [K] at (p,y=z)
c    Dlbub--liquid density [mol/L] at bubble point
c    Dvdew--vapor density [mol/L] at dew point
c     ybub--vapor composition [array of mol frac] at bubble point
c     xdew--liquid composition [array of mol frac] at dew point
c
c  outputs:
c        t--temperature [K]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PHFL2
c     dll_export PHFL2
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension z(ncmax),x(ncmax),y(ncmax),xdew(ncmax),ybub(ncmax)
c
      call ABFL2 (p,h,z,ksat,0,'PH',
     &                  tbub,tdew,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,pp,Dl,Dv,x,y,q,ierr,herr)
c
      RETURN
      end                                              !subroutine PHFL2
c
c ======================================================================
c
      subroutine PSFLSH (p,s,z,t,D,Dl,Dv,x,y,q,e,h,cv,cp,w,ierr,herr)
c
c  flash calculation given pressure, bulk entropy, and bulk composition
c
c  inputs:
c        p--pressure [kPa]
c        s--overall (bulk) entropy [J/mol-K]
c        z--composition array (mol frac)
c
c  outputs:
c        t--temperature [K]
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 of mol frac] for liquid phase
c        y--composition [array of 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 = -998 subcooled liquid, but quality not defined (p > Pc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful (see PBFLSH for others)
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  09-22-95  MM, original version
c  09-25-95  MM, rearrange argument list (outputs in order rho, x, q)
c  10-11-95  MM, RETURN if any error detected
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  12-13-95  MM, set undefined output compositions to zero
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  01-07-97  MM, error message bug (do not concatenate herr with itself)
c  05-01-97  MM, restructure logic to parallel TDFLSH, PHFLSH
c  05-15-97  MM, add q = -998 case for p > Pc, but t < Tc
c  07-15-97  MM, renumber and add detail to error messages
c                get flags for 'not defined' from common /FLAGS/
c  10-01-97  MM, add compiler switch to allow access by DLL
c  10-16-97  MM, evaluate smin at p = 100 kPa, rather than p = pmax
c  02-11-98  MM, check that computed density is within bounds
c  03-31-00 EWL, call limitx at very end to ensure that p<pmelt and d<dmax
c  03-31-00 EWL, if routine for p>pc fails, try another initial t before quitting
c  05-26-00 EWL, check for p<ptrp after calling SATP
c  06-06-00  MM, inputs may have been changed by call to THERM
c  01-18-01 EWL, remove code and call generic PBFLSH
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PSFLSH
c     dll_export PSFLSH
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension x(ncmax),y(ncmax),z(ncmax)
      call PBFLSH (p,s,z,'PS',t,D,Dl,Dv,x,y,q,e,h,ss,cv,cp,w,ierr,herr)
      RETURN
      end                                             !subroutine PSFLSH
c
c ======================================================================
c
      subroutine PSFL1 (p,s,x,kph,t,D,ierr,herr)
c
c  flash calculation given pressure, entropy, and composition.  This routine
c  accepts only single-phase inputs, it is intended primarily for use with
c  the more general flash routine PSFLSH.
c
c  inputs:
c        p--pressure [kPa]
c        s--entropy [J/mol-K]
c        x--composition [array of mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c        t--initial guess for temperature [K]
c        D--initial guess for molar density [mol/L]
c
c  outputs:
c        t--temperature [K]
c        D--molar density [mol/L]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PSFL1
c     dll_export PSFL1
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension x(ncmax)
      call ABFL1 (p,s,x,kph,'PS',0.d0,0.d0,t,pp,D,ierr,herr)
      RETURN
      end                                              !subroutine PSFL1
c
c ======================================================================
c
      subroutine PSFL2 (p,s,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,Dl,Dv,x,y,q,ierr,herr)
c
c  flash calculation given pressure, bulk entropy, and composition
c
c  This routine accepts only two-phase states as input; it is intended
c  primarily for use by the general pressure-entropy flash routine
c  PSFLSH.  It may be called independently if the state is known to be
c  two-phase.  But beware--this routine does not check limits, and it
c  will be significantly faster than PSFLSH only if the bubble and dew
c  point limits can be provided (ksat = 1 option).
c
c  inputs:
c        p--pressure [kPa]
c        s--overall (bulk) molar entropy [J/mol-K]
c        z--overall (bulk) composition [array of mol frac]
c     ksat--flag for bubble and dew point limits
c           0 = dew and bubble point limits computed here
c           1 = must provide values for the following:
c     tbub--bubble point temperature [K] at (p,x=z)
c     tdew--dew point temperature [K] at (p,y=z)
c    Dlbub--liquid density [mol/L] at bubble point
c    Dvdew--vapor density [mol/L] at dew point
c     ybub--vapor composition [array of mol frac] at bubble point
c     xdew--liquid composition [array of mol frac] at dew point
c
c  outputs:
c        t--temperature [K]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PSFL2
c     dll_export PSFL2
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension z(ncmax),x(ncmax),y(ncmax),xdew(ncmax),ybub(ncmax)
c
      call ABFL2 (p,s,z,ksat,0,'PS',
     &                 tbub,tdew,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                 t,pp,Dl,Dv,x,y,q,ierr,herr)
c
      RETURN
      end                                              !subroutine PSFL2
c
c ======================================================================
c
      subroutine PEFLSH (p,e,z,t,D,Dl,Dv,x,y,q,h,s,cv,cp,w,ierr,herr)
c
c  flash calculation given pressure, bulk energy, and bulk composition
c
c  inputs:
c        p--pressure [kPa]
c        e--overall (bulk) internal energy [J/mol]
c        z--composition [array of mol frac]
c
c  outputs:
c        t--temperature [K]
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 of mol frac] for liquid phase
c        y--composition [array of mol frac] for vapor phase
c           if only one phase is present, x = y = z
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c           q < 0 indicates subcooled (compressed) liquid
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q > 1 indicates superheated vapor
c           q = -998 subcooled liquid, but quality not defined (p > Pc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful (see PBFLSH)
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by E.W. Lemmon, NIST Thermophysics Division, Boulder, Colorado
c  01-18-01 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PEFLSH
c     dll_export PEFLSH
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension x(ncmax),y(ncmax),z(ncmax)
      call PBFLSH (p,e,z,'PE',t,D,Dl,Dv,x,y,q,ee,h,s,cv,cp,w,ierr,herr)
      RETURN
      end                                             !subroutine PEFLSH
c
c ======================================================================
c
      subroutine PEFL1 (p,e,x,kph,t,D,ierr,herr)
c
c  flash calculation given pressure, energy, and composition.  This routine
c  accepts only single-phase inputs, it is intended primarily for use with
c  the more general flash routine PEFLSH.
c
c  inputs:
c        p--pressure [kPa]
c        e--energy [J/mol]
c        x--composition [array of mol frac]
c      kph--phase flag:  1 = liquid
c                        2 = vapor
c        t--initial guess for temperature [K]
c        D--initial guess for molar density [mol/L]
c
c  outputs:
c        t--temperature [K]
c        D--molar density [mol/L]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PEFL1
c     dll_export PEFL1
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension x(ncmax)
      call ABFL1 (p,e,x,kph,'PE',0.d0,0.d0,t,pp,D,ierr,herr)
      RETURN
      end                                              !subroutine PEFL1
c
c ======================================================================
c
      subroutine PEFL2 (p,e,z,ksat,tbub,tdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,Dl,Dv,x,y,q,ierr,herr)
c
c  flash calculation given pressure, bulk energy, and composition
c
c  This routine accepts only two-phase states as input; it is intended
c  primarily for use by the general pressure-energy flash routine
c  PEFLSH.  It may be called independently if the state is known to be
c  two-phase.  But beware--this routine does not check limits, and it
c  will be significantly faster than PEFLSH only if the bubble and dew
c  point limits can be provided (ksat = 1 option).
c
c  inputs:
c        p--pressure [kPa]
c        e--overall (bulk) molar energy [J/mol]
c        z--overall (bulk) composition [array of mol frac]
c     ksat--flag for bubble and dew point limits
c           0 = dew and bubble point limits computed here
c           1 = must provide values for the following:
c     tbub--bubble point temperature [K] at (p,x=z)
c     tdew--dew point temperature [K] at (p,y=z)
c    Dlbub--liquid density [mol/L] at bubble point
c    Dvdew--vapor density [mol/L] at dew point
c     ybub--vapor composition [array of mol frac] at bubble point
c     xdew--liquid composition [array of mol frac] at dew point
c
c  outputs:
c        t--temperature [K]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PEFL2
c     dll_export PEFL2
c
      parameter (ncmax=20)        !max number of components in mixture
      character*255 herr
      dimension z(ncmax),x(ncmax),y(ncmax),xdew(ncmax),ybub(ncmax)
c
      call ABFL2 (p,e,z,ksat,0,'PE',
     &                  tbub,tdew,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,pp,Dl,Dv,x,y,q,ierr,herr)
c
      RETURN
      end                                              !subroutine PEFL2
c
c ======================================================================
c
      subroutine PBFLSH (p,b,z,ab,t,D,Dl,Dv,x,y,q,e,h,s,cv,cp,w,
     &                   ierr,herr)
c
c  flash calculation given pressure, bulk composition, and either
c  enthalpy, entropy, or energy
c
c  inputs:
c        p--pressure [kPa]
c        b--second property (energy, enthalpy, or entropy)
c        z--composition array (mol frac)
c       ab--character*2 string defining the inputs: 'PH', 'PS', or 'PE'
c
c  outputs:
c        t--temperature [K]
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 of mol frac] for liquid phase
c        y--composition [array of mol frac] for vapor phase
c           if only one phase is present, x = y = z
c        q--vapor quality on a MOLAR basis [moles vapor/total moles]
c           q < 0 indicates subcooled (compressed) liquid
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q > 1 indicates superheated vapor
c           q = -998 subcooled liquid, but quality not defined (p > Pc)
c           q = 999 indicates supercritical state (t > Tc) and (p > Pc)
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, a flag = -9.99998d6 is returned
c     ierr--error flag:  0 = successful
c                        4 = Pin < 0
c                        8 = x out of range (< 0 or > 1)
c                       12 = x out of range and P < 0
c                      240 = CRITP did not converge
c                      241 = SATP did not converge at bubble point
c                      242 = SATP did not converge at dew point
c                      243 = SATP (bubble pt) did not converge for 2-ph
c                      244 = SATP (dew pt) did not converge for 2-phase
c                      245 = TPFL2 did not converge
c                      246 = 2-phase iteration did not converge
c                      247 = TPRHO did not converge for single-phase
c                      248 = single-phase iteration did not converge
c                      249 = H out of range
c     herr--error string (character*255 variable if ierr<>0)
c
c  rewritten by E.W. Lemmon, NIST Physical & Chemical Properties Div, Boulder, CO
c  01-18-01 EWL, original version, rewritten by combining code of McLinden
c                from PHFLSH and PSFLSH routines to make one unified routine
c  05-16-01 EWL, add checks for h and s in the single phase
c  11-07-02 EWL, check for bmin>bmax
c  12-16-02 EWL, add checks for pseudo-pure fluid, calculate t using q
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
cDEC$ ATTRIBUTES DLLEXPORT :: PBFLSH
c     dll_export PBFLSH
c
      parameter (ncmax=20)        !max number of components in mixture
      character*1 htab,hnull,bt
      character*2 ab
      character*255 herr,herr1,herr2
      dimension x(ncmax),y(ncmax),z(ncmax),
     &          xdew(ncmax),ydew(ncmax),xbub(ncmax),ybub(ncmax)
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
c  flags to GUI indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
      common /FLAGS2/ iamwat,ianc,iwat
      bt=ab(2:2)
c
      icomp=0
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
      endif
      ierr=0
      herr=' '
      D=0.d0
      Dl=0.d0
      Dv=0.d0
c
c  set output liquid and vapor compositions to input values
c  zero output values for undefined components
      do i=1,nc
        x(i)=z(i)
        y(i)=z(i)
      enddo
c     if (nc.lt.ncmax) then
c       do i=nc+1,ncmax
c         x(i)=0.0d0
c         y(i)=0.0d0
c       enddo
c     end if
c
      call CRITP (z,tc,pc,Dc,ierr,herr1)
      t=tc
      if (ierr.ne.0) then
        ierr=240
        write (herr,1002) ab,herr1(1:234),hnull
 1002   format ('[',a2,'FLSH error 240] ',a234,a1)
        call ERRMSG (ierr,herr)
        q=999.d0     !quality undefined
        RETURN
      end if
c
c  check that input conditions are within limits
c
      tdum=-1.0d0
      Ddum=0.0d0
      call LIMITX ('EOS',tdum,Ddum,p,z,tmin,tmax,Dmax,pmax,ierr,herr1)
c  calculate approx limits:
c  lower limit at (tmin,pmax); upper limit at (1.5tmax,rho = 0)
c  !MM, 10-16-97, because of curvature of isotherms, h at pmax is not
c  the minimum, use a lower pressure instead

      p0=100.0d0 !kPa
      call TPRHO (tmin,p0,z,1,0,D0,ierr2,herr2)
      if (bt.eq.'H') call ENTHAL (tmin,D0,z,bmin)
      if (bt.eq.'S') call ENTRO (tmin,D0,z,bmin)
      if (bt.eq.'E') call ENERGY (tmin,D0,z,bmin)
      tmax=1.5d0*tmax
      rho0=1.0d-6
      if (bt.eq.'H') call ENTHAL (tmax,rho0,z,bmax)
      if (bt.eq.'S') call ENTRO (tmax,rho0,z,bmax)
      if (bt.eq.'E') call ENERGY (tmax,rho0,z,bmax)
      if (bmin.gt.bmax) then
c  just in case tmin was so low that the EOS did weird things, raise tmin:
        tmin=tmin+20
        call TPRHO (tmin,p0,z,1,0,D0,ierr2,herr2)
        if (bt.eq.'H') call ENTHAL (tmin,D0,z,bmin)
        if (bt.eq.'S') call ENTRO (tmin,D0,z,bmin)
        if (bt.eq.'E') call ENERGY (tmin,D0,z,bmin)
      endif
c     write (*,1004) tmin,tmax,bmin,bmax
c1004 format (1x,' ABFLSH--limits on t,b:  ',2f8.2,2f12.2)
c  if inputs are outside limits set outputs equal to critical
c  point values and return
      if (ierr.lt.-1) then     !ignore ierr = -1 (t out of range)
c  one or inputs are outside limits--if just a warning proceed w/ calc
        write (herr,1006) ab,ierr,herr1(1:233),hnull
 1006   format ('[',a2,'FLSH warning',i4,'] ',a233,a1)
        call ERRMSG (ierr,herr)
      else
        if (ierr.gt.1) then
c  ignore ierr=1 out of LIMITX (t out of range, but t not an input)
          write (herr,1007) ab,ierr,herr1(1:235),hnull
 1007     format ('[',a2,'FLSH error',i4,'] ',a235,a1)
        else if (b.lt.bmin .or. b.gt.bmax) then
          ierr=249
          write (herr,1008) ab,ierr,bt,b,bmin,bmax,hnull
 1008     format ('[',a2,'FLSH error',i4,'] Input value is ',
     &            'outside limits; ',a1,' =',g11.5,', min,max =',
     &            2(g11.5),a1)
        end if
        if (ierr.gt.1) then
          call ERRMSG (ierr,herr)
          t=tc
          D=Dc
          Dl=Dc
          Dv=Dc
          call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
          q=999.d0     !quality undefined
          RETURN
        end if
      end if
c
      iflag1=0
      if (p.ge.pc .or. p.lt.1.0d-12) then
        iflag1=1     !Single phase
      else
c  if the input enthalpy is greater than h at the critical temperature and zero
c  density, the inputs are single phase.
        if (bt.eq.'H') then
          call enthal(tc,0.d0,x,htc0)
          if (b.gt.htc0) iflag1=1
        endif
c  if the input entropy is greater than S at the critical temperature and given
c  pressure, the inputs are single phase.
        if (bt.eq.'S') then
          call TPRHO (tc,p,x,2,0,dtcp,ierr,herr)
          call entro(tc,dtcp,x,stcp)
          if (b.gt.stcp .and. ierr.eq.0) iflag1=1
        endif
      endif
      if (iflag1.eq.1) then
c  super-critical state or special case for p = 0 (x = y = z as set above)
        kph=2                          !use vapor iteration
        t=tc+10                        !initial guesses
        if (p.ge.pc) then
          D=Dc*2.d0
        else
          D=0.0d0
        end if
        call ENTHAL (tc,dc,z,h)
        call ENERGY (tc,dc,z,e)
        if (bt.eq.'H' .and. b.lt.h/2) t=tc*.8d0
        if (bt.eq.'E' .and. b.lt.e/2) t=tc*.8d0
        call ABFL1 (p,b,z,kph,ab,0.d0,0.d0,t,pp,D,ierr,herr)
        if (ierr.ne.0) then
c  single-phase iteration did not converge--error message written by ABFL1
c  try again with different initial input
          t=tc*1.2d0
          call ABFL1 (p,b,z,kph,ab,0.d0,0.d0,t,pp,D,ierr,herr)
          if (ierr.ne.0) then
            q=999.d0     !quality undefined
            RETURN
          end if
        end if
        call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
        Dl=D
        Dv=D
        if (t.le.tc) then
          q=-998.d0    !subcooled, but cannot define quality
        else
          q=999.d0     !supercritical--quality undefined
        end if
      else
c  sub-critical state--call saturation routine to determine phase
        call SATP (p,z,2,tdew,Dldew,Dvdew,xdew,ydew,ierr,herr1)
        if (ierr.eq.2) then     !p<ptrp
          tdew=tmin
          call SATT (tdew,z,2,pp,Dldew,Dvdew,xdew,ydew,ierr,herr1)
        elseif (ierr.ne.0) then
c  if the dew point iteration fails, make a check to see if the point is
c  in the liquid phase below the liquid saturation state.
          call SATP (p,z,1,tbub,Dlbub,Dvbub,xbub,ybub,ierr2,herr2)
          if (ierr2.eq.0) then
            if (bt.eq.'H') then
              call ENTHAL (tbub,Dlbub,z,blbub)
            elseif (bt.eq.'S') then
              call ENTRO (tbub,Dlbub,z,blbub)
            elseif (bt.eq.'E') then
              call ENERGY (tbub,Dlbub,z,blbub)
            endif
            if (b.le.blbub) then
              kph=1
              t=tbub
              D=Dlbub
              call ABFL1 (p,b,z,kph,ab,0.d0,0.d0,t,pp,D,ierr2,herr2)
              call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
              Dl=D
              Dv=D
              if (ierr2.eq.0) goto 100
            endif
          endif
          ierr=242
          write (herr,1242) ab,herr1(1:192),hnull
          call ERRMSG (ierr,herr)
 1242     format ('[',a2,'FLSH error 242] dew point calculation ',
     &            'did not converge:  ',a192,a1)
          RETURN
        end if
        if (icomp.ne.0 .and. ianc.eq.0) then
c  if pure, bubble point density provided by call to SATT at dew point
          Dlbub=Dldew
          tbub=tdew
        else
c  if mixture, call saturation routine again at bubble point
          call SATP (p,z,1,tbub,Dlbub,Dvbub,xbub,ybub,ierr,herr1)
          if (ierr.ne.0 .and. p.gt.pc*0.9d0) then
            tbub=tc
            Dlbub=dc
          elseif (ierr.ne.0) then
            ierr=241
            write (herr,1241) ab,herr1(1:188),hnull
 1241       format ('[',a2,'FLSH error 241] bubble point calculation ',
     &              'did not converge:  ',a188,a1)
            call ERRMSG (ierr,herr)
            RETURN
          end if
        end if
c  calculate b at dew and bubble points
        if (bt.eq.'H') then
          call ENTHAL (tdew,Dvdew,z,bvdew)
          call ENTHAL (tbub,Dlbub,z,blbub)
        elseif (bt.eq.'S') then
          call ENTRO (tdew,Dvdew,z,bvdew)
          call ENTRO (tbub,Dlbub,z,blbub)
        elseif (bt.eq.'E') then
          call ENERGY (tdew,Dvdew,z,bvdew)
          call ENERGY (tbub,Dlbub,z,blbub)
        endif
        q=(b-blbub)/(bvdew-blbub)
        if (b.le.blbub .or. b.ge.bvdew) then
c  single-phase (liq or vapor) (pure or mixture) (x = y = z as set above)
c  set initial guesses, depending on phase
          if (b.le.blbub) then
            kph=1                        !liquid
            t=tbub                       !initial guesses
            D=Dlbub
          else
            kph=2                        !vapor
            t=tdew                       !initial guesses
            D=Dvdew
          end if
          call ABFL1 (p,b,z,kph,ab,0.d0,0.d0,t,pp,D,ierr,herr)
          if (ierr.ne.0) then
c  single-phase iteration did not converge--error message written by ABFL1
            call ERRMSG (ierr,herr)
            q=999.d0     !quality undefined
            RETURN
          end if
          call THERM (t,D,z,ptherm,e,h,s,cv,cp,w,hjt)
          Dl=D
          Dv=D
        else
          if (icomp.ne.0) then
c  special case:  pure-fluid two-phase (x = y = z as set above)
            t=(1.d0-q)*tbub+q*tdew
            Dl=Dlbub
            Dv=Dvdew
            if (ianc.eq.1) call PTANC (t,p,q,b,bt,Dl,Dv)
          else
c  general case:  mixture 2-phase
            ksat=1   !bubble and dew point data provided to ABFL2
            call ABFL2 (p,b,z,ksat,0,ab,
     &                  tbub,tdew,pbub,pdew,Dlbub,Dvdew,ybub,xdew,
     &                  t,pp,Dl,Dv,x,y,q,ierr,herr)
            if (ierr.ne.0) then
c  two-phase iteration did not converge--error message written by ABFL2
              call ERRMSG (ierr,herr)
              q=999.d0     !quality undefined
              RETURN
            end if
c           write (*,1225) t,q,(x(i),i=1,3),(y(i),i=1,3)
c1225       format (1x,' PBFLSH--ABFL2 return   t,q,x,y = ',8f12.7)
          end if
c  compute remaining properties for 2-phase states
          call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt)
          call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt)
          alpha=1.0d0-q
          D=1.0d0/(alpha/Dl+q/Dv)
          e=alpha*el+q*ev
          h=alpha*hl+q*hv
          s=alpha*sl+q*sv
          w=xnotd       !Cp,w not defined for 2-phase states
          cp=xnotd
          cv=xnotd
        end if
      end if
c
c  check that computed output is within bounds of EOS
 100  continue
      call LIMITX ('EOS',t,d,p,z,tmin,tmax,Dmax,pmax,ierr,herr1)
      if (ierr.ne.0) then
        write (herr,2008) ab,ierr,herr1
 2008   format ('[',a2,'FLSH error',i3,']; ',a235)
        call ERRMSG (ierr,herr)
      end if
c
      RETURN
      end                                             !subroutine PBFLSH
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                    end file flsh_sub.f
c ======================================================================
