c  begin file trns_VIS.f
c
c  This file contains the core routines for viscosity
c
c  contained here are:
c     subroutine SETVS0 (nread,icomp,hcasno,ierr,herr)
c     function ETA0HC (icomp,t,rho,ierr,herr)
c     subroutine SETVS1 (nread,icomp,hcasno,ierr,herr)
c     function ETA1DG (icomp,t)
c     function ETA1B2 (icomp,t)
c     function ETA1RS (icomp,t,rho)
c     subroutine SETVS2 (nread,icomp,hcasno,ierr,herr)
c     function ETA2DG (icomp,t)
c     function ETA2RS (icomp,t,rho)
c     subroutine SETVS3 (nread,icomp,hcasno,ierr,herr)
c     function ETA3DG (icomp,t)
c     function ETA3RS (icomp,t,rho)
c     subroutine SETVS4 (nread,icomp,hcasno,ierr,herr)
c     function ETA4RS (icomp,t,rho)
c     subroutine SETVS6 (nread,icomp,hcasno,ierr,herr)
c     FUNCTION ETAH2(ICOMP,T,D)
c     FUNCTION DELV(icomp,D1,T1,D2,T2)
c     FUNCTION EXVDIL(DD,T)
c     FUNCTION DILV(icomp,T)
c     FUNCTION EXCESV(icomp,DD,T)
c     function ETAHE (icomp,t,rho)
c     function ETAETY (icomp,t,rho)
c     function ETANEO (icomp,t,rho)
c     function ETAR23 (icomp,t,rho)
c     function ETAH2O (icomp,t,rho)
c     function ETAMEO (icomp,t,rho)
c
c =====================================================================
c =====================================================================
c
      subroutine SETVS0 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #0; the model
c  which points to all the hardcoded equations.
c
c  inputs:
c    nread--file to read data from
c           >0 read from logical unit nread (file should have already
c              been opened and pointer set by subroutine SETUP)
c    icomp--component number in mixture (0..nc)
c           1 for pure fluid; 0 for ECS reference fluid
c   hcasno--CAS number of component icomp (not req'd if reading from file)
c
c  outputs:
c     ierr--error flag:  0 = successful
c     herr--error string (character*255 variable if ierr<>0)
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  02-29-00 EWL, original version
c  09-00-00 EWL, change hmdci to hmdeta and hmdtcx to avoid overlapping tcx.
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
      character*3 hetacr,htcxcr
      character*3 hetahc,htcxhc
      character*12 hcasno
      character*255 herr
      character*3 hmdeta,hmdtcx
c
      common /HCHAR/ htab,hnull
c  pointer to critical enhancement auxiliary functions
      common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax)
c  pointer to hardcoded models
      common /HCMOD/ hetahc(nrf0:ncmax),htcxhc(nrf0:ncmax)
c  limits and reducing parameters
      common /WLMETA/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
      common /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
c  commons storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c  Lennard-Jones parameters
      common /WLJETA/ sigma(nrf0:nx),epsk(nrf0:nx)
c  pointer to collision integral model
      common /OMGMOD/ hmdeta(nrf0:nx),hmdtcx(nrf0:nx)
c
c  read data from file (should have been opened by SETUP)
c     write (*,*) ' SETVS0--read component',icomp,' from unit',nread
      write (herr,'(a12)') hcasno  !Use hcasno to avoid warning message
      read (nread,*) tmin(icomp)              !lower temperature limit
      read (nread,*) tmax(icomp)              !upper temperature limit
      read (nread,*) pmax(icomp)              !upper pressure limit
      read (nread,*) rhomax(icomp)            !upper density limit
c
c  read in pointer to the hardcoded model
      read (nread,2003) hetahc(icomp)
c     write (*,*) ' SETVS0--will use model ',hetahc(icomp)
      read (nread,*) ndg(icomp),nB2(icomp),ndel0(icomp),
     &               npoly(icomp),nnum(icomp),nden(icomp),
     &               nexpn(icomp),nexpd(icomp)
c     write (*,*) ' SETVS0--about to read ',ndg(icomp),' dilute terms'
      jterm=0                                 !term counter
      if (ndg(icomp).ge.1) then
        read (nread,2003) hmdeta(icomp)       !pointer to omega model
        read (nread,*) sigma(icomp)           !L-J sigma
        read (nread,*) epsk(icomp)            !L-J epsilon/kappa
        read (nread,*) treddg(icomp),etardg(icomp)  !reducing par
        do j=1,ndg(icomp)                 !read dilute-gas terms
          jterm=jterm+1
          read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
        enddo
      end if
      nrsum=nB2(icomp)+ndel0(icomp)+npoly(icomp)+nnum(icomp)+
     &      nden(icomp)+nexpn(icomp)+nexpd(icomp)
      if (nrsum.ge.1) then
c  read in reducing parameters
        read (nread,*) tred(icomp),Dred(icomp),etared(icomp)
        do j=1,nrsum
          jterm=jterm+1
          read (nread,*) (ceta(icomp,jterm,k),k=1,4),ieta(icomp,jterm)
        enddo
      end if
c
c  read in pointer to critical enhancement model
      read (nread,2003) hetacr(icomp)
c       write (*,*) ' SETVS1--will use critical model ',hetacr(icomp)
      ierr=0
      herr=' '
c
      RETURN
 2003 format (a3)
      end                                             !subroutine SETVS0
c
c ======================================================================
c
      function ETA0HC (icomp,t,rho,ierr,herr)
c
c  model for the hardcoded viscosity models
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c        d--molar density [mol/L]
c  output (as function value):
c   eta0hc--viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  02-29-00 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=20)   !max number of components in mixture
      parameter (nrefmx=10)  !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx)
      parameter (nrf0=n0)    !lower limit for transport ref fluid arrays
      character*3 hetahc,htcxhc
      character*1 htab,hnull
      character*255 herr
      common /HCHAR/ htab,hnull
      common /HCMOD/ hetahc(nrf0:ncmax),htcxhc(nrf0:ncmax)
      ierr=0
      ETA0HC=0.d0
      if (hetahc(icomp).eq.'ETY') THEN
        ETA0HC=ETAETY(icomp,t,rho)
      elseif (hetahc(icomp).eq.'NEO') THEN
        ETA0HC=ETANEO(icomp,t,rho)
      elseif (hetahc(icomp).eq.'R23') THEN
        ETA0HC=ETAR23(icomp,t,rho)
      elseif (hetahc(icomp).eq.'H2O') THEN
        ETA0HC=ETAH2O(icomp,t,rho)
      elseif (hetahc(icomp).eq.'HE') THEN
        ETA0HC=ETAHE(icomp,t,rho)
      elseif (hetahc(icomp).eq.'H2') THEN
        ETA0HC=ETAH2(icomp,t,rho)
      elseif (hetahc(icomp).eq.'MEO') THEN
        ETA0HC=ETAMEO(icomp,t,rho)
      else
        ierr=49
        herr='[ETA0HC error 49] unknown viscosity model specified'
     &        //hnull
        call ERRMSG (ierr,herr)
      endif
      RETURN
      end                                             !subroutine ETA0HC
c
c ======================================================================
c
      subroutine SETVS1 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #1; this, the "composite model,"
c  is written in a general form with terms designed to include several
c  recent correlations including those of Fenghour (1995) for ammonia,
c  Krauss (1996) for R152a, and Laesecke (1997) for R134a.
c
c  inputs:
c    nread--file to read data from
c           <= 0 get data from block data
c           >0 read from logical unit nread (file should have already
c              been opened and pointer set by subroutine SETUP)
c    icomp--component number in mixture (0..nc)
c           1 for pure fluid; 0 for ECS reference fluid
c   hcasno--CAS number of component icomp (not req'd if reading from file)
c
c  outputs:
c     ierr--error flag:  0 = successful
c                      101 = error--block data option not implemented
c     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  01-16-97  MM, original version
c  02-26-97  MM, read pointer for critical enhancement model (future use)
c  08-19-97  MM, change error number for nread<=0; input hcasno is not array
c  09-27-01 MLH, add Laesecke's alternative formulation for del10 term (use neg# terms to flag it)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
      character*3 hetamx,heta,htcxmx,htcx,hmdeta,hmdtcx
      character*3 hetacr,htcxcr
      character*12 hcasno
      character*255 herr
c
      common /HCHAR/ htab,hnull
      common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax)
c  pointer to critical enhancement auxiliary functions
      common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax)
c  pointer to collision integral model
      common /OMGMOD/ hmdeta(nrf0:nx),hmdtcx(nrf0:nx)
c  limits and reducing parameters
      common /WLMETA/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
c  numbers of terms for the various parts of the model:  dilute gas,
c  second viscosity virial (initial density dependence), residual part
      common /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
