c  begin file flash2.f
c
c  This file contains additional iterative routines which call the
c  intermediate level routines
c
      subroutine TQFLSH (t,q,z,kq,p,D,Dl,Dv,x,y,e,h,s,cv,cp,w,ierr,herr)
c
c  flash calculation given temperature, quality, and bulk composition
c
c  This routine accepts saturation or two-phase states as inputs.
c
c  inputs:
c        t--temperature [K]
c        q--vapor quality [basis specified by kq]
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q < 0 or q > 1 are not allowed and will result in warning
c        z--overall (bulk) composition [array of mol frac]
c       kq--flag specifying units for input quality
c           kq = 1 quality on MOLAR basis [moles vapor/total moles]
c           kq = 2 quality on MASS basis [mass vapor/total mass]
c
c  outputs:
c        p--pressure [kPa]
c        D--overall (bulk) molar density [mol/L]
c       Dl--molar density [mol/L] of the liquid phase
c       Dv--molar density [mol/L] of the vapor phase
c        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, -9.992d6 is returned
c     ierr--error flag:   0 = successful
c                         1 = T < Tmin
c                         8 = x out of range
c                         9 = T and x out of range
c                       270 = CRITP did not converge
c                       271 = T > Tcrit
c                       275 = q out of range
c                       278 = TQFLSH did not converge
c                      -279 = 2-phase for mix not implemented yet
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  04-08-96  MM, original version, based on TPFLSH
c  05-16-97  MM, get special-case flags from /FLAGS/
c  07-14-97  MM, add check for q out of range; add call to LIMITS
c  10-01-97  MM, add compiler switches to allow access by DLL
c  12-05-97  MM, change warning -279 to error 279
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: TQFLSH
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
      parameter (ncmax=5)        !max number of components in mixture
      character*1 htab,hnull
      character*255 herr,herr2
      dimension z(ncmax),x(ncmax),y(ncmax)
      dimension xdew(ncmax),ybub(ncmax)
      common /Gcnst/ R
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
c  flags to GUI indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
      ierr=0
      herr=hnull
c
c  zero output liquid and vapor compositions for undefined components
      do 100 i=nc+1,ncmax
      x(i)=0.0d0
      y(i)=0.0d0
  100 continue
c
c  check that input conditions (in this case t and z) are within limits
c
      Ddum=0.0d0
      pdum=0.0d0
      call LIMITX ('EOS',t,Ddum,pdum,z,tmin,tmax,rhomax,pmax,ierr,herr2)
      if (ierr.gt.0) then
