anchor.diff:               i_r = floor((r-r_min)/dr) + 1
                theta = acos( dot_product(unit_r,d)/r )
                i_th = floor(theta/dtheta) + 1
  -             !$omp atomic
  -             this%c(s2, i_th, i_r) = this%c(s2, i_th, i_r) + 1
  +             this_c(s2, i_th, i_r) = this_c(s2, i_th, i_r) + 1
                one_r = d/r
                ! compute v_r, v_theta
                out_of_plane = cross(unit_r, one_r)
  
positive.code:     double precision, dimension(3) :: solvent_v
      integer :: idx, s2, i_r, i_th
  
      integer :: cell_idx, start, n
  
      r_min = this%r_min
      r_min_sq = r_min**2
  
negative.code:              solvent_v = solvent%vel(:,idx) - v
               v_r = dot_product(solvent_v, one_r)
               v_th = dot_product(solvent_v, one_theta)
               this_v(:, s2, i_th, i_r) = this_v(:, s2, i_th, i_r) + [v_r, v_th]
               this_count(s2, i_th, i_r) = this_count(s2, i_th, i_r) + 1
            end if
         end do
      end do
  
      this%v = this%v + this_v
      this%c = this%c + this_c
      this%count = this%count + this_count
  
      deallocate(this_c)
      deallocate(this_v)
      deallocate(this_count)
  
      call this%time_polar_update%tac()
  
    end subroutine update
  