subroutine conductivity(xtemp,xden,massfrac,cond,diff_coeff)

      use dBase, ONLY: ionmax
      
      use ModuleEos, ONLY: eos

! ### NU ADD #### 04.02.02 DJL
      use multifluid_database, ONLY: get_mfluid_property
! ### ### ### 

      implicit none
      save
! 
! this is Frank's routine that approximates the thermal transport 
! coefficients for electron conduction.  
!
! The calling sequence has been modified for use in FLASH.  Instead of 
! passing in the electron pressure and degeneracy pressure, this routine
! makes an eos call to the general FLASH eos driver (eos_fcn).  This 
! makes use of abar and zbar which were computed in this routine anyway.
! The mass fractions are passed in through the array massfrac instead of 
! xmass, to avoid using the common blocks.  The only information which
! comes from the common blocks is the list of isotope proton numbers and
! atomic masses (aion and zion) and the number of isotopes.
!
! input:
!  xtemp      = temperature temp (in K)
!  xden       = density den (in g/cm**3)
!  massfrac   = mass fractions of the composition
!
! output:
!  cond       = conductivity
!  diff_coeff = diffusion coefficient ( = cond/(rho*c_v)) 
!

!
! modified 7-9-99
!
! now returns the conductivity, not opacity
! eliminated zion, aion, and ionmax from the argument list, 
! these come from the network includes
!

! ### NU ADD ### 4.2.02. DJL
! Commented out the following ... seems to be a throwback from FLASH1.6
!      include 'network_common.fh'
! ### ### ###

! ### NU ADD ### 4.2.02. DJL
! Add declaration for zion and aion

     real, save :: zion(ionmax), aion(ionmax)

! ### ### ###


! declare the pass
      real xtemp, xden, massfrac(ionmax), cond, diff_coeff

! declare the variables used by the eos
      real pres, ener, dpt, dpd, det, ded, gammac

! declare the electron contributions -- these are used in this routine
      real pep, xne, eta

      real opac

      real c_p, chid

!..declare the internal variables
      integer iz,i
      real xmu,t6,orad,ocond,ytot1,ymassfrac,abar,zbar,w(6),xh, & 
     &     xhe,xz,xkc,xkap,xkb,xka,dbar,oiben1,d0log,xka1, & 
     &     xkw,xkaz,dbar1log,dbar2log,oiben2,t4,t4r,t44,t45,t46, & 
     &     ck1,ck3,ck2,ck4,ck5,ck6,xkcx,xkcy,xkcz,ochrs,th,fact, & 
     &     facetax,faceta,ocompt,tcut,cutfac,xkf,dlog10,zdel,zdell10, & 
     &     eta0,eta02,thpl,thpla,cfac1,cfac2,oh,pefac,pefacl,pefacal, & 
     &     dnefac,wpar2,walf,walf10,thx,thy,thc,farg,ffac,xmas,ymas, & 
     &     wfac,cint,vie,cie,tpe,yg,xrel,beta2,jy,vee,cee,ov1,ov, & 
     &     drel,drel10,drelim

! various physical and derived constants
!    con2 = con1*sqrt(4*pi*e*e/me)  
!    meff = hbar/(me*c)*(3*pi**2)**(1/3) 
!    weid = (pi*kerg)**2/(3*me)     
!    iec = 4*e**4*me/(3*pi*hbar**3) 
!    xec  = hbar/kerg*e*sqrt(4*pi/me) 
!
      real third,twoth,pi,avo,c,ssol,asol,zbound,t7peek, & 
     &     con2,k2c,meff,weid,iec,xec,rt3
      parameter (third  = 1.0e0/3.0e0,  & 
     &           twoth  = 2.0e0 * third,  & 
     &           pi     = 3.1415926535897932384e0,  & 
     &           avo    = 6.0221367e23,  & 
     &           c      = 2.99792458e10,  & 
     &           ssol   = 5.67050407222e-5,  & 
     &           asol   = 4.0e0*ssol/c,  & 
     &           zbound = 0.1e0,  & 
     &           t7peek = 1.0e20,  & 
     &           k2c    = 4.0e0/3.0e0*asol*c,  & 
     &           meff   = 1.194648642401440e-10,  & 
     &           weid   = 6.884326138694269e-5,  & 
     &           iec    = 1.754582332329132e16,   & 
     &           xec    = 4.309054377592449e-7, & 
     &           rt3    = 1.7320508075688772e0,  & 
     &           con2   = 1.07726359439811217e-7)