c  T and/or x are out of bounds, set error flag and return
        p=0.0d0
        D=0.0d0
        Dl=0.0d0
        Dv=0.0d0
        do 110 i=1,nc
        x(i)=z(i)
        y(i)=z(i)
  110   continue
        write (herr,1110) ierr,herr2(1:236),hnull
 1110   format ('[TQFLSH error',i3,'] ',a236,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
      call CRITP (z,tc,pc,rhoc,ierr,herr2)
      if (ierr.ne.0) then
        ierr=270
        p=0.0d0
        D=0.0d0
        Dl=0.0d0
        Dv=0.0d0
        do 111 i=1,nc
        x(i)=z(i)
        y(i)=z(i)
  111   continue
        write (herr,1111) herr2(1:235),hnull
 1111   format ('[TQFLSH error 270] ',a235,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
      if (t.ge.tc) then
c  supercritical state
        D=0.0d0
        Dl=0.0d0
        Dv=0.0d0
        do 120 i=1,nc
        x(i)=z(i)
        y(i)=z(i)
  120   continue
        ierr=271
        write (herr,1271) t,tc,hnull
 1271   format ('[TQFLSH error 271] ',
     &          'temperature input to temperature-quality flash is ',
     &          'greater than critical temperature; T =',1pe11.3,
     &          ' K, Tcrit =',1pe11.3,' K.',a1)
        call ERRMSG (ierr,herr)
        RETURN
      else if (abs(q).le.1.0d-8) then
c  saturated liquid
        call SATT (t,z,1,p,Dl,Dv,x,y,ierr,herr2)
        if (ierr.ne.0) then
          herr=' ERROR from TQFLSH (sat liquid):  '//herr2
          call ERRMSG (ierr,herr)
          RETURN
        end if
        D=Dl
        call THERM (t,Dl,x,p,e,h,s,cv,cp,w,hjt)
      else if (abs(1.0d0-q).le.1.0d-8) then
c  saturated vapor
        call SATT (t,z,2,p,Dl,Dv,x,y,ierr,herr2)
        if (ierr.ne.0) then
          herr=' ERROR from TQFLSH (sat vapor):  '//herr2
          call ERRMSG (ierr,herr)
          RETURN
        end if
        D=Dv
        call THERM (t,Dv,y,p,e,h,s,cv,cp,w,hjt)
      else if (q.lt.0.0d0 .or. q.gt.1.0d0) then
c  quality out of range
        ierr=275
        write (herr,1275) q,hnull
 1275   format ('[TQFLSH error 275] ',
     &          'input quality is out of range; q =',1pe11.3,
     &          '; quality must be between 0 and 1 ',a1)
        call ERRMSG (ierr,herr)
        RETURN
      else
c  two-phase state
        call SATT (t,z,1,pbub,Dlbub,Dvbub,x,ybub,ierr,herr2)
        call SATT (t,z,2,pdew,Dldew,Dvdew,xdew,y,ierr,herr)
        if (ierr.ne.0) then
          herr=' ERROR from TQFLSH (2-phase):  '//herr2
          call ERRMSG (ierr,herr)
          RETURN
        end if
        if (nc.eq.1) then
c  special case--two-phase state for a pure fluid
          Dl=Dlbub
          Dv=Dvbub
          p=pbub
          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=x2ph                 !Cp, w not defined for 2-phase
          w=x2ph
          cv=x2ph
c  bulk properties are weighted average of liquid and vapor phases
c  note that for a pure fluid quality on mass and molar basis is same
          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
        else
c
c  general 2-phase mixture state
c  generate initial guesses for t, x, y by interpolating sat liq & vap
c
          xsum=0.0d0
          ysum=0.0d0        !sums for normalization of compositions
          do 240 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)
  240     continue
          do 244 i=1,nc
          x(i)=x(i)/xsum
          y(i)=y(i)/ysum
  244     continue
          p=(1.0d0-q)*pbub+q*pdew  !initial guess for pressure
          Dl=Dlbub                 !initial guess for liquid density
          Dv=Dvdew                 !initial guess for vapor density
c  !temp--TQFL2 not yet finished, use initial guesses
c         call TQFL2 (t,z,q,kq,p,Dl,Dv,x,y,ierr,herr)
c  compute 2-phase properties and load output variables
c  call TPRHO to ensure consistency of t, p, rho
          call TPRHO (t,p,x,1,kguess,Dl,ierr,herr)
          call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt)
          call TPRHO (t,p,y,2,kguess,Dv,ierr,herr)
          call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt)
          ierr=279
          herr='[TQFLSH error 279] temperature-quality flash '//
     &         'calculations are not implemented for mixtures'//hnull
          cp=x2ph              !Cp, w not defined for 2-phase
          w=x2ph
          cv=x2ph
c  bulk properties are weighted average of liquid and vapor phases
          if (kq.eq.1) then
            qmole=q                 !qmole is vapor frac [molar basis]
            alpha=1.0d0-q           !alpha is liq fraction
          else
c  convert mass quality to molar quality
            wmliq=WMOL(x)
            wmvap=WMOL(y)
            qmole=q/wmvap/(q/wmvap+(1.0d0-q)/wmliq)
            alpha=1.0d0-q
          end if
          D=1.0d0/(alpha/Dl+qmole/Dv)
          e=alpha*el+qmole*ev
          h=alpha*hl+qmole*hv
          s=alpha*sl+qmole*sv
        end if
      end if
c
      RETURN
      end                                             !subroutine TQFLSH
c
c ======================================================================
c
      subroutine TQFL2 (t,q,z,kq,p,Dl,Dv,x,y,ierr,herr)