c  commons storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c  Lennard-Jones parameters
      common /WLJETA/ sigma(nrf0:nx),epsk(nrf0:nx)
c
      if (nread.le.0) then
c  get coefficients from block data--this option not implemented,
c  place holder to maintain parallel structure with EOS setup routines
        ierr=101
        write (herr,1101) nread,hcasno,hnull
        call ERRMSG (ierr,herr)
 1101   format ('[SETVS1 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
c       write (*,*) ' SETVS1--read component',icomp,' from unit',nread
        read (nread,*) tmin(icomp)              !lower temperature limit
        read (nread,*) tmax(icomp)              !upper temperature limit
        read (nread,*) pmax(icomp)              !upper pressure limit
        read (nread,*) rhomax(icomp)            !upper density limit
        jterm=0                                 !term counter
        read (nread,*) ndg(icomp)               !# dilute-gas terms
c       write (*,*) ' SETVS1--about to read ',ndg(icomp),' dilute terms'
        if (ndg(icomp).ge.1) then
          read (nread,2003) hmdeta(icomp)       !pointer to omega model
          read (nread,*) sigma(icomp)           !L-J sigma
          read (nread,*) epsk(icomp)            !L-J epsilon/kappa
          read (nread,*) treddg(icomp),etardg(icomp)  !reducing par
          do j=1,ndg(icomp)                 !read dilute-gas terms
            jterm=jterm+1
            read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
          enddo
        end if
        read (nread,*) nB2(icomp)               !# visc virial terms
c       write (*,*) ' SETVS1--about to read ',nB2(icomp),' virial terms'
        if (nB2(icomp).ge.1) then
          read (nread,*) tredB2(icomp),etarB2(icomp)  !reducing par
          do j=1,nB2(icomp)        !read viscosity virial terms
            jterm=jterm+1
            read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
          enddo
        end if
c
c  read the number of terms of the various parts of the residual model
c  these are in the order:
c    close-packed density function;
c    simple polynomials in T, rho, rho_0, exp(rho/rhoc);
c    numerator of rational polynomial; denominator of rational polynomial;
c    numerator of exponential term; denominator of exponential term;
c  the coefficients themselves are given in the order:
c    constant multiplier; temperature exponent (all terms);
c    density exponent; close-packed density exponent (all except del0 terms)
c    power of density inside exponential (0 indicates no exponential)
c
        read (nread,*) ndel0(icomp),npoly(icomp),nnum(icomp),nden(icomp)
     &                ,nexpn(icomp),nexpd(icomp)
        nrsum=ndel0(icomp)+npoly(icomp)+nnum(icomp)+nden(icomp)
     &       +nexpn(icomp)+nexpd(icomp)
c       write (*,*) ' SETVS1--about to read ',nrsum,' residual terms'
        if (nrsum.ge.1) then
c  read in reducing parameters
          read (nread,*) tred(icomp),Dred(icomp),etared(icomp)
          if (ndel0(icomp).ge.1) then
            do j=1,ndel0(icomp)     !close-packed density term
              jterm=jterm+1
              read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
            enddo
          elseif(ndel0(icomp).lt.0) then
c  close-packed density term; alternative form
            do j=1,ABS(ndel0(icomp))
              jterm=jterm+1
              read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
            enddo
          end if
c
          if (npoly(icomp).ge.1) then
            do j=1,npoly(icomp)     !simple polynomial terms
              jterm=jterm+1
              read(nread,*)(ceta(icomp,jterm,k),k=1,4),ieta(icomp,jterm)
            enddo
          end if
          if (nnum(icomp).ge.1) then
            do j=1,nnum(icomp)      !numerator of rational polynomial
              jterm=jterm+1
              read(nread,*)(ceta(icomp,jterm,k),k=1,4),ieta(icomp,jterm)
            enddo
          end if
          if (nden(icomp).ge.1) then
            do j=1,nden(icomp)      !denominator of rational poly
              jterm=jterm+1
              read(nread,*)(ceta(icomp,jterm,k),k=1,4),ieta(icomp,jterm)
            enddo
          end if
          if (nexpn(icomp).ge.1) then
            do j=1,nexpn(icomp)     !numerator of exponential term
              jterm=jterm+1
              read(nread,*)(ceta(icomp,jterm,k),k=1,3),ieta(icomp,jterm)
            enddo
          end if
          if (nexpd(icomp).ge.1) then
            do j=1,nexpd(icomp)     !denominator of exponential term
              jterm=jterm+1
              read(nread,*)(ceta(icomp,jterm,k),k=1,3),ieta(icomp,jterm)
            enddo
          end if
        end if
c
c  read in pointer to critical enhancement model
        read (nread,2003) hetacr(icomp)
c       write (*,*) ' SETVS1--will use critical model ',hetacr(icomp)
        ierr=0
        herr=' '
      end if
c
      RETURN
 2003 format (a3)
      end                                             !subroutine SETVS1
c
c ======================================================================
c
      function ETA1DG (icomp,t)
c
c  dilute-gas contribution to the viscosity by the composite model (VS1)
c  and also the friction theory model (VS4)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c  output (as function value):
c   eta1dg--the dilute-gas part of the viscosity [uPa-s]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  01-17-97  MM, original version
c  03-28-97  MM, move calc of tau inside "if" (divide by zero if ndg(i)=0)
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c   1-20-00 EWL, add check for absurd t, coming from the ECS model
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
      character*3 hmdeta,hmdtcx
c
      common /HCHAR/ htab,hnull
c  common storing the fluid constants
      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  reducing parameters
      common /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
c  numbers of terms for the various parts of the model:  dilute gas,
c  second viscosity virial (initial density dependence), residual part
      common /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
c  commons storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c  Lennard-Jones parameters
      common /WLJETA/ sigma(nrf0:nx),epsk(nrf0:nx)
      common /OMGMOD/ hmdeta(nrf0:nx),hmdtcx(nrf0:nx)
c
      i=icomp
      nterm=0                                    !term counter
      eta1dg=0.0d0
c
c  sum the dilute-gas terms
      eta1dg=0.0d0
      if (ndg(i).ge.1 .and. t.lt.1.d8) then
        tau=t/treddg(i)
c  first term is always the Chapman-Enskog term
        eta1dg=ceta(i,1,1)*SQRT(tau)/(sigma(i)**2
     &        *OMEGA(i,t,epsk(i),hmdeta(i)))
        if (ndg(i).ge.2) then
c  possibility for additional, empirical terms
          do j=nterm+1,nterm+ndg(i)
            eta1dg=eta1dg+ceta(i,j,1)*tau**ceta(i,j,2)
c
          enddo
        end if
      end if
c
c  multiply by reducing parameter for viscosity (to convert units, etc.)
      eta1dg=eta1dg*etardg(i)
c      write (*,*) ' ETA1DG--d.g. visc:   ',eta1dg
      RETURN
      end                                               !function ETA1DG
c
c ======================================================================
c
      function ETA1B2 (icomp,t)
c
c  second viscosity "virial coefficient" by the composite model (VS1)
c  This model implements the initial-density dependence of the
c  Rainwater-Friend theory.  It returns a viscosity virial coefficient
c  which must be multiplied by the dilute-gas viscosity and the density
c  to yield the initial-density dependence of viscosity.
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c  output (as function value):
c   eta1B2--the second viscosity virial coefficient [L/(mol-uPa-s)]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  01-17-97  MM, original version
c  03-28-97  MM, move calc of tau inside "if" (divide by zero if nB2(i)=0)
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)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /HCHAR/ htab,hnull
c  common storing the fluid constants
      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  reducing parameters
      common /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
c  numbers of terms for the various parts of the model:  dilute gas,
c  second viscosity virial (initial density dependence), residual part
      common /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
c  commons storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c  Lennard-Jones parameters
      common /WLJETA/ sigma(nrf0:nx),epsk(nrf0:nx)
c
c     write (*,*) ' ETA1B2--tred,etared:  ',tredB2(icomp),etarB2(icomp)
      i=icomp
      nterm=ndg(i)                              !term counter
c
c  sum the terms comprising the viscosity virial
      eta1B2=0.0d0
      if (nB2(i).ge.1) then
        tau=t/tredB2(i)
        do j=nterm+1,nterm+nB2(i)
          eta1B2=eta1B2+ceta(i,j,1)*tau**ceta(i,j,2)
c       write (*,*) ' ETA1B2--j,c1,ti:  ',j,ceta(i,j,1),ceta(i,j,2)
        enddo
c  multiply by reducing parameter (to convert units, etc.)
        eta1B2=eta1B2*etarB2(i)
      end if
c     write (*,*) ' ETA1B2--etaB2:  ',eta1B2
c
      RETURN
      end                                               !function ETA1B2
