c  begin file trns_VIS.f
c
c  This file contains the core routines for viscosity
c
c  contained here are:
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     subroutine SETVS4 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETVS5 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETVS6 (nread,icomp,hcasno,ierr,herr)
c
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                       49 = error--model 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 Thermophysics Division, Boulder, Colorado
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
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)   !max number of components in mixture
c     parameter (n0=-ncmax,nx=ncmax)
      parameter (nx=ncmax)
      parameter (nrf0=0)    !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,hmodci
      character*3 hetacr,htcxcr
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ 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/ hmodci(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
 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) hmodci(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 200 j=1,ndg(icomp)                 !read dilute-gas terms
          jterm=jterm+1
          read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
  200     continue
        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 210 j=1,nB2(icomp)        !read viscosity virial terms
          jterm=jterm+1
          read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
  210     continue
        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 220 j=1,ndel0(icomp)     !close-packed density term
            jterm=jterm+1
            read (nread,*) ceta(icomp,jterm,1),ceta(icomp,jterm,2)
  220       continue
          end if
          if (npoly(icomp).ge.1) then
            do 230 j=1,npoly(icomp)     !simple polynomial terms
            jterm=jterm+1
            read (nread,*) (ceta(icomp,jterm,k),k=1,4),ieta(icomp,jterm)
  230       continue
          end if
          if (nnum(icomp).ge.1) then
            do 240 j=1,nnum(icomp)      !numerator of rational polyonial
            jterm=jterm+1
            read (nread,*) (ceta(icomp,jterm,k),k=1,4),ieta(icomp,jterm)
  240       continue
          end if
          if (nden(icomp).ge.1) then
            do 250 j=1,nden(icomp)      !denominator of rational poly
            jterm=jterm+1
            read (nread,*) (ceta(icomp,jterm,k),k=1,4),ieta(icomp,jterm)
  250       continue
          end if
          if (nexpn(icomp).ge.1) then
            do 260 j=1,nexpn(icomp)     !numerator of exponential term
            jterm=jterm+1
            read (nread,*) (ceta(icomp,jterm,k),k=1,3),ieta(icomp,jterm)
  260       continue
          end if
          if (nexpd(icomp).ge.1) then
            do 270 j=1,nexpd(icomp)     !denominator of exponential term
            jterm=jterm+1
            read (nread,*) (ceta(icomp,jterm,k),k=1,3),ieta(icomp,jterm)
  270       continue
          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=hnull
      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
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
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)   !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /CHAR/ htab,hnull
c  common storing the fluid constants
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(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
      i=icomp
      nterm=0                                    !term counter
      eta1dg=0.0d0
c
c  sum the dilute-gas terms
      eta1dg=0.0d0
      if (ndg(i).ge.1) then
        tau=t/treddg(icomp)
c  first term is always the Chapman-Enskog term
        eta1dg=ceta(i,1,1)*SQRT(tau)/(sigma(i)**2*OMEGA(i,t,epsk(i)))
        if (ndg(i).ge.2) then
c  possibility for additional, empirical terms
          do 200 j=nterm+1,nterm+ndg(i)
          eta1dg=eta1dg+ceta(i,j,1)*tau**ceta(i,j,2)
  200     continue
        end if
      end if
c
c  multiply by reducing parameter for viscosity (to convert units, etc.)
      eta1dg=eta1dg*etardg(i)
c
      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
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)   !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /CHAR/ htab,hnull
c  common storing the fluid constants
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(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 200 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)
  200   continue
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
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)   !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /CHAR/ htab,hnull
c  common storing the fluid constants
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(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
c
      nsum=ndel0(i)+npoly(i)+nnum(i)+nden(i)+nexpn(i)+nexpd(i)
      expdel=1.0d0    !initialize only
      if (nsum.ge.1) then
c  compute tau,del only if terms present (otherwise reducing par not read in)
        tau=t/tred(icomp)
        del=rho/Dred(icomp)
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(icomp)-1.0d0).lt.0.001) then
          expdel=rho/rhoc(icomp)
        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 200 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)
  200   continue
        nterm=nterm+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 220 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 (*,*) ' ETA1RS--j,polynomial term:    ',j,ceta(i,j,1)
  220   continue