c
c  flash calculation given temperature, quality, and bulk composition
c
c  This routine accepts only two-phase states as input and requires good
c  initial guesses for pressure, density, and coexisting compositions.
c  If these are not available, use TQFLSH.
c
c  inputs:
c        t--temperature [K]
c        q--vapor quality [basis specified by kq]
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q < 0 or q > 1 are not allowed and will result in warning
c        z--overall (bulk) composition [array of mol frac]
c       kq--flag specifying units for input quality
c           kq = 1 quality on MOLAR basis [moles vapor/total moles]
c           kq = 2 quality on MASS basis [mass vapor/total mass]
c        p--initial guess for pressure [kPa]
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
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     ierr--error flag:  0 = successful
c                        1 = 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-08-96  MM, original version, based on TPFL2
c  08-19-97  MM, return warning for kq=2--mass quality not yet implemented
c  10-01-97  MM, add compiler switches to allow access by DLL
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: TQFL2
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
      parameter (ncmax=5)        !max number of components in mixture
      character*1 htab,hnull
      character*255 herr
      dimension z(ncmax),x(ncmax),y(ncmax)
      dimension alpha(3),falpha(2),fliq(ncmax),fvap(ncmax),xk(ncmax)
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
c  !temporary--convergence tolerance, put in common block
      data tolr /1.0d-8/
      data itmax /20/
c  !temporary--end
c
c  !temporary--need to implement mass quality option
      if (kq.eq.2) then
         ierr=-279
         herr='[TQFLSH warning 279] mass quality option (kq = 2) '//
     &        'not implemented; input interpreted as molar quality.'//
     &         hnull
      else
        ierr=0
        herr=hnull
      end if
c
c  begin outer iteration loop for the vapor compostion
c
      alpha(1)=1.0d0-q
      kguess=1
      do 800 iy=1,itmax
c  compute densities and fugacities for each phase
      call TPRHO (t,p,x,1,kguess,Dl,ierr,herr)
      call FGCTY (t,Dl,x,fliq)
      call TPRHO (t,p,y,2,kguess,Dv,ierr,herr)
      call FGCTY (t,Dv,y,fvap)
c  compute k-factors (assume constant for a given guess of y)
      do 520 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
  520 continue
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 600 ix=1,itmax
      xsum=0.0d0
      do 540 i=1,nc
      if (z(i).gt.0.0d0) then
c  compute new guess for liquid compostion
        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)
  540 continue
c  normalize liquid composition
      do 560 i=1,nc
      x(i)=x(i)/xsum
  560 continue
      falpha(j)=1.0d0-xsum
