c  begin file trns_TCX.f
c
c  This file contains the core routines for thermal conductivity
c
c  contained here are:
c     subroutine SETTC1 (nread,icomp,hcasno,ierr,herr)
c     function TCX1DG (icomp,t)
c     function TCX1BK (icomp,t,rho)
c     subroutine SETTC2 (nread,icomp,hcasno,ierr,herr)
c     function TCX2DG (icomp,t)
c     function TCX2BK (icomp,t,rho)
c     function TCX2CR (icomp,t,rho)
c     subroutine SETTC3 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETTC4 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETTC5 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETTC6 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETTK1 (nread,icomp,hcasno,ierr,herr)
c     function TCX1CR (icomp,t,rho)
c     subroutine SETTK3 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETTK4 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETTK5 (nread,icomp,hcasno,ierr,herr)
c     subroutine SETTK6 (nread,icomp,hcasno,ierr,herr)
c     function TCCNH3 (icomp,t,rho)
c
c =====================================================================
c =====================================================================
c
      subroutine SETTC1 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid thermal conductivity model #1; this, the
c  "composite model," is written in a general form with terms designed
c  to include several recent correlations including those of Tufeu
c  (1984) for ammonia, Krauss (1996) for R152a, and Laesecke (1996) for
c  R123.
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 Phys & Chem Properties Div, Boulder, CO
c  02-06-97  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 (nx=ncmax)
      parameter (nrf0=0)    !lower limit for transport ref fluid arrays
      parameter (mxtcx=40)  !max no. coefficients for thermal cond
      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  pointer to collision integral model
c     common /OMGMOD/ hmodci(nrf0:nx)
c  limits and reducing parameters
      common /WLMTCX/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCX/ treddg(nrf0:nx),tcxrdg(nrf0:nx),
     &                tredbk(nrf0:nx),Dredbk(nrf0:nx),tcxrbk(nrf0:nx),
     &                tredcr(nrf0:nx),Dredcr(nrf0:nx),tcxrcr(nrf0:nx)
c  numbers of terms for the various parts of the model:  numerator
c  and denominator for dilute gas and background parts
      common /WNTTCX/ ndgnum(nrf0:nx),ndgden(nrf0:nx),
     &                nbknum(nrf0:nx),nbkden(nrf0:nx)
c  commons storing the (real and integer) coefficients to the thermal
c  conductivity model
      common /WCFTCX/ ctcx(nrf0:nx,mxtcx,4)
      common /WIFTCX/ itcx(nrf0:nx,mxtcx)
