diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90
index b81c8c1..694117a 100644
--- a/m_fdm_calc.f90
+++ b/m_fdm_calc.f90
@@ -505,6 +505,21 @@ module m_fdm_calc
 
       f_(:,:,:,1) = 0.0 ! continuity
 
+! diffusivity
+      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)
+
+        ENDDO
+       ENDDO
+      ENDDO
+
+      if (myid.eq.0) write(*,*) 'min(dm)',minval(dm),'max(dm)',maxval(dm)
+
 ! reaction source term
       DO k=1,zz
        DO j=1,yy