diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90
index 13ab743..dfc9709 100644
--- a/m_fdm_calc.f90
+++ b/m_fdm_calc.f90
@@ -22,6 +22,11 @@ module m_fdm_calc
       real*8 :: umax,umin,vmax,vmin,wmax,wmin  ! J. Kwon
       real*8, dimension(3) :: velmax, velmin, velmax1, velmin1
 
+   integer :: vtype1, vtype2
+   integer, dimension(128) :: scnt1, sdisp1, stype1, rcnt1, rdisp1, rtype1
+   integer, dimension(128) :: scnt2, sdisp2, stype2, rcnt2, rdisp2, rtype2
+
+
 !===========================================================================
 !===========================================================================
 
@@ -88,6 +93,34 @@ module m_fdm_calc
 
       CALL ludcmp(nx,ny,nz_all,1,0,0)
 
+      CALL MPI_TYPE_VECTOR (nz, nx*nz, nx*ny, MPI_REAL8, vtype1, mpi_err)
+      CALL MPI_TYPE_COMMIT (vtype1, mpi_err)
+
+      CALL MPI_TYPE_VECTOR (nz, 2*nx*nz, 2*nx*ny, MPI_REAL8, vtype2, mpi_err)
+      CALL MPI_TYPE_COMMIT (vtype2, mpi_err)
+
+
+      do i = 1, (ny/nz)
+         scnt1(i) = 1
+         rcnt1(i) = nx * nz * nz
+
+         sdisp1(i) = (i-1) * nx * nz  * 8
+         rdisp1(i) = (i-1) * nx * nz  * nz * 8
+
+         stype1(i) = vtype1
+         rtype1(i) = MPI_REAL8
+
+         scnt2(i) = 1
+         rcnt2(i) = 2 * nx * nz * nz
+
+         sdisp2(i) = (i-1) * 2 * nx * nz  * 8
+         rdisp2(i) = (i-1) * 2 * nx * nz  * nz * 8
+
+         stype2(i) = vtype2
+         rtype2(i) = MPI_REAL8
+
+      end do
+
       refwr=pre*1.*exp(-ac/(1.+bc*c_ref))  ! Kwon
 
 
@@ -341,6 +374,62 @@ module m_fdm_calc
       real*8 :: uz(4,yy),duz(4,yy),d2uz(yy)
       real*8 :: Ly(yy,zz),Dy
 
+      real*8 :: r1_t(2,xx,zz,yy),f_t(2,xx,zz,yy)
+      real*8 :: ww_t(xx,zz,yy)
+
+
+   !CALL MPI_ALLTOALLW (ww_, scnt1, sdisp1, stype1, ww_t, rcnt1, rdisp1, rtype1, MPI_COMM_TASK, mpi_err)
+
+   !CALL MPI_ALLTOALLW (r1_, scnt2, sdisp2, stype2, r1_t, rcnt2, rdisp2, rtype2, MPI_COMM_TASK, mpi_err)
+
+    do i = 1,numprocs
+        CALL MPI_ISEND(ww_(1,zz*(i-1)+1,1), scnt1(i), stype1(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
+    end do
+    do i = 1,numprocs
+        CALL MPI_RECV(ww_t(1,1,zz*(i-1)+1), rcnt1(i), rtype1(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
+    end do
+
+    do i = 1,numprocs
+        CALL MPI_ISEND(r1_(1,1,zz*(i-1)+1,1), scnt2(i), stype2(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
+    end do
+    do i = 1,numprocs
+        CALL MPI_RECV(r1_t(1,1,1,zz*(i-1)+1), rcnt2(i), rtype2(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
+    end do
+
+!! z-direction
+      DO j=1,zz
+       DO i=1,xx
+        DO k=1,yy
+         uz(1,k)=r1_t(1,i,j,k)                            ! 1:rho
+         uz(2,k)=r1_t(2,i,j,k)/r1_t(1,i,j,k)               ! 2:Y
+         uz(3,k)=uz(1,k)*ww_t(i,j,k)                       ! 3:rho*w
+         uz(4,k)=uz(3,k)*uz(2,k)                         ! 4:rho*w*Y
+        ENDDO
+
+       CALL dfp(yy,hy,uz(1:4,:),duz(1:4,:),4,3)
+       CALL d2fp(yy,hy,uz(2,:),d2uz(:),1,3)
+
+        DO k=1,yy
+! -( d(rho*w)/dz )
+         f_t(1,i,j,k)=-duz(3,k)         ! continuity
+
+! -( 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_t(2,i,j,k)=-duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k))   ! species conserv.
+        ENDDO
+       ENDDO
+      ENDDO
+
+   !CALL MPI_ALLTOALLW (f_t, rcnt2, rdisp2, rtype2, f_, scnt2, sdisp2, stype2, MPI_COMM_TASK, mpi_err)
+
+    do i = 1,numprocs
+        CALL MPI_ISEND(f_t(1,1,1,zz*(i-1)+1), rcnt2(i), rtype2(i), (i-1), myid, MPI_COMM_TASK, mpi_request, mpi_err)
+    end do
+    do i = 1,numprocs
+        CALL MPI_RECV(f_(1,1,zz*(i-1)+1,1), scnt2(i), stype2(i), (i-1), (i-1), MPI_COMM_TASK, mpi_status, mpi_err)
+    end do
+
 ! x-direction
       DO k=1,zz
        DO j=1,yy
@@ -364,12 +453,12 @@ module m_fdm_calc
  
 
 ! -( d(rho*u)/dx )
-         f_(1,i,j,k)=-dux(3,i)   ! continuity
+         f_(1,i,j,k) = f_(1,i,j,k) - dux(3,i)   ! continuity
 
 ! -( 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 )
-         f_(2,i,j,k)=-dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) - wrate  ! species conservation
+         f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) - wrate  ! species conservation
 
         ENDDO
        ENDDO