anchor.diff:      nbinz = floor(BOX_SIZE/h) - 1
   
    ! Need to find good way to estimate this so array not needlessly large
  -  max_part_guess = n
  +  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
  
positive.code:     iter = iter + 1
  
      ! Simulation the fluid flow
      call step
      time = time + dt
  
  
negative.code:     end do
    end do
  
    return
  end subroutine neighbor_find
  
  