diffusivity calculation using sutherland formula

m_fdm_calc.f90

@@ -503,18 +503,25 @@ module m_fdm_calc
       real*8 :: uz(4,zz),duz(4,zz),d2uz(zz)
       real*8 :: uuz(zz)
       real*8 :: y
+      real*8 :: T0, T1, Ts, As
 
 
       f_(:,:,:,1) = 0.0 ! continuity
 
 ! diffusivity
+      T0 = 1.0
+      T1 = (1.0 + bc)
+
+      Ts = (dm_ratio*T1 - (T1**(3./2.))) / (T1**(3./2.) - dm_ratio)
+      As = (T0 + Ts)
+
       DO k=1,zz
        DO j=1,yy
         DO i=1,xx
 
          y=(1.0 - r1_(i,j,k,2)/r1_(i,j,k,1)) * bc + 1.0
 
-         dm(i,j,k) = diff * (y ** 0.76)
+         dm(i,j,k) = diff * As * sqrt(y) / (1. + Ts/y)
 
         ENDDO
        ENDDO

m_parameters.f90

@@ -114,6 +114,7 @@ module m_parameters
   real*8 :: ac,bc,c_cut,c_ref,min_wr,prof_wr,hx,hy,diff,pre
   real*8 :: dt_fullsave,dt_fdmsave,t_fullsave,t_fdmsave
   real*8 :: fdmcfl,startad,dummyu_,schmidt
+  real*8 :: dm_ratio
   real*8 :: bgn_save
 
 ! files for statistics
@@ -581,6 +582,9 @@ contains
     read(in,*,ERR=9000,END=9000) fdmcfl
     write(out,'(a30,3x,e20.10)') 'FDM_CFL     =',fdmcfl
 
+    read(in,*,ERR=9000,END=9000) dm_ratio
+    write(out,'(a30,3x,e20.10)') 'Dm_Ratio    =',dm_ratio
+
    ! initialize ----------------------
       t_fdmsave=dt_fdmsave
       t_fullsave=dt_fullsave