C ***************************************************************
C **                                                           **
C ** routine to return binary diffusion coefficients           ** 
C **                                                           ** 
C ** SNL / Original subroutine 'drfdfk' by phpaul modified by  **
C ** Sophia Lefantzi on 3/28/02 to return the binary diffusion **
C ** coefficients                                              ** 
C **                                                           **
C *************************************************************** 

      subroutine drfbns(ns,xmol,wmix,t,bindc,drtdv,xml2,
     &      sig,epsm,delm,chim,vstm,rstm,error) 
C ****************************************************************
C **                                                            ** 
C                     SNL / phpaul all rights reserved 1997     ** 
C                     last modified Oct/20/97                   ** 
C **************************************************************** 
C inputs: 
C        NS - number of species 
C      XMOL - NS long vector containing molecular weights 
C      WMIX - ns long vector containing mixture fractions 
C    OLMASS - mixture molecular weight 
C         T - temperature (K) 
C       *** - set of molecular parameter arrays from DRINIT 
C returns for primary species: 
C     DRFDM - (NS,NS) array of binary diffusion coefficients 
C             (cm*cm/sec) as reduced to one atmosphere pressure 
C     DRTDV - NS long thermal diffusion ratio  (dimensionless)  
C             per Kihara model and sum rule of Chapman and Cowling  
C ****************************************************************  

      IMPLICIT NONE

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

C passed arrays       
      logical error 
      real*8 sig(KD),epsm(KD,KD),delm(KD,KD),chim(KD,KD) 
      real*8 vstm(KD,KD),rstm(KD,KD),xmol(KD) 
      real*8 wmix(KD),drfdm(KD),drtdv(KD) 
      real*8 xml2(kd,kd), bindc(KD**2)
C internal working arrays  
      real*8 sgx(KD),og2(KD),xs11(KD,KD),fmix(KD) 
      real*8 dcor(KD,KD),alp(KD,KD) 

      integer j, ns, k, jj, kk

      double precision sq2, delta, cndif, o6, zz, t,  zz01, zz02 
      double precision tr, vr, rr, sgr, og1j, og2j, obj, ocj
      double precision oej, og1jk, og2jk, objk, ocjk, oejk, hp1
      double precision hn01, astar, bstar, cstar, ff, y, cjk
      double precision an01, aw, bw

C square root of 2
      data sq2/1.4142135623730951d0/
c delta
      DATA DELTA/1.0d-14/ 
C  DELTA = a small number which is added to the mixture fractions 
C          to prevent divide by zero under pure species conditions   
      DATA CNDIF /2.6280D-3/ 
C  CNDIF = scaling constant for Dij at 1 atm. in cm^2/s 
C for single species 

c ... o6 = 1./6.
      data o6/0.166666666666666667d0/

      do 20 j = 1,ns 

        fmix(j) = wmix(j)+delta 
        ZZ      = 1.D0 + CHIM(j,j) + DELM(j,j)/T  
	zz01    = dexp(o6*dlog(zz))
	zz02    = zz*zz
        TR      = T/(EPSM(j,j)*zz02) 
        VR      = VSTM(j,j)*zz02
        RR      = RSTM(j,j)*zz01
        SGR     = sig(j)/zz01

        CALL DRXIOM12(TR,RR,VR,OG1J,OG2J,OBJ,OCJ,OEJ,ERROR)   

        sgx(j)    = sgr * sgr 
        xs11(j,j) = og1j * sgx(j)

#ifndef FIRST_ORDER_DRFM
        og2(j) = og2j 
	hn01   = 6.d0*ocj-5.d0
        dcor(j,j) = 1.3d0*(hn01*hn01)/58.d0 
        alp(j,j) = 0.d0 
#endif
20      continue           

C for the binary pairs 

      do 40 jj = 2,ns 

        do 30 kk = 1,jj-1        