c
c ======================================================================
c
      function ETA1RS (icomp,t,rho)
c
c  residual contribution to the viscosity by the composite model (VS1)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   eta1rs--the background part of the viscosity [uPa-s]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  01-16-97  MM, original version
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c   2-16-99 EWL, return if rho=0
c  09-27-01 MLH, allow Laesecke's alternative formulation for del10 term
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /HCHAR/ htab,hnull
c  common storing the fluid constants
      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  reducing parameters
      common /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
c  numbers of terms for the various parts of the model (dilute gas,
c  initial density dependence, residual part)
      common /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
c  commons storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c
      i=icomp
      eta1rs=0.0d0
      if (rho.lt.1.0d-10) RETURN
c
      nsum=ndel0(i)+npoly(i)+nnum(i)+nden(i)+nexpn(i)+nexpd(i)
      expdel=1.0d0    !initialize only
      tau=1.d0
      del=1.d0
      if (nsum.ge.1) then
c  compute tau,del only if terms present (otherwise reducing par not read in)
        tau=t/tred(i)
        del=rho/Dred(i)
c  define the density to be used in exponential multipliers
c  if the reducing density is 1.0, must divide by the critical density
        if (abs(Dred(i)-1.0d0).lt.0.001) then
          expdel=rho/rhoc(i)
        else
          expdel=del
        end if
      end if
c
c  compute the various parts of the residual model
c  these are taken in the order:
c    close-packed density function;
c    simple polynomials in T, rho, rho_0, exp(rho/rhoc);
c    numerator of rational polynomial; denominator of rational polynomial;
c    numerator of exponential term; denominator of exponential term;
c  the coefficients themselves are given in the order:
c    constant multiplier; temperature exponent (all terms);
c    density exponent; close-packed density exponent (all except del0 terms)
c    power of density inside exponential (0 indicates no exponential)
c
      nterm=ndg(i)+nB2(i)                        !term counter
c  compute the close-packed density at given temperature, if applicable
      if (ndel0(i).ge.1) then
        del0=0.0d0
        do j=nterm+1,nterm+ndel0(i)
          del0=del0+ceta(i,j,1)*tau**ceta(i,j,2)
c       write (*,*) ' ETA1RS--j,close-packed term:  ',j,ceta(i,j,1)
        enddo
        nterm=nterm+ndel0(i)
      elseif (ndel0(i).lt.0) then   !alternative formulation for del10
        del0=1.0d0
        do j=nterm+2,nterm+ABS(ndel0(i))
          del0=del0+ceta(i,j,1)*tau**ceta(i,j,2)
        enddo
        del0=ceta(i,nterm+1,1)/del0
        nterm=nterm+ABS(ndel0(i))
      else
        del0=1.0d0
      end if
c     write (*,*) ' ETA1RS--icomp, # del0, del0:    ',i,ndel0(i),del0
c
c  sum the simple polynomial terms
      eta1rs=0.0d0
      if (npoly(i).ge.1) then
        do j=nterm+1,nterm+npoly(i)
          visci=ceta(i,j,1)*tau**ceta(i,j,2)*del**ceta(i,j,3)
     &         *del0**ceta(i,j,4)
          if (ieta(i,j).ge.1) then
            visci=visci*exp(-expdel**ieta(i,j))
          end if
          eta1rs=eta1rs+visci
c       write (*,1220) j,ceta(i,j,1),eta1rs
c1220   format (1x,' ETA1RS--j,polynomial coeff,sum:',i3,2e14.6)
        enddo
c       write (*,*) ' ETA1RS--j,polynomial sum:      ',j,eta1rs
        nterm=nterm+npoly(i)
      end if
c
c  calculate the numerator of the rational polynomial
      if (nnum(i).ge.1) then
        xnum=0.0d0
        do j=nterm+1,nterm+nnum(i)
          xnum=xnum+ceta(i,j,1)*tau**ceta(i,j,2)*del**ceta(i,j,3)
     &        *del0**ceta(i,j,4)
          if (ieta(i,j).ge.1) then
            xnum=xnum*exp(-expdel**ieta(i,j))
          end if
        enddo
c       write (*,*) ' ETA1RS--numerator term:       ',xnum
        nterm=nterm+nnum(i)
      else
        xnum=1.0d0
      end if
c
c  calculate the denominator of the rational polynomial
      if (nden(i).ge.1) then
        xden=0.0d0
        do j=nterm+1,nterm+nden(i)
          xden=xden+ceta(i,j,1)*tau**ceta(i,j,2)*del**ceta(i,j,3)
     &        *del0**ceta(i,j,4)
          if (ieta(i,j).ge.1) then
            xden=xden*exp(-expdel**ieta(i,j))
          end if
        enddo
        nterm=nterm+nden(i)
c       write (*,*) ' ETA1RS--denominator term:     ',xden
      else
        xden=1.0d0
      end if
c
c  combine the two parts of the rational polynomial, if applicable
      if (nnum(i).ge.1 .or. nden(i).ge.1) then
c       write (*,*) ' ETA1RS--num,denom:  ',xnum,xden
        eta1rs=eta1rs+xnum/xden
      end if
c
c  calculate the numerator of the exponential term
      if (nexpn(i).ge.1) then
        xnum=0.0d0
        do j=nterm+1,nterm+nexpn(i)
          xnum=xnum+ceta(i,j,1)*tau**ceta(i,j,2)*del**ceta(i,j,3)
     &        *del0**ceta(i,j,4)
        enddo
        nterm=nterm+nexpn(i)
      else
        xnum=1.0d0
      end if
c
c  calculate the denominator of the exponential term
      if (nexpd(i).ge.1) then
        xden=0.0d0
        do j=nterm+1,nterm+nexpd(i)
          xden=xden+ceta(i,j,1)*tau**ceta(i,j,2)*del**ceta(i,j,3)
     &        *del0**ceta(i,j,4)
        enddo
        nterm=nterm+nexpd(i)
      else
        xden=1.0d0
      end if
c
c  combine the two parts of the exponential term, if applicable
      if (nexpn(i).ge.1 .or. nexpd(i).ge.1) then
        eta1rs=eta1rs+EXP(xnum/xden)
      end if
c
c  multiply by reducing parameter for viscosity (to convert units, etc.)
      eta1rs=eta1rs*etared(i)
c
      RETURN
      end                                               !function ETA1RS
c
c ======================================================================
c
      subroutine SETVS2 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #2; the hydrocarbon model of:
c  Younglove, B.A. and Ely, J.F. (1987). Thermophysical properties of
c  fluids. II. Methane, ethane, propane, isobutane and normal butane.
c  J. Phys. Chem. Ref. Data  16: 577-798.
c
c  inputs:
c    nread--file to read data from
c           <= 0 get data from block data
c           >0 read from logical unit nread (file should have already
c              been opened and pointer set by subroutine SETUP)
c    icomp--component number in mixture (0..nc)
c           1 for pure fluid; 0 for ECS reference fluid
c   hcasno--CAS number of component icomp (not req'd if reading from file)
c
c  outputs:
c     ierr--error flag:  0 = successful
c                      101 = error--block data option not implemented
c     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  02-21-97  MM, original version
c  02-26-97  MM, set pointer to critical enhancement to 'NUL' (not used)
c  08-19-97  MM, change error number for nread<=0; input hcasno is not array
c  07-01-98 EWL, read ceta(icomp,1,2) and ceta(icomp,e1,1) if version>6.01
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)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
      character*3 hetamx,heta,htcxmx,htcx,hmdeta,hmdtcx
      character*3 hetacr,htcxcr
      character*12 hcasno
      character*255 herr
c
      common /HCHAR/ htab,hnull
      common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax)
c  pointer to critical enhancement auxiliary functions
      common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax)
c  pointer to collision integral model
      common /OMGMOD/ hmdeta(nrf0:nx),hmdtcx(nrf0:nx)
c  limits
      common /WLMETA/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