d     write (*,1056) ix,alpha(j),xsum,(x(i),i=1,3)
d1056 format (1x,' TQFL2:',i4,35x,2f10.6,3f10.6)
      if (abs(falpha(j)).lt.1.0d-4*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
        j=2
        alpha(2)=alpha(1)+0.05d0*(1.0d0-alpha(1))
      else
c  for subsequent iterations, use secant method
        alpha(3)=alpha(2)-falpha(2)*(alpha(2)-alpha(1))/
     &           (falpha(2)-falpha(1))
c  new guess for alpha must lie within bounds of zero and one
        if (alpha(3).gt.1.0d0) then
          alpha(3)=1.0d0
        else if (alpha(3).lt.0.0d0) then
          alpha(3)=0.0d0
        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
  600 continue
c  inner iteration loop has not converged, issue warning and proceed
      ierr=-1
      herr=' TQFL2 WARNING--inner iteration did not converge'//hnull
d     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 740 i=1,nc
      y(i)=x(i)*xk(i)
      ysum=ysum+y(i)
  740 continue
c  normalize vapor compositions
      do 760 i=1,nc
      y(i)=y(i)/ysum
  760 continue
d     write (*,1076) iy,t,Dl,Dv,ysum,(y(i),i=1,3)
d1076 format (1x,' TQFL2:',i4,f7.2,2e14.6,10x,f10.6,3f10.6)
      if (abs(1.0d0-ysum).lt.tolr) then
c  outer loop has converged
        goto 840
      end if
  800 continue
c  outer iteration loop has not converged, issue error and proceed
      ierr=1
      herr=' TQFL2 ERROR--outer iteration did not converge'//hnull
d     call ERRMSG (ierr,herr)
c
c  end of outer iteration loop for vapor composition
c
  840 continue
c
      RETURN
      end                                              !subroutine TQFL2
c
c ======================================================================
c
      subroutine PQFLSH (p,q,z,kq,t,D,Dl,Dv,x,y,e,h,s,cv,cp,w,ierr,herr)
c
c  flash calculation given pressure, quality, and bulk composition
c
c  This routine accepts saturation or two-phase states as inputs.
c
c  inputs:
c        p--pressure [kPa]
c        q--vapor quality [basis specified by kq]
c           q = 0 indicates saturated liquid
c           q = 1 indicates saturated vapor
c           q < 0 or q > 1 are not allowed and will result in warning
c        z--overall (bulk) composition [array of mol frac]
c       kq--flag specifying units for input quality
c           kq = 1 quality on MOLAR basis [moles vapor/total moles]
c           kq = 2 quality on MASS basis [mass vapor/total mass]
c
c  outputs:
c        t--temperature [K]
c        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        x--composition of liquid phase [array of mol frac]
c        y--composition of vapor phase [array of mol frac]
c        e--overall (bulk) internal energy [J/mol]
c        h--overall (bulk) enthalpy [J/mol]
c        s--overall (bulk) entropy [J/mol-K]
c       Cv--isochoric (constant V) heat capacity [J/mol-K]
c       Cp--isobaric (constant p) heat capacity [J/mol-K]
c        w--speed of sound [m/s]
c           Cp, w are not defined for 2-phase states
c           in such cases, -9.992d6 is returned
c     ierr--error flag:   0 = successful
c                         4 = P < 0
c                         8 = x out of range
c                         9 = P and x out of range
c                       290 = CRITP did not converge
c                       291 = P > Pcrit
c                       295 = q out of range
c                       298 = PQFLSH did not converge
c                      -299 = 2-phase for mix not implemented yet
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  04-08-96  MM, original version, based on TQFLSH
c  04-23-96  MM, alpha = 1 - qmole, not 1 - q when q input as mass
c  05-16-97  MM, get special-case flags from /FLAGS/
c  07-14-97  MM, add check for q out of range; add call to LIMITS
c  10-01-97  MM, add compiler switches to allow access by DLL
c  12-05-97  MM, change warning -299 to error 299
c
c  compiler switches to allow access by DLL; for use with Digital Visual 
c  Fortran; these should be treated as comments by all other compilers
c
!MS$ATTRIBUTES DLLEXPORT :: PQFLSH
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
      parameter (ncmax=5)        !max number of components in mixture
      character*1 htab,hnull
      character*255 herr,herr2
      dimension z(ncmax),x(ncmax),y(ncmax)
      dimension xdew(ncmax),ybub(ncmax)
      common /Gcnst/ R
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
c  flags to GUI indicating 'not applicable', '2-phase', etc.
      common /FLAGS/ xnota,x2ph,xsubc,xsuph,xsupc,xinf,x7,xnotd,xnotc
c
      ierr=0
      herr=hnull
c
c  zero output liquid and vapor compositions for undefined components
      do 100 i=nc+1,ncmax
      x(i)=0.0d0
      y(i)=0.0d0
  100 continue
c
      call CRITP (z,tc,pc,rhoc,ierr,herr2)
      if (ierr.ne.0) then
        ierr=290
        t=300.0d0
        D=0.0d0
        Dl=0.0d0
        Dv=0.0d0
        do 108 i=1,nc
        x(i)=z(i)
        y(i)=z(i)
  108   continue
        write (herr,1108) herr2(1:235),hnull
 1108   format ('[PQFLSH error 290] ',a235,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
c  check that input conditions (in this case p and z) are within limits
c
      Ddum=0.0d0
      tdum=0.8d0*tc
      call LIMITX ('EOS',tdum,Ddum,p,z,tmin,tmax,rhomax,pmax,ierr,herr2)
      if (ierr.gt.0) then
c  p and/or x are out of bounds, set error flag and return
c  set output temperature to "reasonable" value to avoid GUI crash
        t=0.8d0*tc
        D=0.0d0
        Dl=0.0d0
        Dv=0.0d0
        do 110 i=1,nc
        x(i)=z(i)
        y(i)=z(i)
  110   continue
        write (herr,1110) ierr,herr2(1:236),hnull
 1110   format ('[PQFLSH error',i3,'] ',a236,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
      if (p.ge.pc) then
c  supercritical state
c  set output temperature to critical value to avoid GUI crash
        t=tc
        D=0.0d0
        Dl=0.0d0
        Dv=0.0d0
        do 120 i=1,nc
        x(i)=z(i)
        y(i)=z(i)
  120   continue
        ierr=291
        write (herr,1291) p,pc,hnull
 1291   format ('[PQFLSH error 291] ',
     &          'pressure input to pressure-quality flash is ',
     &          'greater than critical pressure; P =',1pe11.3,
     &          ' kPa, Pcrit =',1pe11.3,' kPa.',a1)
        call ERRMSG (ierr,herr)
        RETURN
      else if (abs(q).le.1.0d-8) then
c  saturated liquid
        call SATP (p,z,1,t,Dl,Dv,x,y,ierr,herr2)
        if (ierr.ne.0) then
          herr=' ERROR from PQFLSH (sat liquid):  '//herr2
          call ERRMSG (ierr,herr)
          RETURN
        end if
        D=Dl
        call THERM (t,Dl,x,ptherm,e,h,s,cv,cp,w,hjt)
      else if (abs(1.0d0-q).le.1.0d-8) then
c  saturated vapor
        call SATP (p,z,2,t,Dl,Dv,x,y,ierr,herr2)
        if (ierr.ne.0) then
          herr=' ERROR from PQFLSH (sat vapor):  '//herr2
          call ERRMSG (ierr,herr)
          RETURN
        end if
        D=Dv
        call THERM (t,Dv,y,ptherm,e,h,s,cv,cp,w,hjt)
      else if (q.lt.0.0d0 .or. q.gt.1.0d0) then
c  quality out of range
        ierr=295
        write (herr,1295) q,hnull
 1295   format ('[PQFLSH error 295] ',
     &          'input quality is out of range; q =',1pe11.3,
     &          '; quality must be between 0 and 1 ',a1)
        call ERRMSG (ierr,herr)
        RETURN
      else
c  two-phase state
        call SATP (p,z,1,tbub,Dlbub,Dvbub,x,ybub,ierr,herr2)
        call SATP (p,z,2,tdew,Dldew,Dvdew,xdew,y,ierr,herr2)
        if (ierr.ne.0) then
          herr=' ERROR from PQFLSH (2-phase):  '//herr2
          call ERRMSG (ierr,herr)
          RETURN
        end if
        if (nc.eq.1) then
c  special case--two-phase state for a pure fluid
          t=tbub
          Dl=Dlbub
          Dv=Dvbub
          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=x2ph                 !Cp, w not defined for 2-phase
          w=x2ph
          cv=x2ph
c  bulk properties are weighted average of liquid and vapor phases
c  note that for a pure fluid quality on mass and molar basis is same
          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
        else
c
c  general 2-phase mixture state
c  generate initial guesses for t, x, y by interpolating sat liq & vap
c
          xsum=0.0d0
          ysum=0.0d0        !sums for normalization of compositions
          do 240 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)
  240     continue
          do 244 i=1,nc
          x(i)=x(i)/xsum
          y(i)=y(i)/ysum
  244     continue
          t=(1.0d0-q)*tbub+q*tdew  !initial guess for temperature
          Dl=Dlbub                 !initial guess for liquid density
          Dv=Dvdew                 !initial guess for vapor density
c  !temp--PQFL2 not yet finished, use initial guesses
c         call PQFL2 (p,q,z,kq,t,Dl,Dv,x,y,ierr,herr)
c  compute 2-phase properties and load output variables
c  call TPRHO to ensure consistency of t, p, rho
          call TPRHO (t,p,x,1,kguess,Dl,ierr,herr)
          call THERM (t,Dl,x,ptherm,el,hl,sl,cvl,cp,w,hjt)
          call TPRHO (t,p,y,2,kguess,Dv,ierr,herr)
          call THERM (t,Dv,y,ptherm,ev,hv,sv,cvv,cp,w,hjt)
          ierr=299
          herr='[PQFLSH error 299] pressure-quality flash '//
     &         'calculations are not implemented for mixtures'//hnull
          cp=x2ph              !Cp, w not defined for 2-phase
          w=x2ph
          cv=x2ph
c  bulk properties are weighted average of liquid and vapor phases
          if (kq.eq.1) then
            qmole=q                 !qmole is vapor frac [molar basis]
            alpha=1.0d0-q           !alpha is liq fraction
          else
c  convert mass quality to molar quality
            wmliq=WMOL(x)
            wmvap=WMOL(y)
            qmole=q/wmvap/(q/wmvap+(1.0d0-q)/wmliq)
            alpha=1.0d0-qmole
          end if
          D=1.0d0/(alpha/Dl+qmole/Dv)
          e=alpha*el+qmole*ev
          h=alpha*hl+qmole*hv
          s=alpha*sl+qmole*sv
        end if
      end if
c
      RETURN
      end                                             !subroutine PQFLSH
c
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                      end file flash2.f
c ======================================================================