c ... SNL, Sophia Lefantzi 
c ... Return array with dimensional mixture diffusion coefficients at 
c ... center of cell (i,j)
c ... units in (m^2/s)

      subroutine getdifcoef(temperature,pressure,patm,wref,yspl,ns,
     1                      ospwt,capd,drtdvl)


      IMPLICIT NONE
#include "parameter.com"
#include "drfm_com.com"
#include "defs.h"

      integer ns, k
      double precision vtemp, temperature, wb, wref, rp0, patm, pressure
      
      real*8 capd(ns)
      real*8 ospwt(ns) 
      real*8 yspl(ns), xspl(ns)
      real*8 drfdml(ns), drtdvl(ns), dcma(ns)
      real*8 drfm_bdc(ns,ns), drfm_bdc1(ns*ns), bindkj(ns,ns)


      vtemp = temperature

      call gmlfr(ns,yspl,xspl,wb,ospwt)	! get mole fractions

      wb = wb * wref

c ... DRFM mixture averaged
      call drfdfk(ns,drfm_xmol,xspl,vtemp,wb,
     +            drfdml,drtdvl,drfm_xml2,drfm_sig,drfm_epsm,
     +            drfm_delm,drfm_chim,drfm_vstm,drfm_rstm,
     +            drfm_error)

c ... Chemkin mixture averaged
c      call getmadc(vtemp,pressure,ns,xspl,dcma)
  
      rp0 = patm / pressure

      do k=1,ns
         capd(k) = rp0*drfdml(k)/10000.d0 ! return dimensional D_k_mix
      end do

c ... Find the error between DRFM and chemkin mixtured averaged
      
c ... Chemkin binary diffusion coefficients
c      call getbidij(vtemp,pressure,ns,bindkj)
c      do i=1,ns
c      do j=1,ns
c         bindkj(i,j) = rp0*bindkj(i,j)
c         write(*,*) 'chemkin_bdc(',i,',',j,')=',bindkj(i,j)
c      enddo
c      enddo

c ... DRFM binary diffusion coefficients
c      call drfbns(ns,drfm_xmol,xspl,vtemp,
c     +            drfm_bdc1,drtdvl,drfm_xml2,drfm_sig,drfm_epsm,
c     +            drfm_delm,drfm_chim,drfm_vstm,drfm_rstm,
c     +            drfm_error)

c      do j=1,ns
c      do i=1,ns
c         drfm_bdc(i,j) = rp0*drfm_bdc1( i + (j-1)*ns )
c         write(*,*) 'bdc(',i,',',j,')=',drfm_bdc(i,j)*1.d-4
c         write(*,*) drfm_bdc(i,j)*1.d-4
c      enddo
c      enddo

      return
      end