c  begin file sat_sub.f
c
c  This file contains routines for saturation properties
c
c  contained here are:
c     subroutine SATT (t,x,kph,p,rhol,rhov,xliq,xvap,ierr,herr)
c     subroutine SATP (p,x,kph,t,rhol,rhov,xliq,xvap,ierr,herr)
c     subroutine SURFT (t,rhol,xl,sigma,ierr,herr)
c
c  these routines use the following common blocks from other files
c     common /Gcnst/ R
c     common /NCOMP/ nc
c     common /CHAR/ htab,hnull
c     common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
c    &              tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx),
c    &              accen(n0:nx),dipole(n0:nx)
c
c  various arrays are dimensioned with parameter statements
c     parameter (ncmax=5)        !max number of components in mixture
c     parameter (n0=-ncmax,nx=ncmax)
c
c ======================================================================
c ======================================================================
c
      subroutine SATT (t,x,kph,p,rhol,rhov,xliq,xvap,ierr,herr)
c
c  iterate for saturated liquid and vapor states given temperature
c  and the composition of one phase
c
c  inputs:
c        t--temperature [K]
c        x--composition [array of mol frac] (phase specified by kph)
c      kph--phase flag: 1 = input x is liquid composition (bubble point)
c                       2 = input x is vapor composition (dew point)
c  outputs:
c        p--pressure [kPa]
c     rhol--molar density [mol/L] of saturated liquid
c     rhov--molar density [mol/L] of saturated vapor
c     xliq--liquid phase composition [array of mol frac]
c     xvap--vapor phase composition [array of mol frac]
c     ierr--error flag:   0 = successful
c                         1 = T < Tmin
c                         8 = x out of range
c                         9 = T and x out of range
c                       120 = CRITP did not converge
c                       121 = T > Tcrit
c                       122 = TPRHO-liquid did not converge (pure fluid)
c                       123 = TPRHO-vapor did not converge (pure fluid)
c                       124 = pure fluid iteration did not converge
c           following 3 error codes are advisory--iteration will either
c           converge on later guess or error out (ierr = 128)
c                      -125 = TPRHO did not converge for parent ph (mix)
c                      -126 = TPRHO did not converge for incipient (mix)
c                      -127 = composition iteration did not converge
c                       128 = mixture iteration did not converge
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  07-11-95  MM, original version
c  09-11-95  MM, add error string to argument list
c  09-25-95  MM, rearrange argument list (outputs in order p, rho, x)
c  10-06-95  MM, use stored accentric factor for pure fluids
c  10-11-95  MM, RETURN if any error detected
c  11-26-95  MM, Raoult's law as first guess for mixture
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accomodate ECS reference fluid
c  12-18-95  MM, fill xliq, xvap with zeros for undefined components
c  12-19-20  MM, add full mixture iteration using fugacity
c  12-27-95  MM, pratio for new pressure if no converge for TPRHO for mix
c  01-09-96  MM, move check for supercritical outside nc = 1 block
c  02-27-96  MM, parameter n0=-ncmax to accomodate ECS-thermo model
c                add Zcrit to common /CCON/
c  03-06-96  MM, reset warning from TPRHO if loop eventually converges
c  03-19-96  MM, add dipole moment to /CCON/
c  04-05-96  MM, test for supercritical '.ge. tc' rather than '.gt. tc'
c  05-30-96  MM, check input temperature against limits
c  06-03-96  MM, add 'EOS' to calling list for LIMITX
c  06-05-96  MM, refine error numbers and messages;
c                also ensure that all outputs are set on error condition
c  11-14-96  MM, adjust pratio, etc to get closer to critical
c   2-07-96 EWL, add pressure increment/decrement when TPRHO does not converge
c   2-13-96 EWL, add initial guess for densities near the critical point
c                return critical point values if within delta of Tc
c   6-06-96 EWL, return critical point values for failure to converge if
c                Tc - T < 10 mK
c  10-01-97  MM, add compiler switches to allow access by DLL
c  11-13-97 EWL, initilize fpit(j); potential bomb if no value when writing error message
c  11-14-97 EWL, add line following 422 to improve critical region convergence
c  02-09-98  MM, limit delp step, change delp if TRRHO does not converge
c  02-11-98  MM, check that new guess for mix pressure is < p_crit
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 :: SATT
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
      parameter (n0=-ncmax,nx=ncmax)
      character*1 htab,hnull
      character*255 herr,herr2
      dimension x(ncmax),xliq(ncmax),xvap(ncmax)
      dimension pcomp(ncmax),x2(ncmax),f1(ncmax),f2(ncmax),x2new(ncmax)
      dimension pit(3),fpit(2)
      common /Gcnst/ R
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(n0:nx)
c  !temporary--convergence tolerance, put in common block
      data tolr /1.0d-6/
      data itmax /20/
c  !temporary--end
c     write (*,*) ' SATT--entering with t,kph = ',t,kph
      ierr=0
      herr=hnull
      delp2=0.0d0
      fpit(1)=0.0d0  !initialize to avoid potential problem
      fpit(2)=0.0d0  !when writing error message to GUI
c
c  fill output composition arrays with zeros for undefined components
      if (nc.lt.ncmax) then
        do 100 i=nc+1,ncmax
        xliq(i)=0.0d0
        xvap(i)=0.0d0
  100   continue
      end if
