c  begin file mix_HMX.f
c
c  This file contains the routines implementing the mixture Helmholtz
c  equation of state of Lemmon and Jacobsen and of Kunz and Wagner.
c
c  contained here are:
c     function PHIHMX (itau,idel,tau,del,x)
c     function PHIMIX (jcomp,x,itau,idel,tau,del,x)
c     subroutine CRTHMX (x,tcrit,pcrit,Dcrit,ierr,herr)
c     subroutine REDHMX (x,t0,D0)
c     subroutine SETHMX (hfmix,ierr,herr)
c     subroutine RDBNC (nread,icomp,jcomp,ibin,ierr,herr)
c     subroutine RDMIX (nread,icomp,jcomp,hmodij,lij,ierr,herr)
c     subroutine ESTBNC (i,j,ibin,ierr,herr)
c     subroutine ESTTC (i,j,ierr,herr)
c     subroutine ESTVC (i,j,ierr,herr)
c
c ======================================================================
c ======================================================================
c
      function PHIHMX (itau,idel,tau,del,x)
c
c  compute reduced Helmholtz energy or a derivative as functions
c  of dimensionless temperature and density for the mixture Helmholtz
c  equation of state
c
c  based on notes of E.W. Lemmon and R.T Jacobsen, U. Idaho, 1995
c
c  inputs:
c     itau--flag specifying order of temperature derivative to calc
c     idel--flag specifying order of density derivative to calculate
c           when itau = 0 and idel = 0, compute A/RT
c           when itau = 0 and idel = 1, compute 1st density derivative
c           when itau = 1 and idel = 1, compute cross derivative
c           etc.
c      tau--dimensionless temperature (To/T)
c      del--dimensionless density (D/Do)
c        x--composition array (mol frac)
c  output (as function value):
c      phi--residual (real-gas) part of the Helmholtz energy, or one
c           of its derivatives (as specified by itau and idel),
c           in reduced form (A/RT)
c
c  N.B.  The reducing parameters To and Do are often, but not
c        necessarily, equal to the critical temperature and density.
c
c        The Helmholtz energy consists of ideal gas and residual (real-
c        gas) terms.  The residual term consists of ideal-solution and
c        mixing terms.  This routine calculates only the residual term.
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-02-95  MM, original version
c  11-26-95  MM, rearrange argument list (x at end)
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  02-29-96  MM, replace calls to PHIBWR, PHIFEQ with general PHIK
c  03-21-96  MM, delete reference to /MODEL/, not used
c  07-01-96  MM, move x(i)*x(j) multiplier here (remove from PHIMIX)
c  06-03-97 EWL, add third derivative of Helmholtz energy with respect to density.
c  01-07-98  MM, move x(i)*x(j) multiplier back to PHIMIX
c  09-00-00 EWL, switch order of execution, check x=0 before call to PHIMIX
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
      parameter (ncmax=20)        !max number of components in mixture
      common /NCOMP/ nc,ic
      dimension x(ncmax)
c
c
      if (del.le.1.0d-10) then    !trivial solution at zero density
        phihmx=0.0d0              !for any and all derivatives
        RETURN
      end if
c
      phisum=0.0d0
      do i=1,nc
        if (x(i).gt.1.0d-12) then
c  calculate ideal solution term
          phiii=PHIK(i,itau,idel,tau,del)
          phisum=phisum+x(i)*phiii
c         write (*,1200) i,i,x(1),x(2),x(3),phiii,phisum
c1200     format (1x,' PHIHMX--i,j,x,phi_ij:  ',2i4,3f8.5,1pe14.6)
        end if
      enddo
      do i=1,nc-1
c  calculate mixing term
        do j=i+1,nc
          if (x(i).gt.1.0d-12 .and. x(j).gt.1.0d-12) then
            phiij=PHIMIX(i,j,itau,idel,tau,del,x)
            phisum=phisum+phiij
c           write (*,1200) i,j,x(1),x(2),x(3),phiij
          end if
        enddo
      enddo
      phihmx=phisum
c
      RETURN
      end                                               !function PHIHMX
c
c ======================================================================
c
      function PHIMIX (icomp,jcomp,itau,idel,tau,del,x)
c
c  compute reduced Helmholtz energy of mixing (or a derivative)
c  for the binary interaction of components i and j as a function of
c  composition and dimensionless temperature and density for the
c  mixture Helmholtz equation of state
c
c  inputs:
c    icomp--component i
c    jcomp--component j
c     itau--flag specifying order of temperature derivative to calc
c     idel--flag specifying order of density derivative to calculate
c           when itau = 0 and idel = 0, compute Amix/RT
c           when itau = 0 and idel = 1, compute 1st density derivative
c           when itau = 1 and idel = 1, compute cross derivative
c           etc.
c      tau--dimensionless temperature (To/T)
c      del--dimensionless density (D/Do)
c        x--composition array (mol frac)
c  output (as function value):
c      phi--mixture interaction (excess) part of the Helmholtz energy,
c           or one of its derivatives (as specified by itau and idel),
c           in reduced form (Amix/RT)
c
c        The Helmholtz energy consists of ideal gas and residual (real-
c        gas) terms.  The residual term consists of ideal-solution and
c        mixing terms.  This routine calculates only the mixing term.
c
c  N.B.  This routine incorporates the mixture-specific multiplier
c        (e.g. the Fpq of the Lemmon & Jacobsen model).
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-02-95  MM, original version
c  11-26-95  MM, rearrange argument list (x at end)
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  01-10-96  MM, add commons /CFXHMX/ and /CFIHMX/
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  03-21-96  MM, delete reference to /MODEL/, not used
c  05-06-96  MM, add /MIXMOD/, cases for different mix models
c  06-12-96  MM, add /HMXSAV/, LJi and LMi models
c  07-01-96  MM, move x(i)*x(j) multiplier to PHIHMX; remove x from inputs
c  11-04-96  MM, replace van der Waals rule with linear (LIN),
c                change nmxpar from 4 to 6
c  11-26-97 APP, add mixture model for Ammonia-Water and common block /AMMWAT/
c  12-01-97  MM, add composition to inputs (used by ammonia-water)
c  01-07-98  MM, move x(i)*x(j) back to here from PHIHMX
c  12-22-98 EWL, reorganize ammonia-water code
c  12-22-98 EWL, fix errors in derivatives of ammonia-water equation
c  12-24-98 EWL, add third derivative of tau to LJ model
c  12-24-98 EWL, change AMMWAT common block to FLAGS2 common block
c  09-00-00 EWL, add logic to return old value if inputs are same
c  09-00-00 EWL, change del^d, etc. to LOG format to increase calculation speed
c  09-00-00 EWL, add exponential pieces to mixture terms
c  02-14-01 EWL, increase dimension on amix to 4
c                change use of iexpmx(i,j,k) to amix(i,j,k,4)
c  07-30-01 EWL, move iamwat=0 to beginning of code
c  08-08-01 EWL, change -delli to delli to mimic core_feq routines
c  08-17-06 EWL, update program for GERG-2004 equation of Kunz and Wagner
c  10-19-06 EWL, change calculation of ncode2 to work with "KW"
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nx=ncmax)
      parameter (nmxpar=6)     !number of binary mixture parameters
      parameter (nmxtrm=15)    !number of terms in binary mixing rule
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      parameter (nmsav=25)     !number of mixture models
      character*3 hmodmx
      character*16 drvflg(nmsav)
      dimension x(ncmax)
      common /NCOMP/ nc,ic
      common /MIXMOD/ hmodmx(nx,nx)
c  commons associated with the binary mixing rule(s)
c  the amix(i,j,k,1..3) are the coefs for the (general) mixing term
c  the fmix(i,j,1..nmxpar) are the parameters for the i-j binary
c  the namix(i,j) are the number of a(i,j,k,1..3) terms
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
      common /CFIHMX/ namix(nx,nx),namix2(nx,nx)
      common /HMXSV2/ drvflg
      common /HMXSAV/ phisav(nmsav,0:nmxtrm),tausav(nmsav),delsav(nmsav)
     &                ,drvsav(nmsav,16),delli(nmsav,nmxtrm)
      common /FLAGS2/ iamwat,ianc,iwat
      common /MXTRMS/ phi01(nmxtrm),phi10(nmxtrm),phi02(nmxtrm),
     &                phi20(nmxtrm),phi03(nmxtrm)
c
c
      phimix=0.0d0
      if (del.le.1.0d-10) then    !trivial solution at zero density
        RETURN                    !for any and all derivatives
      end if
      if (tau.le.0.0d0) RETURN
c
      i=icomp
      j=jcomp
      ncode2=0
      if (hmodmx(i,j)(1:2).eq.'LJ') ncode2=ichar(hmodmx(i,j)(3:3))-47
      if (hmodmx(i,j)(1:2).eq.'KW') ncode2=ichar(hmodmx(i,j)(3:3))-37
      if (hmodmx(i,j)(1:2).eq.'LM') ncode2=ichar(hmodmx(i,j)(3:3))-27
      if (ncode2.lt.0 .or. ncode2.gt.nmsav) then
        phimix=1.d20
        RETURN
      endif
      if (hmodmx(i,j).eq.'LIN') then
c
c  linear mixing rules:  all departures from ideal solution are
c  contained in the reducing parameters (i.e. Kt and/or Kv)
c
        phimix=0.0d0
c
      elseif (ncode2.ne.0) then
        ncode=idel*4+itau+1
c
c  Lemmon-Jacobsen or modified Lemmon-Jacobsen mixing function
c
        phisum=0.0d0
        ifij=3
        if (hmodmx(i,j)(1:2).eq.'KW') ifij=5
        if (abs(tau-tausav(ncode2)).lt.1.0d-12 .and.
     &      abs(del-delsav(ncode2)).lt.1.0d-12) then
          if (drvflg(ncode2)(ncode:ncode).eq.'1') then
            PHIMIX=x(i)*x(j)*fmix(i,j,ifij)*drvsav(ncode2,ncode)
            RETURN
          endif
c  retrieve value from previous call
        else
c  otherwise, compute new values and save for possible future use
          drvflg(ncode2)='0000000000000000'
          do k=1,namix(i,j)
            delli(ncode2,k)=0
            if(amix(i,j,k,4).ne.0.d0)delli(ncode2,k)=del**amix(i,j,k,4)
            savijk=amix(i,j,k,1)*EXP(amix(i,j,k,2)*log(tau)
     &            +amix(i,j,k,3)*log(del)-delli(ncode2,k))
            phisav(ncode2,k)=savijk
            phisum=phisum+savijk
