anchor.diff:    part_count = 0
     binpart = 0
   
  -!  call cpu_time(t1)
     do i = 1, n
       binx = NINT((px(i))/(BOX_SIZE)*(real(nbinx,WP)-1.0_WP)) + 1
  -!    print *, px(i), real(nbinx,WP)-1.0_WP,NINT((px(i))/(BOX_SIZE)*(real(nbinx,WP)-1.0_WP)) + 1
       biny = NINT((py(i))/(BOX_SIZE)*(real(nbiny,WP)-1.0_WP)) + 1
       binz = NINT((pz(i))/(BOX_SIZE)*(real(nbinz,WP)-1.0_WP)) + 1
       part_count(binx, biny, binz) = part_count(binx, biny, binz) + 1
       binpart(part_count(binx,biny,binz),binx,biny,binz) = i
     end do
  -!  call cpu_time(t2)
   
     nbs = 0
     nc = 0
  
positive.code:     iter = iter + 1
  
      ! Simulation the fluid flow
      call step
      time = time + dt
  
  
negative.code:     nbinz = floor(BOX_SIZE/h) - 1
  
   ! Need to find good way to estimate this so array not needlessly large
    max_part_guess = 100*n/(nbinx*nbiny*nbinz)
    allocate(part_count(nbinx,nbiny,nbinz))
    allocate(binpart(max_part_guess,nbinx,nbiny,nbinz)); binpart = 0
    allocate(nbs(n,max_part_guess)); nbs = 0
  