c  Lennard-Jones parameters
      common /WLJTCX/ 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 ('[SETTC1 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 (*,*) ' SETTC1--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 demoninator of the
c  dilute-gas function
        read (nread,*) ndgnum(icomp),ndgden(icomp)
c       write (*,*) ' SETTC1--about to read ',ndgnum(icomp),' +',
c    &              ndgden(icomp),' dilute terms'
        if (ndgnum(icomp).ge.1) then
          read (nread,*) treddg(icomp),tcxrdg(icomp)  !reducing par
          do 200 j=1,ndgnum(icomp)  !read dilute-gas terms (numerator)
          jterm=jterm+1
          read (nread,*) ctcx(icomp,jterm,1),ctcx(icomp,jterm,2)
  200     continue
        end if
        if (ndgden(icomp).ge.1) then
          do 210 j=1,ndgden(icomp)  !read dilute-gas terms (denominator)
          jterm=jterm+1
          read (nread,*) ctcx(icomp,jterm,1),ctcx(icomp,jterm,2)
  210     continue
        end if
c
c  read the number of terms in the numerator and demoninator 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 (*,*) ' SETTC1--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),tcxrbk(icomp)
          if (nbknum(icomp).ge.1) then
            do 240 j=1,nbknum(icomp)    !numerator of rational polyonial
            jterm=jterm+1
            read (nread,*) (ctcx(icomp,jterm,k),k=1,4)
  240       continue
          end if
          if (nbkden(icomp).ge.1) then
            do 250 j=1,nbkden(icomp)    !denominator of rational poly
            jterm=jterm+1
            read (nread,*) (ctcx(icomp,jterm,k),k=1,4)
  250       continue
          end if
        end if
c
c  read in pointer to critical enhancement model
        read (nread,2003) htcxcr(icomp)
c       write (*,*) ' SETTC1--will use critical model ',htcxcr(icomp)
c
        ierr=0
        herr=hnull
      end if
c
      RETURN
 2003 format (a3)
      end                                             !subroutine SETTC1
c
c ======================================================================
c
      function TCX1DG (icomp,t)
c
c  dilute-gas contribution to the thermal conductivity by the
c  composite model (TC1)
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   tcx1dg--the dilute-gas part of the thermal conducitivity [W/m-K]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  02-06-97  MM, original version
c  03-28-97  MM, compute tau only if d.g. terms exist
c  06-08-97  MM, add special term; power = -99: mult by (1 + coeff*(Cp0-2.5R))
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 (mxtcx=40)  !max no. coefficients for thermal cond
c
c  reducing parameters
      common /WRDTCX/ treddg(nrf0:nx),tcxrdg(nrf0:nx),
     &                tredbk(nrf0:nx),Dredbk(nrf0:nx),tcxrbk(nrf0:nx),
     &                tredcr(nrf0:nx),Dredcr(nrf0:nx),tcxrcr(nrf0:nx)
c  numbers of terms for the various parts of the model:  numerator
c  and denominator for dilute gas and background parts
      common /WNTTCX/ ndgnum(nrf0:nx),ndgden(nrf0:nx),
     &                nbknum(nrf0:nx),nbkden(nrf0:nx)
c  commons storing the (real and integer) coefficients to the thermal
c  conductivity model
      common /WCFTCX/ ctcx(nrf0:nx,mxtcx,4)
      common /WIFTCX/ itcx(nrf0:nx,mxtcx)
      common /Gcnst/ R
c
      i=icomp
      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(icomp)
      end if
c
      nterm=0                                    !term counter
c  sum the dilute-gas terms, first numerator then denominator
      tcx1dg=0.0d0
      if (ndgnum(i).ge.1) then
        do 200 j=nterm+1,nterm+ndgnum(i)
        if (ABS(ctcx(i,j,2)+99.0d0).lt.1.0d-6) then
c  flag:  exponent -99 indicates: multiply numerator term by
c  [1 + coeff*(Cp0 - 2.5*R)]; this is the Cv_internal as used by
c  Vesovic, et al (1990) for carbon dioxide
          cp01=1.0d0+ctcx(i,j,1)*(CP0K(icomp,t)-2.5d0*R)
c         write (*,*) ' TCX1DG--will mult by Cp0 - 2.5*R: ',cp01
          tcx1dg=tcx1dg*cp01
        else
          tcx1dg=tcx1dg+ctcx(i,j,1)*tau**ctcx(i,j,2)
        end if
  200   continue
      end if
      if (ndgden(i).ge.1) then
        denom=0.0d0
        do 240 j=nterm+1,nterm+ndgden(i)
        denom=denom+ctcx(i,j,1)*tau**ctcx(i,j,2)
  240   continue
c       write (*,*) ' TCX1DG--num,denom: ',tcx1dg,denom
c  divide numerator by denominator
        tcx1dg=tcx1dg/denom
      end if
c
c  multiply by reducing parameter (to convert units, etc.)
      tcx1dg=tcx1dg*tcxrdg(i)
c
c     write (*,1000) icomp,tau,tcx1dg
c1000 format (' TCX1DG--icomp,tau,dilute-gas tc:',i10,d14.6,14x,d14.6)
      RETURN
      end                                               !function TCX1DG
c
c ======================================================================
c
      function TCX1BK (icomp,t,rho)
c
c  background contribution to the thermal conductivity by the
c  composite model (TC1)
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   tcx1bk--the dilute-gas part of the thermal conducitivity [W/m-K]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  02-06-97  MM, original version
c  03-28-97  MM, compute tau,del only if residual terms exist
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 (mxtcx=40)  !max no. coefficients for thermal cond
c
c  reducing parameters
      common /WRDTCX/ treddg(nrf0:nx),tcxrdg(nrf0:nx),
     &                tredbk(nrf0:nx),Dredbk(nrf0:nx),tcxrbk(nrf0:nx),
     &                tredcr(nrf0:nx),Dredcr(nrf0:nx),tcxrcr(nrf0:nx)
c  numbers of terms for the various parts of the model:  numerator
c  and denominator for dilute gas and background parts
      common /WNTTCX/ ndgnum(nrf0:nx),ndgden(nrf0:nx),
     &                nbknum(nrf0:nx),nbkden(nrf0:nx)
c  commons storing the (real and integer) coefficients to the thermal
c  conductivity model
      common /WCFTCX/ ctcx(nrf0:nx,mxtcx,4)
      common /WIFTCX/ itcx(nrf0:nx,mxtcx)
c
      i=icomp
      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 (*,*) ' TCX1BK--tau,del,coeff_1: ',tau,del,ctcx(i,nterm+1,1)
c
c  sum the background terms, first numerator then denominator
      tcx1bk=0.0d0
      if (nbknum(i).ge.1) then
        do 200 j=nterm+1,nterm+nbknum(i)
        tcx1bk=tcx1bk+ctcx(i,j,1)*tau**ctcx(i,j,2)*del**ctcx(i,j,3)
  200   continue
      end if
      if (nbkden(i).ge.1) then
        denom=0.0d0
        do 240 j=nterm+1,nterm+nbkden(i)
        denom=denom+ctcx(i,j,1)*tau**ctcx(i,j,2)*del**ctcx(i,j,3)
  240   continue
c  divide numerator by denominator
        tcx1bk=tcx1bk/denom
      end if
c
c  multiply by reducing parameter (to convert units, etc.)
      tcx1bk=tcx1bk*tcxrbk(i)
c
c     write (*,1000) icomp,tau,del,tcx1bk
c1000 format (' TCX1BK--icomp,tau,del,background tc:  ',i4,3d14.6)
      RETURN
      end                                               !function TCX1BK
c
c ======================================================================
c
      subroutine SETTC2 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid thermal conductivity model #2--the hydrocarbon
c  model of Younglove and Ely, JPCRD 16:577-798 (1987)
c
c  N.B.  Younglove and Ely use a special scaled equation of state to
c        compute derivatives for the critical enhancement; the default
c        EOS is used here
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                       39 = 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 (dummy placeholder)
c  10-16-96  MM, implement model of Younglove & Ely
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 (mxtcx=40)  !max no. coefficients for thermal cond
      character*1 htab,hnull
      character*3 hetamx,heta,htcxmx,htcx,hetacr,htcxcr
      character*3 hmodci
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ htab,hnull
      common /TRNMOD/ hetamx,heta(nrf0:ncmax),htcxmx,htcx(nrf0:ncmax)
c  pointers to critical enhancement and collision int auxiliary functions
      common /CREMOD/ hetacr(nrf0:ncmax),htcxcr(nrf0:ncmax)
      common /OMGMOD/ hmodci(nrf0:nx)
c  commons storing the (real and integer) coefficients to the thermal
c  conductivity model, also the t,p,rho limits
      common /WCFTCX/ ctcx(nrf0:nx,mxtcx,4)
      common /WLMTCX/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WIFTCX/ itcx(nrf0:nx,mxtcx)
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 ('[SETTC2 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 (*,*) ' SETTC2--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
        read (nread,*) ctcx(icomp,1,1)          !L-J sigma
        read (nread,*) ctcx(icomp,2,1)          !L-J epsilon/kappa
c  read constant in Eq 19 = 5/16*(k*MW/1000/pi/Na)**0.5*1.0d12
c  the factor of 1d12 is for sigma in nm and viscosity in micro-Pa-s
        read (nread,*) ctcx(icomp,3,1) !const in Eq 19
        do 200 i=4,5
        read (nread,*) ctcx(icomp,i,1) !dilute-gas terms Gt(1)-Gt(2)
  200   continue
        do 220 i=6,13
        read (nread,*) ctcx(icomp,i,1) !background terms Et(1)-Et(8)
  220   continue
c  read in pointer to critical enhancement model
        read (nread,2003) htcxcr(icomp)
c       write (*,*) ' SETTC2--will use critical model ',htcxcr(icomp)
c  in the case of TC2, the critical enhancement is integral with the
c  model for the dilute gas and background contributuions
        do 240 i=14,17
        read (nread,*) ctcx(icomp,i,1) !critical enhancement X(1)-X(4)
  240   continue
        read (nread,*) ctcx(icomp,18,1)!critical enhancement:  Z
        read (nread,*) ctcx(icomp,19,1)!critical enhancement:  k
c  following coefficients are for the viscosity function of Younglove &
c  Ely, which is also used in the thermal conductivity
        do 260 i=20,23
        read (nread,*) ctcx(icomp,i,1) !initial rho terms:  Fv(1)-Fv(4)
  260   continue
        do 280 i=24,30
        read (nread,*) ctcx(icomp,i,1) !residual viscosity:  Ev(1)-Ev(7)
  280   continue
        ierr=0
        herr=hnull
      end if
c
      RETURN
 2003 format (a3)
      end                                             !subroutine SETTC2
c
c ======================================================================
c
      function TCX2DG (icomp,t)
c
c  dilute-gas contribution to the thermal conductivity by the
c  model of Younglove and Ely, JPCRD 16:577-798 (1987); Eqs 19, 27.
c
c  N.B. there are two terms missing from Eq 27 in the JPCRD article;
c       the 15R/4 is missing and a factor of 1/(mol wt) is needed to
c       convert units
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    tcxdg--the dilute gas part of the thermal conductivity [W/m-K]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  10-16-96 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 (mxtcx=40)  !max no. coefficients for thermal cond
c
c  commons storing the (real and integer) coefficients to the thermal
c  conductivity model, also the t,p,rho limits
      common /WCFTCX/ ctcx(nrf0:nx,mxtcx,4)
      common /WLMTCX/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WIFTCX/ itcx(nrf0:nx,mxtcx)
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)
      common /Gcnst/ R
c
c  compute dilute-gas viscosity
      sigma=ctcx(icomp,1,1)
      epsk=ctcx(icomp,2,1)
      eta0=ctcx(icomp,3,1)*SQRT(t)/(OMEGA(icomp,t,epsk)*sigma**2)
      Cp0=CP0K(icomp,t)    !ideal gas heat capacity by chosen EOS
      Gt1=ctcx(icomp,4,1)
      Gt2=ctcx(icomp,5,1)
      TCX2DG=1.0d-3*eta0/wm(icomp)*(3.75d0*R+
     &       (Cp0-2.5d0*R)*(Gt1+Gt2*epsk/t))
c     write (*,1001) t,Cp0,TCX2ID
c1001 format (1x,' TCX2DG--t,Cp0,tcx_ideal: ',2f10.2,f12.6)
c
      RETURN
      end                                               !function TCX2DG
c
c ======================================================================
c
      function TCX2BK (icomp,t,rho)
c
c  background contribution to the thermal conductivity by the
c  model of Younglove and Ely, JPCRD 16:577-798 (1987); Eqs 26, 28-30.
c
c  N.B. the powers given in the JPCRD article are incorrect for Eqs 29
c       and 30; they should be (4 - n) and (7 - n), respectively; there
c       is an incorrect sign in Eq 26 [(1 + F2*rho), not (1 - F2*rho)]
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    tcxbk--the background part of the thermal conductivity [W/m-K]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  10-16-96 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 (mxtcx=40)  !max no. coefficients for thermal cond
c
c  commons storing the (real and integer) coefficients to the thermal
c  conductivity model, also the t,p,rho limits
      common /WCFTCX/ ctcx(nrf0:nx,mxtcx,4)
      common /WLMTCX/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WIFTCX/ itcx(nrf0:nx,mxtcx)
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
c  compute functions given as Eqs 28-30
      tinv=1.0d0/t
      tinv2=tinv*tinv
      F0=ctcx(icomp,6,1)+ctcx(icomp,7,1)*tinv+ctcx(icomp,8,1)*tinv2
      F1=ctcx(icomp,9,1)+ctcx(icomp,10,1)*tinv+ctcx(icomp,11,1)*tinv2
      F2=ctcx(icomp,12,1)+ctcx(icomp,13,1)*tinv
c  background term is given by Eq 26 in Younglove & Ely
      TCX2BK=(F0+F1*rho)*rho/(1.0d0+F2*rho)
c     write (*,1001) t,rho,TCX2RS
c1001 format (1x,' TCX2RS--t,rho,tcx_resid: ',f10.2,f10.4,f12.6)
c
      RETURN
      end                                               !function TCX2BK
c
c ======================================================================
c
      function TCX2CR (icomp,t,rho)
c
c  critical enhancement to the thermal conductivity by the model of
c  Younglove and Ely, JPCRD 16:577-798 (1987); Eqs D1-D4
c
c  N.B.  there are numerous errors in the equations presented in the
c        Younglove & Ely paper; the present code is derived from the
c        code of Younglove used to generate the tables in JPCRD
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    tcxcr--the critical enhancement part of the thermal conductivity [W/m-K]
c
c  written by M. McLinden, NIST Phys & Chem Properties Div, Boulder, CO
c  10-28-96  MM, original version
c  11-14-97  MM, return zero if rho = 0 (avoid division by zero)
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 (mxtcx=40)  !max no. coefficients for thermal cond
c
c  commons storing the (real and integer) coefficients to the thermal
c  conductivity model, also the t,p,rho limits
c  the array ctcx contains, in order:
c       1-2:  Lennard-Jones sigma and epsilon/kappa
c         3:  constant in Eq 19 (5/16*(k*MW/1000/pi/Na)**0.5*1.0d12)
c       4-5:  dilute gas terms, Gt(1), Gt(2)
c      6-13:  background terms, Et(1) - Et(8)
c     14-17:  critical enhancement terms, X(1) - X(4)
c        18:  critical enhancement term  Z
c        19:  Boltzmann's constant, k
c  and the following terms from the Younglove & Ely viscosity model:
c     20-23:  initial density dependence terms, Fv(1) - Fv(4)
c     24-30:  residual viscosity terms, Ev(1) - Ev(7)
      common /WCFTCX/ ctcx(nrf0:nx,mxtcx,4)
      common /WLMTCX/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WIFTCX/ itcx(nrf0:nx,mxtcx)
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)
      common /Gcnst/ R