c           write (*,1010) k,phisav(ncode2,k),phisum
c1010       format (1x,i3,2d30.20)   !write out each term for debugging
          enddo
          if (namix2(i,j).ne.0) then
            do k=namix(i,j)+1,namix(i,j)+namix2(i,j)
              delli(ncode2,k)=amix(i,j,k,4)*(del-amix(i,j,k,5))**2
              delli2         =amix(i,j,k,6)*(del-amix(i,j,k,7))
              savijk=amix(i,j,k,1)*EXP(amix(i,j,k,2)*log(tau)
     &              +amix(i,j,k,3)*log(del)-delli(ncode2,k)
     &                                     -delli2)
              phisav(ncode2,k)=savijk
              phisum=phisum+savijk
            enddo
          endif
          tausav(ncode2)=tau
          delsav(ncode2)=del
          drvsav(ncode2,1)=phisum
          drvflg(ncode2)(1:1)='1'
        end if
c
c  check if derivatives are requested, calculations make use of fact
c  that terms in derivative summations are very similar to A/RT terms
c
        if (idel.eq.1) then
c  compute derivative w.r.t. del (dimensionless density)
c  save individual terms for possible use in cross derivative
          phisum=0.0d0
          do k=1,namix(i,j)
            phi01(k)=phisav(ncode2,k)
     &              *(amix(i,j,k,3)-amix(i,j,k,4)*delli(ncode2,k))
            phisum=phisum+phi01(k)
          enddo
          if (namix2(i,j).ne.0) then
            do k=namix(i,j)+1,namix(i,j)+namix2(i,j)
              phi01(k)=phisav(ncode2,k)
     &                *(amix(i,j,k,3)
     &   -2.d0*del*amix(i,j,k,4)*(del-amix(i,j,k,5))-del*amix(i,j,k,6))
              phisum=phisum+phi01(k)
            enddo
          endif
        else if (idel.eq.2) then
c  compute 2nd derivative w.r.t. del (dimensionless density)
          phisum=0.0d0
          do k=1,namix(i,j)
            ax3=amix(i,j,k,3)
            dell=amix(i,j,k,4)*delli(ncode2,k)
            phi02(k)=phisav(ncode2,k)
     &           *(ax3**2-ax3-dell*(2.0d0*ax3+amix(i,j,k,4)-1.0d0-dell))
            phisum=phisum+phi02(k)
c         end if
          enddo
          if (namix2(i,j).ne.0) then
            do k=namix(i,j)+1,namix(i,j)+namix2(i,j)
              aa3=amix(i,j,k,3)
              aa4=amix(i,j,k,4)
              phi02(k)=phisav(ncode2,k)
     &  *(del**2*(aa3/del-2.d0*aa4*(del-amix(i,j,k,5))-amix(i,j,k,6))**2
     &  -aa3-2.d0*aa4*del**2)
              phisum=phisum+phi02(k)
            enddo
          endif
        else if (idel.eq.3) then
c  compute 3nd derivative w.r.t. del (dimensionless density)
          phisum=0.0d0
          do k=1,namix(i,j)
            ax3=amix(i,j,k,3)
            dell=amix(i,j,k,4)*delli(ncode2,k)
            dl=amix(i,j,k,4)
            phi03(k)=phisav(ncode2,k)
     &      *(2.0d0*ax3-3.0d0*ax3*ax3+ax3**3+dell*(-2.0d0+3.0d0*dl-dl**2
     &       +6.0d0*ax3-3.0d0*ax3*dl-3.0d0*ax3**2
     &       +dell*(-3.0d0+3.0d0*dl+3.0d0*ax3-dell)))
            phisum=phisum+phi03(k)
          enddo
          if (namix2(i,j).ne.0) then
            do k=namix(i,j)+1,namix(i,j)+namix2(i,j)
c  3rd derivative has not been implemented yet
c             phi03(k)=...
c             phisum=phisum+phi03(k)
            enddo
          endif
        end if
c
        if (itau.eq.1) then
c  compute derivative w.r.t. tau (dimensionless temperature)
c  save individual terms for possible use in cross derivative
          phisum=0.0d0
          do k=1,namix(i,j)+namix2(i,j)
            phi10(k)=phisav(ncode2,k)*amix(i,j,k,2)
            phisum=phisum+phi10(k)
          enddo
        else if (itau.eq.2) then
c  compute 2nd derivative w.r.t. tau (dimensionless temperature)
          phisum=0.0d0
          do k=1,namix(i,j)+namix2(i,j)
            ax2=amix(i,j,k,2)
            phi20(k)=phisav(ncode2,k)*ax2*(ax2-1.0d0)
            phisum=phisum+phi20(k)
          enddo
        else if (itau.eq.3) then
c  compute 3rd derivative w.r.t. tau (dimensionless temperature)
          phisum=0.0d0
          do k=1,namix(i,j)+namix2(i,j)
            ax2=amix(i,j,k,2)
            phisum=phisum+phisav(ncode2,k)*ax2*(ax2-1.0d0)*(ax2-2.0d0)
          enddo
        end if
c
        if (itau.eq.1 .and. idel.eq.1) then
c  compute cross derivative using terms from 1st derivatives
          phisum=0.0d0
          do k=1,namix(i,j)+namix2(i,j)
            if (phisav(ncode2,k).ne.0.d0)
     &          phisum=phisum+phi10(k)*phi01(k)/phisav(ncode2,k)
          enddo
        end if
        drvsav(ncode2,ncode)=phisum
        drvflg(ncode2)(ncode:ncode)='1'
        PHIMIX=x(i)*x(j)*fmix(i,j,ifij)*phisum    !fmix(i,j,3) is Fpq
c
      else if (hmodmx(i,j).eq.'TR1') then
c
c  ammonia-water model of Tillner-Roth
c
c  flip compositions; fmix(i,j,5) is gamma
         ij=iamwat
         phisum=0.0d0
c
         do k=1,14
           ak=amix(i,j,k,1)
           tk=amix(i,j,k,2)
           dk=amix(i,j,k,3)
           ik=amix(i,j,k,4)
           xk=1.0d0
           if (k.ge.7) xk=x(ij)
           if (k.eq.14) xk=xk**2
           delik=-del**ik
           bk=1.0d0
           if (ik.gt.0) bk=exp(delik)
           aka=xk*ak*tau**tk*del**dk*bk
c
           if (itau.eq.0 .and. idel.eq.0) then           !no deriv
             phisum=phisum+aka
           else if (itau.eq.0 .and. idel.eq.1) then    !1st deriv of del
             phisum=phisum+aka*(dk+ik*delik)
           else if (itau.eq.0 .and. idel.eq.2) then    !2nd deriv of del
             phisum=phisum+aka*(dk*(dk-1.0d0)
     &             +ik*delik*(ik-1.0d0+2.0d0*dk+ik*delik))
           else if (itau.eq.0 .and. idel.eq.3) then    !3rd deriv of del
             dki = ik*delik
             phisum=phisum+(dk*(dk-1.0d0)*(dk-2.0d0)
     &         +dki*(3.0d0*dk*((dk-2.0d0)+ik)+(ik-1.0d0)*(ik-2.0d0)
     &         +3.0d0*dki*(dk+(ik-1.0d0))+dki*dki))*aka
           else if (itau.eq.1 .and. idel.eq.0) then    !1st deriv of tau
             phisum=phisum+tk*aka
           else if (itau.eq.2 .and. idel.eq.0) then    !2nd deriv of tau
             phisum=phisum+tk*(tk-1.0d0)*aka
           else if (itau.eq.3 .and. idel.eq.0) then    !3rd deriv of tau
             phisum=phisum+tk*(tk-1.0d0)*(tk-2.0d0)*aka
           else if (itau.eq.1 .and. idel.eq.1) then!deriv of tau and del
             phisum=phisum+aka*tk*(dk+ik*delik)
           end if
         enddo
         PHIMIX=phisum*x(ij)*(1.0d0-x(ij)**fmix(i,j,5))
      else
c
c  space for additional mixing rule
c
        PHIMIX=0.0d0                !not currently implemented
*       write (*,*) ' PHIMIX ERROR--unknown mixing rule'
      end if
c
      RETURN
      end                                               !function PHIMIX
c
c ======================================================================
c
      subroutine CRTHMX (x,tcrit,pcrit,Dcrit,ierr,herr)
c
c  returns critical parameters associated with mixture Helmholtz EOS
c
c  N.B.  these critical parameters are estimates based on polynomial
c        fits to the binary critical lines and (for 3 or more components)
c        combining rules applied to the constituent binaries
c
c  input:
c        x--composition array (mol frac)
c  outputs:
c    tcrit--critical temperature (K)
c    pcrit--critical pressure (kPa)
c    Dcrit--molar density (mol/L) at critical point
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  11-02-95  MM, original version
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  01-29-96  MM, critical parameters based on average of pure fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                add Zcrit to common /CCON/
c  03-19-19  MM, add dipole moment to /CCON/
c  05-06-96  MM, replace temporary func for crit pars with reducing pars
c  11-19-96  MM, add empirical function for critical line(s)
c  11-14-97 EWL, fix bug in critical line correlation
c  11-24-97  MM, Tsum = x(i)*x(j)*Tij, not xij*tcij, ditto for Vsum
c  01-21-98 APP, add Rainwater's ammonia-water critical line equations
c                for hmodmx(i,j) = 'TR1'
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  12-22-98 EWL, if x(i)=1 for any component, return with pure fluid parameters
c  12-23-98 EWL, rearrange code to eliminate extra calculations with amm-water
c  03-04-02 EWL, check sum(x)=0
c  08-17-06 EWL, update program for GERG-2004 equation of Kunz and Wagner
c  08-21-06 EWL, if TC1 or VC1 are not available, return reducing parameters instead of an error message
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines
      character*1 htab,hnull
      character*3 hmodtc,hmodvc
      character*255 herr
      common /NCOMP/ nc,ic
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(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 /HCHAR/ htab,hnull
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c  commons associated with the binary critical lines
      common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx)
      common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc)
c
      dimension x(ncmax)
c     double precision term1,term2,term3,term4,term5,term6
c     double precision xt
      ierr=0
      herr=' '
c
      if (nc.eq.1 .or. ic.ne.0) then
        icomp=1
        if (ic.ne.0) icomp=ic
        tcrit=tc(icomp)
        Dcrit=rhoc(icomp)
        pcrit=pc(icomp)
        RETURN
      end if
c  check for x(i)=1 in a mixture setup
      do i=1,nc
        if (ABS(x(i)-1.0d0).lt.1.0d-8) then
          tcrit=tc(i)
          Dcrit=rhoc(i)
          pcrit=pc(i)
          RETURN
        end if
      enddo
