anchor.diff:      cell_idx = 1
       count = 0
       te = 0
  -    !$omp parallel do private(start, n, local_v, local_kin, i, cell) &
  +    !$omp parallel do private(start, n, local_v, local_kin, i, cell, n_effective) &
       !$omp& reduction(+:count) reduction(+:te)
       do cell_idx = 1, cells% N
          if (cells% cell_count(cell_idx) <= 1) cycle
  
positive.code:     call particles%time_step%tic()
      !$omp parallel private(thread_id)
      thread_id = omp_get_thread_num() + 1
      !$omp do private(start, n, local_v, i, vec, omega)
      do cell_idx = 1, cells% N
         if (cells% cell_count(cell_idx) <= 1) cycle
  
         start = cells% cell_start(cell_idx)
         n = cells% cell_count(cell_idx)
  
         local_v = 0
         do i = start, start + n - 1
            local_v = local_v + particles% vel(:, i)
         end do
         local_v = local_v / n
  
         vec = rand_sphere(state(thread_id))
         omega = &
  
negative.code:        n = cells% cell_count(cell_idx)
         count = count + 1
  
         n_effective = 0
         local_v = 0
         local_kin = 0
         do i = start, start + n - 1
            if (particles%species(i) > 0) then
               local_v = local_v + particles% vel(:, i)
               n_effective = n_effective + 1
            end if
         end do
         local_v = local_v / n
         do i = start, start + n - 1
            if (particles%species(i) > 0) then
               local_kin = local_kin + sum((particles% vel(:, i)-local_v)**2)
            end if
         end do
         local_kin = local_kin / dble(3*(n_effective-1))
         te = te + local_kin
  
         if (do_tz) then
  