c  commons storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c  Lennard-Jones parameters
      common /WLJETA/ sigma(nrf0:nx),epsk(nrf0: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 /VERS/ verfl(n0:nx),vermx    !fluid & mix file version nos.
c
      if (nread.le.0) then
c  get coefficients from block data--this option not implemented,
c  place holder to maintain parallel structure with EOS setup routines
        ierr=101
        write (herr,1101) nread,hcasno,hnull
        call ERRMSG (ierr,herr)
 1101   format ('[SETVS2 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
c       write (*,*) ' SETVS2--read component',icomp,' from unit',nread
        read (nread,*) tmin(icomp)              !lower temperature limit
        read (nread,*) tmax(icomp)              !upper temperature limit
        read (nread,*) pmax(icomp)              !upper pressure limit
        read (nread,*) rhomax(icomp)            !upper density limit
        read (nread,2003) hmdeta(icomp)         !pointer to omega model
 2003   format (a3)
        read (nread,*) sigma(icomp)             !L-J sigma
        read (nread,*) epsk(icomp)              !L-J epsilon/kappa
        read (nread,*) ceta(icomp,1,1)          !Chapman-Enskog term
        ceta(icomp,1,2)=0.50d0                  !Chapman-Enskog exponent
        if (verfl(icomp).ge.6.099d0) read (nread,*) ceta(icomp,1,2)
        do j=2,12                  !read residual viscosity terms
          read (nread,*) ceta(icomp,j,1)
        enddo
        ceta(icomp,13,1)=rhoc(icomp)
        if (verfl(icomp).ge.6.099d0) read (nread,*) ceta(icomp,13,1)
        hetacr(icomp)='NUL'       !no critical enhancement in this model
        ierr=0
        herr=' '
      end if
c
      RETURN
      end                                             !subroutine SETVS2
c
c ======================================================================
c
      function ETA2DG (icomp,t)
c
c  dilute-gas contribution to the viscosity by the hydrocarbon model of:
c  Younglove, B.A. and Ely, J.F. (1987). Thermophysical properties of
c  fluids. II. Methane, ethane, propane, isobutane and normal butane.
c  J. Phys. Chem. Ref. Data  16: 577-798.
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c  output (as function value):
c   eta1dg--the dilute-gas part of the viscosity [uPa-s]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  02-21-97  MM, original version
c  07-01-98 EWL, replace SQRT(tau) with tau**ceta(i,1,2)
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)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
      character*3 hmdeta,hmdtcx
c
      common /HCHAR/ htab,hnull
c  common storing the fluid constants
      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 storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c  Lennard-Jones parameters
      common /WLJETA/ sigma(nrf0:nx),epsk(nrf0:nx)
      common /OMGMOD/ hmdeta(nrf0:nx),hmdtcx(nrf0:nx)
c
      i=icomp
      tau=t
c  in this case, the dilute gas is simply the Chapman-Enskog term
      eta2dg=ceta(i,1,1)*tau**ceta(i,1,2)/
     &       (sigma(i)**2*OMEGA(i,t,epsk(i),hmdeta(i)))
c     write (*,*) ' ETA2DG--dilute-gas viscosity:  ',eta2dg
c
      RETURN
      end                                               !function ETA2DG
c
c ======================================================================
c
      function ETA2RS (icomp,t,rho)
c
c  residual contribution to the viscosity by the hydrocarbon model of:
c  Younglove, B.A. and Ely, J.F. (1987). Thermophysical properties of
c  fluids. II. Methane, ethane, propane, isobutane and normal butane.
c  J. Phys. Chem. Ref. Data  16: 577-798.
c
c  Although this correlation has a separate initial density term, it is
c  not of the form required by ETAK1; thus the initial density term is
c  combined with the residual term.
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   eta2rs--the background part of the viscosity [uPa-s]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  02-21-97  MM, original version
c  07-01-98 EWL, replace rhoc with ceta(icomp,13,1)
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)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /HCHAR/ htab,hnull
c  common storing the fluid constants
      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 storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c
c  initial density term for viscosity
      eta1=rho*(ceta(icomp,2,1)+ceta(icomp,3,1)
     &    *(ceta(icomp,4,1)-LOG(t/ceta(icomp,5,1)))**2)           !Eq 21
c  now compute the residual viscosity (viscosity minus the dilute gas
c  and initial density terms)
      G=ceta(icomp,6,1)+ceta(icomp,7,1)/t                         !Eq 23
      H=SQRT(rho)*(rho-ceta(icomp,13,1))/ceta(icomp,13,1)         !Eq 25
      F=G+(ceta(icomp,8,1)+ceta(icomp,9,1)*t**(-1.5d0))*rho**0.1d0+
     &  (ceta(icomp,10,1)+ceta(icomp,11,1)/t+ceta(icomp,12,1)/(t*t))*H
      eta2=EXP(F)-EXP(G)                                          !Eq 22
      ETA2RS=eta1+eta2
c     write (*,*) ' ETA2RS--residual viscosity:    ',eta2rs
c
      RETURN
      end                                               !function ETA2RS
c
c ======================================================================
c
      subroutine SETVS3 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #3, minim subroutine SETTC1
c
c  inputs:
c    nread--file to read data from
c           <= 0 get data from block data
c           >0 read from logical unit nread (file should have already
c              been opened and pointer set by subroutine SETUP)
c    icomp--component number in mixture (0..nc)
c           1 for pure fluid; 0 for ECS reference fluid
c   hcasno--CAS number of component icomp (not req'd if reading from file)
c
c  outputs:
c     ierr--error flag:  0 = successful
c                       49 = error--model not implemented
c                      101 = error--block data option not implemented
c     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  02-16-99 EWL, mimic routines from SETTC1
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients
      character*1 htab,hnull
      character*3 hetamx,heta,htcxmx,htcx
      character*3 hetacr,htcxcr
      character*12 hcasno
      character*255 herr
c
      common /HCHAR/ htab,hnull
      common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax)
c  pointer to critical enhancement auxiliary functions
      common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax)
c  limits and reducing parameters
      common /WLMETA/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WR3ETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredbk(nrf0:nx),Dredbk(nrf0:nx),etarbk(nrf0:nx),
     &                tredcr(nrf0:nx),Dredcr(nrf0:nx),etarcr(nrf0:nx)
c  numbers of terms for the various parts of the model:  numerator
c  and denominator for dilute gas and background parts
      common /WN3ETA/ ndgnum(nrf0:nx),ndgden(nrf0:nx),
     &                nbknum(nrf0:nx),nbkden(nrf0:nx)
c  commons storing the (real and integer) coefficients to the viscosity
c  model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c
      if (nread.le.0) then