c
c      if (hmodmx(1,2).eq.'TR1') then
cc  Moldover and Rainwater, J. Chem. Phys. 88(12):7772, 1988.
c        call Rmix (x)
c        term1=x(1)/(R*tc(1))
c        term2=x(2)/(R*tc(2))
c        term3=x(1)*x(2)
c        term4=1.0d0-2.0d0*x(2)
cc       ftcij(1,2,k),k=1,3 are T1,T2, and T3 constants
c        term5=ftcij(1,2,1)+term4*ftcij(1,2,2)+term4**2*ftcij(1,2,3)
c        term6=term1+term2+term3*term5
c        tcrit=1.0d0/(term6*R)
cc
c        xt=(1.0d0/tc(2)-1.0d0/tcrit)/(1.0d0/tc(2)-1.0d0/tc(1))
c        term1=(1.0d0-xt)*pc(1)/(R*tc(1))
c        term2=xt*pc(2)/(R*tc(2))
c        term3=xt*(1.0d0-xt)
c        term4=1.0d0-2.0d0*xt
c        term5=ftcij(1,2,4)+term4*ftcij(1,2,5)+term4**2*ftcij(1,2,6)
c        term6=term1+term2+term3*term5
c        pcrit=term6*R*tcrit
cc
c        term1=xt*rhoc(1)
c        term2=(1.0d0-xt)*rhoc(2)
c        term5=fvcij(1,2,1)+term4*fvcij(1,2,2)+term4**2*fvcij(1,2,3)
c        term6=term1+term2+term3*term5
c        dcrit=term6
cc       write (*,*) 'new tc, pc, rhoc for NH3=',x(1),tcrit,pcrit,dcrit
c        RETURN
c      end if
c
      tsum=0.0d0
      Vsum=0.0d0
      xsum=0.0d0
c     write (*,1010) (x(k),k=1,nc),(tc(k),k=1,nc)
c1010 format (' CRTHMX input x,Tc:  ',10f10.4)
      do i=1,nc-1
        do j=i+1,nc
          xsum=xsum+x(i)*x(j)
c  the tcij and Vcij are the binary critical line correlations, and
c  these terms (ONLY) are evaluated at the pseudo-binary compositions
c  defined by xi and xj:
          xi=0.5d0*(x(i)+(1.0d0-x(j)))
          xj=0.5d0*(x(j)+(1.0d0-x(i)))
          xij=xi*xj
          if (hmodtc(i,j).eq.'TC1') then
            tcij=xi*tc(i)+xj*tc(j)+xij*ftcij(i,j,1)
c  terms assuming mixture specified in same order as in hmx.bnc file
     &          +xj*xi*(ftcij(i,j,2)+(ftcij(i,j,3)+(ftcij(i,j,4)
     &          +(ftcij(i,j,5)+ftcij(i,j,6)*xi)*xi)*xi)*xi)*xi
c  extra terms to allow specifying mix in opposite order (j,i)
c  note:  either terms 2-6 or 7-11 should be set to zero in RDBNC
     &          +xi*xj*((ftcij(i,j,7)+(ftcij(i,j,8)+(ftcij(i,j,9)
     &          +(ftcij(i,j,10)+ftcij(i,j,11)*xj)*xj)*xj)*xj)*xj)
c           write (*,*) 'Tcij using mixture model TC1:  ',i,j,tcij
          else
c           ierr=1
c           herr='[CRTHMX error] Tcrit model not found'//hnull
c           call ERRMSG (ierr,herr)
c           tcij=xi*tc(i)+xj*tc(j)
            call REDHMX (x,tcrit,Dcrit)
            goto 100
          end if
          tsum=tsum+x(i)*x(j)*tcij
          if (hmodvc(i,j).eq.'VC1') then
            Vcij=xi/rhoc(i)+xj/rhoc(j)+xij*fvcij(i,j,1)
c  terms assuming mixture specified in same order as in hmx.bnc file
     &          +xj*xi*((fvcij(i,j,2)+(fvcij(i,j,3)+(fvcij(i,j,4)
     &          +(fvcij(i,j,5)+fvcij(i,j,6)*xi)*xi)*xi)*xi)*xi)
c  extra terms to allow specifying mix in opposite order (j,i)
c  note:  either terms 2-6 or 7-11 should be set to zero in RDBNC
     &          +xi*xj*((fvcij(i,j,7)+(fvcij(i,j,8)+(fvcij(i,j,9)
     &          +(fvcij(i,j,10)+fvcij(i,j,11)*xj)*xj)*xj)*xj)*xj)
          else
c           ierr=1
c           herr='[CRTHMX error] Vcrit model not found'//hnull
c           call ERRMSG (ierr,herr)
c           Vcij=xi/rhoc(i)+xj/rhoc(j)
            call REDHMX (x,tcrit,Dcrit)
            goto 100
          end if
          Vsum=Vsum+x(i)*x(j)*Vcij
c         write (*,1016) i,j,tcij,Vcij,xsum
c1016     format (' CRTHMX:  i,j,Tc_ij,Vc_ij,xsum:  ',2i3,3f12.6)
c         write (*,1017) i,j,(ftcij(i,j,k),k=1,nmxtc)
c1017     format (' CRTHMX:  i,j,ftcij:  ',2i3,8f10.6)
        enddo
      enddo
      if (xsum.gt.1.0d-10) then
        tcrit=tsum/xsum
        Dcrit=xsum/Vsum
      else
        tcrit=300
        Dcrit=10
        Pcrit=10
        ierr=1
        herr='[CRTHMX error] sum(x)=0'//hnull
        call ERRMSG (ierr,herr)
      end if
c
c  base Pc on Tc, Dc
 100  continue
      call PRESS (tcrit,Dcrit,x,pcrit)
      if (pcrit.lt.0) pcrit=pc(1)
c     write (*,1024) tcrit,1.0/Dcrit,pcrit
c1024 format (' CRTHMX:  Tcrit,Vcrit,Pcrit:  ',15x,3f12.6)
c
      RETURN
      end                                             !subroutine CRTHMX
c
c ======================================================================
c
      subroutine REDHMX (x,tred,Dred)
