diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90
index b3ee762..b81c8c1 100644
--- a/m_fdm_calc.f90
+++ b/m_fdm_calc.f90
@@ -9,6 +9,7 @@ module m_fdm_calc
 !variables
 
       real*8, dimension(:,:,:), allocatable :: u_,v_,w_
+      real*8, dimension(:,:,:), allocatable :: dm
       real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf
       real*8, dimension(:,:), allocatable :: fzzl, fzzu
 
@@ -233,6 +234,10 @@ module m_fdm_calc
       v_=0.0
       w_=0.0
 
+      allocate(dm(nx,ny,nz))
+
+      dm = diff
+
 ! DQ initializing
       fdmsavecount=1        !FDM save count
       sum_wrate=0.
@@ -529,7 +534,7 @@ module m_fdm_calc
       DO j=1,yy
        DO i=1,xx
         DO k=1,zz
-         uz(1,k)=r1_(i,j,k,1)                             ! 1:rho
+         uz(1,k)=r1_(i,j,k,1)*dm(i,j,k)                   ! 1:rho*D
          uz(2,k)=r1_(i,j,k,2)/r1_(i,j,k,1)                ! 2:Y
          uz(3,k)=uz(1,k)*ww_(i,j,k)                       ! 3:rho*w
          uz(4,k)=uz(3,k)*uz(2,k)                          ! 4:rho*w*Y
@@ -545,8 +550,8 @@ module m_fdm_calc
 
 ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
 ! = -( d(rho*w*Yr)/dz )
-!    + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
-         f_(i,j,k,2) = f_(i,j,k,2) - duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k))   ! species conserv.
+!    + (rho*D * d2(Yr)/dz2 + d(rho*D)/dz * d(Yr)/dz )
+         f_(i,j,k,2) = f_(i,j,k,2) - duz(4,k) + (uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k))   ! species conserv.
         ENDDO
        ENDDO
       ENDDO
@@ -563,7 +568,7 @@ module m_fdm_calc
 
 ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
 ! = -( d(rho*w*Yr)/dz )
-!    + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
+!    + (rho*D * d2(Yr)/dz2 + d(rho*D)/dz * d(Yr)/dz )
 
       fbuf1(:,:,:) = r1_(:,:,:,2)*ww_(:,:,:)        ! rho*w*Y
       CALL pdfp (fbuf1,  fzzl,  fzzu,  fbuf2,  hy, xx*yy, zz, yy, 3)
@@ -574,13 +579,13 @@ module m_fdm_calc
       CALL pdfp (fbuf1,  fzzl,  fzzu,  fbuf2,  hy, xx*yy, zz, yy, 3)
       CALL pd2fp(fbuf1,  fzzl,  fzzu,  fbuf3,  hy, xx*yy, zz, yy, 3)
 
-      fbuf1(:,:,:) = r1_(:,:,:,1)                   ! rho
+      fbuf1(:,:,:) = r1_(:,:,:,1)                   ! rho*D
       CALL pdfp (fbuf1,  fzzl,  fzzu,  fbuf4,  hy, xx*yy, zz, yy, 3)
 
-      fbuf2 = fbuf2 * fbuf4
-      fbuf1 = fbuf1 * fbuf3 + fbuf2
+      fbuf2 = fbuf2 * fbuf4                         ! dY/dz * d(rho*D)/dz
+      fbuf1 = fbuf1 * fbuf3 + fbuf2                 ! rho*D * d2(rho*D)/dz2 + ...
 
-      f_(:,:,:,2) = f_(:,:,:,2) + diff*fbuf1(:,:,:)  ! species conserv.
+      f_(:,:,:,2) = f_(:,:,:,2) + fbuf1(:,:,:)      ! species conserv.
 
     endif
 
@@ -595,7 +600,7 @@ module m_fdm_calc
 
 ! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx
 ! = -( d(rho*u*Yr)/dx )
-!    + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx )
+!    + (rho*D * d2(Yr)/dx2 + d(rho*D)/dx * d(Yr)/dx )
 
         CALL tp2mul (yxbuf1, r1_(:,:,k,2), uu_(:,:,k), yy, xx)  ! rho*u*Y
         CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1)                   ! d/dx(rho*u*Y)
@@ -605,11 +610,11 @@ module m_fdm_calc
         CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1)                   ! d/dx(Y)
         CALL d2fnonp(xx,hx,yxbuf1,yxbuf3,yy,1)                  ! d2/dx2(Y)
 
-        CALL tp2 (yxbuf1, r1_(:,:,k,1), yy, xx)                 ! rho
-        CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1)                   ! d/dx(rho)
+        CALL tp2mul (yxbuf1, r1_(:,:,k,1), dm(:,:,k), yy, xx)   ! rho*D
+        CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1)                   ! d/dx(rho*D)
 
-        yxbuf2 = yxbuf2 * yxbuf4                                ! d/dx(Y) * d/dx(rho)
-        yxbuf1 = -diff*(yxbuf2 + yxbuf1 * yxbuf3)               ! -D( ... + rho * d2/dx2(Y))
+        yxbuf2 = yxbuf2 * yxbuf4                                ! d/dx(Y) * d/dx(rho*D)
+        yxbuf1 = -(yxbuf2 + yxbuf1 * yxbuf3)                    ! -( ... + rho*D * d2/dx2(Y))
 
         CALL tp2subasgn (f_(:,:,k,2), yxbuf1, xx, yy)           ! species conservation
 
@@ -628,7 +633,7 @@ module m_fdm_calc
 
 ! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy
 ! = -( d(rho*v*Yr)/dy ) 
-!    + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy )
+!    + (rho*D* d2(Yr)/dyy2 + d(rho*D)/dy * d(Yr)/dy )
 
         xybuf1(:,:) = r1_(:,:,k,2)*vv_(:,:,k)   ! rho*v*Y
         CALL dfp(yy,hy,xybuf1,xybuf2,xx,2)      ! d/dy(rho*v*Y)
@@ -639,13 +644,13 @@ module m_fdm_calc
         CALL dfp(yy,hy,xybuf1,xybuf2,xx,2)      ! d/dy(Y)
         CALL d2fp(yy,hy,xybuf1,xybuf3,xx,2)     ! d2/dy2(Y)
 
-        xybuf1(:,:) = r1_(:,:,k,1)              ! rho
-        CALL dfp(yy,hy,xybuf1,xybuf4,xx,2)      ! d/dy(rho)
+        xybuf1(:,:) = r1_(:,:,k,1) * dm(:,:,k)  ! rho*D
+        CALL dfp(yy,hy,xybuf1,xybuf4,xx,2)      ! d/dy(rho*D)
 
-        xybuf2 = xybuf2 * xybuf4                ! d/dy(Y) * d/dy(rho)
-        xybuf1 = xybuf2 + xybuf1 * xybuf3       ! ... + rho * d2/dy2(Y)
+        xybuf2 = xybuf2 * xybuf4                ! d/dy(Y) * d/dy(rho*D)
+        xybuf1 = xybuf2 + xybuf1 * xybuf3       ! ... + rho*D * d2/dy2(Y)
 
-        f_(:,:,k,2) = f_(:,:,k,2) + diff*xybuf1(:,:)   ! species conserv.
+        f_(:,:,k,2) = f_(:,:,k,2) + xybuf1(:,:)   ! species conserv.
 
       ENDDO