! conversion coefficient for opacity to conductivity
      real cond_factor
      parameter (cond_factor = 4.e0*asol*c/3.e0)


! ### NU ADD ### 04.02.02 DJL
! First call - fill in zion and aion
 
      logical, save :: firstCall = .TRUE.


!===========================================================================
 
      if (firstCall) then

! Define zion and aion via multifluid database
 
      call get_mfluid_property ("num positive", zion)
      call get_mfluid_property ("num total", aion)

      firstCall = .FALSE.

      endif

!===========================================================================

! ### ### ###



!..initialize
      opac      = 0.0e0
      orad      = 0.0e0
      ocond     = 0.0e0
      oiben1    = 0.0e0
      oiben2    = 0.0e0
      ochrs     = 0.0e0
      oh        = 0.0e0
      ov        = 0.0e0
      zbar      = 0.0e0
      ytot1     = 0.0e0

!..set the composition variables 
      do i=1,6
         w(i) = 0.0e0
      enddo

      do i = 1,ionmax 
         iz        = min(3,max(1,int(zion(i))))
         ymassfrac = massfrac(i)/aion(i)
         w(iz)     = w(iz) + massfrac(i) 
         w(iz+3)   = w(iz+3) + zion(i) * zion(i) * ymassfrac
         zbar      = zbar + zion(i) * ymassfrac
         ytot1     = ytot1 + ymassfrac
      enddo
      abar = 1.0e0/ytot1
      zbar = zbar * abar
      t6   = xtemp * 1.0e-6 
      xh   = w(1)
      xhe  = w(2)
      xz   = w(3)

! now that we have abar and zbar, call the eos and get the electron
! contributions
      call eos(xden,xtemp,pres,ener,massfrac,abar,zbar, & 
     &     dpt,dpd,det,ded,gammac,pep,xne,eta,1)

! now generate the specific heat at constant pressure from these quantities
! see C&G 9.98
      chid = dpd*xden/pres
      c_p = det*gammac/chid

!..radiative section:
!..from iben apj 196 525 1975
      if (xh .lt. 1.0e-5) then 
         xmu      = max(1.0e-99, w(4)+w(5)+w(6)-1.0e0)    
         xkc      = (2.019e-4*xden/t6**1.7e0)**(2.425e0) 
         xkap     = 1.0e0 + xkc * (1.0e0 + xkc/24.55e0)
         xkb      = 3.86e0 + 0.252e0*sqrt(xmu) + 0.018e0*xmu 
         xka      = 3.437e0 * (1.25e0 + 0.488e0*sqrt(xmu) + 0.092e0*xmu)
         dbar     = exp(-xka + xkb*log(t6))  
         oiben1   = xkap * (xden/dbar)**(0.67e0)
      end if 

      if ( .not.((xh.ge.1.0e-5) .and. (t6.lt.1.0)) .and.  & 
     &     .not.((xh.lt.1.0e-5) .and. (xz.gt.zbound)) ) then 
         if (t6 .gt. 1.0) then
            d0log = -(3.868e0 + 0.806e0*xh) + 1.8e0*log(t6)
         else 
            d0log = -(3.868e0 + 0.806e0*xh) +  & 
     &           (3.42e0 - 0.52e0*xh)*log(t6)
         endif
         xka1 = 2.809e0 * exp(-(1.74e0  - 0.755e0*xh)  & 
     &        * (log10(t6) - 0.22e0 + 0.1375e0*xh)**2)  
         xkw  = 4.05e0 * exp(-(0.306e0  - 0.04125e0*xh) & 
     &        * (log10(t6) - 0.18e0 + 0.1625e0*xh)**2) 
         xkaz = 50.0e0*xz*xka1 * exp(-0.5206e0* & 
     &        ((log(xden)-d0log)/xkw)**2)
         dbar2log = -(4.283e0 + 0.7196e0*xh) + 3.86e0*log(t6)
         dbar1log = -5.296e0 + 4.833e0*log(t6)  
         if (dbar2log .lt. dbar1log) dbar1log = dbar2log
         oiben2   = (xden/exp(dbar1log))**(0.67e0) * exp(xkaz) 
      end if 