c
c  returns reducing parameters associated with mixture Helmholtz EOS;
c  used to calculate the 'tau' and 'del' which are the independent
c  variables in the EOS
c
c  input:
c        x--composition array [mol frac]
c  outputs:
c     tred--reducing temperature [K]
c     Dred--reducing molar density [mol/L]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-02-95  MM, original version
c  11-29-95  MM, variable lower limit on coefficient/constant arrays
c                to accommodate ECS reference fluid
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                add Zcrit to common /CCON/
c  03-19-19  MM, add dipole moment to /CCON/
c  03-21-96  MM, replace call to REDFEQ with general REDX, delete /MODEL/
c  05-06-96  MM, add /MIXMOD/,/CFXHMX/,/CFIHMX/; add Kt, Kv dependence;
c                branch for arithmetic or geometric mean of temperature
c  05-14-96  MM, modify volume mixing rule for VDW model; use arithmetic
c                average of critical volumes, rather than Reiner's rule
c  06-10-96  MM, modify volume mixing rule for VDW model (again); use
c                average of cube root of critical volumes
c                add Lemmon-Jacobsen mixing rules
c  06-12-96  MM, add LMx mixing rules (Lemmon-Jacobsen modified to be of
c                form compatible with van der Waal's (VDW) model
c  07-01-96  MM, fix bug--missing **3 in VDW volume term
c  11-04-96  MM, replace van der Waals rule with linear (LIN),
c                change nmxpar from 4 to 6
c  11-08-96  MM, fix wrong array reference for Kv
c  11-14-96  MM, fix wrong array reference for zeta, xi (LJi model)
c  05-22-97  MM, bug fix:  beta, gamma were applied to V, should not be
c  11-26-97 APP, add Ammonia-Water model and common block /AMMWAT/
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  12-24-98 EWL, change AMMWAT common block to FLAGS2 common block
c  12-24-98 EWL, fix ammonia-water code to fit within do-loop
c  12-24-98 EWL, avoid reducing values of 0 when a fluid has not been loaded
c  02-14-01 EWL, increase dimension on amix to 4
c                change use of iexpmx(i,j,k) to amix(i,j,k,4)
c  08-08-01 EWL, remove the 0.5* in t12= and V12= used in the linear and LM models
c  08-17-06 EWL, update program for GERG-2004 equation of Kunz and Wagner
c  08-21-06 EWL, add Kunz and Wagner mixing rules
c  03-05-07 EWL, add fmix2 array in case fluids are entered in opposite order as required by Kunz and Wagner
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      parameter (nmxpar=6)     !number of binary mixture parameters
      parameter (nmxtrm=15)    !number of terms in binary mixing func
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      character*3 hmodmx
      character*255 herr
      common /NCOMP/ nc,ic
      common /MIXMOD/ hmodmx(nx,nx)
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
      common /CF2HMX/ fmix2(nx,nx,nmxpar)
      common /FLAGS2/ iamwat,ianc,iwat
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
      dimension x(ncmax),treda(ncmax),Dreda(ncmax)
c
      tsum=0.0d0
      Vsum=0.0d0
      do i=1,nc
        treda(i)=tz(i)
        Dreda(i)=rhoz(i)
      enddo
c
      do i=1,nc
        do j=i,nc
          if (i.eq.j) then
c  special case when i = j
c           write (*,*) ' REDHMX--special case for i = j',i,j
            tsum=tsum+x(i)*treda(i)
            Vsum=Vsum+x(i)/Dreda(i)
          else if (hmodmx(i,j)(1:2).eq.'KW') then
c  Kunz-Wagner rule
            if (x(i).gt.1.0d-10 .and. x(j).gt.1.0d-10) then
             xij=2.d0*x(i)*x(j)
             xijt=(x(i)+x(j))/(fmix(i,j,1)**2*x(i)+fmix2(i,j,1)**2*x(j))
             xijV=(x(i)+x(j))/(fmix(i,j,3)**2*x(i)+fmix2(i,j,3)**2*x(j))
             t12=(treda(i)*treda(j))**0.5d0
             V12=(Dreda(i)**(-1.d0/3.d0)+Dreda(j)**(-1.d0/3.d0))**3/8.d0
             tsum=tsum+xij*fmix(i,j,1)*fmix2(i,j,1)*fmix(i,j,2)*xijt*t12
             Vsum=Vsum+xij*fmix(i,j,3)*fmix2(i,j,3)*fmix(i,j,4)*xijV*V12
c  the following addition allows the use of linear mixture rules x(i)*Tc(i)
c  when the specified rule (as in this case for the GERG-2004 equation) is
c  not linear:  x(i)^2*Tc(i)
             tsum=tsum-x(i)*x(j)*(treda(i)+treda(j))
             Vsum=Vsum-x(i)*x(j)*(1.d0/Dreda(i)+1.d0/Dreda(j))
            else
c         zero composition for one or both components, no change in sums
            end if
          else if (hmodmx(i,j).eq.'LIN') then
c  simple linear mixing rules
            xij=x(i)*x(j)
            t12=(treda(i)+treda(j))
            V12=(1.0d0/Dreda(i)+1.0d0/Dreda(j))
            tsum=tsum+xij*(fmix(i,j,1)-1.0d0)*t12     !fmix(i,j,1) = Kt(i,j)
            Vsum=Vsum+xij*(fmix(i,j,2)-1.0d0)*V12     !fmix(i,j,2) = Kv(i,j)
          else if (hmodmx(i,j)(1:2).eq.'LM') then
c  Lemmon-Jacobsen rule modified to use Kt, Kv, rather than zeta, xi
            if (x(i).gt.1.0d-10 .and. x(j).gt.1.0d-10) then
              xij=x(i)*x(j)
              xijt=x(i)**fmix(i,j,4)*x(j)**fmix(i,j,5)!fmix(i,j,4) = beta
              t12=(treda(i)+treda(j))           !fmix(i,j,5) = gamma
              V12=(1.0d0/Dreda(i)+1.0d0/Dreda(j))
              tsum=tsum+xijt*(fmix(i,j,1)-1.0d0)*t12  !fmix(i,j,1) = Kt(i,j)
              Vsum=Vsum+xij*(fmix(i,j,2)-1.0d0)*V12   !fmix(i,j,2) = Kv(i,j)
            else
c         zero composition for one or both components, no change in sums
            end if
          else if (hmodmx(i,j)(1:2).eq.'LJ') then
c  original Lemmon-Jacobsen rule using zeta, xi
            if (x(i).gt.1.0d-12 .and. x(j).gt.1.0d-12) then
              xij=x(i)*x(j)
              xijt=x(i)**fmix(i,j,4)*x(j)**fmix(i,j,5)!fmix(i,j,5) = gamma
              tsum=tsum+xijt*fmix(i,j,1)              !fmix(i,j,1) = zeta(T)
              Vsum=Vsum+xij*fmix(i,j,2)               !fmix(i,j,2) = xi  (x)
            else
c             zero composition for one or both components, no change in sums
            end if
c  mixing parameters for the Ammonia-Water model
c  fmix(i,j,1)=Kt, fmix(i,j,2)=Kv, fmix(i,j,3)=alpha, fmix(i,j,4)=beta
          elseif (hmodmx(i,j).eq.'TR1') then
             ij=i
             if (iamwat.eq.2) ij=j
             t12=treda(i)+treda(j)
             v12=1.0d0/Dreda(i)+1.0d0/Dreda(j)
             tsum=tsum+x(ij)*(1.0d0-x(ij)**fmix(i,j,3))*fmix(i,j,1)*t12
     &                -x(i)*x(j)*t12
             Vsum=Vsum+x(ij)*(1.0d0-x(ij)**fmix(i,j,4))*fmix(i,j,2)*v12
     &                -x(i)*x(j)*v12
          else
c  mixing rule not specified or not found--should not get here
c  call ERRMSG anyways in case error messages are being written to the screen
            ierr=1
            herr='[REDHMX error] Mixing rule not found for i,j'
            call ERRMSG (ierr,herr)
*           write (*,*) herr
          end if
 1    format (2i3,a5,f12.5)
        enddo
      enddo
      if (ABS(tsum).gt.0.5d-12 .and. ABS(Vsum).gt.0.5d-12) then
        tred=tsum
        Dred=1.0d0/Vsum
      else
        tred=100
        Dred=10
      endif
c
c      write (*,*) "tred and dred for X:",x(1),x(2),tred,Dred
c1120 format ('  REDHMX--x,t0,D0: ',f10.7,f10.3,f10.5)
c
      RETURN
      end                                             !subroutine REDHMX
c
c ======================================================================
c
      subroutine SETHMX (hfmix,ierr,herr)
c
c  read parameters for the HMX mixture model from a file
c
c  input:
c    hfmix--file name containing mixture parameters [character*255]
c  outputs:
c     ierr--error flag:  0 = successful
c                        1 = error (e.g. fluid not found)
c     herr--error string (character*255 variable if ierr<>0)
c     [fluid parameters, etc. returned via various common blocks]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  01-09-96  MM, original version
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                add Zcrit to common /CCON/
c  03-19-96  MM, add dipole moment to /CCON/
c  03-22-96  MM, replace /MODEL/ with /EOSMOD/
c  05-06-96  MM, add /MIXMOD/ to specify model for binary pairs;
c                initialize j,i same as i,j; initialize i = j to 1.0
c  05-07-96  MM, move read of binary model(s) to RDMIX
c  05-08-96  MM, add /MXINFO/ and read corresponding info from file
c  11-19-96  MM, move initialization of common blocks to RDBNC, ESTBNC
c  03-27-97  MM, if file not found, search in likely directories
c  07-03-97  MM, pass any error from ESTBNC up the chain
c  08-22-97  MM, use double backslash on file open to avoid problem with Unix machines
c  11-26-97 APP, additions for Ammonia-Water calculations, new common block /AMMWAT/
c  04-15-98 EWL, add check for nist14 directory
c  07-08-98 EWL, change character strings from *80 to *255
c  12-24-98 EWL, change AMMWAT common block to FLAGS2 common block
c  12-24-98 EWL, restructure check for ammonia-water mixtures
c  01-10-01 EWL, only call checks for estimation scheme if ierr1<>0, not ierr
c  02-02-01 EWL, remove double backslash, which caused problems in Windows
c  02-14-01 EWL, increase dimension on amix to 4
c                change use of iexpmx(i,j,k) to amix(i,j,k,4)
c  01-23-02 EWL, split common block MXINFO into 2 pieces for some compilers
c  07-13-04 EWL, read in a version number on the second line of hmx.bnc
c  09-28-05  DT, change occurrences of a backslash to 'char(92)' for compatibility
c                with some compilers using F77
c  08-17-06 EWL, update program for GERG-2004 equation of Kunz and Wagner
c  09-06-06 EWL, change nbrule to 20
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      parameter (nbrule=20)    !number of binary mixing rules
      parameter (nbin=ncmax*(ncmax-1)/2)   !# possible binary pairs
      parameter (nmxpar=6)     !number of binary mixture parameters
      parameter (nmxtrm=15)    !number of terms in binary mixing rule
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines
      parameter (nmsav=25)     !number of mixture models
      character*1 htab,hnull
      character*3 hpheq,heos,hmxeos,hmodcp
      character*3 hmodmx,hmodtc,hmodvc,hprkij
      character*3 htype
      character*12 hcas
      character*255 hfmix
      character*255 herr,herr1
c  next 3 declarations associated with /MXINFO/
      character*8 hbpar
      character*255 hmfile,hrule
c
      logical lmix,lij(nx,nx),lfbnc,lftc,lfvc
c
      common /NCOMP/ nc,ic
      common /EOSMOD/ hpheq,heos,hmxeos(n0:nx),hmodcp(n0:nx)
      common /CCAS/ hcas(n0:nx)
      common /HCHAR/ htab,hnull
      common /FLAGS2/ iamwat,ianc,iwat
c  commons associated with the binary mixing rule(s)
c  the amix(i,j,k,1..3) are the coefs for the (general) mixing term
c  the fmix(i,j,1..nmxpar) are the parameters for the i-j binary
c  the namix(i,j) are the number of a(i,j,k,1..3) terms
      common /MIXMOD/ hmodmx(nx,nx)
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
      common /CF2HMX/ fmix2(nx,nx,nmxpar)
      common /CFIHMX/ namix(nx,nx),namix2(nx,nx)
c  /MXINFO/ contains information on the mixing rules and parameters
c  hmfile specifies the files from which mixing rules originate
c  (stored in order 1,2  1,3  2,3; element zero is file called in SETUP)
c  hbin provides documentation for the current binary parameters
c  hrule contains descriptions of the currently available mixing rules
c  hbpar contains descriptions of the binary parameters (e.g. Kt, Kv)
c  associated with the currently available mixing rules
      common /MXINFO/ hmfile(0:nbin),hrule(nbrule),hbpar(nbrule,nmxpar)
c  commons associated with the binary critical lines
      common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx)
      common /HMXSAV/ phisav(nmsav,0:nmxtrm),tausav(nmsav),delsav(nmsav)
     &                ,drvsav(nmsav,16),delli(nmsav,nmxtrm)
      common /HMXFLG/ ivrsnm
      common /PRMOD/ hprkij(nx,nx)
      common /CFXPR/ fprkij(nx,nx)
      common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc)
c
      hmfile(0)=hfmix
      ierr=0
      ierr1=0
      herr=' '
      do i=1,nx
        do j=1,nx
          fprkij(i,j)=0.d0
        enddo
      enddo
c  (re)initialize contents of /HMXSAV/ when a new mixture is read in
      do i=1,nmsav
        delsav(i)=0.0d0
        tausav(i)=0.0d0
        do k=0,nmxtrm
          phisav(i,k)=0.0d0
        enddo
      enddo
      iamwat=0
c
c  initialize the diagonal terms (i = j)
c
      do i=1,nc
        hmodmx(i,i)='i=j'
        hmodtc(i,i)='i=j'
        hmodvc(i,i)='i=j'
        hprkij(i,i)='   '
c  the i = j terms in the fmix summation should not be called (neither
c  the i > j), but set the number of terms to zero as precaution
        namix(i,i)=0
        namix2(i,i)=0
        do k=1,nmxpar
          fmix(i,i,k)=0.0d0
          fmix2(i,i,k)=0.0d0
        enddo
        do k=1,nmxtc
          ftcij(i,i,k)=0.0d0
        enddo
        do k=1,nmxvc
          fvcij(i,i,k)=0.0d0
        enddo
      enddo
c
      if (hfmix(1:5).eq.'BDATA' .or. hfmix(1:5).eq.'bdata') then
c  get coefficients from common blocks--not currently implemented
        ierr=1
        herr='[SETHMX error] Block data option not available'//hnull
        call ERRMSG (ierr,herr)
        RETURN
      else
c
c  read parameters for the binary pairs from file
c
c        write (*,1105) hfmix
c1105   format (/1x,'reading mixture coefficients from file: (',a40,')')
        nread=12             !logical unit for file reads
        call OPENFL (nread,hfmix,1,ierr,herr)
        if (ierr.ne.0) goto 999
        rewind (nread)
        read (nread,2003) htype
        if (htype.ne.heos) then
          ierr=1
          herr='[SETHMX error] Mixture file does not match model'//hnull
          call ERRMSG (ierr,herr)
          RETURN
        end if
        read (nread,'(i2)') ivrsnm
        ibin=0     !initialize counter for binary pair number
        do i=1,nc-1
          do j=i+1,nc
            ibin=ibin+1
            hmfile(ibin)=hfmix
            call RDBNC (nread,i,j,ibin,ierr1,herr1)
            if (ierr1.ne.0 .and. ierr.le.0) then
              ierr=ierr1
              herr=herr1
            endif
