diff --git a/EZFIO.tar.gz b/EZFIO.tar.gz
index 2b3689d..af469b9 100644
Binary files a/EZFIO.tar.gz and b/EZFIO.tar.gz differ
diff --git a/Makefile b/Makefile
index 5f6ea60..7edb4fc 100644
--- a/Makefile
+++ b/Makefile
@@ -1,3 +1,5 @@
+include version
+
 
 all: bin/eplf
 
@@ -6,6 +8,7 @@ EZFIO/config/eplf.config: EZFIO.tar.gz
 	cd EZFIO ; ./configure
 	if [ -e $@ ] ; then rm $@ ; fi
 	ln -s $$PWD/eplf.config $@
+	touch EZFIO/config/eplf.config
 
 EZFIO/lib/libezfio.so: EZFIO/config/eplf.config
 	make -C EZFIO/
@@ -24,5 +27,8 @@ clean:
 	- rm bin/ezfio.py
 	- make -C src veryclean
 
+
 archive:
-	git archive --format=tar HEAD | gzip > eplf.tar.gz
+	git archive --format=tar HEAD | gzip > eplf.$(VERSION).tar.gz
+	cp bin/to_ezfio.exe to_ezfio.$(VERSION).exe
+
diff --git a/configure b/configure
index 890bfc6..34f7873 100755
--- a/configure
+++ b/configure
@@ -2587,7 +2587,7 @@ echo "$as_me: error: Please download EZFIO at
 echo "$as_me: error: Download EZFIO.tar.gz into $PWD and run this configure script again." >&2;}
    { (exit 1); exit 1; }; }
   fi
-  mv Downloads/EZFIO*.tar.gz .
+  mv Downloads/EZFIO*.tar.gz EZFIO.tar.gz
 fi
 
 ######################################################
diff --git a/configure.ac b/configure.ac
index f7cefc3..1a6d30c 100644
--- a/configure.ac
+++ b/configure.ac
@@ -122,7 +122,7 @@ if [[ -z $EZFIO ]] ; then
   if [[ -z $EZFIO ]] ; then
     AC_MSG_ERROR([Download EZFIO.tar.gz into $PWD and run this configure script again.])
   fi
-  mv Downloads/EZFIO*.tar.gz .
+  mv Downloads/EZFIO*.tar.gz EZFIO.tar.gz
 fi
 
 ######################################################
diff --git a/src/aim.irp.f b/src/aim.irp.f
new file mode 100644
index 0000000..574cf1f
--- /dev/null
+++ b/src/aim.irp.f
@@ -0,0 +1,45 @@
+program aim
+ implicit none
+ provide mpi_rank
+
+ integer :: isize
+ double precision, allocatable :: population(:), variance(:)
+ isize = int(maxval(grid_density_partition))
+ allocate ( population(0:isize) )
+ allocate ( variance(0:isize) )
+ 
+ integer :: i,j,k,l
+ do i=1,isize
+  population(i) = 0.d0
+  variance(i) = 0.d0
+ enddo
+
+ do i=1,grid_z_num
+  do j=1,grid_y_num
+   do k=1,grid_x_num
+    l = int(grid_density_partition(k,j,i))
+    population(l) += grid_density(k,j,i)
+    variance(l) += grid_density(k,j,i)*grid_density(k,j,i)
+   enddo
+  enddo
+ enddo
+
+ real :: factor
+ factor = grid_step(1)*grid_step(2)*grid_step(3)
+ do i=1,isize
+  variance(i) -= population(i)*population(i)
+  population(i) *= factor
+  variance(i) *= factor*factor
+ enddo
+
+ print *, 'Basin  Population '
+ do i=0,isize
+  print *, i, population(i) !, abs(variance(i)), sqrt(abs(variance(i)))
+ enddo
+ print *, 'Total : ', sum(population)
+
+ deallocate (population, variance)
+
+ call finish()
+end
+
diff --git a/src/grid.irp.f b/src/grid.irp.f
index 16b833b..046ec1e 100644
--- a/src/grid.irp.f
+++ b/src/grid.irp.f
@@ -284,10 +284,12 @@ BEGIN_PROVIDER [ real, grid_$X_partition, (grid_x_num,grid_y_num,grid_z_num) ]
     enddo
    enddo
 
+   integer :: icount, imax
+   imax = 100
    do iz=1,grid_z_num 
     do iy=1,grid_y_num 
      do ix=1,grid_x_num 
-      call find_attribution(ix,iy,iz,grid_$X,grid_$X_partition)
+      call find_attribution(ix,iy,iz,grid_$X,grid_$X_partition,imax)
      enddo
     enddo
    enddo
@@ -313,36 +315,40 @@ BEGIN_PROVIDER [ real, current_partition_index ]
  current_partition_index = 0.
 END_PROVIDER
 
-recursive subroutine find_attribution(ix,iy,iz,g,a)
+recursive subroutine find_attribution(ix,iy,iz,g,a,depth)
  implicit none
 
  integer, intent(in) :: ix, iy, iz
  real, intent(in) :: g(grid_x_num,grid_y_num,grid_z_num)
  real, intent(inout) :: a(grid_x_num,grid_y_num,grid_z_num)
+ integer, intent(in) :: depth
 
  real :: buffer(-1:1,-1:1,-1:1)
 
  integer :: i,j,k
 
+ if (depth == 0) then
+   return
+ endif
  if (a(ix,iy,iz) /= 0.) then
   return
  else if (ix == 1) then
-  call find_attribution(2,iy,iz,g,a)
+  call find_attribution(2,iy,iz,g,a,depth-1)
   a(ix,iy,iz) = a(2,iy,iz)
  else if (ix == grid_x_num) then
-  call find_attribution(ix-1,iy,iz,g,a)
+  call find_attribution(ix-1,iy,iz,g,a,depth-1)
   a(ix,iy,iz) = a(ix-1,iy,iz)
  else if (iy == 1) then
-  call find_attribution(ix,2,iz,g,a)
+  call find_attribution(ix,2,iz,g,a,depth-1)
   a(ix,iy,iz) = a(ix,2,iz)
  else if (iy == grid_y_num) then
-  call find_attribution(ix,iy-1,iz,g,a)
+  call find_attribution(ix,iy-1,iz,g,a,depth-1)
   a(ix,iy,iz) = a(ix,iy-1,iz)
  else if (iz == 1) then
-  call find_attribution(ix,iy,2,g,a)
+  call find_attribution(ix,iy,2,g,a,depth-1)
   a(ix,iy,iz) = a(ix,iy,2)
  else if (iz == grid_z_num) then
-  call find_attribution(ix,iy,iz-1,g,a)
+  call find_attribution(ix,iy,iz-1,g,a,depth-1)
   a(ix,iy,iz) = a(ix,iy,iz-1)
 
  else
@@ -399,7 +405,7 @@ recursive subroutine find_attribution(ix,iy,iz,g,a)
      endif
    else
      ! Otherwise, get the partition index of the max value
-     call find_attribution(ix+im(1), iy+im(2), iz+im(3),g,a)
+     call find_attribution(ix+im(1), iy+im(2), iz+im(3),g,a,depth-1)
      a(ix,iy,iz) = a(ix+im(1), iy+im(2), iz+im(3))
    endif