subroutine drfstc(ns,xmol,wmix,t,drftc,drfvs,xml1,
     &      sig,epsm,delm,chim,vstm,rstm,tvta,itox,error)
C ** routine to calculate mixture viscosity per WILKE '50 rule  **
C ** & mixture thermal conductivity per MASON-SAXENA '58 rule.  **  
C                                                               **
C                     SNL / phpaul all rights reserved 1997     **
C                     last modified Oct/22/97                   **
C ****************************************************************
C inputs:
C      NS - number of species
C       T - temperature (K)
C    WMIX - NS long vector containing mixture fractions
C    XMOL - NS long vector containing molecular weights
C     *** - set of molecular parameter arrays from DRINIT
C returns:
C   DRFTC - mixture thermal conductivity (milliwatts/meter Kelvin)
C   DRFVS - mixture dynamic viscosity (microPascal seconds)
C ****************************************************************

      IMPLICIT NONE

#include "defs.h"
      INCLUDE 'drfmcm.for'

C passed arrays      
      logical error,itox(kd)
      real*8 sig(kd),epsm(kd,kd),delm(kd,kd),chim(kd,kd)
      real*8 vstm(kd,kd),rstm(kd,kd),tvta(kd,6),xmol(kd),wmix(kd)
      real*8 xml1(kd,kd)
C internal working arrays      
      real*8 osgx(kd),sqosgx(kd),eta(kd)


      integer j, ns, k
      double precision  CNVSC, RGAS, ds8, hpcon 


      double precision zz, t, zz01, zz02, tr, vr, rr, sgr 
      double precision sqog2j, og2j, obj, ocj, oej, hp1, sumt 
      double precision drfvs, u, sumg, osqj, gn01, tcmult 
      double precision sds8, hp2, drftc

C   
C  CNVSC = scaling constant for viscosity in micro_Pa sec
      DATA CNVSC /2.6695D0/

C  RGAS = universal gas constant in joules / mol K 
      DATA RGAS /8.31441D0/

c  ds8 = 1.d0/dsqrt(8.d0)
      data ds8 /0.35355339059327373D0/

c  hpcon = 3.d0/196.d0
      data hpcon /0.15306122448979591D-01/

C generate property primitives
 
      do 10 j = 1,ns

        ZZ   = 1.D0 + CHIM(j,j) + DELM(j,j)/T
	zz01 = zz**0.166666667D0
	zz02 = zz*zz
        TR   = T/(EPSM(j,j)*zz02)
        VR   = VSTM(j,j)*zz02
        RR   = RSTM(j,j)*zz01
        SGR  = SIG(j)/zz01

        CALL DRXIOM2(TR,RR,VR,sqog2j,OG2J,OBJ,OCJ,OEJ,ERROR)

        osgx(j) = og2j * sgr * sgr 
	sqosgx(j) = sqog2j * sgr

#ifdef FIRST_ORDER_DRFM
        hp1 = 1.d0
#else
        hn01 = 8.d0*oej-7.d0
        hp1 = 1.d0+hpcon*(hn01*hn01)
#endif

        eta(j) = cnvsc*DSQRT(xmol(j)*T)*hp1/osgx(j)

10      continue

C compute mixture properties

      sumt = 0.d0
      drfvs = 0.d0
      u = dlog(T)

      do 100 j = 1,ns

        sumg = 0.d0
	osqj = 1.0d0 / sqosgx(j)

	if(j.gt.1)then
	   do k=1,j-1
             gn01 = 1.d0 + sqosgx(k) * osqj
             sumg = sumg + wmix(k) * gn01 * gn01 * xml1(k,j)
	   end do
        endif

	if(j.lt.ns)then
	   do k=j+1,ns
             gn01 = 1.d0 + sqosgx(k) * osqj
             sumg = sumg + wmix(k) * gn01 * gn01 * xml1(k,j)
	   end do
	endif

        if(itox(j))then 
          tcmult=(tvta(j,1)+u*(tvta(j,3)+u*tvta(j,5)))/
     &      (1.d0+u*(tvta(j,2)+u*(tvta(j,4)+u*tvta(j,6))))
        else
            tcmult = 1.d0
        endif

	sds8  = sumg * ds8
        hp2   = wmix(j)*eta(j)/(wmix(j)+1.065d0*sds8)
        drfvs = drfvs + wmix(j)*eta(j)/(wmix(j)+sds8)
        sumt  = sumt + tcmult*hp2/xmol(j)

100     continue

      drftc = 3.75D0 * rgas * sumt

      return
      end
c*******************************************************************************

       SUBROUTINE drxiom2(TR,RS,VS,sqog22,OG22,OBX,OCX,OEX,ERROR) 

