diff --git a/code/m_calculate.f90 b/code/m_calculate.f90
new file mode 100644
index 0000000..d5a2a38
--- /dev/null
+++ b/code/m_calculate.f90
@@ -0,0 +1,279 @@
+module m_calculate
+
+  use m_parameters
+  use m_compact
+
+  implicit none
+
+    real*8, allocatable, dimension(:,:) :: xsrc
+    real*8, allocatable, dimension(:,:) :: xdst
+
+    real*8, allocatable, dimension(:,:) :: rsrc
+    real*8, allocatable, dimension(:,:) :: rdst
+
+    integer, parameter :: nb = BLOCKSIZE
+
+    private :: nb
+    private :: xsrc, xdst, rsrc, rdst
+
+contains
+
+     subroutine m_calculate_init
+
+     integer :: ierr
+
+     call ludcmp(nxp,nyp,nzp,1,0,0)      ! 1,1,0
+
+     allocate(xsrc(nb, nxp), stat=ierr)
+     allocate(xdst(nb, nxp), stat=ierr)
+
+     allocate(rsrc(nxp, nzp), stat=ierr)
+     allocate(rdst(nxp, nzp), stat=ierr)
+
+     end subroutine m_calculate_init
+
+
+     subroutine m_calculate_finalize
+
+     deallocate(xsrc)
+     deallocate(xdst)
+
+     deallocate(rsrc)
+     deallocate(rdst)
+
+     end subroutine m_calculate_finalize
+
+
+     subroutine ddx(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp), intent(in) :: src
+        real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
+
+        integer :: i, j ,k
+        integer :: ju
+
+
+        do k = 1,nzp
+        do j = 1,nyp,nb
+
+        ju = min(j+nb-1,nyp)
+
+        call tp2(xsrc, src(:,j:ju,k), nb, nxp)
+
+        call dfnonp(nxp, hxp, xsrc, xdst, nb, 1)
+
+        call tp2(dst(:,j:ju,k), xdst, nxp, nb)
+
+        end do
+        end do
+
+     end subroutine ddx
+
+
+     subroutine ddy(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp), intent(in) :: src
+        real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
+
+        integer :: i, j ,k
+
+
+        do k = 1,nzp
+
+        call dfp(nyp, hyp, src(:,:,k), dst(:,:,k), nxp, 2)
+
+        end do
+
+     end subroutine ddy
+
+
+     subroutine ddz(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp), intent(in) :: src
+        real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
+
+        integer :: i, j ,k
+
+        do j = 1,nyp
+
+        do k = 1,nzp
+        rsrc(:,k) = src(:,j,k)
+        end do
+
+        call dfp(nzp, hzp, rsrc, rdst, nxp, 3)
+
+        do k = 1,nzp
+        dst(:,j,k) = rdst(:,k)
+        end do
+
+        end do
+
+     end subroutine ddz
+
+
+     subroutine d2dx(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp), intent(in) :: src
+        real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
+
+        integer :: i, j ,k
+        integer :: ju
+
+
+        do k = 1,nzp
+        do j = 1,nyp,nb
+
+        ju = min(j+nb-1,nyp)
+
+        call tp2(xsrc, src(:,j:ju,k), nb, nxp)
+
+        call d2fnonp(nxp, hxp, xsrc, xdst, nb, 1)
+
+        call tp2(dst(:,j:ju,k), xdst, nxp, nb)
+
+        end do
+        end do
+
+     end subroutine d2dx
+
+
+     subroutine d2dy(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp), intent(in) :: src
+        real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
+
+        integer :: i, j ,k
+
+
+        do k = 1,nzp
+
+        call d2fp(nyp, hyp, src(:,:,k), dst(:,:,k), nxp, 2)
+
+        end do
+
+     end subroutine d2dy
+
+
+     subroutine d2dz(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp), intent(in) :: src
+        real*8, dimension(nxp,nyp,nzp), intent(out) :: dst
+
+        integer :: i, j ,k
+
+        do j = 1,nyp
+
+        do k = 1,nzp
+        rsrc(:,k) = src(:,j,k)
+        end do
+
+        call d2fp(nzp, hzp, rsrc, rdst, nxp, 3)
+
+        do k = 1,nzp
+        dst(:,j,k) = rdst(:,k)
+        end do
+
+        end do
+
+     end subroutine d2dz
+
+
+     subroutine tp(a, b, nx)
+       ! a(nb,nx) = transpose(b(nx,nb))
+
+        integer,intent(in) :: nx
+
+        real*8,intent(out) :: a(nb,nx)
+        real*8,intent(in) :: b(nx,nb)
+
+        call tp2(a, b, nb, nx)
+
+     end subroutine tp
+
+     subroutine tp2 (a, b, n1, n2)
+       ! a = transpose(b)
+
+         implicit none
+
+         integer,intent(in) :: n1, n2
+         real*8,intent(out) :: a(n1,n2)
+         real*8,intent(in) :: b(n2,n1)
+         integer :: i,j,ii,jj
+
+            DO jj=1,n2,nb
+             DO ii=1,n1,nb
+
+            DO j=jj,min(jj+nb-1,n2)
+             DO i=ii,min(ii+nb-1,n1)
+               a(i,j) = b(j,i)
+             ENDDO
+            ENDDO
+
+             ENDDO
+            ENDDO
+
+     end subroutine tp2
+
+     real function rxn_rate (c)
+
+      real :: c
+
+        if(c.lt.0.) c=0.
+        if(c.gt.1.) c=1.
+
+        if (c.le.c_cut) then
+
+          rxn_rate = min_wr
+
+        else if (c.gt.c_ref) then
+
+          rxn_rate = pre*(1.-c)*exp(-ac/(1.+bc*c))
+
+        else
+
+          rxn_rate = &
+            ((refwr-min_wr)*exp(prof_wr*(c-c_ref)) + min_wr - refwr*exp(prof_wr*(c_cut-c_ref))) &
+            / (1.-exp(prof_wr*(c_cut-c_ref)))
+
+        endif
+
+     end function rxn_rate
+
+     real function threshold_min_max (c, minc, maxc)
+
+      real :: c
+      real :: minc, maxc
+
+        if ((c.lt.minc) .or. (c.gt.maxc)) then
+            threshold_min_max = 0.
+        else
+            threshold_min_max = 1.0
+        end if
+
+     end function threshold_min_max
+
+     real function positive (c)
+
+      real :: c
+
+        if (c > 0.0d0) then
+            positive = c
+        else
+            positive = 0.
+        end if
+
+     end function positive
+
+     real function negative (c)
+
+      real :: c
+
+        if (c < 0.0d0) then
+            negative = c
+        else
+            negative = 0.
+        end if
+
+     end function negative
+
+end module m_calculate
diff --git a/code/m_compact.f90 b/code/m_compact.f90
index 231aefd..1d534e7 100644
--- a/code/m_compact.f90
+++ b/code/m_compact.f90
@@ -1,14 +1,13 @@
     MODULE m_compact
      IMPLICIT NONE