c
c  check that input conditions (in this case t and x) are within limits
c
      Ddum=0.0d0
      pdum=0.0d0
      call LIMITX ('EOS',t,Ddum,pdum,x,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
        rhol=0.0d0
        rhov=0.0d0
        do 110 i=1,nc
        xliq(i)=x(i)
        xvap(i)=x(i)
  110   continue
        write (herr,1000) ierr,herr2(1:238),hnull
 1000   format ('[SATT error',i3,'] ',a238,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
      call CRITP (x,tc,pc,rhoc,ierr,herr2)
      if (ierr.ne.0) then
c  error condition--set outputs, issue warning, and return
        ierr=120
        p=0.0d0
        rhol=0.0d0
        rhov=0.0d0
        do 111 i=1,nc
        xliq(i)=x(i)
        xvap(i)=x(i)
  111   continue
        write (herr,1120) herr2(1:237),hnull
 1120   format ('[SATT error 120] ',a237,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
      if (t.gt.tc-1.0d-8) then
c  input temperature is equal to or greater than critical,
c  return critical parameters for output pressure and densities
        p=pc
        rhol=rhoc
        rhov=rhoc
        do 120 i=1,nc
        xliq(i)=x(i)
        xvap(i)=x(i)
  120   continue
        if (t.gt.tc+1.0d-8) then
c  supercritical temperature as input, set error flag and return
c  critical parameters for output pressure and densities
          ierr=121
          write (herr,1121) t,tc,hnull
 1121     format ('[SATT error 121] ',
     &            'temperature input to saturation routine is ',
     &            'greater than critical temperature; T =',1pe11.3,
     &            ' K, Tcrit =',1pe11.3,' K.',a1)
        end if
        RETURN
      endif
c
      if (nc.eq.1) then
c  pure fluid iteration
        xliq(1)=x(1)
        xvap(1)=x(1)
c
c  Iterate for saturated liquid and vapor states given temperature using
c  a simple successive substitution method.  The independent variable
c  in the iteration is the vapor pressure.  The convergence criteria
c  is equality of Gibbs free energy in both phases.
c
c  generate initial guess for pressure using accentric factor
c
        p=pc*10.0**(-2.333333*(1.0+accen(1))*(tc/t-1.0))
c       write (*,*) 'SATT:  P(0): ',p
c
        if (t.gt.0.99*tc) then
          kguess=1
          lcrit=.true.
c  for near-critical states, generate initial guesses for density;
c  using correlation developed by E.W. Lemmon, NIST
c         theta=(1.0d0-t/tc)**0.5d0*(2.5d0+accen(1))
c         p=pc*(1.0d0-1.000959d0*theta**2+0.308538d0*theta**3)
          theta=(1.0d0-t/tc)**(1.0d0/3.0d0)*(1.5d0+accen(1))
          rholi=rhoc*(1.0d0+1.113614d0*theta+0.080400d0*theta**2)
          rhol=rholi
          rhovi=rhoc*(1.0d0-1.078683d0*theta+5.014057d-2*theta**2)
          rhov=rhovi
          call PRESS (t,rhol,x,p1)
          call PRESS (t,rhov,x,p2)
          p = (p1 + p2)/2.0d0
c         write (*,1109) rhol,rhov,p1,p2,p
c1109     format (1x,' SATT--near-critical guesses for rhol,rhov,',
c    &               'pliq,pvap,pavg: ',5e14.6)
c  check that above guesses yield valid roots on calls to TPRHO
          call TPRHO (t,p,x,1,kguess,rhol,ierr1,herr2) !find liq density
          call TPRHO (t,p,x,2,kguess,rhov,ierr2,herr2) !find vap density
          if (ABS(rhol-rhov).lt.1.0d-6) then
c  try alternate guess for the pressure and try again
            p=pc*10.0**(-2.333333*(1.0+accen(1))*(tc/t-1.0))
c           write (*,*) ' SATT--pc, p_guess_alt:  ',pc,p
            rhol=rholi
            rhov=rhovi
c           call TPRHO (t,p,x,1,kguess,rhol,ierr1,herr2)
c           call TPRHO (t,p,x,2,kguess,rhov,ierr2,herr2)
          end if
        else
c  assume nothing about densities on initial calls to TPRHO
          kguess=0
          lcrit=.false.
        endif
        do 200 it=1,itmax
c       write (*,*) 'SATT--t,p input to TPRHO: ',t,p
        call TPRHO (t,p,x,1,kguess,rhol,ierr,herr2) !find liquid density
        if (ierr.ne.0) then
c         write (*,1111) 'liq',t,p
c1111     format (1x,' SATT--TPRHO (',a3,') did not conv for t,p = ',
c    &               2e14.6)
          p=p*1.005d0
          goto 200
c         ierr=122
c         write (herr,1122) it,herr2(1:146),hnull
c         call ERRMSG (ierr,herr)
c  set all outputs to some value (to keep GUI happy)
c         rhov=rhol
c         RETURN
        end if
        call TPRHO (t,p,x,2,kguess,rhov,ierr,herr2) !find vapor density
        if (ierr.ne.0) then
c         write (*,1111) 'vap',t,p
          p=p*0.95d0
          goto 200
c         ierr=123
c         write (herr,1123) it,herr2(1:147),hnull
c         call ERRMSG (ierr,herr)
c         RETURN
        end if
c  use previous densities as initial guesses for calls to TPRHO after
c  first iteration
        kguess=1
c
        call GIBBS (t,rhol,x,Aliq,Gliq)
        call GIBBS (t,rhov,x,Avap,Gvap)
        ZG=Gliq-Gvap
c       write (*,1014) it,t,rhol,rhov,p,Gliq,Gvap,ZG
c1014   format (1x,'SATT:  it,t,rhol,rhov,p,Gliq,Gvap,ZG: ',
c    &              i4,f8.3,2f12.8,e14.6,2f12.4,e14.6)
c
c  check convergence
c
c  check that liquid and vapor densities are different
        if (ABS(rhol-rhov).lt.1.0d-8) then
          ierr=124
          herr = '[SATT error 124] density roots equal'
          call ERRMSG (ierr,herr)
          RETURN
        end if
        delp=ZG/(1.0/rhol-1.0/rhov)
c  the delp-delp2 check is only important for very low pressures on the
c  liquid surface (propane or R124).  See comments in TPRHO.
        if (abs(delp/p).lt.tolr .or. abs(delp-delp2).lt.1.0d-11) then
          p=p-delp
          call TPRHO (t,p,x,1,kguess,rhol,ierr,herr2) !find liq density
          if (ierr.ne.0 .or. p.gt.pc .or. rhol.lt.rhoc) then
            ierr=122
            write (herr,1122) it,herr2(1:146),hnull
 1122       format ('[SATT error 122] liquid density iteration in ',
     &            'saturation routine did not converge for pressure ',
     &            'iteration',i3,'; ',a146,a1)
            call ERRMSG (ierr,herr)
c  return critical parameters if not converged and very close to Tc
            if (t.gt.0.999975*tc) then
              p=pc
              rhol=rhoc
              rhov=rhoc
            end if
            RETURN
          end if
          call TPRHO (t,p,x,2,kguess,rhov,ierr,herr2) !find vap density
          if (ierr.ne.0 .or. p.gt.pc .or. rhov.gt.rhoc) then
            ierr=123
            write (herr,1123) it,herr2(1:147),hnull
 1123       format ('[SATT error 123] vapor density iteration in ',
     &            'saturation routine did not converge for pressure ',
     &            'iteration',i3,'; ',a147,a1)
            call ERRMSG (ierr,herr)
c  return critical parameters if not converged and very close to Tc
            if (t.gt.0.999975*tc) then
              p=pc
              rhol=rhoc
              rhov=rhoc
            end if
            RETURN
          end if
c           ierr=124
c           herr = '[SATT error 124] density roots equal'
c           call ERRMSG (ierr,herr)
c           RETURN
c  !debug--next six lines for debug only
c         call GIBBS (t,rhol,x,Aliq,Gliq)
c         call GIBBS (t,rhov,x,Avap,Gvap)
c         ZG=Gliq-Gvap
c         write (*,1015) it,t,rhol,rhov,p,Gliq,Gvap,ZG
c1015   format (1x,'SATT:  it,t,rhol,rhov,p,Gliq,Gvap,ZG: ',
c    &              i4,f8.3,2f12.8,e14.6,2f12.4,e14.6)
          ierr=0
          herr=hnull
          RETURN           !normal termination for pure fluid
        end if
        delp2=delp
c
c  continue iteration, define next guess (check that delp
c  will not result in negative [or very small] pressure)
c  02-09-98 MM:  limit step size
c       if (delp.gt.0.9d0*p) then
        if (ABS(delp).gt.0.4d0*p) then
          do 198 j=1,10
c         write (*,1198) p,delp
c1198     format (1x,'% SATT advisory; delp gives p<0; p,delp:',2e12.4)
          delp=0.5*delp
          if (ABS(delp).lt.0.4d0*p) goto 199
  198     continue
        end if
  199   p=p-delp
c
  200   continue
c  iteration has not converged
        ierr=124
        write (herr,1124) t,Gliq,Gvap,ZG,delp,hnull
 1124   format ('[SATT error 124] ',
     &          'iteration for saturation state did not converge; ',
     &          'T =',1pe11.3,' K; Gibbs-liq & vap =',2(1pe11.3),
     &          ' J/mol; deltaG =',1pe11.3,' J/mol; deltaP =',1pe11.3,
     &          ' kPa.',a1)
        call ERRMSG (ierr,herr)
c  return critical parameters if not converged and very close to Tc
        if (t.gt.0.999975*tc) then
          p=pc
          rhol=rhoc
          rhov=rhoc
        end if
        RETURN
c
c  end of pure fluid iteratin
c
      else
c
c  begin mixture iteration
c
c  Iterate for the pressure and the composition of the incipient phase
c  (vapor phase for a bubble point calculation, liquid for dew point)
c  given temperature and the composition of the parent phase.  Iteration
c  is generally based on the algorithm given by Smith & Van Ness (Intro
c  to Chem Engr Thermo, McGraw-Hill, 1975); convergence criteria is the
c  equality of fugacity for each component in both phases.
c
        x2sum=0.0d0   !initialize only
c
c  generate initial guesses for mixture pressure and compositions
c  using Raoult's law
c
c  estimates for vapor pressures of pure components
        do 401 i=1,nc
        pcomp(i)=pcrit(i)
     &          *10.0**(-2.333333*(1.0+accen(i))*(tcrit(i)/t-1.0))
  401   continue
c
        psum=0.0d0
        if (kph.eq.1) then
c  bubble point
          kph2=2
c  variable kph2 specifies the state of the incipient phase (x2):
c     1 = liq,  2 = vap
c  it is used in calls to TPRHO
          pratio=1.05d0
c  pratio is multiplier for pressure to use when TPRHO does not converge
c  different values for liquid and vapor phases, such that new guess for
c  pressure is further into corresponding single-phase region
          do 411 i=1,nc
          psum=psum+x(i)*pcomp(i)
  411     continue
          do 412 i=1,nc
          x2(i)=x(i)*pcomp(i)/psum
  412     continue
        else
c  dew point
          kph2=1
          pratio=0.95d0
          ysum=0.0d0
          do 421 i=1,nc
          ysum=ysum+x(i)/pcomp(i)
  421     continue
          do 422 i=1,nc
          x2(i)=x(i)/pcomp(i)/ysum
          psum=psum+x2(i)*pcomp(i)
  422     continue
        end if
        if (psum.gt.pc) psum = 0.99d0*pc  !helps critical region conv.
c
c  begin main mixture iteration--outer loop for pressure,
c  using Raoult's Law result from above as first guess
c
        lppos=.false.         !flags for reguli-falsi iteration
        lpneg=.false.
        pneg=0.0d0
        ppos=0.0d0
        fpneg=0.0d0
        fppos=0.0d0
        j=1
        pit(1)=psum           !first guess for pressure = sum (x1*Pi)
        kguess=0              !for first calls to TPRHO
c       write (*,*) ' SATT--begin outer iteration loop for pressure'
        do 800 itp=1,itmax
        p=pit(j)
c       write (*,*) ' SATT--pressure iteration',itp,' w/ p =',pit(j)
        lx2con=.false.        !flag for convergence of inner loop
c  compute density and fugacities for parent phase
        call TPRHO (t,p,x,kph,kguess,rho1,ierr,herr2)   !parent phase
        if (ierr.gt.0) then
          ierr=-125
          write (herr,1125) itp,herr2(1:149),hnull
 1125     format ('[SATT advisory -125] density iteration in ',
     &            'saturation routine did not converge for pressure ',
     &            'iteration',i3,'; ',a149,a1)
          call ERRMSG (ierr,herr)
          pit(j)=pit(j)*pratio !try another pressure and use up one
          goto 800             !iteration (to prevent infinite loop)
        end if
        call FGCTY (t,rho1,x,f1)
c       write (*,1082) itp,t,p,rho1,(x(i),i=1,2),(f1(i),i=1,2)
c1082   format (1x,' SATT--phase 1:  ',i3,f8.2,2e14.6,2e16.8,2e18.10)
c
c  begin inner iteration loop for composition of phase 2
c
c       write (*,*) 'SATT--begin inner loop for composition of phase 2'
        do 600 itx=1,itmax
c  compute density and fugacities for phase 2
        call TPRHO (t,p,x2,kph2,kguess,rho2,ierr,herr2) !incipient phase
        if (ierr.gt.0) then
          ierr=-126
          write (herr,1126) itx,herr2(1:146),hnull
 1126     format ('[SATT advisory -126] density iteration in ',
     &            'saturation routine did not converge for composition',
     &            ' iteration',i3,'; ',a146,a1)
          call ERRMSG (ierr,herr)
          pit(j)=pit(j)/pratio !try another pressure and use up one
          goto 800             !iteration (to prevent infinite loop)
        end if
        call FGCTY (t,rho2,x2,f2)
c       write (*,1086) itx,rho2,(x2(i),i=1,2),(f2(i),i=1,2)
c1086   format (1x,'       phase 2:  ',i3,22x,e14.6,2e16.8,2e18.10)
c  calculate new x2's by ratio of fugacities; inner loop has converged
c  when x2's change by less than a convergence tolerance
        x2sum=0.0d0
        do 540 i=1,nc
        if (f2(i).gt.0.0d0) then
          x2new(i)=x2(i)*f1(i)/f2(i)
        else
c  in case fugacity is zero; e.g., if x(i)=0 x2(i) does not change
          x2new(i)=x2(i)
        end if
        x2sum=x2sum+x2new(i)
  540   continue
c  normalize the x2 compositions; this yields next guess for x2 and
c  ensures that the x2 always sum to one
        sumdel=0.0d0
        do 560 i=1,nc
        x2new(i)=x2new(i)/x2sum
        sumdel=sumdel+abs(x2(i)-x2new(i))  !change in compositions
        x2(i)=x2new(i)
  560   continue
        if (sumdel.lt.tolr) then
c  inner iteration loop has converged
          lx2con=.true.
          ierr=0
          goto 640
        end if
c  if not, continue inner iteration loop
  600   continue
c  inner iteration loop has not converged
        ierr=-127
        write (herr,1127) t,sumdel,hnull
 1127   format ('[SATT advisory -127] ',
     &          'iteration for composition in saturation routine ',
     &          'did not converge; T =',1pe11.3,
     &          ' K; deltaX =',1pe11.3,' mol frac.',a1)
        call ERRMSG (ierr,herr)
c
c  end of inner (x2) iteration loop
c
  640   continue
        fpit(j)=1.0d0-x2sum
c  outer (pressure) loop has converged when the x2's sum to one, i.e.,
c  when the fugacities of each component in each phase are equal
c       write (*,*) ' SATT--check conv, p, fp: ',pit(j),fpit(j) 
        if (ABS(fpit(j)).lt.tolr) then
          ierr=0
          herr=hnull
          goto 840
        else
c  provided that the inner loop has converged, update positive and
c  negative bounds on pressure for possible use in reguli-falsi iteration
          if (lx2con) then
            if (fpit(j).lt.0.0d0) then
              lpneg=.true.
              pneg=pit(j)
              fpneg=fpit(j)
            else
              lppos=.true.
              ppos=pit(j)
              fppos=fpit(j)
            end if
          end if
        end if
c
c  compute new guess for saturation pressure
c
        if (j.eq.1) then
c  for first iteration, new pressure is ratio of old
          j=2
          kguess=1    !use previous density as initial guess to TPRHO
          if (kph.eq.1) then
c  bubble point
            pit(2)=pit(1)*x2sum
          else
c  dew point
            pit(2)=pit(1)/x2sum
          end if
        else
c  subsequent iterations--use secant method
          pit(3)=pit(2)-fpit(2)*(pit(2)-pit(1))/(fpit(2)-fpit(1))
c  check that new pressure is not outside bounds, if so use reguli-falsi
          if (lpneg .and. lppos .and. (pit(3).gt.MAX(ppos,pneg)
     &        .or. pit(3).lt.MIN(ppos,pneg))) then
            pit(3)=ppos-fppos*(ppos-pneg)/(fppos-fpneg)
          end if
c  check that new pressure is below critical pressure
          if (pit(3).gt.pc) then
            pit(3)=0.5d0*(pit(2)+pc)
c           write (*,*) ' SATT--new p > p_crit; will reset'
          end if
c  discard oldest iteration
          pit(1)=pit(2)
          pit(2)=pit(3)
          fpit(1)=fpit(2)
        end if
c       write (*,1799) itp,j,pit(1),pit(2),fpit(2)
c1799   format (1x,' SATT--itp,j,p1,p2,fp2:  ',2i4,3e14.6)
  800   continue
c  outer iteration loop has not converged
        ierr=128
        write (herr,1128) t,sumdel,fpit(j),hnull
 1128   format ('[SATT error 128] ',
     &          'iteration for saturation state did not converge; T =',
     &          1pe11.3,' K; deltaX (inner & outer loops) =',2(1pe11.3),
     &          ' mol frac; ',a1)
        call ERRMSG (ierr,herr)
c
c  end of outer (pressure) iteration loop
c
  840   continue
c
c  assign final compositions and densities for parent and incipient
c  phases (x and x2, rho1 and rho2, respectively) to outputs
c
        if (kph.eq.1) then
c  bubble point
          rhol=rho1
          rhov=rho2
          do 481 i=1,nc
          xliq(i)=x(i)
          xvap(i)=x2(i)
  481     continue
        else
c  dew point
          rhol=rho2
          rhov=rho1
          do 482 i=1,nc
          xliq(i)=x2(i)
          xvap(i)=x(i)
  482     continue
        end if
c
c  end of mixture iteration
      end if
      RETURN
c
      end                                               !subroutine SATT
c
c ======================================================================
c
      subroutine SATP (p,x,kph,t,rhol,rhov,xliq,xvap,ierr,herr)
c
c  iterate for saturated liquid and vapor states given temperature
c  and the composition of one phase
c
c  inputs:
c        p--pressure [kPa]
c        x--composition [array of mol frac] (phase specified by kph)
c      kph--phase flag:  1 = input x is liquid composition
c                        2 = input x is vapor composition
c
c  outputs:
c        t--temperature [K]
c     rhol--molar density [mol/L] of saturated liquid
c     rhov--molar density [mol/L] of saturated vapor
c     xliq--liquid phase composition [array of mol frac]
c     xvap--vapor phase composition [array of mol frac]
c     ierr--error flag:  0 = successful
c                        4 = P < 0
c                        8 = x out of range
c                       12 = P and x out of range
c                      140 = CRITP did not converge
c                      141 = P > Pcrit
c                      142 = TPRHO-liquid did not converge (pure fluid)
c                      143 = TPRHO-vapor did not converge (pure fluid)
c                      144 = pure fluid iteration did not converge
c           following 3 error codes are advisory--iteration will either
c           converge on later guess or error out (ierr = 148)
c                     -144 = Raoult's law (mixture intial guess) did
c                            not converge
c                     -145 = TPRHO did not converge for parent ph (mix)
c                     -146 = TPRHO did not converge for incipient (mix)
c                     -147 = composition iteration did not converge
c                      148 = mixture iteration did not converge
c     herr--error string if ierr<>0 (character*255)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  07-13-95  MM, original version
c  09-11-95  MM, add error string to argument list
c  09-25-95  MM, rearrange argument list (outputs in order t, rho, x)
c  10-11-95  MM, RETURN if any error detected
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accomodate ECS reference fluid
c  12-11-95  MM, Raoult's law as first guess for mixture
c  12-18-95  MM, fill xliq, xvap with zeros for undefined components
c  12-27-95  MM, add full mixture iteration using fugacity, based on SATT
c  02-27-96  MM, parameter n0=-ncmax to accomodate ECS-thermo model
c                add Zcrit to common /CCON/
c  03-19-96  MM, add dipole moment to /CCON/
c  04-05-96  MM, test for supercritical '.ge. pc' rather than '.gt. pc'
c  11-14-96  MM, adjust initial guesses, tratio, etc to get closer to critical
c   2-12-96 EWL, special initial guess for temperature near the critical point
c  07-15-97  MM, add errors/warnings to parallel SATT
c  10-01-97  MM, add compiler switches to allow access by DLL
c  11-13-97 EWL, initilize ft(j); potential bomb if no value when writing error message
c  12-05-97  MM, check that TPRHO gives density within bounds
c                if Raoult's law iteration D.N.C., revert to initial guess
c  02-10-98  MM, add reguli-falsi, quadratic interpolation and bisection
c                to Raoult's law iteration
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 :: SATP
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
      parameter (n0=-ncmax,nx=ncmax)
      character*1 htab,hnull
      character*255 herr,herr2
      dimension x(ncmax),xliq(ncmax),xvap(ncmax)
      dimension pcomp(ncmax),x2(ncmax),f1(ncmax),f2(ncmax),x2new(ncmax)
      dimension tk(4),ft(3)                  !used for T iteration
      common /Gcnst/ R
      common /NCOMP/ nc
      common /CHAR/ htab,hnull
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tcrit(n0:nx),pcrit(n0:nx),Dcrit(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(n0:nx)
c  !temporary--convergence tolerance, put in common block
      data tolr /1.0d-6/
      data itmax /40/
c  !temporary--end
      ierr=0
      herr=hnull
      ft(1)=0.0d0  !initialize to avoid potential problem
      ft(2)=0.0d0  !when writing error message to GUI
c
c  initialize outputs in event of failure of routines
      t=300.0d0
      rhol=0.0d0
      rhov=0.0d0
      do 100 i=1,nc
      xliq(i)=x(i)
      xvap(i)=x(i)
  100 continue
c  fill output composition arrays with zeros for undefined components
      if (nc.lt.ncmax) then
        do 104 i=nc+1,ncmax
        xliq(i)=0.0d0
        xvap(i)=0.0d0
  104   continue
      end if
c
      call CRITP (x,tc,pc,rhoc,ierr,herr2)
      if (ierr.ne.0) then
c  error condition--set outputs, issue warning, and return
        ierr=140
        write (herr,1140) ierr,herr2(1:237),hnull
 1140   format ('[SATP error 140] ',a237,a1)
        call ERRMSG (ierr,herr2)
        RETURN
      end if
c
c  check that input conditions (in this case p and x) are within limits
c
      Ddum=0.0d0
      tdum=0.8d0*tc
      call LIMITX ('EOS',tdum,Ddum,p,x,tmin,tmax,rhomax,pmax,ierr,herr2)
c     write (*,*) ' SATP--density limit:  ',rhomax
      if (ierr.gt.1) then     !ignore ierr = 1 (t out of range)
c  p and/or x are out of bounds, set error flag and return
        t=0.8d0*tc
        write (herr,1000) ierr,herr2(1:238),hnull
 1000   format ('[SATP error',i3,'] ',a238,a1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
c
      if (p/pc.ge.0.999999d0) then
c  input pressure is equal to or greater than critical point value,
c  return critical parameters for output temperature and densities
c  output compositions initialized above
        t=tc
        rhol=rhoc
        rhov=rhoc
        if (p/pc.gt.1.000001d0) then
c  supercritical pressure as input, set error flag and return
c  critical parameters for output pressure and densities
          ierr=141
          write (herr,1141) p,pc,hnull
 1141     format ('[SATP error 141] ',
     &            'pressure input to saturation routine is ',
     &            'greater than critical pressure; P =',1pe11.3,
     &            ' kPa, Pcrit =',1pe11.3,' kPa.',a1)
        end if
        RETURN
      end if
c
      if (nc.eq.1) then
c
c  pure fluid iteration
        xliq(1)=x(1)
        xvap(1)=x(1)
c
c  Iterate for saturated liquid and vapor states given pressure using
c  a simple successive substitution method.  The independent variable
c  in the iteration is the temperature.  The convergence criteria is
c  equality of Gibbs free energy in both phases.
c
c  generate initial guess using accentric factor
c
        t=tc/(1.0-0.428571*log10(p/pc)/(1.0+accen(1)))
c       write (*,*) 'SATP:  T(0): ',t
c
c  assume nothing about densities on initial calls to TPRHO
        kguess=0
        if (p.gt.0.98*pc) then
c  for near-critical states, generate initial guesses for density;
c  using correlation developed by E.W. Lemmon, NIST
          theta=(1-p/pc)*(1.6d0-accen(1))
          t=tc*(1.0d0-0.103947d0*theta-4.108265d-2*theta**2)
          theta=(1-p/pc)**(1.0d0/3.0d0)*(3.0d0+accen(1))
          rhov=rhoc*(1.0d0-0.290039d0*theta-7.120197d-3*theta**2)
          rhol=rhoc*(1.0d0+0.298544d0*theta+1.870808d-2*theta**2)
          kguess = 1
        endif
        do 200 it=1,itmax
        call TPRHO (t,p,x,1,kguess,rhol,ierr,herr2) !find liquid density
        if (ierr.ne.0 .or. rhol.lt.rhoc) then
          t=t/1.0002d0
          goto 200
c         herr=' ERROR from SATP:  '//herr2
c         call ERRMSG (ierr,herr)
c         RETURN
        end if
        call TPRHO (t,p,x,2,kguess,rhov,ierr,herr2)  !find vapor density
        if (ierr.ne.0 .or. rhov.gt.rhoc) then
          t=t*1.001d0
          goto 200
c         herr=' ERROR from SATP:  '//herr2
c         call ERRMSG (ierr,herr)
c         RETURN
        end if
c  use previous densities as initial guesses for calls to TPRHO after
c  first iteration
        kguess=1
        call GIBBS (t,rhol,x,Aliq,Gliq)
        call GIBBS (t,rhov,x,Avap,Gvap)
        call ENTRO (t,rhol,x,sliq)
        call ENTRO (t,rhov,x,svap)
        ZG=Gliq-Gvap
c       write (*,1014) it,p,rhol,rhov,t,Gliq,Gvap,ZG
c1014   format (1x,'SATP:',i4,e14.6,2f12.8,f10.5,2f12.4,e14.6)
c
c  check convergence
c
        delt=ZG/(sliq-svap)
        if (abs(delt) .lt. tolr) then
c  pure component iteration is done (make use of current delt)
          t=t+delt
          call TPRHO (t,p,x,1,kguess,rhol,ierr,herr2) !find liq density
          if (ierr.ne.0 .or. t.gt.tc .or. rhol.lt.rhoc) then
            ierr=142
            write (herr,1142) it,herr2(1:142),hnull
 1142       format ('[SATP error 142] liquid density iteration in ',
     &            'saturation routine did not converge for temperature',
     &            ' iteration',i3,'; ',a142,a1)
            call ERRMSG (ierr,herr)
c  return critical parameters if not converged and very close to Pc
            if (p.gt.0.9999*pc) then
              t=tc
              rhol=rhoc
              rhov=rhoc
            end if
            RETURN
          end if
          call TPRHO (t,p,x,2,kguess,rhov,ierr,herr2) !find vap density
          if (ierr.ne.0 .or. p.gt.pc .or. rhov.gt.rhoc) then
            ierr=143
            write (herr,1143) it,herr2(1:143),hnull
 1143       format ('[SATP error 143] vapor density iteration in ',
     &            'saturation routine did not converge for temperature',
     &            ' iteration',i3,'; ',a143,a1)
            call ERRMSG (ierr,herr)
c  return critical parameters if not converged and very close to Pc
            if (p.gt.0.9999*pc) then
              t=tc
              rhol=rhoc
              rhov=rhoc
            end if
            RETURN
          end if
c  !debug--next four lines for debug only
c         call GIBBS (t,rhol,x,Aliq,Gliq)
c         call GIBBS (t,rhov,x,Avap,Gvap)
c         ZG=Gliq-Gvap
c         write (*,1015) it,p,rhol,rhov,t,Gliq,Gvap,ZG
c1015   format (1x,'SATP:',i4,e14.6,2f12.8,f10.5,2f12.4,e14.6)
          xliq(1)=x(1)
          xvap(1)=x(1)
          ierr=0
          herr=hnull
          RETURN           !normal termination for pure fluid
        end if
c
c  continue iteration, define next guess (check that delt
c  will not result in too large a change in temperature)
        if (delt.gt.0.5*t) then
          do 198 j=1,100
c         write (*,1198) t,delt
c1198     format (1x,'% SATP advisory; delt > 0.5*t; t,delt:',2f12.6)
          delt=0.25*delt
          if (delt.lt.0.5*t) goto 199
  198     continue
        end if
  199   t=t+delt
c
  200   continue
c  iteration has not converged
        ierr=144
        write (herr,1144) p,Gliq,Gvap,ZG,delt,hnull
 1144   format ('[SATP error 144] ',
     &          'iteration for saturation state did not converge; ',
     &          'P =',1pe11.3,' kPa; Gibbs-liq & vap =',2(1pe11.3),
     &          ' J/mol; deltaG =',1pe11.3,' J/mol; deltaT =',1pe11.3,
     &          ' K.',a1)
        call ERRMSG (ierr,herr)
c  return critical parameters if not converged and very close to Tc
        if (p.gt.0.9999*pc) then
          t=tc
          rhol=rhoc
          rhov=rhoc
        end if
        RETURN
c
c  end of pure fluid iteratin
c
      else
c
c  begin mixture iteration
c
c  Iterate for the temperature and the composition of the incipient
c  phase (vapor phase for a bubble point calculation, liquid for dew
c  point) given pressure and the composition of the parent phase.
c  Iteration is generally based on the algorithm given by Van Ness
c  & Abbott (Classical Thermodynamics of Nonelectrolyte Solutions with
c  Applications to Phase Equilibria, McGraw-Hill, 1982); convergence
c  criteria is the equality of fugacity for each component in both
c  phases.
c
        ltpos=.false. !flags for reguli-falsi iteration
        ltneg=.false.
        tneg=0.0d0
        tpos=0.0d0
        ftneg=0.0d0
        ftpos=0.0d0
        x2sum=0.0d0   !initialize only
c
c  generate initial guess using Raoult's law,
c  but first generate initial guess for temperature by simple ratio of
c  pure component saturation temperatures
        tk(1)=0.0d0
        acf=0.0d0              !average accentric factor
        do 401 i=1,nc
        ti=tcrit(i)/(1.0-0.428571*log10(p/pcrit(i))/(1.0+accen(i)))
c       write (*,*) ' SATP:  tsat-guess for component ',i,ti
        tk(1)=tk(1)+x(i)*ti
        acf=acf+x(i)*accen(i)
  401   continue
        tzero=tk(1)  !save in case Raoult iteration does not converge
c       write (*,*) ' SATP--initial guess for tsat: ',tk(1)
c
c  secant method iteration to find t which satisfies Raoult's law
c
        j=1                     !flag for first iteration
        do 420 it=1,itmax
c  approximate pure component vapor pressures with accentric factor
        do 410 i=1,nc
        pcomp(i)=pcrit(i)
     &          *10.0**(-2.333333*(1.0+accen(i))*(tcrit(i)/tk(j)-1.0))
  410   continue
c
        if (kph.eq.1) then
c  bubble point
          psum=0.0d0
          do 411 i=1,nc
          psum=psum+x(i)*pcomp(i)
  411     continue
          ft(j)=1.0d0-p/psum
          if (abs(ft(j)).lt.tolr) then
c  iteration has converged
            goto 422
          else
c  update + & - bounds on temperature for possible use in reguli-falsi
            if (ft(j).lt.0.0d0) then
              ltneg=.true.
              tneg=tk(j)
              ftneg=ft(j)
            else
              ltpos=.true.
              tpos=tk(j)
              ftpos=ft(j)
            end if
c  generate next guess
            if (it.eq.1) then
              tratio=1.0-0.42857*log10(psum/p)/(1.0d0+acf)
              tk(2)=tk(1)*tratio
              j=2
            else if (it.ge.itmax/3 .and. ltpos .and. ltneg) then
c  if iteration has not converged after many iterations, use bisection
c  (provided that guesses bounding the root are available)
              tk(2)=0.5d0*(tpos+tneg)
              tk(1)=tk(2)
              ft(1)=ft(j)
              j=2
            else
c  use secant method
              if (ABS(ft(2)-ft(1)).gt.1.0d-10) then
                tk(3)=tk(2)-ft(2)*(tk(2)-tk(1))/(ft(2)-ft(1))
              else
                tk(3)=0.5d0*(tk(1)+tk(2))
              end if
c  check that new temperature is not outside bounds, if so use reguli-falsi
              if (ltneg .and. ltpos .and. (tk(j+1).gt.MAX(tpos,tneg)
     &            .or. tk(j+1).lt.MIN(tpos,tneg))) then
                tk(j+1)=tpos-ftpos*(tpos-tneg)/(ftpos-ftneg)
              end if
              tk(1)=tk(2)
              tk(2)=tk(3)
              ft(1)=ft(2)
            end if
          end if
        else
c  dew point
          xdamp=1.0d0     !damping ratio for secant method
          ypsum=0.0d0
          do 412 i=1,nc
          ypsum=ypsum+x(i)*p/pcomp(i)
          ft(j)=1.0-ypsum
  412     continue
          if (ABS(ft(j)).lt.1.0d3*tolr) then
c  iteration has converged (initial guess only, tolerance is relaxed)
            goto 422
          else
c  update + & - bounds on pressure for possible use in reguli-falsi
            if (ft(j).lt.0.0d0) then
              ltneg=.true.
              tneg=tk(j)
              ftneg=ft(j)
            else
              ltpos=.true.
              tpos=tk(j)
              ftpos=ft(j)
            end if
            if (it.eq.1) then
              tratio=1.0d0/(1.0-0.42857*log10(ypsum)/(1.0d0+acf))
              tk(2)=tk(1)*tratio
              j=2
            else if (it.eq.2) then
c  secant method for 2nd guesse
                if (ABS(ft(2)-ft(1)).gt.1.0d-10) then
                  tk(3)=tk(2)-xdamp*ft(2)*(tk(2)-tk(1))/(ft(2)-ft(1))
                else
                  tk(3)=0.5d0*(tk(1)+tk(2))
                end if
c  check that new temperature is not outside bounds, if so use reguli-falsi
                if (ltneg .and. ltpos .and. (tk(j+1).gt.MAX(tpos,tneg)
     &              .or. tk(j+1).lt.MIN(tpos,tneg))) then
                  tk(j+1)=tpos-ftpos*(tpos-tneg)/(ftpos-ftneg)
                end if
                j=3
            else if (it.le.itmax/3 .or. .not.(ltpos.and.ltneg)) then
c  2nd order secant (inverse quadratic interpolation) for subseqent guesses
c  see Numerical Recipes, p 252
              rr=ft(3)/ft(2)  !these are the R,S,T used in Num. Rec.
              rs=ft(3)/ft(1)
              rt=ft(1)/ft(2)
              tk(4)=tk(3)
     &             +rs*(rt*(rr-rt)*(tk(2)-tk(3))
     &             -(1.0d0-rr)*(tk(3)-tk(1)))
     &             /((rt-1.0d0)*(rr-1.0d0)*(rs-1.0d0))
c  check that new temperature is not outside bounds, if so use reguli-falsi
              if (ltneg .and. ltpos .and. (tk(j+1).gt.MAX(tpos,tneg)
     &            .or. tk(j+1).lt.MIN(tpos,tneg))) then
                tk(j+1)=tpos-ftpos*(tpos-tneg)/(ftpos-ftneg)
              end if
c  discard oldest iteration
              tk(1)=tk(2)
              tk(2)=tk(3)
              tk(3)=tk(4)
              ft(1)=ft(2)
              ft(2)=ft(3)
            else
c  if iteration has not converged by now, use bisection
              tk(2)=0.5d0*(tpos+tneg)
              tk(1)=tk(2)
              ft(1)=ft(j)
              j=2
            end if
          end if
        end if
c       write (*,1006) it,j,tk(j-1),tk(j),ft(j-1)
c1006   format (1x,' SATP Raoult''s: it,j,t1,t2,ft:',2i4,2f12.4,e16.7)
  420   continue                                  !next trial for t
c
c  iteration has not converged, issue warning and proceed
        ierr=-144
        herr='[SATP advisory -144] Raoult''s law iteration (to '//
     &       'generate mixture initial guess) has not converged.'//
     &       hnull
        call ERRMSG (ierr,herr)
        if (ABS(ft(j)).lt.1.0d4*tolr .or. 
     &      ABS(tk(j)-tk(j-1)).lt.1.0d0) then
c         write (*,*) ' SATP--Raoult iteration D.N.C.; tk(j) = ',tk(j)
        else
c  if current guess is not even close go back to initial guess
c         write (*,*) ' SATP--Raoult iteration D.N.C.; tk(j), tk_new = '
c    &                ,tk(j),tzero
          tk(j)=tzero
        end if
c
  422   continue
c
c  initial temperature (satisfying Raoult's law) has been found,
c  check that this temperature is not above critical
        tmax=0.998d0*tc
        if (tk(j).gt.tmax) then
c         write (*,1042) tk(j),tmax
c1042     format (1x,' SATP:  mix temp by Raoults law > Tc;',
c    &               ' Raoults law guess, new guess:  ',2f10.3)
          tk(j)=tmax
        end if
c
c  generate initial guesses for densities & incipient phase composition;
c  the do loop allows for the possibility that TPRHO does not converge,
c  it should normally exit with just one pass
c
        do 460 it=1,itmax
        t=tk(j)
        psum=0.0d0
        if (kph.eq.1) then
c  bubble point
          kph2=2
c  variable kph2 specifies the state of the incipient phase (x2):
c     1 = liq,  2 = vap
c  it is used in calls to TPRHO
          tratio=0.995d0
c  tratio is temperature multiplier to use when TPRHO does not converge;
c  different values for liquid and vapor phases, such that new guess
c  is further into corresponding single-phase region
          do 431 i=1,nc
          psum=psum+x(i)*pcomp(i)
  431     continue
          do 432 i=1,nc
          x2(i)=x(i)*pcomp(i)/psum
  432     continue
        else
c  dew point
          kph2=1
          tratio=1.005d0
          ysum=0.0d0
          do 441 i=1,nc
          ysum=ysum+x(i)/pcomp(i)
  441     continue
          do 442 i=1,nc
          x2(i)=x(i)/pcomp(i)/ysum
          psum=psum+x2(i)*pcomp(i)
  442     continue
        end if
c  first guess for densities
        kguess=0
c       write (*,*) ' SATP call TPRHO (parent) for it,T = ',it,t
        call TPRHO (t,p,x,kph,kguess,rho1,ierr,herr2)     !parent phase
        if (ierr.gt.0) then
          ierr=-145
          write (herr,1145) herr2(1:146),hnull
 1145     format ('[SATP advisory -145] density iteration in ',
     &            'saturation routine did not converge for the parent',
     &            ' phase; ',a146,a1)
          call ERRMSG (ierr,herr)
          if (t.gt.0.8*tc .and. kph.eq.1) then
c   non-convergence probably because too close to critical
            tnew=0.995*t
          else
            tnew=1.005*t
          end if
          tk(j)=tnew
        else
c         write (*,*) ' SATP call TPRHO (incip) for it,T =  ',it,t
          call TPRHO (t,p,x2,kph2,kguess,rho2,ierr,herr2) !incipient ph
          if (ierr.gt.0) then
            ierr=-144
            write (herr,1146) herr2(1:146),hnull
 1146       format ('[SATP advisory -146] density iteration in ',
     &            'saturation routine did not converge for the ',
     &            'incipient phase; ',a146,a1)
            call ERRMSG (ierr,herr)
            if (t.gt.0.8*tc .and. kph2.eq.1) then
c   non-convergence probably because too close to critical
              tnew=0.995*t
            else
              tnew=1.005*t
            end if
c           write (*,1047) kph2,t,tnew
c1047       format (1x,' SATP--TPRHO (initial incipient-ph density)',
c    &                 ' did not converge; kph,T,Tnew:  ',i3,2f10.3)
            tk(j)=tnew
          else
            goto 461   !both parent and incipient phases have converged
          end if
        end if
  460   continue
c
c  now, finally, proceed with main iteration loop for mixtures, outer
c  loop for temperature, using Raoult's Law result (above) as first guess
c
  461   ltpos=.false.         !flags for reguli-falsi iteration
        ltneg=.false.
        tneg=0.0d0
        tpos=0.0d0
        ftneg=0.0d0
        ftpos=0.0d0
c       write (*,*) ' SATP--start main iteration; j,tk(j):  ',j,tk(j)
        tk(1)=tk(j)           !first guess for temperature from above
        j=1                   !reset iteration flag
        kguess=0              !for first calls to TPRHO
        do 800 itt=1,itmax
        t=tk(j)
        lx2con=.false.        !flag for convergence of inner loop
c  compute density and fugacities for parent phase
c       write (*,*) ' SATP call TPRHO (parent) for it,T = ',itt,t
        call TPRHO (t,p,x,kph,kguess,rho1,ierr,herr2)   !parent phase
        if (ierr.gt.0) then
          ierr=-145
          write (herr,1245) itt,herr2(1:146),hnull
 1245     format ('[SATP advisory -145] density iteration in ',
     &            'saturation routine did not converge for temperature',
     &            ' iteration',i3,'; ',a146,a1)
          call ERRMSG (ierr,herr)
          tk(j)=tk(j)*tratio   !try another temperature and use up one
          goto 800             !iteration (to prevent infinite loop)
        else if (rho1.gt.1.1d0*rhomax) then
c  density from TPRHO is out of range, reset
          rho1=rhomax
          ierr=-145
          write (herr,1245) itt,herr2(1:146),hnull
        end if
        call FGCTY (t,rho1,x,f1)
c       write (*,1082) itt,t,p,rho1,(x(i),i=1,2),(f1(i),i=1,2)
c1082   format (1x,' SATP--phase 1:  ',i3,f11.5,2e14.6,2e16.8,2e18.10)
c
c  begin inner iteration loop for composition of phase 2
c
        do 600 itx=1,itmax
c  compute density and fugacities for phase 2
        call TPRHO (t,p,x2,kph2,kguess,rho2,ierr,herr2) !incipient phase
        if (ierr.gt.0) then
          ierr=-146
          write (herr,1246) itx,herr2(1:146),hnull
 1246     format ('[SATP advisory -146] density iteration in ',
     &            'saturation routine did not converge for composition',
     &            ' iteration',i3,'; ',a146,a1)
          call ERRMSG (ierr,herr)
          tk(j)=tk(j)/tratio   !try another pressure and use up one
          goto 800             !iteration (to prevent infinite loop)
        else if (rho2.gt.1.2*rhomax) then
c  density from TPRHO is out of range, reset
          rho2=rhomax
          ierr=-146
          write (herr,1246) itt,herr2(1:146),hnull
        end if
        call FGCTY (t,rho2,x2,f2)
c       write (*,1086) itx,rho2,(x2(i),i=1,2),(f2(i),i=1,2)
c1086   format (1x,'       phase 2:  ',i3,25x,e14.6,2e16.8,2e18.10)
c  calculate new x2's by ratio of fugacities; inner loop has converged
c  when x2's change by less than a convergence tolerance
        x2sum=0.0d0
        do 540 i=1,nc
        if (f2(i).gt.1.0d-20) then
          x2new(i)=x2(i)*f1(i)/f2(i)
        else
c  in case fugacity is zero; e.g., if x(i)=0 x2(i) does not change
          x2new(i)=x2(i)
        end if
        x2sum=x2sum+x2new(i)
  540   continue
c  normalize the x2 compositions; this yields next guess for x2 and
c  ensures that the x2 always sum to one
        sumdel=0.0d0
        do 560 i=1,nc
        x2new(i)=x2new(i)/x2sum
        sumdel=sumdel+abs(x2(i)-x2new(i))  !change in compositions
        x2(i)=x2new(i)
  560   continue
        if (sumdel.lt.tolr) then
c  inner iteration loop has converged
          lx2con=.true.
          goto 640
        end if
c  if not, continue inner iteration loop
  600   continue
c  inner iteration loop has not converged
        ierr=-147
        write (herr,1147) p,sumdel,hnull
 1147   format ('[SATP advisory -147] ',
     &          'iteration for composition in saturation routine ',
     &          'did not converge; P =',1pe11.3,
     &          ' kPa; deltaX =',1pe11.3,' mol frac.',a1)
        call ERRMSG (ierr,herr)
c
c  end of inner (x2) iteration loop
c
  640   continue
        ft(j)=1.0d0-x2sum
c  outer (temperature) loop has converged when the x2's sum to one, i.e.,
c  when the fugacities of each component in each phase are equal
        if (abs(ft(j)).lt.tolr) then
          goto 840
        else
c  provided that the inner loop has converged, update positive and
c  negative bounds on pressure for possible use in reguli-falsi iteration
          if (lx2con) then
            if (ft(j).lt.0.0d0) then
              ltneg=.true.
              tneg=tk(j)
              ftneg=ft(j)
            else
              ltpos=.true.
              tpos=tk(j)
              ftpos=ft(j)
            end if
          end if
        end if
c
c  compute new guess for saturation temperature
c
        if (j.eq.1) then
c  for first iteration, new temperature is ratio of old
          j=2
          kguess=1    !use previous density as initial guess to TPRHO
c  ratio for next guess of temperature; the 0.10d0 is adjustable
          tnew=1.0d0+0.10d0*(x2sum-1.0d0)
c         write (*,*) ' PSAT--x2sum,tnew:  ',x2sum,tnew
          if (kph.eq.1) then
c  bubble point
            tk(2)=tk(1)/tnew
          else
c  dew point
            tk(2)=tk(1)*tnew
          end if
        else
c  subsequent iterations--use secant method, check for divide by zero
          if (ABS(ft(2)-ft(1)).lt.1.0d-10) then
            tk(3)=0.5*(tk(1)+tk(2))
          else
            tk(3)=tk(2)-ft(2)*(tk(2)-tk(1))/(ft(2)-ft(1))
          end if
c  check that new temperature is not outside bounds, if so use reguli-falsi
          if (ltneg .and. ltpos .and. (tk(3).gt.MAX(tpos,tneg)
     &        .or. tk(3).lt.MIN(tpos,tneg))) then
            tk(3)=tpos-ftpos*(tpos-tneg)/(ftpos-ftneg)
          end if
c         write (*,1137) tk(1),tk(2),tk(3),ft(1),ft(2)
c1137     format (1x,' SATP--tguess_1,2,3; ft_1,2:  ',5e14.6)
c  discard oldest iteration
          tk(1)=tk(2)
          tk(2)=tk(3)
          ft(1)=ft(2)
        end if
  800   continue
c  outer iteration loop has not converged
        ierr=148
        write (herr,1148) p,sumdel,ft(j),hnull
 1148   format ('[SATP error 148] ',
     &          'iteration for saturation state did not converge; P =',
     &          1pe11.3,' kPa; deltaX (inner & outer loops) =',
     &          2(1pe11.3),' mol frac; ',a1)
        call ERRMSG (ierr,herr)
c
c  end of outer (temperature) iteration loop
c
  840   continue
c
c  assign final compositions and densities for parent and incipient
c  phases (x and x2, rho1 and rho2, respectively) to outputs
c
        t=tk(j)
        if (kph.eq.1) then
c  bubble point
          rhol=rho1
          rhov=rho2
          do 881 i=1,nc
          xliq(i)=x(i)
          xvap(i)=x2(i)
  881     continue
        else
c  dew point
          rhol=rho2
          rhov=rho1
          do 882 i=1,nc
          xliq(i)=x2(i)
          xvap(i)=x(i)
  882     continue
        end if
c
c  end of mixture iteration
      end if
      RETURN
c
      end                                               !subroutine SATP
c
c ======================================================================
c
      subroutine SURFT (t,rhol,xl,sigma,ierr,herr)
c
c  compute surface tension
c
c  inputs:
c        t--temperature [K]
c     rhol--molar density of liquid phase [mol/L]
c           if rho > 0 use as input value
c                  < 0 call SATT to find density
c       xl--composition of liquid phase [array of mol frac]
c  outputs:
c    sigma--surface tension [N/m]
c     ierr--error flag:   0 = successful
c                         1 = T < Tmin
c                         8 = x out of range
c                         9 = T and x out of range
c                       120 = CRITP did not converge
c                       121 = T > Tcrit
c                       122 = TPRHO-liquid did not converge in SATT
c                       123 = TPRHO-vapor did not converge in SATT
c                       124 = SATT pure fluid iteration did not converge
c                       128 = SATT mixture iteration did not converge
c     herr--error string if ierr<>0 (character*255)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  03-24-96  MM, original version
c  03-27-96  MM, add error checks; move calculations to STN (in core_STN)
c  04-05-96  MM, test for supercritical '.ge. tc' rather than '.gt. tc'
c  05-31-96  MM, if error on call to SATT modify herr and return
c  06-03-96  MM, check input temperature against limits
c  06-07-96  MM, fix loss of LIMITX warnings on calls to SATT, STN
c  04-21-97  MM, delete tcrit from call to STN (t > tcrit check moved to STN)
c  10-01-97  MM, add compiler switches to allow access by DLL
c  12-15-97  MM, pass any ierr,herr from STN as outputs
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 :: SURFT
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      parameter (ncmax=5)        !max number of components in mixture
      character*1 htab,hnull
      character*255 herr,herr2
      common /CHAR/ htab,hnull
      dimension xl(ncmax),xliq(ncmax),xvap(ncmax)
c
      ierr=0
      herr=hnull
c
c  check that input conditions (in this case t and x) are within limits
c  (check that t < tcrit done in STN)
c
      Ddum=0.0d0
      pdum=0.0d0
      call LIMITX ('STN',t,Ddum,pdum,xl,tmin,tmax,Dmx,pmx,ierr,herr2)
      if (ierr.lt.0) then
c  temperature is outside limits, but in region where extrapolation is
c  usually reliable, set warning flag
        ierr=ierr-20
        write (herr,1002) ierr,herr2(1:234),hnull
 1002   format ('[SURFT warning',i4,'] ',a234,a1)
        call ERRMSG (ierr,herr)
      else if (ierr.gt.0) then
c  temperature and/or x are outside limits, set error flag
        write (herr,1003) ierr,herr2(1:237),hnull
 1003   format ('[SURFT error',i3,'] ',a237,a1)
        call ERRMSG (ierr,herr)
      end if
c
c  calculate density of saturated liquid, if required
      if (rhol.le.0.0d0) then
        kph=1
        call SATT (t,xl,kph,p,rhol,rhov,xliq,xvap,ierr2,herr2)
        if (ierr2.ne.0) then
          write (herr,1005) ierr2,herr2(1:237),hnull
 1005     format ('[SURFT error',i3,'] ',a237,a1)
          call ERRMSG (ierr,herr)
          if (ierr2.gt.0) then
            sigma=0.0d0
            RETURN
          end if
        end if
      end if
      call STN (t,rhol,xl,sigma,ierr,herr)
c
      RETURN
      end                                              !subroutine SURFT
c
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                     end file sat_sub.f
c ======================================================================