c  get coefficients from block data--this option not implemented,
c  place holder to maintain parallel structure with EOS setup routines
        ierr=101
        write (herr,1101) nread,hcasno,hnull
        call ERRMSG (ierr,herr)
 1101   format ('[SETVS3 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
c       write (*,*) ' SETVS3--read component',icomp,' from unit',nread
        read (nread,*) tmin(icomp)              !lower temperature limit
        read (nread,*) tmax(icomp)              !upper temperature limit
        read (nread,*) pmax(icomp)              !upper pressure limit
        read (nread,*) rhomax(icomp)            !upper density limit
        jterm=0                                 !term counter
c  read the number of terms in the numerator and denominator of the
c  dilute-gas function
        read (nread,*) ndgnum(icomp),ndgden(icomp)
c       write (*,*) ' SETVS3--about to read ',ndgnum(icomp),' +',
c    &              ndgden(icomp),' dilute terms'
        if (ndgnum(icomp).ge.1) then
          read (nread,*) treddg(icomp),etardg(icomp)  !reducing par
          do j=1,ndgnum(icomp)  !read dilute-gas terms (numerator)
            jterm=jterm+1
            read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
          enddo
        end if
        if (ndgden(icomp).ge.1) then
          do j=1,ndgden(icomp)  !read dilute-gas terms (denominator)
            jterm=jterm+1
            read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
          enddo
        end if
c
c  read the number of terms in the numerator and denominator of the
c  background model; the coefficients themselves are given in the order:
c    constant multiplier; temperature exponent; density exponent; spare
c
        read (nread,*) nbknum(icomp),nbkden(icomp)
        nbksum=nbknum(icomp)+nbkden(icomp)
c       write (*,*) ' SETVS3--about to read ',nbknum(icomp),' +',
c    &              nbkden(icomp),' background terms'
        if (nbksum.ge.1) then
c  read in reducing parameters
          read (nread,*) tredbk(icomp),Dredbk(icomp),etarbk(icomp)
          if (nbknum(icomp).ge.1) then
            do j=1,nbknum(icomp)    !numerator of rational polynomial
              jterm=jterm+1
              read (nread,*) (ceta(icomp,jterm,k),k=1,4)
            enddo
          end if
          if (nbkden(icomp).ge.1) then
            do j=1,nbkden(icomp)    !denominator of rational poly
              jterm=jterm+1
             read (nread,*) (ceta(icomp,jterm,k),k=1,4)
            enddo
          end if
        end if
c
        ierr=0
        herr=' '
      end if
c
      RETURN
      end                                             !subroutine SETVS3
c
c ======================================================================
c
      function ETA3DG (icomp,t)
c
c  dilute-gas contribution to the viscosity by the composite model (VS3)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c  output (as function value):
c   eta3dg--the dilute-gas part of the viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  02-16-99 EWL, mimic routines from TCX1DG
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients
c
c  reducing parameters
      common /WR3ETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredbk(nrf0:nx),Dredbk(nrf0:nx),etarbk(nrf0:nx),
     &                tredcr(nrf0:nx),Dredcr(nrf0:nx),etarcr(nrf0:nx)
c  numbers of terms for the various parts of the model:  numerator
c  and denominator for dilute gas and background parts
      common /WN3ETA/ ndgnum(nrf0:nx),ndgden(nrf0:nx),
     &                nbknum(nrf0:nx),nbkden(nrf0:nx)
c  commons storing the (real and integer) coefficients to the viscosity
c  model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c
      i=icomp
      tau=1.d0
      if ((ndgnum(i)+ndgden(i)).ge.1) then
c  compute tau only if dilute-gas terms exist, otherwise treddg may
c  not be defined
        tau=t/treddg(i)
c       write (*,*) 'ETA3DG--tred,tau:  ',treddg(i),tau
      end if
c
      nterm=0                                    !term counter
c  sum the dilute-gas terms, first numerator then denominator
      eta3dg=0.0d0
      if (ndgnum(i).ge.1) then
        do j=nterm+1,nterm+ndgnum(i)
          eta3dg=eta3dg+ceta(i,j,1)*tau**ceta(i,j,2)
c         write (*,*) ' ETA3DG--j,eta_dg(j): ',j,eta3dg
        enddo
        nterm=nterm+ndgnum(i)
      end if
      if (ndgden(i).ge.1) then
        denom=0.0d0
        do j=nterm+1,nterm+ndgden(i)
c       write (*,*) 'denom coeff:  ',ceta(i,j,1),ceta(i,j,2)
          denom=denom+ceta(i,j,1)*tau**ceta(i,j,2)
        enddo
c       write (*,*) ' ETA3DG--num,denom: ',eta3dg,denom
c  divide numerator by denominator
        eta3dg=eta3dg/denom
      end if
c
c  multiply by reducing parameter (to convert units, etc.)
      eta3dg=eta3dg*etardg(i)
c
c     write (*,1000) icomp,tau,eta3dg
c1000 format (' ETA3DG--icomp,tau,dilute-gas tc:',i10,d14.6,14x,d14.6)
      RETURN
      end                                               !function ETA3DG
c
c ======================================================================
c
      function ETA3RS (icomp,t,rho)
c
c  background contribution to the viscosity by the composite model (VS3)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   eta3rs--the background part of the viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  02-16-99 EWL, mimic routines from TCX1BK
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients
c
c  reducing parameters
      common /WR3ETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredbk(nrf0:nx),Dredbk(nrf0:nx),etarbk(nrf0:nx),
     &                tredcr(nrf0:nx),Dredcr(nrf0:nx),etarcr(nrf0:nx)
c  numbers of terms for the various parts of the model:  numerator
c  and denominator for dilute gas and background parts
      common /WN3ETA/ ndgnum(nrf0:nx),ndgden(nrf0:nx),
     &                nbknum(nrf0:nx),nbkden(nrf0:nx)
c  commons storing the (real and integer) coefficients to the viscosity
c  model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c
      i=icomp
      tau=1.d0
      del=1.d0
      if ((nbknum(i)+nbkden(i)).ge.1) then
c  compute tau only if residual terms exist, otherwise tredbk, Dredbk may
c  not be defined
        tau=t/tredbk(i)
        del=rho/Dredbk(i)
      end if
      nterm=ndgnum(i)+ndgden(i)                        !term counter
c     write (*,*) ' ETA3RS--tau,del,coeff_1: ',tau,del,ceta(i,nterm+1,1)
c
c  sum the background terms, first numerator then denominator
      eta3rs=0.0d0
      if (nbknum(i).ge.1) then
        do j=nterm+1,nterm+nbknum(i)
          eta3rs=eta3rs+ceta(i,j,1)*tau**ceta(i,j,2)*del**ceta(i,j,3)
        enddo
      end if
      if (nbkden(i).ge.1) then
        denom=0.0d0
        if (del.gt.0) then
          do j=nterm+nbknum(i)+1,nterm+nbknum(i)+nbkden(i)
            denom=denom+ceta(i,j,1)*tau**ceta(i,j,2)*del**ceta(i,j,3)
          enddo
        endif
c  divide numerator by denominator
        if (denom.ne.0) eta3rs=eta3rs/denom
      end if
c
c  multiply by reducing parameter (to convert units, etc.)
      eta3rs=eta3rs*etarbk(i)
c
c     write (*,1000) icomp,tau,del,eta3rs
c1000 format (' ETA3RS--icomp,tau,del,background tc:  ',i4,3d14.6)
      RETURN
      end                                               !function ETA3RS
c
c ======================================================================
c
      subroutine SETVS4 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #4; this, the "friction theory model,"
c  S.E. Quinones-Cisneros and U.K. Deiters,
c  "Generalization of the Friction Theory for Viscosity Modeling,"
c  J. Phys. Chem. B 2006, 110,12820-12834.
c
c  inputs:
c    nread--file to read data from
c           <= 0 get data from block data
c           >0 read from logical unit nread (file should have already
c              been opened and pointer set by subroutine SETUP)
c    icomp--component number in mixture (0..nc)
c           1 for pure fluid; 0 for ECS reference fluid
c   hcasno--CAS number of component icomp (not req'd if reading from file)
c
c  outputs:
c     ierr--error flag:  0 = successful
c                      101 = error--block data option not implemented
c     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  12-26-06  MLH, original version
c  01-23-07  MLH, fix mixture bug
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
      character*3 hetamx,heta,htcxmx,htcx,hmdeta,hmdtcx
      character*3 hetacr,htcxcr
      character*12 hcasno
      character*255 herr
c
      common /HCHAR/ htab,hnull
      common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax)
c  pointer to critical enhancement auxiliary functions
      common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax)
c  pointer to collision integral model
      common /OMGMOD/ hmdeta(nrf0:nx),hmdtcx(nrf0:nx)
c  limits and reducing parameters
      common /WLMETA/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
c  Lennard-Jones parameters
      common /WLJETA/ sigma(nrf0:nx),epsk(nrf0:nx)
c  numbers of terms for the various parts of the model (dilute gas,
c  initial density dependence, residual part)
      common /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
c  commons storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)

c  common storing residual coefficients to the FT visc model
      common /TRNFT/ as(nrf0:nx,0:2),bs(nrf0:nx,0:2),cs(nrf0:nx,0:2),
     &               AB(nrf0:nx,0:2),BB(nrf0:nx,0:2),CB(nrf0:nx,0:2),
     &               DB(nrf0:nx,0:2)
c
      if (nread.le.0) then
c  get coefficients from block data--this option not implemented,
c  place holder to maintain parallel structure with EOS setup routines
        ierr=101
        write (herr,1101) nread,hcasno,hnull
        call ERRMSG (ierr,herr)
 1101   format ('[SETVS4 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
c       write (*,*) ' SETVS4--read component',icomp,' from unit',nread
        read (nread,*) tmin(icomp)              !lower temperature limit
        read (nread,*) tmax(icomp)              !upper temperature limit
        read (nread,*) pmax(icomp)              !upper pressure limit
        read (nread,*) rhomax(icomp)            !upper density limit
        jterm=0                                 !term counter
        read (nread,*) ndg(icomp)               !# dilute-gas terms
c       write (*,*) ' SETVS4--about to read ',ndg(icomp),' dilute terms'
        if (ndg(icomp).ge.1) then
          read (nread,2003) hmdeta(icomp)       !pointer to omega model
          read (nread,*) sigma(icomp)           !L-J sigma
          read (nread,*) epsk(icomp)            !L-J epsilon/kappa
          read (nread,*) treddg(icomp),etardg(icomp)  !reducing par
          do j=1,ndg(icomp)                 !read dilute-gas terms
            jterm=jterm+1
            read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
          enddo
        end if
        read (nread,*) nB2(icomp)               !# visc virial terms bot presently used
c       write (*,*) ' SETVS4--about to read ',nB2(icomp),' virial terms'
        if (nB2(icomp).ge.1) then
          read (nread,*) tredB2(icomp),etarB2(icomp)  !reducing par
          do j=1,nB2(icomp)        !read viscosity virial terms
            jterm=jterm+1
            read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
          enddo
        end if
c
c  read the 18 parameters in the original generalized model
c  these are in the order: a(0:2),b(0:2),c(0:2),A(0:2),B(0:2),C(0:2)
c       write (*,*) ' SETVS4--about to read residual terms'
        READ(nread,*)(as(icomp,j),j=0,2)
        READ(nread,*)(bs(icomp,j),j=0,2)
        READ(nread,*)(cs(icomp,j),j=0,2)
        READ(nread,*)(AB(icomp,j),j=0,2)
        READ(nread,*)(BB(icomp,j),j=0,2)
        READ(nread,*)(CB(icomp,j),j=0,2)
c  one additional higher order term needed for high pressures (>500 MPa)
        READ(nread,*)(DB(icomp,j),j=0,2)
c
c  read in pointer to critical enhancement model
        read (nread,2003) hetacr(icomp)
c       write (*,*) ' SETVS4--will use critical model ',hetacr(icomp)
        ierr=0
        herr=' '
      end if
c
      RETURN
 2003 format (a3)
      end                                             !subroutine SETVS4
c
c ======================================================================
c
      function ETA4RS (icomp,t,rho)
c
c  residual contribution to the viscosity for the friction theory model (VS4)
c  S.E. Quinones-Cisneros and U.K. Deiters,
c  "Generalization of the Friction Theory for Viscosity Modeling,"
c  J. Phys. Chem. B 2006, 110,12820-12834.
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   eta4rs--the background part of the viscosity [uPa-s]
c
c  12-28-06 MLH, original version
c  01-18-07 MLH, put fgam =psi2 to accommodate new methane equation
c  01-23-07 MLH, fix mix bug
c  02-21-07 MLH, allow for higher order repulsive terms, check for 2-phase
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)   !max no. coefficients for viscosity
      character*1 htab,hnull
      CHARACTER*255 herr
      character*3 hdl,hdlk,hdv,hdvk