c
c  if density is approx zero, return zero for the critical enhancement
c  (avoid division by zero)
      if (rho.lt.1.0d-6) then
        TCX2CR=0.0d0
        RETURN
      end if
c
c  find derivatives dP/dD and dP/dT
      call REDK (icomp,t0,rho0)
      tau=t0/t
      del=rho/rho0
      phi01=PHIK(icomp,0,1,tau,del)  !real-gas terms
      phi02=PHIK(icomp,0,2,tau,del)
      phi11=PHIK(icomp,1,1,tau,del)
c  factor 1.0d3 in next 3 lines converts kPa -> Pa
      dpdrho=R*t*(1.0d0+2.0d0*del*phi01+del*del*phi02)*1.0d3
      dpt=R*rho*(1.0d0+del*phi01-del*tau*phi11)*1.0d3
      pcrit=pc(icomp)*1.0d3
      xi=(pcrit*rho/(rhoc(icomp)**2*dpdrho))**ctcx(icomp,16,1)    !Eq D3
      dellam=ctcx(icomp,17,1)*ctcx(icomp,19,1)/pcrit
     &      *(t*dpt*rhoc(icomp)/rho)**2*xi                        !Eq D2
      delt=ABS(t-tc(icomp))/tc(icomp)
      delD=ABS(rho-rhoc(icomp))/rhoc(icomp)
      eterm=ctcx(icomp,14,1)*delt**4+ctcx(icomp,15,1)*delD**4