-     PRIVATE
-     REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & 
-                                        lyf,lys,wyf,wys, &
-                                        lzf,lzs,wzf,wzs
-!                                           lyzf,lyzs,wyzf,wyzs
+
+     REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
+                                                lyf,lys,wyf,wys, &
+                                                lzf,lzs,wzf,wzs
+
      INTEGER :: nxc,nyc,nzc
-     REAL, PARAMETER :: ezero = 1.0e-14
-      
-     PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
+
+     REAL(KIND=8), PARAMETER :: ezero = 1.0e-14
 
     CONTAINS
 
@@ -20,6 +19,22 @@
        nxc=nx
        nyc=ny
        nzc=nz
+
+       CALL ludcmp_allocate(nx,ny,nz,xp,yp,zp)
+
+       CALL ludcmp_calculate(nx,ny,nz,xp,yp,zp)
+
+     END SUBROUTINE ludcmp
+
+     SUBROUTINE ludcmp_allocate(nx,ny,nz,xp,yp,zp)
+      INTEGER, INTENT(IN) :: nx,ny,nz
+      INTEGER, INTENT(IN) :: xp,yp,zp
+      INTEGER :: ierr
+
+       nxc=nx
+       nyc=ny
+       nzc=nz
+
 !       IF(nyc /= nzc) PRINT*,'ny should be equal nz'
 
 !      xp, yp, zp = 0 : periodic
@@ -32,9 +47,6 @@
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
          ALLOCATE(wxs(nxc),STAT=ierr)
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
-         CALL p_lud(1,nxc)
-       ELSE
-         CALL nonp_lud(1,nxc)
        ENDIF
 
        ALLOCATE(lyf(nyc),STAT=ierr)
@@ -46,9 +58,6 @@
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
          ALLOCATE(wys(nyc),STAT=ierr)
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
-         CALL p_lud(2,nyc)
-       ELSE
-         call nonp_lud(2,nyc)
        ENDIF
 
        ALLOCATE(lzf(nzc),STAT=ierr)