C generate primatives           

          ZZ   = 1.0D0 + CHIM(kk,jj) + DELM(kk,jj)/T  
	  zz01 = dexp(o6*dlog(zz))
	  zz02 = zz*zz
          TR   = T/(EPSM(kk,jj)*zz02)
          VR   = VSTM(kk,jj)*zz02
          RR   = RSTM(kk,jj)*zz01
          SGR  = 0.5d0*(sig(jj)+sig(kk))/zz01

          CALL DRXIOM12(TR,RR,VR,OG1JK,OG2JK,OBJK,OCJK,OEJK,ERROR)            

          xs11(kk,jj) = og1jk*sgr*sgr  
          xs11(jj,kk) = xs11(kk,jj)          

#ifndef FIRST_ORDER_DRFM
          astar = og2jk/og1jk 
          bstar = 5.d0-2.4d0*objk 
          cstar = 6.d0*ocjk - 5.d0 
C order according to moelcular weights Mj > Mk         
          if(xmol(jj).gt.xmol(kk))then 
            j = jj 
            k = kk 
            ff = 1.0d0 
          else 
            j = kk 
            k = jj 
            ff = -1.0d0 
          endif           
C Kihara correction to the binary diffusion coefficient           
          y = xmol(k)/xmol(j) 
          cjk = fmix(j)/(fmix(j)+fmix(k)) 
	  an01 = 1.0d0+1.8d0*y
          aw = sq2*og1jk/(8.0d0*og2(k)*(an01*an01)) 
          bw = 10.0d0*(1.0d0+1.8d0*y+3.0d0*y*y) 
          dcor(jj,kk) = 1.3d0*(cstar*cstar)* 
     &      aw*cjk/(1.d0+(aw*bw-1.d0)*cjk)          
          dcor(kk,jj) = dcor(jj,kk)          
#endif

C thermal diffusion factor per Kihara           
C         z = 1.d0/y 
C         zp1 = z + 1.d0
C         zm1 = z - 1.d0
C         zr  = zm1/zp1
C         ckj =fmix(k)/(fmix(j)+fmix(k))        
C         hpj = sgx(j)*og2(j)/xs11(jj,kk) 
C         hpk = sgx(k)*og2(k)/xs11(jj,kk) 
C         s1 = z*DSQRT(2.d0/zp1)*hpj +  
C    &      (7.5d0*zm1 - 4.d0*z*astar) / (zp1*zp1) 
C         q1 = (2.d0/zp1)*DSQRT(2.d0/zp1)*hpj* 
C    &      (0.5d0*bstar*z*z+3.d0+1.6d0*z*astar) 
C         q12=15.d0*(zr*zr)*(0.5d0*bstar) 
C    &        + 4.d0*z*astar*(6.d0+bstar)/(zp1*zp1)           
C    &        + 1.6d0*zp1*hpj*hpk/DSQRT(z) 
C         s2 = DSQRT(2.d0/(z+z*z))*hpk +  
C    &      (7.5d0*(z-z*z) - 4.d0*z*astar) / (zp1*zp1)       
C         q2 = (2.d0/zp1)*DSQRT(2.d0/(z+z*z))*hpk* 
C    &      ((0.5d0*bstar)+3.d0*z*z+1.6d0*z*astar) 
C         alp(jj,kk) = ff*cstar*(cjk*s1 - ckj*s2) /  
C    &            (cjk*cjk*q1 + ckj*ckj*q2 + cjk*ckj*q12) 
C         alp(kk,jj) = (-1.0d0)*alp(jj,kk)       

30        continue 

40      continue         

        do 200 j = 1,ns 
           do 100 k = 1,ns
              hp1 = 2.d0*xmol(k)*xmol(j)/(xmol(k)+xmol(j)) 
              bindc(k+(j-1)*ns) = cndif*t*dsqrt(t/hp1)
     &                            * (1.d0+dcor(k,j)) / xs11(k,j)
 100       continue
 200    continue

        return 
        end