c  check that exponential term will not result in underflow
      if (eterm.gt.500.0d0) then
        eterm=EXP(-500.0d0)
      else
        eterm=EXP(-eterm)
      end if
c  now compute the viscosity, first the dilute gas contribution
      sigma=ctcx(icomp,1,1)
      epsk=ctcx(icomp,2,1)
      eta0=ctcx(icomp,3,1)*SQRT(t)/(OMEGA(icomp,t,epsk)*sigma**2)
c  initial density term for viscosity
      eta1=rho*(ctcx(icomp,20,1)+ctcx(icomp,21,1)
     &    *(ctcx(icomp,22,1)-LOG(t/ctcx(icomp,23,1)))**2)         !Eq 21
c  now compute the residual viscosity (viscosity minus the dilute gas
c  and initial density terms)
      G=ctcx(icomp,24,1)+ctcx(icomp,25,1)/t                       !Eq 23
      H=SQRT(rho)*(rho-rhoc(icomp))/rhoc(icomp)                   !Eq 25
      F=G+(ctcx(icomp,26,1)+ctcx(icomp,27,1)*t**(-1.5d0))*rho**0.1d0+
     &  (ctcx(icomp,28,1)+ctcx(icomp,29,1)/t+ctcx(icomp,30,1)/(t*t))*H
      eta2=EXP(F)-EXP(G)                                          !Eq 22
      visc=(eta0+eta1+eta2)*1.0d-6       !factor of d-6 converts to Pa-s