@@ -60,19 +69,140 @@
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
          ALLOCATE(wzs(nzc),STAT=ierr)
          IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
+       ENDIF
+
+     END SUBROUTINE ludcmp_allocate
+
+     SUBROUTINE ludcmp_deallocate(xp,yp,zp)
+      INTEGER, INTENT(IN) :: xp,yp,zp
+
+!       IF(nyc /= nzc) PRINT*,'ny should be equal nz'
+
+!      xp, yp, zp = 0 : periodic
+       DEALLOCATE(lxf)
+       DEALLOCATE(lxs)
+       IF(xp.eq.0) THEN
+         DEALLOCATE(wxf)
+         DEALLOCATE(wxs)
+       ENDIF
+
+       DEALLOCATE(lyf)
+       DEALLOCATE(lys)
+       IF(yp.eq.0) THEN
+         DEALLOCATE(wyf)
+         DEALLOCATE(wys)
+       ENDIF
+
+       DEALLOCATE(lzf)
+       DEALLOCATE(lzs)
+       IF(zp.eq.0) THEN
+         DEALLOCATE(wzf)
+         DEALLOCATE(wzs)
+       ENDIF
+
+     END SUBROUTINE ludcmp_deallocate
+
+     SUBROUTINE ludcmp_testalloc
+
+       IF (.not. ALLOCATED(lxf)) print *, "lxf not allocated"
+       IF (.not. ALLOCATED(lxs)) print *, "lxs not allocated"
+       IF (.not. ALLOCATED(wxf)) print *, "wxf not allocated"
+       IF (.not. ALLOCATED(wxs)) print *, "wxs not allocated"
+
+       IF (.not. ALLOCATED(lyf)) print *, "lyf not allocated"
+       IF (.not. ALLOCATED(lys)) print *, "lys not allocated"
+       IF (.not. ALLOCATED(wyf)) print *, "wyf not allocated"
+       IF (.not. ALLOCATED(wys)) print *, "wys not allocated"
+
+       IF (.not. ALLOCATED(lzf)) print *, "lzf not allocated"
+       IF (.not. ALLOCATED(lzs)) print *, "lzs not allocated"
+       IF (.not. ALLOCATED(wzf)) print *, "wzf not allocated"
+       IF (.not. ALLOCATED(wzs)) print *, "wzs not allocated"
+
+     END SUBROUTINE ludcmp_testalloc
+
+     SUBROUTINE ludcmp_calculate(nx,ny,nz,xp,yp,zp)
+      INTEGER, INTENT(IN) :: nx,ny,nz
+      INTEGER, INTENT(IN) :: xp,yp,zp
+      INTEGER :: ierr
+
+       nxc=nx
+       nyc=ny
+       nzc=nz
+
+       ! CALL ludcmp_testalloc
+
+!       IF(nyc /= nzc) PRINT*,'ny should be equal nz'
+
+!      xp, yp, zp = 0 : periodic
+       IF(xp.eq.0) THEN
+         CALL p_lud(1,nxc)
+       ELSE
+         CALL nonp_lud(1,nxc)
+       ENDIF
+
+       IF(yp.eq.0) THEN
+         CALL p_lud(2,nyc)
+       ELSE
+         call nonp_lud(2,nyc)
+       ENDIF
+
+       IF(zp.eq.0) THEN
          CALL p_lud(3,nzc)
        ELSE 
          call nonp_lud(3,nzc)
        ENDIF
 
-!       CALL x_lud
-!       CALL yz_lud
+     END SUBROUTINE ludcmp_calculate
+
+
+     SUBROUTINE test_nonp_lud1(xx, coef)
+      INTEGER :: xx
+      REAL(KIND=8), DIMENSION(xx) :: aa
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
+       aa=3.
+       aa(1)=0.5 ; aa(2)=4.
+       aa(xx-1)=4. ; aa(xx)=0.5
+
+       CALL stdlu(aa,xx,coef)
+     END SUBROUTINE test_nonp_lud1
+
+
+     SUBROUTINE test_nonp_lud2(xx, coef)
+      INTEGER :: xx
+      REAL(KIND=8), DIMENSION(xx) :: aa
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
+       aa=5.5
+       aa(1)=2./11. ; aa(2)=10.
+       aa(xx-1)=10. ; aa(xx)=2./11.
+
+       CALL stdlu(aa,xx,coef)
+     END SUBROUTINE test_nonp_lud2
+
+
+     SUBROUTINE test_p_lud1(xx, coef1, coef2)
+      INTEGER :: xx
+      REAL(KIND=8) :: a
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
+      a=3. ! first derivative
+
+      CALL ptdlu(a,xx,coef1,coef2) ! x-direction
+     END SUBROUTINE test_p_lud1
+
+
+     SUBROUTINE test_p_lud2(xx, coef1, coef2)
+      INTEGER :: xx
+      REAL(KIND=8) :: a
+      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
+      a=11./2. ! second derivative
+
+      CALL ptdlu(a,xx,coef1,coef2) ! x-direction
+     END SUBROUTINE test_p_lud2
 