c  check if any binary parameters need to be estimated
            lfbnc=.true.
            lftc=.true.
            lfvc=.true.
            if (ierr1.ne.0) then
              if (ierr.lt.0) then
                ierr2=-ierr-80
                if (mod(ierr2,4).eq.0) then
                  lfvc=.true.
                  ierr2=ierr2-4
                end if
                if (mod(ierr2,2).eq.0) then
                  lftc=.true.
                  ierr2=ierr2-2
                end if
                if (ierr2.eq.1) then
                  lfbnc=.true.
                end if
              end if
              if (lfbnc) then
                call ESTBNC (i,j,ibin,ierr1,herr1)
                hmfile(ibin)='estimated by ESTBNC'
                if (ierr1.ne.0 .and. ierr.le.0) then
                  ierr=ierr1
                  herr=herr1
                endif
              end if
              if (lftc) call ESTTC (i,j,ierr1,herr1)
              if (lfvc) call ESTVC (i,j,ierr1,herr1)
c  any error from ESTBNC takes precedence over those from ESTTC or ESTVC
              if (ierr1.ne.0 .and. ierr.le.0) then
                ierr=ierr1
                herr=herr1
              end if
            end if
c  check to see if Peng-Robinson parameters need to be estimated
            if (fprkij(i,j).eq.0.d0) call ESTPR (i,j,ibin,ierr1,herr1)
          enddo
        enddo
c
c  read coefficients for general mixing rule(s) from file
c
        call RDMIX (nread,icomp,jcomp,'ALL',lij,ierr1,herr1)
        if (ierr1.ne.0 .and. ierr.le.0) then
          ierr=ierr1
          herr=herr1
        end if
c  check that mixing term for each binary has been found
        lmix=.true.
        do i=1,nc-1
          do j=i+1,nc
            lmix=lmix.and.lij(i,j)
*           write (*,1024) i,j,lij(i,j),lmix
*1024       format (1x,' SETHMX  i,j,lij,lmix:        ',2i3,1x,2i2)
          enddo
        enddo
c
c  check for the presence of water in the mixture
        do i=1,nc
          if (hcas(i).eq.'7732-18-5') iwat=1
        enddo

c
c  check for Ammonia-Water model
        if (hmodmx(1,2).eq.'TR1') then
c  set whether mixture is water-ammonia or ammonia-water
          if (hcas(1).eq.'7664-41-7'.and.hcas(2).eq.'7732-18-5')iamwat=1
          if (hcas(1).eq.'7732-18-5'.and.hcas(2).eq.'7664-41-7')iamwat=2
          if (iamwat.eq.0) then
            ierr=106
            herr='[SETMHX error 106] Wrong components; '//
     &           'model TR1 valid for ammonia + water only'//hnull
            call ERRMSG (ierr,herr)
          endif
          if (nc.ne.2) then
            ierr=107
            herr='[SETMHX error 107] Wrong number of components;'//
     &           ' model TR1 valid for binary mixtures only'//hnull
            call ERRMSG (ierr,herr)
          end if
        end if
c
        if (.not.lmix .and. abs(ierr).ne.117) then
          ierr=-117
          herr='[SETHMX error] Mixing term not found for one or '//
     &         'more binary pair(s)'//hnull
          call ERRMSG (ierr,herr)
        end if
c  rewind and close the file
        rewind (nread)
        close (nread)
c
c  write all parameters for debugging
*       write (*,*)
*       write (*,*) ' SETHMX:  mixing parameters'
*       write (*,1281)
*1281   format ('   i  j  mod    par_1     par_2     par_3',
*    &           '     par_4     par_5     par_6'/)
*       do i=1,nc
*       do j=1,nc
*       write (*,1280) i,j,hmodmx(i,j),(fmix(i,j,k),k=1,nmxpar)
*       enddo
*       enddo
*1280   format (1x,2i3,2x,a3,11f10.5)
*       write (*,*)
*       write (*,*) ' SETHMX:  critical line (Tc) parameters'
*       write (*,1283)
*1283   format ('   i  j  mod    par_1     par_2     par_3',
*    &           '     par_4     par_5     par_6     par_7',
*    &           '     par_8     par_9    par_10    par_11'/)
*       do i=1,nc
*       do j=1,nc
*       write (*,1280) i,j,hmodtc(i,j),(ftcij(i,j,k),k=1,nmxtc)
*       enddo
*       enddo
*       write (*,*)
*       write (*,*) ' SETHMX:  critical line (Vc) parameters'
*       write (*,1283)
*       do i=1,nc
*       do j=1,nc
*       write (*,1280) i,j,hmodvc(i,j),(fvcij(i,j,k),k=1,nmxvc)
*       enddo
*       enddo
        RETURN
      end if
c
c  error terminations
c
 999  continue
      ierr=111
      herr='[SETHMX error] Error in opening mixture file'//hnull
      call ERRMSG (ierr,herr)
      RETURN
c
 2003 format (a3)
c
      end                                             !subroutine SETHMX
c
c ======================================================================
c
      subroutine RDBNC (nread,icomp,jcomp,ibin,ierr,herr)
c
c  read binary parameters for the HMX mixture model and also the
c  critical line model from a file
c
c  inputs:
c    nread--file to read data from
c           <= 0 get data from block data (not currently implemented)
c           >0 read from logical unit n (file should have already
c              been opened and pointer set by subroutine SETHMX)
c    icomp--component i
c    jcomp--component j
c     ibin--counter for binary pair (1,2 = 1; 1,3 = 2; 2,3 = 3; etc.)
c  outputs:
c     ierr--error flag:  0 = successful
c                     -114 = critical temperature line will be estimated
c                     -115 = critical volume line will be estimated
c                     -116 = critical temperature and volume lines will be estimated
c                     -117 = binary pair not found, all parameters will be estimated
c                     -118 = binary pair found, but no thermo interaction parameters
c                      118 = error in reading mixture file
c     herr--error string (character*255 variable if ierr<>0)
c     other outputs written to common
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  01-10-96  MM, original version
c  01-12-96  MM, call ESTBNC if binary parameter not found in file
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c  05-07-96  MM, change inputs from CAS numbers to icomp, jcomp
c  05-08-96  MM, add /MXINFO/ and read corresponding info from file
c  11-04-96  MM, add gamma to LJi and LMi models (i,j not eq to j,i),
c                change nmxpar from 4 to 6
c  11-19-96  MM, move initialization of common blocks to RDBNC, ESTBNC
c                read critical line coefficients
c  02-19-97  MM, fix bug: parameter name is ncvc, not ncvv
c  05-22-97  MM, change format to read coeff from f10.x to f12.x
c  06-04-97 EWL, change value of nctc and ncvc to read in six coefficients
c  11-28-97 APP, add code for Ammonia-Water model and common /AMMWAT/
c  12-22-98 EWL, remove common /AMMWAT/
c  05-26-00 EWL, make the i,j and j,i entries in the fvcij and ftcij identical
c  12-05-00  MM, set either 7-11 or 2-6 of ftcij and fvcij to zero's (allow for
c                specified mixture order to be opposite that in coeff file)
c  01-23-02 EWL, split common block MXINFO into 2 pieces for some compilers
c  07-15-04 EWL, read in kij for Peng-Robinson equation of state
c  10-08-04 MLH, initialize PR coefficients to zero
c  08-17-06 EWL, rodo logic to make more efficient and allow nonformatted inputs
c  08-17-06 EWL, update program for GERG-2004 equation of Kunz and Wagner
c  12-24-06 MLH, allow for 4 binary interaction parameters for viscosity
c  01-09-07 EWL, modify the code so that other entries than just 'KW#' or 'LJ#'
c                can be on the same line as the CAS numbers (thus, the KW or LJ
c                entries do not have to exist)
c  01-28-07 MLH, allow for 3 binary interaction parameters for thermal conductivity
c  03-05-07 EWL, add fmix2 array in case fluids are entered in opposite order as required by Kunz and Wagner
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      parameter (nbin=ncmax*(ncmax-1)/2)   !# possible binary pairs
      parameter (nmxpar=6)       !number of binary mixture parameters
      parameter (nmxtrm=15)      !number of terms in binary mixing func
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines
      character*1 htab,hnull,hstar,h1
      character*3 hmodmx,hmodij,hmodt,hmodtc,hmodv,hmodvc,hprkij
      character*3 hmodc,hmodz
      character*3 hflag,htype
      character*12 hcas,hcas1,hcas2,hcasi,hcasj
      character*255 herr,hbinp,hbin,h255,hmdijx,hmodtx,hmodvx
      character*255 hmodcx,hmodzx
c
      common /HCHAR/ htab,hnull
      common /CCAS/ hcas(n0:nx)
      common /MXINF1/ hbin(nbin)
      common /MIXMOD/ hmodmx(nx,nx)
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
      common /CF2HMX/ fmix2(nx,nx,nmxpar)
      common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx)
      common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc)
      common /PRMOD/ hprkij(nx,nx)
      common /CFXPR/ fprkij(nx,nx)
      common /TRNBIN/ xljs(nx,nx),xlje(nx,nx),xkij(nx,nx),xlij(nx,nx),
     &                xaji(nx,nx),xkijk(nx,nx),xlijk(nx,nx)
c
      ierr=0
      herr=' '
      i=icomp
      j=jcomp
      hcasi=hcas(icomp)
      hcasj=hcas(jcomp)
      ifmt=27   !Starting location of all LJ1, KW1, etc., identifiers
c
c  initialize arrays storing coeff for critical lines
      fprkij(i,j)=0.0d0
      fprkij(j,i)=0.0d0
      xljs(i,j)=0.0d0
      xljs(j,i)=0.0d0
      xlje(i,j)=0.0d0
      xlje(j,i)=0.0d0
      xkij(i,j)=0.0d0
      xkij(j,i)=0.0d0
      xlij(i,j)=0.0d0
      xlij(j,i)=0.0d0
      xaji(i,j)=0.0d0
      xaji(j,i)=0.0d0
      xkijk(i,j)=0.0d0
      xkijk(j,i)=0.0d0
      xlijk(i,j)=0.0d0
      xlijk(j,i)=0.0d0
      do k=1,nmxpar
        fmix(i,j,k)=0.0d0
        fmix(j,i,k)=0.0d0
        fmix2(i,j,k)=0.0d0
        fmix2(j,i,k)=0.0d0
      enddo
      do k=1,nmxtc
        ftcij(i,j,k)=0.0d0
        ftcij(j,i,k)=0.0d0
      enddo
      do k=1,nmxvc
        fvcij(i,j,k)=0.0d0
        fvcij(j,i,k)=0.0d0
      enddo
c
      if (nread.le.0) then
c  get coefficients from common blocks--not currently implemented
        ierr=1
        herr='[RDBNC error] Block data option not available'//hnull
        call ERRMSG (ierr,herr)
        RETURN
      else
c
c  rewind file and read coefficients
c
        rewind (nread)