c
      common /HCHAR/ htab,hnull
c  common storing the fluid constants
      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  numbers of terms for the various parts of the model (dilute gas,
c  initial density dependence, residual part)
      common /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
c  common storing residual coefficients to the FT visc model
      common /TRNFT/ as(nrf0:nx,0:2),bs(nrf0:nx,0:2),cs(nrf0:nx,0:2),
     &               AB(nrf0:nx,0:2),BB(nrf0:nx,0:2),CB(nrf0:nx,0:2),
     &               DB(nrf0:nx,0:2)
      common /DLMOD/ hdl,hdlk(n0:nx)
      common /DVMOD/ hdv,hdvk(n0:nx)
      COMMON /extraeta /sumkapi
      DIMENSION FEED (ncmax),xliq(ncmax),xvap(ncmax)
c
      i=icomp
      eta4rs=0.0d0
      if (rho.lt.1.0d-10) RETURN

c  compute the residual term
      kph=1
      feed(1)= 1.0d0
      tr=tc(i)/t
      psi1= EXP(tr)-1.0d0
      psi2= EXP(tr**2)-1.0d0
      fgam = psi2
      sumkapa=(as(i,0)+as(i,1)*psi1 +as(i,2)*psi2)*tr
      sumkapaa=(AB(i,0)+AB(i,1)*psi1 +AB(i,2)*psi2)*tr**3
      sumkapr=(bs(i,0)+bs(i,1)*psi1 +bs(i,2)*psi2)*tr
      sumkaprr=(BB(i,0)+BB(i,1)*psi1 +BB(i,2)*psi2)*tr**3
      sumkapi=(cs(i,0)+cs(i,1)*psi1 +cs(i,2)*fgam)*tr
      sumkapii=(CB(i,0)+CB(i,1)*psi1 +CB(i,2)*psi2)*tr**3

      call PUREFLD(i)             !obtain properties of the pure
      ! next compute p at actual t, rho and dpdt
      call PRESS (t,rho,feed,pkpa)
      p=pkpa/100.0d0              !convert to bar
      call DPDT (t,rho,feed,dpt)


      IF(t.gt.tc(i))then
        prep = t*dpt/100.0d0      !bar
        patt = p - prep           !bar
        pid = rho * 0.08314472 * t
        delpr = prep - pid
        eta4rs = sumkapr*delpr +sumkapa*patt+sumkaprr*delpr**2 +
     &         sumkapaa*patt**2 +sumkapi*pid + sumkapii*pid**2
     &         + DB(i,0)*(prep**3)*(tr**2)
      ELSE   !could be in two-phase region
        !find sat rhol, rhov for this t and dpdt on the boundary
        !check if short sat equations are available
        if ((hdlk(i).ne.' ' .and. hdlk(i).ne.'NBS').and.
     &      (hdvk(i).ne.' ' .and. hdvk(i).ne.'NBS'))then
          call DLSATK (i,t,rhol,ierr,herr)
          call DVSATK (i,t,rhov,ierr,herr)
          call PSATK (i,t,psatkpa,ierr,herr)
        else  !no short forms present; call full
          call SATT (t,feed,kph,psatkpa,rhol,rhov,xliq,xvap,ierr,herr)
        endif
        psat=psatkpa/100.0d0      !convert to bar

      ! if in the two-phase region, interpolate liquid and vapor endpoints to avoid
      ! bizarre behavior caused by loops from the EOS
        IF((rho.gt.rhov).AND.(rho.lt.rhol)) then   ! it is two-phase region
          p=psat
          ! get viscosity on liquid side
          call DPDT (t,rhol,feed,dptl)
          prep = t*dptl/100.0d0   !bar
          patt = p - prep         !bar
          pid = rhol * 0.08314472 * t
          delpr = prep - pid
          eta4rsl = sumkapr*delpr +sumkapa*patt+sumkaprr*delpr**2 +
     &         sumkapaa*patt**2 +sumkapi*pid + sumkapii*pid**2
     &         + DB(i,0)*(prep**3)*(tr**2)
          ! get viscosity on vapor side
          call DPDT (t,rhov,feed,dptv)
          prep = t*dptv/100.0d0   !bar
          patt = p - prep         !bar
          pid = rhov * 0.08314472 * t
          delpr = prep - pid
          eta4rsv = sumkapr*delpr +sumkapa*patt+sumkaprr*delpr**2 +
     &         sumkapaa*patt**2 +sumkapi*pid + sumkapii*pid**2
     &         + DB(i,0)*(prep**3)*(tr**2)
          drat=(rhol-rho)/(rhol-rhov)
          eta4rs=eta4rsl-drat*(eta4rsl-eta4rsv)
        ELSE !not in two-phase
          prep = t*dpt/100.0d0    !bar
          patt = p - prep         !bar
          pid = rho * 0.08314472 * t
          delpr = prep - pid
          eta4rs = sumkapr*delpr +sumkapa*patt+sumkaprr*delpr**2 +
     &         sumkapaa*patt**2 +sumkapi*pid + sumkapii*pid**2
     &         + DB(i,0)*(prep**3)*(tr**2)
        ENDIF
      ENDIF

      eta4rs=eta4rs*1000.0d0 ! convert to uPa.s
      call purefld(0)        ! reset
      RETURN
      end                                               !function ETA4RS
c
c ======================================================================
c
      subroutine SETVS6 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #6
c
c  temporary place holder--this model is not yet implemented
c
c  inputs:
c    nread--file to read data from
c           <= 0 get data from block data
c           >0 read from logical unit nread (file should have already
c              been opened and pointer set by subroutine SETUP)
c    icomp--component number in mixture (0..nc)
c           1 for pure fluid; 0 for ECS reference fluid
c   hcasno--CAS number of component icomp (not req'd if reading from file)
c
c  outputs:
c     ierr--error flag:  0 = successful
c                       49 = error--model not implemented
c                      101 = error--block data option not implemented
c     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  06-18-96  MM, original version
c  08-19-97  MM, change error number for nread<=0; input hcasno is not array
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=20)   !max number of components in mixture
      parameter (nrefmx=10)  !max number of fluids for transport ECS
      parameter (n0=-ncmax-nrefmx)
      parameter (nrf0=n0)    !lower limit for transport ref fluid arrays
      character*1 htab,hnull
      character*3 hetamx,heta,htcxmx,htcx
      character*3 hetacr,htcxcr
      character*12 hcasno
      character*255 herr
c
      common /HCHAR/ htab,hnull
      common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax)
c  pointer to critical enhancement auxiliary functions
      common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax)
c
      if (nread.le.0) then