-     END SUBROUTINE ludcmp
 
      SUBROUTINE nonp_lud(xyz,xx)
       INTEGER :: i,xyz,xx
-      REAL, DIMENSION(xx) :: aa
+      REAL(KIND=8), DIMENSION(xx) :: aa
        aa=3.
        aa(1)=0.5 ; aa(2)=4.
        aa(xx-1)=4. ; aa(xx)=0.5
@@ -92,7 +222,7 @@
 
      SUBROUTINE p_lud(xyz,xx)
       INTEGER :: i,xyz,xx
-      REAL :: a
+      REAL(KIND=8) :: a
       a=3. ! first derivative
       IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction
       IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction
@@ -105,21 +235,23 @@
 
      SUBROUTINE stdlu(a,n,l)
       INTEGER :: n
-      REAL :: a(n),l(n)
-      REAL :: d
+      REAL(KIND=8), INTENT(IN)  :: a(n)
+      REAL(KIND=8), INTENT(OUT) :: l(n)
+      REAL(KIND=8) :: d
       INTEGER :: i
-      l(1)=1.0/a(1)
+      l(1)=1.0d0/a(1)
       DO i=2,n
         d=a(i)-l(i-1)
-        l(i)=1.0/d
+        l(i)=1.0d0/d
       ENDDO
      END SUBROUTINE stdlu
 
      SUBROUTINE ptdlu(a,n,l,w)
       INTEGER :: n
-      REAL :: a,l(n),w(n)
+      REAL(KIND=8), INTENT(OUT) :: a
+      REAL(KIND=8), INTENT(OUT) :: l(n),w(n)
       INTEGER :: i
-      REAL :: aa(n),d
+      REAL(KIND=8) :: aa(n),d
 
       DO i=1,n-1
         aa(i)=a
@@ -141,18 +273,19 @@
       l(n)=1./d
      END SUBROUTINE ptdlu
 
-     SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
-      INTEGER,INTENT(IN) :: n,nd,dir
-      REAL,INTENT(IN) :: h
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+
+     SUBROUTINE rhs1np(n,h,x,dx,nd)
+      INTEGER,INTENT(IN) :: n,nd
+      REAL(KIND=8),INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
+      REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
 
-      h1=1./h
+      h1=1.d0/h
 
-      r1=7./3.
-      r2=1./12.
+      r1=7.d0/3.d0
+      r2=1.d0/12.d0
       r3=3.
       a=-1.25
       b=1.
@@ -178,6 +311,20 @@
         dx(j,i)=h1*(r1*t1+r2*t2)
        ENDDO
       ENDDO
+
+     END SUBROUTINE rhs1np
+
+
+     SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
+      INTEGER,INTENT(IN) :: n,nd,dir
+      REAL(KIND=8),INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
+      INTEGER :: i,j
+      REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
+
+      CALL rhs1np (n,h,x,dx,nd)
+
       IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction
       IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction
       IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction
@@ -185,11 +332,14 @@
 
      SUBROUTINE dfp(n,h,x,dx,nd,dir)
       INTEGER,INTENT(IN) :: n,nd,dir
-      REAL,INTENT(IN) :: h
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+      REAL(KIND=8),INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: r1,r2,h1
+      REAL(KIND=8) :: r1,r2,h1
+
+
+      ! print *, "dfnonp received (nd,n)", nd, n
 
       h1=1./h
       r1=7./3.
@@ -221,10 +371,10 @@
 
      SUBROUTINE ptdslv(r,n,l,w,nd)
       INTEGER,INTENT(IN) :: n,nd
-      REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
-      REAL,INTENT(IN),DIMENSION(:) :: l,w
+      REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
+      REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w
       INTEGER i,j