c     write (*,1060) t,rho,eta0,eta1,eta2,visc
c1060 format (1x,' TCX2CR: t,rho,eta0,1,2,visc:  ',f8.3,f10.6,4e14.6)
c
c  combine all of the above to arrive at the critical enhancement
c  (the denominator is missing from the Younglove & Ely paper,
      TCX2CR=dellam*eterm/(6.0d0*3.141592654d0*ctcx(icomp,18,1)*visc)
c
      RETURN
      end                                               !function TCX2CR
c
c ======================================================================
c
      subroutine SETTC3 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid thermal conductivity 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                       39 = 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, error 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,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 ('[SETTC3 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        htcxcr(icomp)='NUL'
        ierr=39
        herr='[SETUP error 39] thermal conductivity model #3 specified '
     &     //'in fluid file but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETTC3
c
c ======================================================================
c
      subroutine SETTC4 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid thermal conductivity 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                       39 = 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, error 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,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 ('[SETTC3 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        htcxcr(icomp)='NUL'
        ierr=39
        herr='[SETUP error 39] thermal conductivity model #4 specified '
     &     //'in fluid file but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETTC4
c
c ======================================================================
c
      subroutine SETTC5 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid thermal conductivity 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                       39 = 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  10-30-96  MM, original version
c  08-19-97  MM, error 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,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 ('[SETTC3 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        htcxcr(icomp)='NUL'
        ierr=39
        herr='[SETUP error 39] thermal conductivity model #5 specified '
     &     //'in fluid file but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETTC5
c
c ======================================================================
c
      subroutine SETTC6 (nread,icomp,hcasno,ierr,herr)
c
c  initialize pure fluid thermal conductivity 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                       39 = 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  10-30-96  MM, original version
c  08-19-97  MM, error 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,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 ('[SETTC3 error 101] illegal file specified; nread = ',
     &          i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        htcxcr(icomp)='NUL'
        ierr=39
        herr='[SETUP error 39] thermal conductivity model #6 specified '
     &     //'in fluid file but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETTC6
c
c ======================================================================
c
      subroutine SETTK1 (nread,icomp,hcasno,ierr,herr)
c
c  initialize model #1 for the thermal conductivity critical enhancement--
c  the empirical model used by Perkins and Laesecke
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                        1 = error
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-24-97  MM, original version
c  03-27-97  MM, change powers in exp term to integer
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 (mxtck=40)  !max no. coefficients for t.c. crit
      character*1 htab,hnull
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ htab,hnull
c  limits and reducing parameters (separate reducing par for poly & exp)
      common /WLMTCK/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCK/ tred(nrf0:nx),Dred(nrf0:nx),tcxred(nrf0:nx),
     &                tredex(nrf0:nx),Dredex(nrf0:nx)
c  numbers of terms for the various parts of the model:
c  polynomial (numerator & denominator), exponential, spare
      common /WNTTCK/ nnum(nrf0:nx),nden(nrf0:nx),nexp(nrf0:nx),
     &                nspare(nrf0:nx)
c  commons storing the (real and integer) coefficients to the model
      common /WCFTCK/ ctck(nrf0:nx,mxtck,5)
      common /WIFTCK/ itck(nrf0:nx,mxtck,0:5)
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=39
        write (herr,1039) hcasno,hnull
 1039   format ('[SETUP error 39] block data option for t.c. crit ',
     &          'model #1 specified for CAS # ',a12,' but not ',
     &          'implemented in code.',a1)
      else
c  read data from file (should have been opened by SETUP)
c       write (*,*) ' SETTK1--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  # terms in numerator & denominator of polynomial multiplier,
c  exponential term, and spare for future use
        read (nread,*) nnum(icomp),nden(icomp),nexp(icomp),nspare(icomp)
        if (nnum(icomp)+nden(icomp).ge.1) then  !read reducing pars
          read (nread,*) tred(icomp),Dred(icomp),tcxred(icomp)
          if (nnum(icomp).ge.1) then
            do 200 j=1,nnum(icomp)              !read numerator terms
            jterm=jterm+1
            read (nread,*) (ctck(icomp,jterm,k),k=1,5),
     &                     itck(icomp,jterm,0)
  200       continue
          end if
          if (nden(icomp).ge.1) then
            do 220 j=1,nden(icomp)              !read denominator terms
            jterm=jterm+1
            read (nread,*) (ctck(icomp,jterm,k),k=1,5),
     &                     itck(icomp,jterm,0)
  220       continue
          end if
        end if
        if (nexp(icomp).ge.1) then
          read (nread,*) tredex(icomp),Dredex(icomp)
          do 240 j=1,nexp(icomp)                !read exponential terms
          jterm=jterm+1
          read (nread,*) (ctck(icomp,jterm,k),k=1,2),itck(icomp,jterm,2)
     &                  ,ctck(icomp,jterm,3),itck(icomp,jterm,3)
     &                  ,itck(icomp,jterm,0)
  240     continue
        end if
c       if (nspare(icomp).ge.1) then
c         do 260 j=1,nspare(icomp)              !read spare terms
c         jterm=jterm+1
c         read (nread,*) (ctck(icomp,jterm,k),k=1,1),itck(icomp,jterm,1)
c 260     continue
c       end if
      end if
c
      RETURN
      end                                             !subroutine SETTK1
c
c ======================================================================
c
      function TCX1CR (icomp,t,rho)
c
c  model #1 for the thermal conductivity critical enhancement--
c  the empirical model used by Perkins and Laesecke
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   tcx1cr--the critical enhancement to the thermal conducitivity [W/m-K]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  02-24-97  MM, original version
c  03-27-97  MM, change powers in exp term to integer
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 (mxtck=40)  !max no. coefficients for t.c. crit
      character*1 htab,hnull
c
      common /CHAR/ htab,hnull
c  limits and reducing parameters (separate reducing par for poly & exp)
      common /WLMTCK/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCK/ tred(nrf0:nx),Dred(nrf0:nx),tcxred(nrf0:nx),
     &                tredex(nrf0:nx),Dredex(nrf0:nx)
c  numbers of terms for the various parts of the model:
c  polynomial (numerator & denominator), exponential, spare
      common /WNTTCK/ nnum(nrf0:nx),nden(nrf0:nx),nexp(nrf0:nx),
     &                nspare(nrf0:nx)
c  commons storing the (real and integer) coefficients to the model
      common /WCFTCK/ ctck(nrf0:nx,mxtck,5)
      common /WIFTCK/ itck(nrf0:nx,mxtck,0:5)
c
c  compute the various parts of the critical enhancement
c  these are taken in the order:
c    rational polynomial in T, rho (first numerator, then denominator)
c    exponential term
c    spare for future use
c  the coefficients themselves are given in the order:
c    constant multiplier;
c    additive term to temperature, exponent for (T + Tadd)
c    additive term to density, exponent for (D + Dadd)
c    spare 1 [if = 99 for denominator take MAX(T, T+Tadd]
c    spare 2
c
      i=icomp
      nterm=0                                 !term counter
      tcx1cr=0.0d0
      if (nnum(i)+nden(i).ge.1) then
        if (tred(i).le.0.0d0) then
          tau=-tred(i)/t     !negative tred indicates reverse order
        else
          tau=t/tred(i)
        end if
        del=rho/Dred(i)
        if (nnum(i).ge.1) then
          do 200 j=nterm+1,nterm+nnum(i)              !numerator terms
          tcx1cr=tcx1cr+ctck(i,j,1)*(tau+ctck(i,j,2))**ctck(i,j,3)
     &          *(del+ctck(i,j,4))**ctck(i,j,5)
  200     continue
          nterm=nterm+nnum(i)
        end if
        if (nden(i).ge.1) then
          xden=0.0d0
          do 220 j=nterm+1,nterm+nden(i)              !denominator terms
          if (itck(i,j,0).eq.99) then
c  flag for special form used by Perkins (MAX(T,2Tc-T))
            tcmax=MAX(tau,ctck(i,j,2)-tau)
            xden=xden+ctck(i,j,1)*tcmax**ctck(i,j,3)
     &          *(del+ctck(i,j,4))**ctck(i,j,5)
          else
            xden=xden+ctck(i,j,1)*(tau+ctck(i,j,2))**ctck(i,j,3)
     &          *(del+ctck(i,j,4))**ctck(i,j,5)
          end if
  220     continue
          nterm=nterm+nden(i)
          tcx1cr=tcx1cr/xden
        end if
      end if
      if (nexp(i).ge.1) then
        if (tredex(i).le.0.0d0) then
          tau=-tredex(i)/t   !negative tred indicates reverse order
        else
          tau=t/tredex(i)
        end if
        del=rho/Dredex(i)
        xexp=0.0d0
        do 240 j=nterm+1,nterm+nexp(i)              !exponential terms
        xexp=xexp+ctck(i,j,1)*(tau+ctck(i,j,2))**itck(i,j,2)
     &      *(del+ctck(i,j,3))**itck(i,j,3)
  240   continue
        nterm=nterm+nexp(i)
        if (xexp.lt.-500.0d0) then
c         tcx1cr=tcx1cr   !avoid underflow if exponential is almost one
        else
          tcx1cr=tcx1cr*EXP(xexp)
        end if
      end if
c     if (nspare(i).ge.1) then
c       xspare=0.0d0
c       do 260 j=nterm+1,nterm+nspare(i)            !spare terms
c 260   continue
c       tcx1cr=tcx1cr+xspare
c     end if
c     write (*,*) ' TCX1CR--xnum,xden,xexp: ',xnum,xden,xexp
c     write (*,*) ' TCX1CR--t,tcx_crit: ',t,tcx1cr
c
c  multiply by reducing parameter (to convert units, etc.)
      TCX1CR=tcx1cr*tcxred(i)
c
      RETURN
      end                                             !subroutine TCX1CR
c
c ======================================================================
c
      subroutine SETTK3 (nread,icomp,hcasno,ierr,herr)
c
c  initialize model #3 for the thermal conductivity critical enhancement--
c  the simplified critical enhancement of Vesovic, et al for CO2:
c  Vesovic, V., Wakeham, W.A., Olchowy, G.A., Sengers, J.V., Watson, J.T.R.
c  and Millat, J. (1990). The transport properties of carbon dioxide. 
c  J. Phys. Chem. Ref. Data  19: 763-808.
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                        1 = error
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-08-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
c     parameter (n0=-ncmax,nx=ncmax)
      parameter (nx=ncmax)
      parameter (nrf0=0)    !lower limit for transport ref fluid arrays
      parameter (mxtck=40)  !max no. coefficients for t.c. crit
      character*1 htab,hnull
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ htab,hnull
c  limits and reducing parameters (separate reducing par for poly & exp)
      common /WLMTCK/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCK/ tred(nrf0:nx),Dred(nrf0:nx),tcxred(nrf0:nx),
     &                tredex(nrf0:nx),Dredex(nrf0:nx)
c  numbers of terms for the various parts of the model:
c  polynomial (numerator & denominator), exponential, spare
      common /WNTTCK/ nnum(nrf0:nx),nden(nrf0:nx),nexp(nrf0:nx),
     &                nspare(nrf0:nx)
c  commons storing the (real and integer) coefficients to the model
      common /WCFTCK/ ctck(nrf0:nx,mxtck,5)
      common /WIFTCK/ itck(nrf0:nx,mxtck,0:5)
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=39
        write (herr,1039) hcasno,hnull
 1039   format ('[SETUP error 39] block data option for t.c. crit ',
     &          'model #3 specified for CAS # ',a12,' but not ',
     &          'implemented in code.',a1)
      else
c  read data from file (should have been opened by SETUP)
c       write (*,*) ' SETTK3--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  # terms in numerator & denominator of polynomial multiplier,
c  exponential term, and spare for future use
c  put CO2 terms into arrays used for "numerator," others for future use
        read (nread,*) nnum(icomp),nden(icomp),nexp(icomp),nspare(icomp)
        if (nnum(icomp)+nden(icomp).ge.1) then
c  read reducing pars
          read (nread,*) tred(icomp),Dred(icomp),tcxred(icomp)
          if (nnum(icomp).ge.1) then
            do 200 j=1,nnum(icomp)              !read "numerator" terms
c  as originally implemented for CO2, these terms are:
c  ctck(icomp,1,1) = gnu (universal exponent, approx 0.63)
c  ctck(icomp,2,1) = gamma (universal exponent, approx 1.2145)
c  ctck(icomp,3,1) = R0 (universal amplitude, 1.01 +/- 0.04)
c  ctck(icomp,4,1) = z (universal exponent, 0.065 +/- 0.005)
c  ctck(icomp,5,1) = c (visc const, approx 1.075, but often set to 1)
c  ctck(icomp,6,1) = xi0 (amplitude, order 1d-10 m)
c  ctck(icomp,7,1) = gam0 (amplitude, order 0.05 - 0.06)
c  ctck(icomp,8,1) = qd_inverse (cutoff diameter, order 10d-9 m)
c  ctck(icomp,9,1) = tref (reference temperature, 1.5 - 2.0 * Tc)
            jterm=jterm+1
            read (nread,*) ctck(icomp,jterm,1)
  200       continue
          end if
        end if
      end if
c
      RETURN
      end                                             !subroutine SETTK3
c
c ======================================================================
c
      function TCX3CR (icomp,t,rho)
c
c  model #3 for the thermal conductivity critical enhancement--
c  the simplified critical enhancement of:
c  Olchowy, G.A. and Sengers, J.V. (1989). A simplified representation for
c  the thermal conductivity of fluids in the critical region.
c  Int. J. Thermophysics  10: 417-426.
c
c  also applied to CO2 by:
c  Vesovic, V., Wakeham, W.A., Olchowy, G.A., Sengers, J.V., Watson, J.T.R.
c  and Millat, J. (1990). The transport properties of carbon dioxide.
c  J. Phys. Chem. Ref. Data  19: 763-808.
c  equation numbers in comments refer to Vesovic paper
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   tcx3cr--the critical enhancement to the thermal conducitivity [W/m-K]
c
c  written by M. McLinden, NIST Thermophysics Division, Boulder, Colorado
c  06-08-97  MM, original version
c  06-16-97  MM, change DPDD, CVCP calls to DPDDK, CVCPK i.e. (t,x) to (icomp,t)
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 (mxtck=40)  !max no. coefficients for t.c. crit
      character*1 htab,hnull
      character*255 herr
c
      common /CHAR/ htab,hnull
c  limits and reducing parameters (separate reducing par for poly & exp)
      common /WLMTCK/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCK/ tred(nrf0:nx),Dred(nrf0:nx),tcxred(nrf0:nx),
     &                tredex(nrf0:nx),Dredex(nrf0:nx)
c  numbers of terms for the various parts of the model:
c  polynomial (numerator & denominator), exponential, spare
c  the "CO2" terms are stored in the "numerator" area
      common /WNTTCK/ nnum(nrf0:nx),nden(nrf0:nx),nexp(nrf0:nx),
     &                nspare(nrf0:nx)
c  commons storing the (real and integer) coefficients to the model
      common /WCFTCK/ ctck(nrf0:nx,mxtck,5)
      common /WIFTCK/ itck(nrf0:nx,mxtck,0:5)
c
      if (rho.lt.1.0d-6) then
c  critical enhancement is zero for low densities (avoid divide by zero)
        tcx3cr=0.0d0
        RETURN
      end if
c  recover the parameters from the storage array:
      gnu=ctck(icomp,1,1)      !nu (universal exponent, approx 0.63)
      gamma=ctck(icomp,2,1)    !gamma (universal exponent, approx 1.24)
      R0=ctck(icomp,3,1)       !R0 (universal amplitude, 1.01 +/- 0.04)
c  following two parameters not used here, but in file for future use
c     z=ctck(icomp,4,1)        !z (universal exponent, 0.065 +/- 0.005)
c     cvisc=ctck(icomp,5,1)    !visc const, approx 1.075, but often 1)
      xi0=ctck(icomp,6,1)      !xi0 (amplitude, order 1d-10 m)
      gam0=ctck(icomp,7,1)     !gam0 (amplitude, order 0.05 - 0.06)
      qd=1.0d0/ctck(icomp,8,1) !qd_inverse (cutoff dia, order 10d-9 m)
      tref=ctck(icomp,9,1)     !tref (reference temp, 1.5 - 2.0 * Tc)
c
      call INFO (icomp,wm,ttp,tnbp,tc,pc,Dc,Zc,acf,dip,Rgas)
      call DPDDK (icomp,t,rho,dpdrho)
c     write (*,*) ' TCX3CR--t,dpdrho: ',t,dpdrho
c  function chi (Eq 40 in Vesovic) evaluated at t,rho and tref,rho
c  Vesovic introduces a t/tc term which is absent in other papers
c     chi=pc/(Dc*Dc*tc)*rho*t/dpdrho               !Vesovic form
      chi=pc/(Dc*Dc)*rho/dpdrho
      call DPDDK (icomp,tref,rho,dpdrho)
c     chiref=pc/(Dc*Dc*tc)*rho*tref/dpdrho*tref/t  !Vesovic form
      chiref=pc/(Dc*Dc)*rho/dpdrho*tref/t
      delchi=chi-chiref
      if (delchi.le.0.0d0) then
c  delchi can go negative far from critical
c       write (*,*) ' TCX3CR--chi < chiref:  ',chi,chiref
        tcx3cr=0.0d0
        RETURN
      end if
c  function xi (Eq 46)
c  Olchowy, Vesovic put gam0 inside exponent, Krauss (R134a) puts this
c  term outside, but this yields incorrect values
      xi=xi0*(delchi/gam0)**(gnu/gamma)            !Vesovic form
c     xi=xi0/gam0*delchi**(gnu/gamma)              !Krauss' form
c     write (*,1150) chi,chiref,xi
c1150 format (1x,' TCX3CR--chi,chiref,xi:  ',3e14.6)
c  functions omega and omega_zero (Eqs 59 & 60)
c     write (*,1154) icomp,t,rho
c1154 format (1x,' TCX3CR--call CVCPK for icomp,t,rho =  ',i2,2e14.6)
      call CVCPK (icomp,t,rho,cv,cp)
c     write (*,1156) cv,cp
c1156 format (1x,' TCX3CR--cv,cp returned from CVCPK =  ',2e14.6)
      piinv=1.0d0/3.141592654d0
c     write (*,1160) cv,cp,qd,xi
c1160 format (1x,' TCX3CR--cv,cp,qd,xi:  ',4e14.6)
      xomg=2.0d0*piinv*((cp-cv)/cp*ATAN(qd*xi)+cv/cp*qd*xi)
      xomg0=2.0d0*piinv*(1.0d0-EXP(-1.0d0/(1.0d0/(qd*xi)
     &     +((qd*xi*Dc/rho)**2)/3.0d0)))
      call ETAK (icomp,t,rho,eta,ierr,herr)
c     write (*,1162) xomg,xomg0,eta
c1162 format (1x,' TCX3CR--omega,omega_0,eta:  ',3e14.6)
      boltz=Rgas/6.0221367d23       !Boltzman's const
c  factor of 1d9 in next equation to convert from mol/L --> mol/m**3
c  and from micro-Pa-s to Pa-s
      tcx3cr=rho*1.0d9*cp*R0*boltz*t*piinv/(6.0d0*eta*xi)*(xomg-xomg0)
c  multiply by reducing parameter (to convert units, etc.)
      TCX3CR=tcx3cr*tcxred(icomp)
c     write (*,*) ' TCX3CR--tcx3cr: ',TCX3CR
c
      RETURN
      end                                             !subroutine TCX3CR
c
c ======================================================================
c
      subroutine SETTK4 (nread,icomp,hcasno,ierr,herr)
c
c  initialize model #4 for the thermal conductivity critical enhancement--
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                       39 = 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-24-97  MM, original version
c  08-19-97  MM, error 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 (mxtck=40)  !max no. coefficients for t.c. crit
      character*1 htab,hnull
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ htab,hnull
c  limits and reducing parameters (separate reducing par for poly & exp)
      common /WLMTCK/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCK/ tred(nrf0:nx),Dred(nrf0:nx),tcxred(nrf0:nx),
     &                tredex(nrf0:nx),Dredex(nrf0:nx)
c  numbers of terms for the various parts of the model:
c  polynomial (numerator & denominator), exponential, spare
      common /WNTTCK/ nnum(nrf0:nx),nden(nrf0:nx),nexp(nrf0:nx),
     &                nspare(nrf0:nx)
c  commons storing the (real and integer) coefficients to the model
      common /WCFTCK/ ctck(nrf0:nx,mxtck,5)
      common /WIFTCK/ itck(nrf0:nx,mxtck,0:5)
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) icomp,nread,hcasno,hnull
 1101   format ('[SETTK4 error 101] illegal file specified for icomp =',
     &          i3,'; nread = ',i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        ierr=39
        herr='[SETUP error 39] model TK4 for the thermal conductivity '
     &     //'critical enhancement specified in fluid file '
     &     //'but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETTK4
c
c ======================================================================
c
      subroutine SETTK5 (nread,icomp,hcasno,ierr,herr)
c
c  initialize model #5 for the thermal conductivity critical enhancement--
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                       39 = 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-24-97  MM, original version
c  08-19-97  MM, error 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 (mxtck=40)  !max no. coefficients for t.c. crit
      character*1 htab,hnull
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ htab,hnull
c  limits and reducing parameters (separate reducing par for poly & exp)
      common /WLMTCK/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCK/ tred(nrf0:nx),Dred(nrf0:nx),tcxred(nrf0:nx),
     &                tredex(nrf0:nx),Dredex(nrf0:nx)
c  numbers of terms for the various parts of the model:
c  polynomial (numerator & denominator), exponential, spare
      common /WNTTCK/ nnum(nrf0:nx),nden(nrf0:nx),nexp(nrf0:nx),
     &                nspare(nrf0:nx)
c  commons storing the (real and integer) coefficients to the model
      common /WCFTCK/ ctck(nrf0:nx,mxtck,5)
      common /WIFTCK/ itck(nrf0:nx,mxtck,0:5)
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) icomp,nread,hcasno,hnull
 1101   format ('[SETTK5 error 101] illegal file specified for icomp =',
     &          i3,'; nread = ',i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        ierr=39
        herr='[SETUP error 39] model TK5 for the thermal conductivity '
     &     //'critical enhancement specified in fluid file '
     &     //'but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETTK5
c
c ======================================================================
c
      subroutine SETTK6 (nread,icomp,hcasno,ierr,herr)
c
c  initialize model #6 for the thermal conductivity critical enhancement--
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                       39 = 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-24-97  MM, original version
c  08-19-97  MM, error 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 (mxtck=40)  !max no. coefficients for t.c. crit
      character*1 htab,hnull
      character*12 hcasno
      character*255 herr
c
      common /CHAR/ htab,hnull
c  limits and reducing parameters (separate reducing par for poly & exp)
      common /WLMTCK/ tmin(nrf0:nx),tmax(nrf0:nx),pmax(nrf0:nx),
     &                rhomax(nrf0:nx)
      common /WRDTCK/ tred(nrf0:nx),Dred(nrf0:nx),tcxred(nrf0:nx),
     &                tredex(nrf0:nx),Dredex(nrf0:nx)
c  numbers of terms for the various parts of the model:
c  polynomial (numerator & denominator), exponential, spare
      common /WNTTCK/ nnum(nrf0:nx),nden(nrf0:nx),nexp(nrf0:nx),
     &                nspare(nrf0:nx)
c  commons storing the (real and integer) coefficients to the model
      common /WCFTCK/ ctck(nrf0:nx,mxtck,5)
      common /WIFTCK/ itck(nrf0:nx,mxtck,0:5)
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) icomp,nread,hcasno,hnull
 1101   format ('[SETTK6 error 101] illegal file specified for icomp =',
     &          i3,'; nread = ',i4,'; CAS no. = ',a12,a1)
      else
c  read data from file (should have been opened by SETUP)
        ierr=39
        herr='[SETUP error 39] model TK6 for the thermal conductivity '
     &     //'critical enhancement specified in fluid file '
     &     //'but not implemented in code.'//hnull
      end if
c
      RETURN
      end                                             !subroutine SETTK6
c
c ======================================================================
c
      function TCCNH3 (icomp,t,rho)
c
c  model for the thermal conductivity critical enhancement of ammonia
c  by the empirical model of Tufeu et al.
c
c  this is a special model only for ammonia with all constants "hardwired"
c
c  inputs:
c    icomp--component number in mixture (1..nc); 1 for pure fluid
c           [not used--included to maintain parallel structure]
c        t--temperature [K]
c      rho--molar density [mol/L]
c  output (as function value):
c   tccnh3--the critical enhancement to the thermal conducitivity [W/m-K]
c
c  written by S.A. Klein & M. McLinden,
c  NIST Thermophysics Division, Boulder, Colorado
c  02-26-97  MM, original version (adapted from TCENHC of SAK)
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)
      character*1 htab,hnull
c
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)
      common /CHAR/ htab,hnull
c
      rhokg=rho*wm(icomp)                 !correlation in mass units
      tr=abs(t-405.4d0)/405.4d0           !405.4 is Tc
      etab=1.0d-5*(2.6d0+1.6d0*tr)        !viscosity--Eq 9
      dPdT=1.0d5*(2.18d0-0.12d0/EXP(17.8d0*tr))  !dP/dT for rho = rhoc
c  Eq 9 of Tufeu for conductivity along critical isochore
      tcrhoc=1.2d0*1.38066d-23*t**2*dPdT**2*0.423d-8/(tr**1.24d0)*
     &   (1.0d0+1.429d0*tr**0.50d0)/(6.0d0*3.14159*etab*(1.34d-10/tr**
     &   (0.63d0)*(1.0d0+1.0d0*tr**0.50d0)))
      dtcid=tcrhoc*EXP(-36.0d0*tr**2)     !Eq 10
      xcon=0.61d0*253d0+16.5d0*log(tr)    !Eq 12
      if (rho/rhoc(icomp) .lt. 0.6d0) then
c  Eq 14 for rho < 0.6*rhoc (141 = 0.6*rhoc)
        tcc_sw=dtcid*xcon**2/(xcon**2+(141.0d0-0.96d0*141.0d0)**2)
        TCCNH3=tcc_sw*rhokg**2/141.0d0**2 !SAK had 253 in denominator
      else
c  Eq 11 for rho > 0.6*rhoc
        TCCNH3=dtcid*xcon**2/(xcon**2+(rhokg-0.96d0*253.0d0)**2)
      end if
c
      RETURN
      end                                             !subroutine TCCNH3

c
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c
c ======================================================================
c                                                    end file trns_TCX.f
c ======================================================================