c  get coefficients from block data--this option not implemented,
c  place holder to maintain parallel structure with EOS setup routines
        ierr=101
        write (herr,1101) nread,hcasno,hnull
        call ERRMSG (ierr,herr)
 1101   format ('[SETVS6 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        hetacr(icomp)='NUL'
        ierr=49
        herr='[SETUP error 49] viscosity model #6 specified in fluid '//
     &       'file but not implemented in code.'//hnull
        call ERRMSG (ierr,herr)
      end if
c
      RETURN
      end                                             !subroutine SETVS6
c
c ======================================================================
c
c     The following functions (through EXCESV) were taken from NIST12,
c     Version 3.1, and modified to work with the current version.
c
      FUNCTION ETAH2(ICOMP,T,D)
c
c  model for the viscosity of para and normal hydrogen
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c        d--molar density [mol/L]
c  output (as function value):
c     eta5--viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  10-20-99 EWL, original version
c  09-13-02 EWL, removed code for d>27 and d>dtest
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c  The code for d>27 and d>dtest caused vis to be constant at any pressure
c  above 27 or dtest.  These checks were once used because the viscosity
c  starts to drop as the temperature drops in the liquid at high densities.
c  By changing the maximum density of the formulation to 44 mol/L, the drops
c  are not significant, and the tests can be removed.
      IF (T.GT.100.D0) THEN
c       IF (D.GT.27.D0) THEN
c         ETAH2=DILV(icomp,100.D0)+EXVDIL(27.D0,100.D0)
c    &        +DELV(icomp,D,100.D0,27.D0,100.D0)
c    &        +DELV(icomp,D,T,D,100.D0)
c       ELSE
          ETAH2=DILV(icomp,100.D0)+EXVDIL(D,100.D0)
     &          +DELV(icomp,D,T,D,100.D0)
c       ENDIF
      ELSE
c       DTEST=49.D0-0.2204D0*T
c       IF (D.LE.DTEST) THEN
          ETAH2=DILV(icomp,T)+EXVDIL(D,T)
c       ELSE
c         ETAH2=DILV(icomp,T)+EXVDIL(DTEST,T)+DELV(icomp,D,T,DTEST,T)
c       ENDIF
      ENDIF
      END

      FUNCTION DELV(icomp,D1,T1,D2,T2)
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      DELV=DILV(icomp,T1)+EXCESV(icomp,D1,T2)
     &    -DILV(icomp,T2)-EXCESV(icomp,D2,T2)
      END

      FUNCTION EXVDIL(DD,T)
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      EXVDIL=0.D0
      IF (DD.LE.0) RETURN
      D=DD*2.0159D0/1000.D0
      A=5.7694D0+DLOG(D)+.65D2*D**1.5D0-6.D-6*DEXP(127.2D0*D)
      A=DEXP(A)
      B=1.D1+7.2D0*((D/.07D0)**6-(D/.07D0)**(3.D0/2.D0))
     &-17.63D0*DEXP(-58.75D0*(D/.07D0)**3)
      EXVDIL=A*DEXP(B/T)*0.1D0
      END

      FUNCTION DILV(icomp,T)
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      SUM=0.0D0
      TF=T**(1.D0/3.D0)
      TFF=T**(-4.D0/3.D0)
      DO I=1,9
        TFF=TFF*TF
        SUM=SUM+ceta(icomp,i,1)*TFF
      ENDDO
      DILV=SUM*100.D0
      END

      FUNCTION EXCESV(icomp,DD,T)
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      DOUBLE PRECISION EV(8)
      EV(1)=ceta(icomp,10,1)
      EV(2)=ceta(icomp,11,1)
      EV(3)=ceta(icomp,12,1)
      EV(4)=ceta(icomp,13,1)
      EV(5)=ceta(icomp,14,1)
      EV(6)=ceta(icomp,15,1)
      EV(7)=ceta(icomp,16,1)
      EV(8)=ceta(icomp,17,1)
      D=DD*2.0159D0/1000.D0
      R2=D**0.5D0*(D-EV(8))/EV(8)
      RR=D**0.1D0
      X=EV(1)+EV(2)*R2+EV(3)*RR+EV(4)*R2/(T*T)+EV(5)*RR/T**1.5D0+EV(6)/T
     &+EV(7)*R2/T
      X1=EV(1)+EV(6)/T
      EXCESV=(DEXP(X)-DEXP(X1))*0.1D0
      END
c
c ======================================================================
c
      function ETAHE (icomp,t,rho)
c
c  viscosity of helium taken from:
c    Arp, V.D., McCarty, R.D., and Friend, D.G.,
c    "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with
c     Pressures ot 2000 MPa,"
c    NIST Technical Note 1334, Boulder, CO 1998.
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   etahe --viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  07-06-98 EWL, original version
c  12-01-98 EWL, add Reos and triple point pressure and density to /CCON/
c  02-16-99 EWL, rename to ETA4RS
c  11-06-00 EWL, change to ETAHE using the hardcoded model
c  01-25-07 EWL, check for large value in the calculation of expe
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /HCHAR/ htab,hnull
c  common storing the fluid constants
      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 storing the (real and integer) coefficients to the visc model
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /WIFETA/ ieta(nrf0:nx,mxeta)
c

      x=5.7037825d0
      if (t.le.300.0d0) x=dlog(t)
      dd=rho*wm(icomp)/1000.0d0
      b=-47.5295259d0/x+87.6799309d0-42.0741589d0*x
     &  +8.33128289d0*x**2-0.589252385d0*x**3
      c= 547.309267d0/x-904.870586d0+431.404928d0*x
     &  -81.4504854d0*x**2+5.37008433d0*x**3
      d=-1684.39324d0/x+3331.08630d0-1632.19172d0*x
     &  +308.804413d0*x**2-20.2936367d0*x**3

      eta0a=exp(-0.135311743d0/x+1.00347841d0+1.20654649d0*x
     &     -0.149564551d0*x**2+0.0125208416d0*x**3)
      ee=b*dd+c*dd**2+d*dd**3
      if (ee.gt.100.d0) ee=100.d0
      etae=exp(ee)
      if (t.gt.100.0d0) then
        eta0b=196.d0*t**0.71938d0
     &       *exp(12.451d0/t-295.67d0/t**2-4.1249d0)
        if (t.lt.110.0d0) then
          eta0=eta0a+(eta0b-eta0a)*(t-100.0d0)/10.0d0
        else
          eta0=eta0b
        endif
        viscos=eta0a*etae+eta0-eta0a
      else
        viscos=eta0a*etae
      endif

      etahe=viscos/10.0d0
c     write (*,*) ' ETAHE--residual viscosity:    ',etahe
c
      RETURN
      end                                                !function ETAHE
c
c ======================================================================
c
      function ETAETY (icomp,t,rho)
c
c  viscosity model for ethylene by Holland et al. (1983)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   etaety--viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  02-29-00 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      double precision gv(9)
      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)
      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
      etaety=0.0d0
      gv(1)=-3.5098225018d6
      gv(2)= 2.5008406184d6
      gv(3)=-5.8365540744d5
      gv(4)= 4.5549146583d3
      gv(5)= 2.2881683403d4
      gv(6)=-4.7318682077d3
      gv(7)= 4.5022249258d2
      gv(8)=-2.1490688088d1
      gv(9)= 4.1649263233d-1
      v1=-4.8544486732d0
      v2= 1.3033585236d1
      v3= 2.7808928908d4
      v4=-1.8241971308d3
      v5= 1.5913024509d0
      v6=-2.0513573927d2
      v7=-3.9478454708d4
      d=rho*wm(icomp)/1000.0d0
      dc=.221d0
      th=(d-dc)/dc
      tt=t**(1.0d0/3.0d0)
      eta0=gv(1)/t+gv(2)/tt**2+gv(3)/tt+gv(4)+gv(5)*tt
     &    +gv(6)*tt**2+gv(7)*t+gv(8)*tt**4+gv(9)*tt**5
      etapr=exp(v1+v4/t)*(exp(d**0.1d0*(v2+v3/t**1.5d0)
     &     +th*d**0.5d0*(v5+v6/t+v7/t**2))-1.0d0)
      etaety=(eta0+etapr)/10.0d0
c
      RETURN
      end                                               !function ETAETY
c
c ======================================================================
c
      function ETANEO (icomp,t,rho)
c
c  viscosity model for neon of:
c
c  Rabinovich, V.A., Vasserman, A.A., Nedostup, V.I. and Veksler, L.S.
c   "Thermophysical Properties of Neon, Argon, Krypton, and Xenon,"
c   Hemisphere Publishing Corp., 1988.
c
c  The numbers calculated here do not exactly match those given by Rabinovich.
c  The ECS model is currently the preferred model.
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   etaneo--viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  03-27-00 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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)
      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
      etaneo=0.0d0
      a0= 17.67484d0
      a1=-2.78751d0
      a2= 311498.7d0
      a3=-48826500d0
      a4= 3938774000d0
      a5=-1.654629d11
      a6= 2.86561d12
      tred=t/0.29944
      y=a0+a1*log10(tred)+a2/tred**2+a3/tred**3+a4/tred**4
     & +a5/tred**5+a6/tred**6
      y=y*0.68321d0
      etat=266.93d0*SQRT(t*wm(icomp))/y
      a0= 1.03010d0
      a1=-0.99175d0
      a2= 2.47127d0
      a3=-3.11864d0
      a4= 1.57066d0
      b0= 0.48148d0
      b1=-1.18732d0
      b2= 2.80277d0
      b3=-5.41058d0
      b4= 7.04779d0
      b5=-3.76608d0
      om=rho/(1673.0d0/wm(icomp))
