anchor.diff:    real(WP) :: y_height, distx, disty, distz, tdist
   
     deallocate(nbs)
  -  y_height = maxval(py(:))
  -  nbinx = ceiling(BOX_SIZE/h) - 1
  -  nbiny = ceiling(y_height/h) - 1
  -!  nbinz = ceiling(BOX_SIZE/h) - 1
  -nbinz = 1
  -
  +
  +  nbinx = floor(BOX_SIZE/h) - 1
  +  nbiny = floor(BOX_SIZE/h) - 1
  +  nbinz = floor(BOX_SIZE/h) - 1
   
     max_part_guess = 10000*n/(nbinx*nbiny*nbinz)
     allocate(part_count(nbinx,nbiny,nbinz))
  
positive.code:     do j = 1, l
        nb = nbs(i,j)
        rhosum = rhosum + densitykernel(i,nb)
  
      end do
    rho(i) = rho(i) + rhosum
  
    end do
  
    return
  end subroutine compute_density
  
negative.code:       dz = pz(i) - pz(nb)
  
        k = pressurekernel(i,nb)
  !      print*, i, nb, k, dx, k*dx, dy, k*dy, dz, k*dz
        sx = sx + k*dx
        sy = sy + k*dy
        sz = sz + k*dz
  