c
c  search for key characters in cols 1 and 2-4 and branch accordingly
c  hstar = '#' or '@' indicates start of model specification
c  hflag indicates type of model:
c       'BNC' = binary mixture coefficients
c       'MXM' = mixture model specification
c
        do iline=1,20000
        read (nread,2013,end=199,err=199) hstar,hflag
        if (hstar.eq.'@' .and. (hflag.eq.'END'.or.hflag.eq.'end')) then
          ierr=-117
          herr='[RDBNC error] Binary pair not found (end of file):  '//
     &         hcasi//hcasj//hnull
          call ERRMSG (ierr,herr)
          RETURN
        end if
        if (hstar.eq.'#' .or. hstar.eq.'@') then
          read (nread,2003) htype
 100      continue
          read (nread,2001) h1
          if (h1.eq.'?') goto 100   !comment block terminated by '!'
        end if
        if (hflag.eq.'BNC' .and. htype.eq.'BNC') then
 110      continue
            hmodij=' '
            hmdijx=' '
            hmodt =' '
            hmodtx=' '
            hmodv =' '
            hmodvx=' '
            hmodc =' '
            hmodcx=' '
            hmodz =' '
            hmodzx=' '
            hbinp =' '
            hcas1 =' '
            hcas2 =' '
c  read in one block of interaction parameters at a time
 115        continue
            read (nread,2255,err=199,end=199) h255
            if (h255(1:1).eq.'?') then
              if (hbinp.eq.' ') hbinp=h255(2:255)
              goto 115
            elseif (h255(1:1).eq.'!') then
            elseif (h255(ifmt:ifmt+1).ne.' ') then
              if (h255(1:ifmt-1).ne.' ') then
                ii=index(h255(1:ifmt-1),'/')
                if (ii.gt.0) then
                  hcas1=h255(1:ii-1)
                  hcas2=h255(ii+1:ifmt-1)
                else
                  goto 199
                endif
              endif
              if (h255(ifmt:ifmt+1).eq.'KW' .or.
     &            h255(ifmt:ifmt+1).eq.'LJ' .or.
     &            h255(ifmt:ifmt+1).eq.'LM' .or.
     &            h255(ifmt:ifmt+2).eq.'TR1') then
                hmodij=h255(ifmt:ifmt+2)
                hmdijx=h255(ifmt+3:255)
              elseif (h255(ifmt:ifmt+1).eq.'TC') then
                hmodt =h255(ifmt:ifmt+2)
                hmodtx=h255(ifmt+3:255)
              elseif  (h255(ifmt:ifmt+1).eq.'VC') then
                hmodv =h255(ifmt:ifmt+2)
                hmodvx=h255(ifmt+3:255)
              elseif  (h255(ifmt:ifmt+1).eq.'PR') then
                hmodc =h255(ifmt:ifmt+2)
                hmodcx=h255(ifmt+3:255)
              elseif  (h255(ifmt:ifmt+2).eq.'TRN') then
                hmodz =h255(ifmt:ifmt+2)
                hmodzx=h255(ifmt+3:255)
              else
                goto 199
              endif
              goto 115
            else
              goto 120
            endif

c  check for match with CAS numbers
            if ((hcasi.eq.hcas1 .and. hcasj.eq.hcas2) .or.
     &          (hcasi.eq.hcas2 .and. hcasj.eq.hcas1)) then

              hbin(ibin)=hbinp
              hmodmx(i,j)=hmodij
              hmodmx(j,i)=hmodij      !pair j,i is same as i,j
              hmodtc(i,j)=hmodt
              hmodtc(j,i)=hmodt       !pair j,i is same as i,j
              hmodvc(i,j)=hmodv
              hmodvc(j,i)=hmodv       !pair j,i is same as i,j
              hprkij(i,j)=hmodc
              hprkij(j,i)=hmodc
              if (hmodij.ne.' ') then
                read (hmdijx,*) (fmix(i,j,k),k=1,nmxpar)
                read (hmdijx,*) (fmix(j,i,k),k=1,nmxpar)
                do k=1,nmxpar
                  fmix2(i,j,k)=1.d0  !These are not in the hmx.bnc file, manually set here
                  fmix2(j,i,k)=1.d0
                enddo
              else
               ierr=-118
               herr='[RDBNC error] Binary pair found, but equilibrium'//
     &              ' interaction parameters are not available:  '//
     &              hcasi//hcasj//hnull
               call ERRMSG (ierr,herr)
              endif
c             write (*,1016) icomp,jcomp,hmodij,(a(k),k=1,nmxpar)
c1016         format (1x,' RDBNC--i,j,mod,fij: ',9x,2i3,2x,a3,8f8.4)
c
c  check for specification that critical line(s) are to be estimated,
c  if so, set warning (estimation routine will be called from SETHMX)
              ierr2=0
              if (hmodt.eq.'EST') then
                ierr2=ierr2-1
              end if
              if (hmodv.eq.'EST') then
                ierr2=ierr2-2
              end if
              if (ierr2.lt.0) then
                ierr=-113+ierr2
              end if
c
c  for model TC1, the order i,j vs j,i is significant; first coeff is
c  symmetric, elements 2-6 store coeff for order i,j; 7-11 for j,i
              if (hmodt.eq.'TC1') then
                read (hmodtx,*) (ftcij(i,j,k),k=1,6)
                ftcij(j,i,1)=ftcij(i,j,1)
                do k=2,6
                  ftcij(j,i,k)=ftcij(i,j,k)
                  ftcij(i,j,k+5)=0.0d0
                  ftcij(j,i,k+5)=0.0d0
                enddo
              end if
c
c  for model VC1, the order i,j vs j,i is significant; first coeff is
c  symmetric, elements 2-6 store coeff for order i,j; 7-11 for j,i
              if (hmodv.eq.'VC1') then
                read (hmodvx,*) (fvcij(i,j,k),k=1,6)
                fvcij(j,i,1)=fvcij(i,j,1)
                do k=2,6
                  fvcij(j,i,k)=fvcij(i,j,k)
                  fvcij(i,j,k+5)=0.0d0
                  fvcij(j,i,k+5)=0.0d0
                enddo
              end if
c
c  read in Peng-Robinson interaction parameters
              if (hmodc.eq.'PR') then
                read (hmodcx,*) fprkij(i,j)   !only uses the first number for now
                fprkij(j,i)=fprkij(i,j)
              endif
c
c  read in transport interaction parameters
              if (hmodz.eq.'TRN') then
                read (hmodzx,*) xljs(i,j),xlje(i,j),xkij(i,j),xlij(i,j),
     &                          xaji(i,j),xkijk(i,j),xlijk(i,j)
                xljs(j,i)=xljs(i,j)
                xlje(j,i)=xlje(i,j)
                xkij(j,i)=xkij(i,j)
                xlij(j,i)=xlij(i,j)
                xaji(j,i)=xaji(i,j)
                xkijk(j,i)=xkijk(i,j)
                xlijk(j,i)=xlijk(i,j)
              end if


              if (hcasi.eq.hcas2 .and. hcasj.eq.hcas1) then
c  specified order is opposite that in file hfmix
c  pair j,i is same as i,j except for beta and gamma parameters in
c  the LJi and LMi models are associated with a particular fluid
                if (hmodij(1:2).eq.'LM' .or. hmodij(1:2).eq.'LJ') then
                  a4=fmix(i,j,4)
                  a5=fmix(i,j,5)
                  fmix(j,i,4)=a5
                  fmix(j,i,5)=a4
                  fmix(i,j,4)=a5
                  fmix(i,j,5)=a4
                end if
                if (hmodij(1:2).eq.'KW') then
                  fmix2(j,i,1)=fmix(j,i,1)
                  fmix2(j,i,3)=fmix(j,i,3)
                  fmix2(i,j,1)=fmix(i,j,1)
                  fmix2(i,j,3)=fmix(i,j,3)
                  fmix(j,i,1)=1.d0
                  fmix(j,i,3)=1.d0
                  fmix(i,j,1)=1.d0
                  fmix(i,j,3)=1.d0
                end if
                if (hmodt.eq.'TC1') then       !Similar for TC1
                  do k=2,6
                    ftcij(i,j,k+5)=ftcij(i,j,k)
                    ftcij(j,i,k+5)=ftcij(j,i,k)
                    ftcij(i,j,k)=0.d0
                    ftcij(j,i,k)=0.d0
                  enddo
                end if
                if (hmodv.eq.'VC1') then       !Similar for VC1
                  do k=2,6
                    fvcij(i,j,k+5)=fvcij(i,j,k)
                    fvcij(j,i,k+5)=fvcij(j,i,k)
                    fvcij(i,j,k)=0.0d0
                    fvcij(j,i,k)=0.0d0
                  enddo
                end if
              end if
              RETURN
            endif
          goto 110
        end if
 120    continue
        enddo
      end if
      RETURN
c
 199  continue       !abnormal termination
*     write (*,*) ' RDBNC line 199, hstar,hflag: ',hstar,hflag
      ierr=118
      herr='[RDBNC error 118] Error in reading mixture file in RDBNC'
     &    //hnull
      call ERRMSG (ierr,herr)
      RETURN
c
 2001 format (a1)
 2003 format (a3)
 2013 format (a1,a3)
 2255 format (a255)
c
      end                                              !subroutine RDBNC
c
c ======================================================================
c
      subroutine RDMIX (nread,icomp,jcomp,hmodij,lij,ierr,herr)
c
c  read coefficients for mixing rule(s) from file
c
c  inputs:
c    nread--file to read data from
c           <= 0 get data from block data (not currently implemented)
c           >0 read from logical unit nread (file should have already
c              been opened by subroutine SETHMX)
c    icomp--component i; ignored if hmodij='ALL'
c    jcomp--component j; ignored if hmodij='ALL'
c   hmodij--model for binary interaction of i and j [character*3]
c           'ALL' indicates read model(s) for all binary pair(s)
c                 (used when this routine is called from SETHMX)
c           when a specific model is specified (e.g. 'LJ1' or 'MKT')
c                 read only that model (used when called from SETKTV)
c  outputs:
c      lij--flags indicating if mixing term found for binary pair (i,j)
c           [logical array of dimension (nc,nc)]; used only in
c           conjunction with hmodij='ALL'
c     ierr--error flag:  0 = successful
c                      118 = error in reading mixture file
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  05-07-96  MM, original version (extracted from SETHMX)
c  05-08-96  MM, add /MXINFO/ and read corresponding info from file
c  05-10-96  MM, read default binary parameters (used only by interface)
c  11-04-96  MM, change nmxpar from 4 to 6
c  12-00-00 EWL, add bideal variable, storing fideal
c  02-14-01 EWL, increase dimension on aterm to 4, change use of iexp to aterm
c  01-23-02 EWL, split common block MXINFO into 2 pieces for some compilers
c  07-08-02 EWL, initialize 4th term in amix to 0
c  08-17-06 EWL, update program for GERG-2004 equation of Kunz and Wagner
c  01-03-07 EWL, change the 2000 to 20000 in 'do iline=1,2000' to read in more lines
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nx=ncmax)
      parameter (nbrule=20)    !number of binary mixing rules
      parameter (nbin=ncmax*(ncmax-1)/2)   !# possible binary pairs
      parameter (nmxpar=6)     !number of binary mixture parameters
      parameter (nmxtrm=15)    !number of terms in binary mixing func
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      character*1 htab,hnull
      character*1 hstar,h1
      character*3 hmodmx,hmodij
      character*3 hflag,htype
      character*8 hpar(nmxpar)
      character*251 hinfo
      character*255 herr
