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