C**  Collision integrals using LJ for TR<10 and BM for TR>10  ** 
C**  coefs. from Bzowski et al. JPCRD v19 p1179 (1990).       ** 
C**  Range for 0.2<TR<1.0 has been fit to tables from         ** 
C**  Mason et al. with f, f' and f'' continuous at TR=1.      ** 
C                                                             ** 
C                  SNL / PHPaul  all rights reserved 1997     ** 
C                  last modified Sept/6/97                    ** 
C*************************************************************** 
C inputs: 
C  TR = reduced temperature 
C  RS = reduced Born-Mayer radius 
C  VS = reduced Born-Mayer energy 
C returns: 
C  ERROR = .FALSE. for normal end 
C        = .TRUE. for TR < 0.2 user should treat as fatal 
C  sqOG22 = sqrt(og22)
C  OG22 = normalized collision integral. type (2,2) 
C  OEX = higher collision integral functional Estar 
C  OCX = higher collision integral functional Cstar 
C  OBX = higher collision integral functional Bstar 
C*************************************************************** 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      LOGICAL ERROR 
      DIMENSION BT(4),G1(3),G2(3),G3(3),G(5),GL(6) 
      DIMENSION DT(5),F1(3),F2(3),F3(3),F(5),FL(6)        
      DATA BT /0.89D0,-267.0D0,26700.0D0,-8.9D5/ 
      DATA G1 /201.570D0,174.672D0,7.36916D0/ 
      DATA G2 /19.2265D0,27.6938D0,3.29559D0/ 
      DATA G3 /6.31013D0,10.2266D0,2.33033D0/ 
      DATA G /0.295402D0,-0.510069D0,0.189395D0, 
     &         -0.045427D0,0.0037928D0/ 
      DATA GL /0.295402D0,-0.510069D0,0.189395D0, 
     &          0.484463D0,0.417806D0,0.122148D0/  
      DATA DT /1.04D0,0.0D0,-33.0838,101.571,-87.7036/ 
      DATA F1 /20.0862D0,72.1059D0,8.27648D0/ 
      DATA F2 /56.4472D0,286.393D0,17.7610D0/ 
      DATA F3 /46.3130D0,277.146D0,19.0573D0/ 
      DATA F /0.46641D0,-0.56991D0,0.19591D0, 
     &         -0.03879D0,0.00259D0/ 
      DATA FL /0.46641D0,-0.56991D0,0.19591D0, 
     &          0.747363D0,0.662153D0,0.188447D0/      
c---------------------------------------------------------------------------------
      IF(TR.LT.0.2D0) THEN 
        ERROR = .TRUE. 
      ELSE 
        ERROR = .FALSE. 
      ENDIF 
c---------------------------------------------------------------------------------
C 0.2 <= TR < 1.0 

      IF (TR.LT.1.0D0) THEN        

        AT = DLOG(TR)        

C  type II           

        sqog22 = DEXP( 0.5d0 * (
     &	   FL(1)+AT*(FL(2)+AT*(FL(3)+AT*(FL(4)+ AT*(FL(5)+AT*FL(6)))))
     &                          )) 
        og22   = sqog22 * sqog22

#ifndef FIRST_ORDER_DRFM
        OEX = 1.D0 +( FL(2)+AT*(2.D0*FL(3)+AT* 
     &   (3.D0*FL(4)+AT*(4.D0*FL(5)+5.D0*AT*FL(6)))) )/4.D0           
#endif

        RETURN 
        ENDIF 

c---------------------------------------------------------------------------------
C 1 <= TR <= 10.0 

      IF (TR.LE.10.D0) THEN 

        AT = DLOG(TR) 

C  type II         

        sqog22=DEXP(0.5d0*(F(1)+AT*(F(2)+AT*(F(3)+AT*(F(4)+AT*F(5)))))) 
	og22  = sqog22 * sqog22

#ifndef FIRST_ORDER_DRFM
        OEX = 1.D0 + ( F(2)+AT*(2.D0*F(3)+AT* 
     &           (3.D0*F(4)+AT*4.D0*F(5))) )/4.D0              
#endif
        RETURN 
        ENDIF 

c---------------------------------------------------------------------------------
C TR > 10.0     

      A10 = DLOG(0.1D0*VS) 
      AA  = DLOG(VS/TR) 

C  type II         

      AT   = 1.D0/DLOG(TR) 
      fn01 = f1(3)/a10
      fn02 = f2(3)/a10
      fn03 = f3(3)/a10
      B1   = F1(1)+F1(2)/A10+fn01*fn01
      B2   = (F2(1)+F2(2)/A10+fn02*fn02)*(-1.D0) 
      B3   = (F3(1)+F3(2)/A10+fn03*fn03) 
      SUM1 = RS*RS*(DT(1)+AT*AT*(DT(3)+AT*(DT(4)+AT*DT(5)))) 
      SUM2 = AT*AT*(B1+AT*(B2+AT*B3))/(A10*A10) 
      OG22 = (SUM1 + SUM2) * (AA*AA) 

      sqog22 = dsqrt(og22)

#ifndef FIRST_ORDER_DRFM
      SUM3 = RS*RS*AT*AT*AT*(2.D0*DT(3)+AT*(3.D0*DT(4)+ 
     &          4.D0*AT*DT(5))) 
      SUM4 = AT*AT*AT*(2.D0*B1+AT*(3.D0*B2+AT*4.D0*B3))/(A10*A10) 
      OEX = 1.D0-( 2.D0/AA + (SUM3+SUM4)/(SUM1+SUM2) )/4.D0         
#endif

      RETURN  
      END