!..from christy apj 144 108 1966
      if ((t6.lt.1.5) .and. (xh .ge. 1.0e-5)) then 
         t4    = xtemp * 1.0e-4 
         t4r   = sqrt(t4)  
         t44   = t4**4
         t45   = t44 * t4 
         t46   = t45 * t4 
         ck1   = 2.0e6/t44 + 2.1e0*t46  
         ck3   = 4.0e-3/t44 + 2.0e-4/xden**(0.25e0) 
         ck2   = 4.5e0*t46 + 1.0e0/(t4*ck3)    
         ck4   = 1.4e3*t4 + t46  
         ck5   = 1.0e6 + 0.1e0*t46  
         ck6   = 20.0e0*t4 + 5.0e0*t44 + t45  
         xkcx  = xh*(t4r/ck1 + 1.0e0/ck2)    
         xkcy  = xhe*(1.0e0/ck4 + 1.5e0/ck5)   
         xkcz  = xz*(t4r/ck6)
         ochrs = pep * (xkcx + xkcy + xkcz)
      end if 

!..opacity in presence of hydrogen  
      if (xh .ge. 1.0e-5) then 
         if (t6 .lt. 1.0) then 
            orad   = ochrs  
         else if (t6 .le. 1.5) then 
            orad   = 2.0e0*(ochrs*(1.5e0 - t6) + oiben2*(t6 - 1.0e0)) 
         else  
            orad   = oiben2    
         end if 

!..opacity in absence of hydrogen   
      else  
         if (xz .gt. zbound) then 
            orad   = oiben1   
         else 
            orad   = oiben1*(xz/zbound) + oiben2*((zbound-xz)/zbound)
         end if 
      end if 

!..add in the compton scattering opacity, weaver et al. apj 1978 225 1021
      th      = min(511.0e0, xtemp * 8.617e-8)
      fact    = 1.0e0 + 2.75e-2*th - 4.88e-5*th*th   
      facetax = 1.0e100    
      if (eta .le. 500.0) facetax = exp(0.522e0*eta - 1.563e0)  
      faceta  = 1.0e0 + facetax  
      ocompt  = 6.65205e-25/(fact * faceta) * xne/xden
      orad    = orad   + ocompt 

!..cutoff radiative opacity when 4kt/hbar is less than the plasma frequency  
      tcut = con2 * sqrt(xne) 
      if (xtemp .lt. tcut) then 
         if (tcut .gt. 200.0*xtemp) then 
            orad   = orad * 2.658e86  
         else 
            cutfac   = exp(tcut/xtemp - 1.0e0)  
            orad     = orad * cutfac
         end if 
      end if 

!..fudge molecular opacity for low temps
      xkf    = t7peek * xden * (xtemp * 1.0e-7)**4
      orad   = xkf * orad/(xkf + orad)


!..conductivity section:
!..drel is the dividing line between nondegenerate and degenerate regions,
!..taken from clayton eq. 2-34. if the density is larger than drel, then
!..use the degenerate expressions. if the density is smaller than
!..drelim, use the non-degenerate formulas. in between drel and drelim,
!..apply a smooth blending of the two.
      drel     =  2.4e-8 * zbar/abar * xtemp * sqrt(xtemp)
      if (xtemp .le. 1.0e5) drel = drel * 15.0e0
      drel10   = log10(drel)
      drelim   = drel10 + 0.3e0
      dlog10   = log10(xden)    


!..from iben apj 196 525 1975 for non-degenerate regimes
      if (dlog10 .lt. drelim) then 
         zdel    = xne/(avo*t6*sqrt(t6))
         zdell10 = log10(zdel)    
         eta0    = exp(-1.20322e0 + twoth * log(zdel))  
         eta02   = eta0*eta0

!..thpl factor 
         if (zdell10 .lt. 0.645) then 
            thpl    = -7.5668e0 + log(zdel * (1.0e0 + 0.024417e0*zdel))  
         else 
            if (zdell10 .lt. 2.5) then 
               thpl   = -7.58110e0 + log(zdel*(1.0e0 + 0.02804e0*zdel))  
               if (zdell10 .ge. 2.0) then 
                  thpla = thpl 
                  thpl  = -11.0742e0 +  & 
     &                 log(zdel**2 * (1.0e0 + 9.376e0/eta02))
                  thpl  = 2.0e0*((2.5e0-zdell10)*thpla +  & 
     &                 (zdell10-2.0e0)*thpl)   
               end if 
            else 
               thpl   = -11.0742e0 +  & 
     &              log(zdel**2 * (1.0e0 + 9.376e0/eta02))
            end if 
         end if 