c  The sign has been changed in Eq. 5.23 to negative
      s=a0+a1*om+a2*om**2+a3*om**3+a4*om**4
     &-(b0+b1*om+b2*om**2+b3*om**3+b4*om**4+b5*om**5)*log10(t/122.1d0)
      s=s*0.000000000305d0
      b1= 0.27676d0
      b2= 0.014355d0
      b3= 2.6480d0
      b4=-1.9643d0
      b5= 0.89161d0
      b=2.0d0/3.0d0*3.1415927d0*6.02221367d23*s**3
      brho=rho*1000.d0*b
      etad=1.0d0+b1*brho+b2*brho**2+b3*brho**3+b4*brho**4+b5*brho**5
      etaneo=etad*etat/100.d0
c
      RETURN
      end                                               !function ETANEO
c
c ======================================================================
c
      function ETAR23 (icomp,t,rho)
c
c  viscosity model for R23
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   etaR23--viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  11-01-00 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      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 /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c
      etaR23=0.0d0
      rhoL=ceta(icomp,2,1)
      C1=ceta(icomp,3,1)
      C2=ceta(icomp,4,1)
      DG=ceta(icomp,5,1)
      etamax=ceta(icomp,6,1)
      drho=rhoL-rho
      del=rho-dred(icomp)
      tau=T-tred(icomp)
      etadg=ETA1DG(icomp,t)*(drho/rhoL)**C1
      etars=(rho/rhoL)**C1*C2*rhoL**2/drho*t**0.5d0
     &     *EXP(rho/drho*DG/R/t)
      etacrt=4.d0*etamax/(exp(del)+exp(-del))/(exp(tau)+exp(-tau))
      etaR23=etadg+etars+etacrt
c
      RETURN
      end                                               !function ETAR23
c
c ======================================================================
c
      function ETAH2O (icomp,t,rho)
c
c  viscosity model for water and heavy water
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   etaH2O--viscosity [uPa-s]
c
c  written by E.W. Lemmon, NIST Phys & Chem Properties Div, Boulder, CO
c  11-06-00 EWL, original version
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 (nrf0=n0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      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 /WNTETA/ ndg(nrf0:nx),nB2(nrf0:nx),ndel0(nrf0:nx),
     &                npoly(nrf0:nx),nnum(nrf0:nx),nden(nrf0:nx),
     &                nexpn(nrf0:nx),nexpd(nrf0:nx)
      common /WRDETA/ treddg(nrf0:nx),etardg(nrf0:nx),
     &                tredB2(nrf0:nx),etarB2(nrf0:nx),
     &                tred(nrf0:nx),Dred(nrf0:nx),etared(nrf0:nx)
      common /WCFETA/ ceta(nrf0:nx,mxeta,4)
      common /Gcnst/ R,tz(n0:nx),rhoz(n0:nx)
c
      etaH2O=0.0d0
      tr=t/tred(icomp)
      dr=rho/dred(icomp)
      tau=1.d0/tr-1.d0
      del=dr-1.d0
      if (ABS(tau).lt.1.d-12 ) tau=1.d-12
      if (ABS(del).lt.1.d-12 ) del=1.d-12
      s=0.d0
      do i=1,ndel0(icomp)
        s=s+ceta(icomp,i,1)/tr**ceta(icomp,i,2)
      enddo
      eta0=SQRT(tr)/s
      s=0.d0
      do i=1,npoly(icomp)
        j=i+ndel0(icomp)
        s=s+ceta(icomp,j,1)
     &        *tau**INT(ceta(icomp,j,2))*del**INT(ceta(icomp,j,3))
      enddo
      eta1=EXP(dr*s)
      eta2=1.d0
      if (nnum(icomp).gt.0) then
        if (tr.gt.0.997d0 .and. tr.lt.1.0082d0) then
          if (dr.gt.0.755d0 .and. dr.lt.1.29d0) then
            call DPDDK (icomp,t,rho,dpdrho)
            i=ndel0(icomp)+npoly(icomp)+1
            x=dr/dpdrho*ceta(icomp,i,1)/dred(icomp)
            if (x.ge.ceta(icomp,i+1,3))
     &          eta2=ceta(icomp,i+1,1)*x**ceta(icomp,i+1,2)
          endif
        endif
      endif
      etaH2O=etared(icomp)*eta0*eta1*eta2
c
      RETURN
      end                                               !function ETAH2O
c
c ======================================================================

      function etaMEO(icomp,t,rhom)
c
c  viscosity model for methanol
c  Xiang, H.W., Huber, M.L. and Laesecke, A., "A New Reference
c  Correlation for the Viscosity of Methanol", submitted to
c  J. Phys. Chem. Ref. Data (2006)
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c        t--temperature [K]
c      rhom--molar density [mol/L]
c  output (as function value):
c   etaMEOH--viscosity [uPa-s]
c
c  11.09.05 version based on code from HWX
c
      implicit double PRECISION (a-h,o-z)
      DIMENSION h(9),d(7),e(9)

       data h/-19.572881d0,219.73999d0,-1015.3226d0,2471.01251d0,
     &  -3375.1717d0,2491.6597d0,-787.26086d0,14.085455d0,-.34664158d0/
       data d/-1.181909d0,0.5031030d0,-0.6268461d0,0.5169312d0,
     &  -0.2351349d0,5.3980235d-02,-4.9069617d-03/
       data e/4.018368d0,-4.239180d0,2.245110d0,-0.5750698d0,
     &  2.3021026d-02,2.5696775d-02,-6.8372749d-03,7.2707189d-04,
     & -2.9255711d-05/
      etaMEO=0.0d0
C     constants used in the correlation
      av=6.0221415d23
      ak=8.314472d0/av
      amm=32.04216d0 !value used in the methanol EOS
      am=amm/1000.0d0/av
c
      tc=512.6d0
      rc=273.0d0

c     convert from mol/l to kg/m3
      rho=rhom*amm

c     0-density parameters
      da0p=.3408d-9
      tcep=577.87d0
      delt=0.4575d0

c     B-coefficient parameters, the same as those of 0-density
      da0s=.3408d-9
      tces=577.87d0

c     C-coefficient parameters
      da0t=.43d-9
      tcet=440.0d0
      da0t=da0p
      tcet=tcep

      vmc=da0p
c     vmc is defined as in the critical-like viscosity
      vmc=(6.0d0/3.14159265d0*amm/rc/av/1000.0d0)**(1.0d0/3.0d0)

      r=rho  !in kg/m3
      rr=rho/rc
      tr=t/tc
      tres=t/tces
      trep=t/tcep
      tret=t/tcet
      cf=1.0d0/(1.0d0+exp(5.0d0*(rr-1.0d0)))

C     calculate the hard sphere diameter, da, eq.16
      da=0.0d0
      do kk=1,7
        da=da+d(kk)/tr**(kk-1)
      enddo
      do mm=1,9
        da=da+e(mm)*(rr)**(mm)
      enddo
      da=da*vmc

      r=amm/rho/1.0d3
       ele0=1.16145d0/trep**.14874d0+.52487d0/exp(.77320d0*trep)
     & +2.16178d0/exp(2.43787d0*trep)
       eled=0.10225d0/trep**.97346d0+.10657d0/exp(.34528d0*trep)
     & -.44557d0/exp(2.58055d0*trep)
       ele0=ele0*(1.0d0+delt**2/(1.0d0+0.95976d-3*delt**6)*eled)
       u0=5.0d0/16.0d0/da0p**2*am**.5d0*ak**.5d0*(T/3.14159265d0)**.5d0
     &  /ele0

c      B-coefficient,subcript,s
       bnr=h(1)+h(2)/tres**.25d0+h(3)/tres**.5d0+h(4)/tres**.75d0
     &  +h(5)/tres +h(6)/tres**1.25d0+h(7)/tres**1.5d0+h(8)/tres**2.5d0
     &  +h(9)/tres**5.5d0
       bn=bnr*av*da0s**3

c      C-coefficient,subcript,t
      cetar=18.6222085d-4*EXP(9.990338d0/tret**0.5d0)*tret**3
      ceta=cetar*(av*da0t**3)**2
      paa0=u0*(1.0d0+bn*rho/amm*1000.0d0+ceta*(rho/amm*1000.0d0)**2)
      b=2.0d0*3.14159265d0*av*da**3/3.0d0
      ay=b/4.0d0/r
      gd=(1.0d0 -.5d0*ay)/(1.0d0-ay)**3
      ue0=1.0d0/gd+0.8d0*b/r+.761d0*gd*(b/r)**2
      paa=ue0*u0
      paf=cf*paa0+(1.0d0-cf)*paa
      etaMEO=paf *(1.0d6)

      return
      end                                               !function ETAMEO
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                    end file trns_VIS.f
c ======================================================================