c  next 3 declarations associated with /MXINFO/
      character*8 hbpar
      character*255 hmfile,hrule
c
      dimension aterm(nmxtrm,nmxcof),fideal(nmxpar)
      logical lij(nx,nx)
c
      common /NCOMP/ nc,ic
      common /HCHAR/ htab,hnull
c  commons associated with the binary mixing rule(s)
c  the amix(i,j,k,1..3) are the coefs for the (general) mixing term
c  the fmix(i,j,1..nmxpar) are the parameters for the i-j binary
c  the namix(i,j) are the number of a(i,j,k,1..3) terms
      common /MIXMOD/ hmodmx(nx,nx)
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
      common /CFIHMX/ namix(nx,nx),namix2(nx,nx)
c  /MXINFO/ contains information on the mixing rules and parameters
c  hmfile specifies the files from which mixing rules originate
c  (stored in order 1,2  1,3  2,3; element zero is file called in SETUP)
c  hbin provides documentation for the current binary parameters
c  hrule contains descriptions of the currently available mixing rules
c  hbpar contains descriptions of the binary parameters (e.g. Kt, Kv)
c  associated with the currently available mixing rules
      common /MXINFO/ hmfile(0:nbin),hrule(nbrule),hbpar(nbrule,nmxpar)
      common /MXINF2/ bideal(nbrule,nmxpar)
      common /MXRULE/ nrule
c
c  search for key characters in cols 1 and 2-4 and branch accordingly
c  hstar = '#' or '@' indicates start of model specification
c  hflag indicates type of model:
c       'BNC' = binary mixture coefficients
c       'MXM' = mixture model specification
c
      ierr=0
      do i=1,nmxtrm
        do j=1,nmxcof
          aterm(i,j)=0
        enddo
      enddo
      if (hmodij.eq.'ALL') then
c  initialize flags indicating if mixing terms have been found
        irule=0
        do i=1,nc-1
          do j=i+1,nc
            lij(i,j)=.false.
          enddo
        enddo
      else
        do k=1,nrule
          if (hmodij.eq.hrule(k)(1:3)) then
            irule=k
            goto 100
          end if
        enddo
        irule=nrule+1
c        write (*,*) ' RDHMX--new mixing rule found',hmodij
 100  continue
      end if
c  rewind file (RDBNC may have left it at the end)
      rewind (nread)
      do iline=1,20000
      read (nread,2013,end=998,err=998) hstar,hflag
      if ((hstar.eq.'#' .or. hstar.eq.'@').and. hflag.eq.'MXM') then
        if (hmodij.eq.'ALL') then
          irule=irule+1
        end if
        read (nread,2053) htype,hinfo
        do k=1,1000             !skip over source comments
          read (nread,2001) h1
          if (h1.ne.'?') goto 110     !comment block terminated by '!'
        enddo
 110    continue
c  read descriptors for binary-specific parameters
        read (nread,2048) (hpar(k),k=1,nmxpar)
c  read default (ideal-solution) parameters
        read (nread,*) (fideal(k),k=1,nmxpar)
c  read coefficients for mixing rule
        nterm=0
        ncoef=0
        nspare=0
        nterm2=0
        ncoef2=0
        read (nread,*,err=120) nterm,ncoef,nspare,nterm2,ncoef2
 120    continue
        if (nspare.eq.1) ncoef=ncoef+1  !Compatibility for version 7.0
c  nterm is number of terms
c  ncoef is number of coefficients per term
        do k=1,nterm
          aterm(k,4)=0.d0
          read (nread,*) (aterm(k,n),n=1,ncoef)
        enddo
        if (nterm2.ne.0) then
          do k=1,nterm2
            read (nread,*) (aterm(k+nterm,n),n=1,ncoef2)
          enddo
        endif
        if (hmodij.eq.'ALL') then
          ipair=0               !counter for binary pair
c         write (*,1342) irule,htype,hinfo
c1342     format (1x,' SETHMX--mixing function',i3,' (',a3,'): ',a251)
          hrule(irule)=htype//hinfo
          do k=1,nmxpar
            hbpar(irule,k)=hpar(k)
            bideal(irule,k)=fideal(k)
          enddo
          do i=1,nc-1
            do j=i+1,nc
              if (htype.eq.hmodmx(i,j)) then
c  mixing term found; load coefficients into working arrays
                lij(i,j)=.true.
                ipair=ipair+1
c               write (*,1344) irule,htype,i,j,hrule(irule)
c1344           format (1x,' SETHMX:  use mixing rule #',i3,' (',a3,
c    &                 ') for i,j =',2i3,':  ',a255)
                namix(i,j)=nterm
                namix2(i,j)=nterm2
                do k=1,nterm+nterm2
                  do n=1,nmxcof
                    amix(i,j,k,n)=aterm(k,n)
                  enddo
                enddo
              end if
            enddo
          enddo
          nrule=irule
        else if (htype.eq.hmodij) then
c  mixing term found; load coefficients into working arrays
          lij(icomp,jcomp)=.true.
c          write (*,1370) htype,icomp,jcomp
c1370     format (1x,' SETKTV:  use ',a3,' mixing term for i,j =',2i3)
c          write (*,*) ' SETHMX--mixing term irule = ',irule,htype,hinfo
          hrule(irule)=htype//hinfo
          do k=1,nmxpar
            hbpar(irule,k)=hpar(k)
            bideal(irule,k)=fideal(k)
          enddo
          namix(icomp,jcomp)=nterm
          namix(jcomp,icomp)=nterm          !pair (j,i) is same as (i,j)
          namix2(icomp,jcomp)=nterm2
          namix2(jcomp,icomp)=nterm2        !pair (j,i) is same as (i,j)
          do k=1,nterm+nterm2
            do n=1,nmxcof
              amix(icomp,jcomp,k,n)=aterm(k,n)
              amix(jcomp,icomp,k,n)=aterm(k,n)!pair (j,i) is same as (i,j)
            enddo
          enddo
        end if
      else if (hflag.eq.'END') then
c       write (*,*) ' RDMIX end of file encountered'
        RETURN
      end if
      enddo
      RETURN
c
c  error terminations
c
 998  continue
*     write (*,*) ' SETHMX line 998, hstar,hflag: ',hstar,hflag
      ierr=118
      herr='[RDMIX error 118] Error in reading mixture file in RDMIX'
     &    //hnull
      call ERRMSG (ierr,herr)
      RETURN
c
 2001 format (a1)
 2013 format (a1,a3)
c  multiplier in format should match "parameter (nmxpar=6)"
 2048 format (6(a8,1x))     !descriptors for binary-specific parameters
 2053 format (a3,a251)
c
      end                                              !subroutine RDMIX
c
c ======================================================================
c
      subroutine ESTBNC (i,j,ibin,ierr,herr)
c
c  estimate binary parameters for the HMX mixture model
c
c  inputs:
c        i--component i
c        j--component j
c     ibin--counter for binary pair (1,2 = 1; 1,3 = 2; 2,3 = 3; etc.)
c  outputs:
c     ierr--error flag:  0 = successful
c                     -117 = no binary parameters found (warning)
c                      117 = no binary parameters found (critical error)
c     herr--error string (character*255 variable if ierr<>0)
c     other outputs written to common
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  01-12-96  MM, original version
c  02-27-96  MM, parameter n0=-ncmax to accommodate ECS-thermo model
c                add Zcrit to common /CCON/
c  03-19-19  MM, add dipole moment to /CCON/
c  05-06-96  MM, change default model to VDW
c  05-08-96  MM, add /MXINFO/ and load corresponding info
c  11-04-96  MM, replace van der Waals rule with linear (LIN),
c                change nmxpar from 4 to 6
c  11-19-96  MM, move initialization of common blocks to RDBNC, ESTBNC
c                change argument list
c  07-03-97  MM, return ierr=117 if this routine is called (causes error
c                alert to appear in GUI); !temporary until prediction in place
c  09-08-97  MM, change error 117 to warning--GUI displayed alert, then quit
c  11-06-97  MM, implement preliminary predictive model of Lemmon
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  04-19-99 EWL, do not use estimation scheme with hydrogen sulfide
c  06-14-99 EWL, do not use estimation scheme if estimate value > 200
c  01-23-02 EWL, split common block MXINFO into 2 pieces for some compilers
c  02-05-02 EWL, allow estimation scheme to return values, even if out of range
c  05-09-05 EWL, add check for very small acentric factor ratio
c  09-07-06 EWL, change LJ1 to LJ6
c  09-12-06 EWL, keep fmix(1) bounds within 500 to avoid crashes
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      parameter (nmxpar=6)       !number of binary mixture parameters
      parameter (nbin=ncmax*(ncmax-1)/2)   !# possible binary pairs
      parameter (nmxtrm=15)      !number of terms in binary mixing rule
      parameter (nmxcof=7)     !number of coefficients in binary mixing rule
      character*1 htab,hnull
      character*3 hmodmx
      character*12 hcas
      character*255 herr
c  next 3 declarations associated with /MXINFO/
      character*255 hbin
      logical lest