!..pefac and walf factors
         if (zdell10 .lt. 2.0) then 
            pefac   = 1.0e0 + 0.021876e0*zdel 
            if (zdell10 .gt. 1.5) then 
               pefacal   = log(pefac) 
               pefacl    = log(0.4e0 * eta0 + 1.64496e0/eta0)  
               cfac1     = 2.0e0 - zdell10  
               cfac2     = zdell10 - 1.5e0  
               pefac     = exp(2.0e0 * (cfac1*pefacal + cfac2*pefacl))    
            end if 
         else 
            pefac   = 0.4e0 * eta0 + 1.64496e0/eta0
         end if 
         
         if (zdel.lt.40.0) then
            dnefac = 1.0e0 + zdel * (3.4838e-4 * zdel - 2.8966e-2)
         else
            dnefac = 1.5e0/eta0 * (1.0e0 - 0.8225e0/eta02)
         endif
         wpar2  = 9.24735e-3 * zdel *  & 
     &        (xden*avo*(w(4)+w(5)+w(6))/xne + dnefac)/(sqrt(t6)*pefac)
         walf   = 0.5e0 * log(wpar2)    
         walf10 = 0.5e0 * log10(wpar2)    
         
!..thx, thy and thc factors
         if (walf10 .le. -3.0) then 
            thx   = exp(2.413e0 - 0.124e0*walf) 
         else if (walf10 .le. -1.0) then 
            thx   = exp(0.299e0 - walf*(0.745e0 + 0.0456e0*walf))   
         else 
            thx   = exp(0.426e0 - 0.558e0*walf) 
         end if 
         
         if (walf10 .le. -3.0) then 
            thy   = exp(2.158e0 - 0.111e0*walf)  
         else if (walf10 .le. 0.0) then 
            thy   = exp(0.553e0 - walf*(0.55e0 + 0.0299e0*walf))  
         else 
            thy   = exp(0.553e0 - 0.6e0*walf)   
         end if 
         
         if (walf10 .le. -2.5) then 
            thc   = exp(2.924e0 - 0.1e0*walf)   
         else if (walf10 .le. 0.5) then 
            thc   = exp(1.6740e0 - walf*(0.511e0 + 0.0338e0*walf))  
         else 
            thc   = exp(1.941e0 - 0.785e0*walf) 
         end if 
         
         oh   = (xh*thx + xhe*thy + w(6)*third*thc) / (t6*exp(thpl)) 
      end if 

!..from yakovlev & urpin soviet astro 1980 24 303 and 
!..potekhin et al. 1997 aa 323 415 for degenerate regimes
      if (dlog10 .gt. drel10) then 
         xmas   = meff * xne**third 
         ymas   = sqrt(1.0e0 + xmas*xmas) 
         wfac   = weid * xtemp/ymas * xne 
         cint   = 1.0e0 

!..ion-electron collision frequency and the thermal conductivity 
         vie   = iec * zbar * ymas * cint
         cie   = wfac/vie 

!..electron-electron collision frequency and thermal conductivity 
         tpe  = xec * sqrt(xne/ymas)
         yg   = rt3 * tpe/xtemp 
         xrel = 1.009e0 * (zbar/abar * xden * 1.0e-6)**third
         beta2 = xrel**2/(1.0e0 + xrel**2)
         jy   = (1.0e0 + 6.0e0/(5.0e0*xrel**2) + 2.0e0/(5.0e0*xrel**4)) & 
     &        * ( yg**3 / (3.0e0 * (1.0e0 + 0.07414 * yg)**3)  & 
     &        * log((2.81e0 - 0.810*beta2 + yg)/yg)     & 
     &        + pi**5/6.0e0 * (yg/(13.91e0 + yg))**4 )
         vee = 0.511e0 * xtemp**2 * xmas/ymas**2 * sqrt(xmas/ymas) * jy 
         cee = wfac/vee

!..total electron thermal conductivity and conversion to an opacity 
         ov1   = cie * cee/(cee + cie)
         ov    = k2c/(ov1*xden) * xtemp**3 
      end if 

!..blend the opacities in the intermediate region
      if (dlog10 .le. drel10) then 
         ocond   = oh 
      else if (dlog10 .gt. drel10  .and. dlog10 .lt. drelim) then 
         farg    = pi * (dlog10 - drel10) / 0.3e0
         ffac    = 0.5e0 * (1.0e0 - cos(farg))
         ocond   = exp((1.0e0-ffac)*log(oh) + ffac*log(ov))
      else if (dlog10 .ge. drelim) then 
         ocond   = ov 
      end if 

!..total opacity
      opac    = orad * ocond / (ocond + orad)

! now return the conductivity
      cond = cond_factor*xtemp**3/(opac*xden)

      diff_coeff = cond/(xden*c_p)

      return 
      end