anchor.diff:          do stz = -1, 1
             cbinz = binz + stz
             if(cbinz.lt.1 .or. cbinz.gt.nbinz) cycle
  -          do j = 1, max_part_guess
  -            if (binpart(j,cbinx,cbiny,cbinz) .eq. 0) exit
  +          do j = 1, part_count(cbinx, cbiny, cbinz)
               p = binpart(j,cbinx,cbiny,cbinz)
               distx = px(i) - px(p)
               disty = py(i) - py(p)
               distz = pz(i) - pz(p)
  -            tdist = sqrt(distx*distx + disty*disty + distz*distz)
  -            if(tdist .lt. h .and. p.ne.i) then
  +            tdist2 = distx*distx + disty*disty + distz*distz
  +            if(tdist2 .lt. h2 .and. p.ne.i) then
                 q = q + 1
                 nbs(i,q) = binpart(j,cbinx,cbiny,cbinz)
               end if
  
positive.code:   integer :: binx, biny, binz
    integer :: cbinx, cbiny, cbinz
    integer :: stx, sty, stz
    real(WP) :: distx, disty, distz, tdist
  
    part_count = 0
    binpart = 0
  
negative.code:   call cpu_time(t4)
    call compute_visc
    call cpu_time(t5)
    if(USE_SURFACE_TENSION) call compute_SF
    call ext_forces
    call cpu_time(t6)
    call posvel_update
  