c
      common /HCHAR/ htab,hnull
      common /CCAS/ hcas(n0:nx)
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(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  commons associated with the binary mixing rule(s)
c  the amix(i,j,k,1..3) are the coefs for the (general) mixing term
c  the fmix(i,j,1..nmxpar) are the parameters for the i-j binary
      common /MIXMOD/ hmodmx(nx,nx)
      common /CFXHMX/ amix(nx,nx,nmxtrm,nmxcof),fmix(nx,nx,nmxpar)
      common /MXINF1/ hbin(nbin)
c
*     write (*,1002) i,j,hcas(i),hcas(j)
*1002 format (1x,' ESTBNC will estimate i,j =   ',2i3,2x,a12,1x,a12)
      lest=.true.   !temporary--enable estimation scheme
      if (.not.lest) then
c  return default values for now
        hmodmx(i,j)='LIN'
        hmodmx(j,i)='LIN'
        hbin(ibin)='No mixture data are available for this binary pair.'
     &         //' Ideal solution behavior will be assumed.'
     &         //hnull
        do k=1,2
          fmix(i,j,k)=1.0d0           !fmix(i,j,1) = Kt
          fmix(j,i,k)=1.0d0           !fmix(i,j,2) = Kv
        enddo
        fmix(i,j,3)=0.0d0           !Fpq(i,j)
        fmix(j,i,3)=0.0d0           !Fpq(k,i)
        do k=4,5
          fmix(i,j,k)=1.0d0           !fmix(i,j,4) = beta
          fmix(j,i,k)=1.0d0           !fmix(i,j,5) = gamma
        enddo
        do k=6,nmxpar
          fmix(i,j,k)=0.0d0           !not used
          fmix(j,i,k)=0.0d0
        enddo
c  return warning message to GUI
        ierr=117
        write (herr,1000) hnull
 1000   format('[SETUP error 117] No mixture data are available for',
     &   ' one or more binary pairs in the specified mixture. The',
     &   ' mixture is outside the range of the model and',
     &   ' calculations will not be made.',A1)
        call ERRMSG (ierr,herr)
        RETURN
      end if
      iflg=0
      do k=1,2
        n=0
        if (k.eq.1) n=i
        if (k.eq.2) n=j
        if (hcas(n).eq.'7732-18-5' ) iflg=1    !water
        if (hcas(n).eq.'7789-20-0' ) iflg=1    !heavy water
        if (hcas(n).eq.'7664-41-7' ) iflg=1    !ammonia
        if (hcas(n).eq.'7440-59-7' ) iflg=1    !helium
        if (hcas(n).eq.'1333-74-0' ) iflg=1    !hydrogen
        if (hcas(n).eq.'7782-39-0' ) iflg=1    !deuterium
        if (hcas(n).eq.'1333-74-0p') iflg=1    !parahydrogen
        if (hcas(n).eq.'67-56-1'   ) iflg=1    !methanol
        if (hcas(n).eq.'64-17-5'   ) iflg=1    !ethanol
        if (hcas(n).eq.'67-64-1'   ) iflg=1    !acetone
        if (hcas(n).eq.'7782-41-4' ) iflg=1    !fluorine
        if (hcas(n).eq.'7440-01-9' ) iflg=1    !neon
        if (hcas(n).eq.'7783-54-2' ) iflg=1    !NF3
      enddo
      if (iflg.eq.1) then
        ierr=117
        herr='[SETUP error 117] Estimation of mixing parameters is not '
     &    //'possible for most mixtures containing water, heavy water, '
     &    //'helium, hydrogen, parahydrogen, deuterium, neon, ammonia, '
     &    //'acetone, fluorine, methanol, or ethanol.'//hnull
        hbin(ibin)='Estimation of mixing parameters is not possible '
     &    //'for this mixture.'//hnull
        hmodmx(i,j)='LJ6'
        hmodmx(j,i)='LJ6'
        call ERRMSG (ierr,herr)
        RETURN
      endif
c
c  return warning message to GUI
      ierr=-117
      herr='[SETUP warning -117] No mixture data are available for one '
     &   //'or more binary pairs in the specified mixture. The mixing '
     &   //'parameters have been estimated.'//hnull
      call ERRMSG (ierr,herr)
c
c  ordering of the components depends on ratio of the dipole moments
c  interim predictive scheme of Lemmon
      hmodmx(i,j)='LJ6'
      hmodmx(j,i)='LJ6'
      hbin(ibin)='No mixture data are available for this binary pair. '
     &         //'The mixing parameters have been estimated.'
     &         //hnull
      k=0
      if (ABS(dipole(i)-dipole(j)).lt.1.0d-5) then
c  special case for equal dipole moments
        if (tc(i)/(pc(i)*accen(i)).gt.tc(j)/(pc(j)*accen(j))) k=1
      else
        if (dipole(i).lt.dipole(j)) k=1
      end if
      if (k.eq.1) then
        tratio=tc(i)/tc(j)
        pratio=pc(i)/pc(j)
        accenr=accen(i)/accen(j)
      else
        tratio=tc(j)/tc(i)
        pratio=pc(j)/pc(i)
        accenr=accen(j)/accen(i)
      endif
      if (tratio/pratio/accenr.gt.5 .and. abs(accenr).lt.0.1d0) then
        if (tratio .gt. 0.75d0 .and. tratio.lt.1.25d0) then
          tratio=1.d0/tratio
          pratio=1.d0/pratio
          accenr=1.d0/accenr
        endif
      endif
c
c  define the zeta parameter (the one applied to the reducing temperature)
      fmix(i,j,1)=(40.4d0-25.03d0*2.0d0**(tratio/pratio/accenr))/tratio
c
c
c     c1=-4.47546077875894
c     c2=6.95801916004273
c     c3=0.301565763059889
c     c4=0.429854120106733
c     c5=0.869600175257425
c     xx=(tc(j)-tc(i))*(tc(i)/100.d0)**c4/100.d0
c     if (tc(i).gt.tc(j))
c    &  xx=(tc(i)-tc(j))*(tc(j)/100.d0)**c4/100.d0
c     fmix(i,j,1)=c1+c2*xx+c3*xx**2+c5*xx**3
c
c
c
c
      fmix(j,i,1)=fmix(i,j,1)
c  set remaining parameters to defaults; 2 = xi; 3 = Fpq; 4 = beta; 5 = gamma
      do k=2,3
        fmix(i,j,k)=0.0d0  !k = 2:  xi (term in reducing volume)
        fmix(j,i,k)=0.0d0  !k = 3:  Fpq (multiplier for excess function)
      enddo
      do k=4,5
        fmix(i,j,k)=1.0d0  !k = 4: beta (exponent for first component)
        fmix(j,i,k)=1.0d0  !k = 5: gamma (exponent for second component)
      enddo
      do k=6,nmxpar
        fmix(i,j,k)=0.0d0  !not used
        fmix(j,i,k)=0.0d0
      enddo
      if (ABS(fmix(i,j,1)).gt.200.0d0) then
        if (fmix(i,j,1).gt.500.d0) fmix(i,j,1)=500.d0
        if (fmix(i,j,1).lt.-500.d0) fmix(i,j,1)=-500.d0
        ierr=117
        write (herr,1000) hnull
        call ERRMSG (ierr,herr)
      endif
      RETURN
c
      end                                             !subroutine ESTBNC
c
c ======================================================================
c
      subroutine ESTTC (i,j,ierr,herr)
c
c  estimate parameters for binary critical temperature line
c
c  inputs:
c        i--component i
c        j--component j
c  outputs:
c     ierr--error flag:  0 = successful
c                        1 = error
c     herr--error string (character*255 variable if ierr<>0)
c     other outputs written to common
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-19-96  MM, original version
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nx=ncmax)
      parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines
      character*1 htab,hnull
      character*3 hmodtc,hmodvc,hprkij
      character*255 herr
c
      common /HCHAR/ htab,hnull
c  commons associated with the binary critical lines
      common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx)
      common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc)
      common /PRMOD/ hprkij(nx,nx)
c
*     write (*,1002) i,j
*1002 format (1x,' ESTTC will estimate i,j =    ',2i3)
c  return default values for now
      hmodtc(i,j)='TC1'
      hmodtc(j,i)='TC1'
      hprkij(j,i)='   '
      do k=1,nmxtc
        ftcij(i,j,k)=0.0d0
        ftcij(j,i,k)=0.0d0
      enddo
      ierr=0
      herr=' '
      RETURN
c
      end                                              !subroutine ESTTC
c
c ======================================================================
c
      subroutine ESTVC (i,j,ierr,herr)
c
c  estimate parameters for binary critical volume line
c
c  inputs:
c        i--component i
c        j--component j
c  outputs:
c     ierr--error flag:  0 = successful
c                        1 = error
c     herr--error string (character*255 variable if ierr<>0)
c     other outputs written to common
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  11-19-96  MM, original version
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nx=ncmax)
      parameter (nmxtc=11,nmxvc=11) !number of coeff for critical lines
      character*1 htab,hnull
      character*3 hmodtc,hmodvc
      character*255 herr
c
      common /HCHAR/ htab,hnull
c  commons associated with the binary critical lines
      common /CRTMOD/ hmodtc(nx,nx),hmodvc(nx,nx)
      common /CFXCRT/ ftcij(nx,nx,nmxtc),fvcij(nx,nx,nmxvc)
c
*     write (*,1002) i,j
*1002 format (1x,' ESTVC will estimate i,j =    ',2i3)
c  return default values for now
      hmodvc(i,j)='VC1'
      hmodvc(j,i)='VC1'
      do k=1,nmxvc
        fvcij(i,j,k)=0.0d0
        fvcij(j,i,k)=0.0d0
      enddo
      ierr=0
      herr=' '
      RETURN
c
      end                                              !subroutine ESTVC
c ======================================================================
c
      subroutine AMH2OR
c
c  check for ammonia-water mixture, if set, use international standard settings
c  check that the user has not modified the reference state for their own use
c
      implicit double precision (a-h,o-z)
      implicit integer (i-k,m,n)
      implicit logical (l)
      parameter (ncmax=20)        !max number of components in mixture
      parameter (nrefmx=10)       !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx,nx=ncmax)
      character*3 hrf
      character*12 hcas
      character*255 herr
      dimension x0(ncmax)
      common /NCOMP/ nc,ic
      common /CREF/ tref(n0:nx),rhoref(n0:nx),href(n0:nx),sref(n0:nx)
      common /CREFDF/ tdef(n0:nx),pdef(n0:nx),hdef(n0:nx),sdef(n0:nx)
      common /WRDCPP/ tred(n0:nx),Cred(n0:nx)
      common /CCON/ tc(n0:nx),pc(n0:nx),rhoc(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 /CCAS/ hcas(n0:nx)
      if (nc.eq.2) then
        i=0
        if (hcas(1).eq.'7732-18-5'.and. hcas(2).eq.'7664-41-7') i=1
        if (hcas(2).eq.'7732-18-5'.and. hcas(1).eq.'7664-41-7') i=2
        if (i.ne.0) then
          if (abs(tref(i)-300.) .lt. 0.001 .and.
     &    abs(hdef(i)-45957.).lt.1. .and. abs(sdef(i)-164.).lt.1. .and.
     &    abs(hdef(3-i)-28945.).lt.1. .and. abs(sdef(3-i)-159.).lt.1.)
     &    then
            reos(i)=8.314471d0
            cred(i)=8.314471d0
            hdef(i)=45957.744019720d0
            sdef(i)=164.007259009d0
            hdef(3-i)=28945.5688500127d0
            sdef(3-i)=159.336684302313d0
c           hsvrfs=' '
            ixflag=1
            h0=0.d0
            s0=0.d0
            t0=0.d0
            p0=0.d0
            hrf='DEF'
            call SETREF (hrf,ixflag,x0,h0,s0,t0,p0,ierr,herr)
          endif
        endif
      endif
      end                                             !subroutine AMH2OR
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                     end file mix_HMX.f
c ======================================================================