-      REAL, DIMENSION(nd) :: sum
+      REAL(KIND=8), DIMENSION(nd) :: sum
       DO j=1,nd
         sum(j)=w(1)*r(j,1)
         r(j,1)=r(j,1)*l(1)
@@ -249,11 +399,13 @@
 
      SUBROUTINE d2fp(n,h,x,dx,nd,dir)
       INTEGER,INTENT(IN) :: n,nd,dir
-      REAL,INTENT(IN) :: h
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+      REAL(KIND=8),INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: h2,r1,r2,t1,t2
+      REAL(KIND=8) :: h2,r1,r2,t1,t2
+
+
       h2=1./(h*h)
       r1=6.
       r2=3./8.
@@ -312,10 +464,10 @@
 
      SUBROUTINE tdslv(r,n,l,nd)
       INTEGER,INTENT(IN) :: n,nd
-      REAL,INTENT(INOUT),DIMENSION(nd,n) :: r
-      REAL,INTENT(IN),DIMENSION(:) :: l
+      REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
+      REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l
       INTEGER i,j
-      REAL t1
+      REAL(KIND=8) t1
       DO j=1,nd
         r(j,1)=r(j,1)*l(1)
       ENDDO
@@ -334,11 +486,12 @@
 
      SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
       INTEGER,INTENT(IN) :: n,nd,dir
-      REAL,INTENT(IN) :: h
-      REAL,INTENT(IN),DIMENSION(nd,n) :: x
-      REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
+      REAL(KIND=8),INTENT(IN) :: h
+      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
+      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
       INTEGER :: i,j
-      REAL :: h2,r1,r2,r3,a,b,c,e,t1,t2
+      REAL(KIND=8) :: h2,r1,r2,r3,a,b,c,e,t1,t2
+
 
       h2=1./(h*h)
       r1=6.
diff --git a/code/m_parameters.f90 b/code/m_parameters.f90
new file mode 100644
index 0000000..34c7dd5
--- /dev/null
+++ b/code/m_parameters.f90
@@ -0,0 +1,36 @@
+module m_parameters
+
+  implicit none
+
+     integer :: nxp,nyp,nzp
+
+     real :: hxp,hyp,hzp
+
+     real :: l_0
+
+     integer :: startnum,endnum,skipnum,shiftnum
+
+     real :: scp,prp,lep,vis0p,rod
+
+     real :: prof_wr,min_wr,min_fsd,min_c,refwr
+
+     real :: pre,ac,bc,c_cut,c_ref
+
+     integer :: syp,eyp,twod
+
+     real :: SL_u
+
+     integer :: omitnum
+
+     integer, allocatable :: omit_t(:,:)
+
+     integer, allocatable :: file_dist(:)
+
+     integer, allocatable :: export_offset(:)
+
+     real, parameter :: pi=3.14159265358979323846
+     real, parameter :: me=1.00e-20
+
+contains
+
+end module m_parameters
diff --git a/code/makefile b/code/makefile
index 0507951..b49f9ce 100644
--- a/code/makefile
+++ b/code/makefile
@@ -1,4 +1,6 @@
-flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8  -march=native 
+BLOCKSIZE?=32
+
+flags = -Wall -O3 -cpp -DBLOCKSIZE=$(BLOCKSIZE) -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8  -march=native 
 
 ifdef CHECK
 flags += -fcheck=all
@@ -10,18 +12,24 @@ endif
 
 compiler = gfortran
 
-ex  : test.o ysolve.o m_compact.o
+ex  : test.o ysolve.o m_parameters.o m_calculate.o m_compact.o
 	${compiler} -o ex test.o ysolve.o m_compact.o 
 
-test.o  : test.f90 ysolve.mod compact.mod
+test.o  : test.f90 ysolve.mod m_compact.mod 
 	${compiler} -c ${flags} test.f90
 
-ysolve.o ysolve.mod : ysolve.f90 compact.mod
+ysolve.o ysolve.mod : ysolve.f90 m_compact.mod
 	${compiler} -c ${flags} ysolve.f90
 
-m_compact.o compact.mod : m_compact.f90
+m_compact.o m_compact.mod : m_compact.f90
 	${compiler} -c ${flags} m_compact.f90
 
+m_calculate.o m_calculate.mod : m_calculate.f90
+	${compiler} -c ${flags} m_calculate.f90
+
+m_parameters.o m_parameters.mod : m_parameters.f90
+	${compiler} -c ${flags} m_parameters.f90
+
 test: ex
 	sh test.sh