anchor.diff:  contains
   
     subroutine init(this, L, a, has_walls)
  +    use omp_lib
       class(cell_system_t), intent(out) :: this
       integer, intent(in) :: L(3)
       double precision, intent(in) :: a
  
positive.code:      double precision :: a
       double precision :: origin(3)
       integer, allocatable :: cell_count(:)
       integer, allocatable :: cell_start(:)
       integer :: M(3)
       logical :: has_walls
  
negative.code:     end do
      this%N = 2**this%M(1)*2**this%M(2)*2**this%M(3)
      allocate(this%cell_count(this%N))
      allocate(this%cell_count_tmp(this%N, omp_get_max_threads()))
      allocate(this%cell_start(this%N))
  
    end subroutine init
  
    subroutine count_particles(this, position)
      use omp_lib
      class(cell_system_t), intent(inout) :: this
      double precision, intent(in) :: position(:, :)
  
      integer :: i, idx, N_particles, t_id, t_max
      integer :: p(3)
      integer :: L(3)
      logical :: nowalls
  
      N_particles = size(position, 2)
      t_max = omp_get_max_threads()
  
      L = this% L
      nowalls = .not. this% has_walls
  
      this%cell_count_tmp = 0
  
      !$omp parallel private(i, p, idx, t_id)
      t_id = omp_get_thread_num() + 1
      !$omp do
      do i=1, N_particles
         p = floor( (position(:, i) / this% a ) - this%origin )
         if ( p(1) == L(1) ) p(1) = 0
  