c       write (*,*) ' ETA1RS--j,polynomial term:    ',j,eta1rs
        nterm=nterm+npoly(i)
      end if
c
c  calculate the numerator of the rational polyonial
      if (nnum(i).ge.1) then
        xnum=0.0d0
        do 240 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
  240   continue
c       write (*,*) ' ETA1RS--numerator term:       ',xnum
        nterm=nterm+nnum(i)
      else
        xnum=1.0d0
      end if
c
c  calculate the denominator of the rational polyonial
      if (nden(i).ge.1) then
        xden=0.0d0
        do 250 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
  250   continue
        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
        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 260 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)
  260   continue
        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 270 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)
  270   continue
        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                       49 = error--model 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 Thermophysics Division, Boulder, Colorado
c  02-21-97  MM, original version
c  02-26-97  MM, set pointer to critical ehancement to 'NUL' (not used)
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=5)   !max number of components in mixture
c     parameter (n0=-ncmax,nx=ncmax)
      parameter (nx=ncmax)
      parameter (nrf0=0)    !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,hmodci
      character*3 hetacr,htcxcr
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ 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/ hmodci(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)
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
 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) hmodci(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
        do 210 j=2,12                  !read residual viscosity terms
        read (nread,*) ceta(icomp,j,1)
  210   continue
        hetacr(icomp)='NUL'       !no critical enhancement in this model
        ierr=0
        herr=hnull
      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
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)   !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /CHAR/ htab,hnull
c  common storing the fluid constants
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(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)
c
      i=icomp
      tau=t
c  in this case, the dilute gas is simply the Chapman-Enskog term
      eta2dg=ceta(i,1,1)*SQRT(tau)/(sigma(i)**2*OMEGA(i,t,epsk(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
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      parameter (ncmax=5)   !max number of components in mixture
      parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !lower limit for transport ref fluid arrays
      parameter (mxeta=40)  !max no. coefficients for viscosity
      character*1 htab,hnull
c
      common /CHAR/ htab,hnull
c  common storing the fluid constants
      common /CCON/ wm(n0:nx),ttp(n0:nx),tnbp(n0:nx),
     &              tc(n0:nx),pc(n0:nx),rhoc(n0:nx),Zcrit(n0:nx),
     &              accen(n0:nx),dipole(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-rhoc(icomp))/rhoc(icomp)                   !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
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     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
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=5)   !max number of components in mixture
c     parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !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 /CHAR/ 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
 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)
        hetacr(icomp)='NUL'
        ierr=49
        herr='[SETUP error 49] viscosity model #3 specified in fluid '//
     &       'file but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETVS3
c
c ======================================================================
c
      subroutine SETVS4 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #4
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     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
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=5)   !max number of components in mixture
c     parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !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 /CHAR/ 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
 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)
        hetacr(icomp)='NUL'
        ierr=49
        herr='[SETUP error 49] viscosity model #4 specified in fluid '//
     &       'file but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETVS4
c
c ======================================================================
c
      subroutine SETVS5 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid viscosity model #5
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     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
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=5)   !max number of components in mixture
c     parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !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 /CHAR/ 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
 1101   format ('[SETVS5 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 #5 specified in fluid '//
     &       'file but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETVS5
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     herr--error string (character*255 variable if ierr<>0)
c     other quantities returned via arrays in commons
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
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=5)   !max number of components in mixture
c     parameter (n0=-ncmax,nx=ncmax)
      parameter (nrf0=0)    !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 /CHAR/ 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
 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
      end if
c
      RETURN
      end                                             !subroutine SETVS6
c
c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                    end file trns_VIS.f
c ======================================================================