anchor.diff:  
       rho(i) = rho(i) + rhosum
     end do
  +  !$OMP END DO
   
     return
   end subroutine compute_density
  
positive.code:     part_count(binx, biny, binz) = part_count(binx, biny, binz) + 1
      binpart(part_count(binx,biny,binz),binx,biny,binz) = i
    end do
  
    nbs = 0
    nc = 0
  
    !$OMP PARALLEL DO default(none) shared(nbs,nc,nbinx,nbiny,nbinz,px,py,pz,h2,part_count,binpart) &
    !$OMP private(q,i,binx,biny,binz,stx,sty,stz,cbinx,cbiny,cbinz,nb,distx,disty,distz,tdist2) &
    !$OMP NUM_THREADS(procs)
    do i = 1, n
      binx = 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
  
negative.code:     fz(i) = sz
  
    end do
    !$OMP END DO
  
    return
  end subroutine compute_forces
  