anchor.diff:              stop "Unknown value of flag_compute_mpi"
         end select
         !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ FIX BELOW
  -
         !--------------------------------
         ! reduce accelerations
         !--------------------------------
  
positive.code: 
        call integrate(v, a, dt/2)  ! Compute v(t+dt/2)
        call integrate(r, v, dt)    ! Compute r(t+dt)
  
  
        select case(flag_compute_mpi)
           case(0)
  
negative.code:                     MPI_COMM_WORLD, stat, err)
  
                 if (rank == i) then
                    r_i = r
                 else if (rank == nprocs - i)
                    r_np_i = r
                 end if
  
                 ! Broadcast i-th data
                 call mpi_bcast(r_i, 3*N, MPI_REAL_XP, i, MPI_COMM_WORLD, err)
                 ! Broadcast n-i-th data
                 call mpi_bcast(r_np_i, 3*N, MPI_REAL_XP, nprocs - i, MPI_COMM_WORLD, err)
  
                 ! If i == rank, compute both diagonals
                 if (rank == i) then
                    r_right = r_i
                    call compute_force_diag(m, r, r, a)
                    call compute_force_diag(m, r_right, r_right, a_right)
                    a_comm = a
  
                    ! compute on right side but communicate nothing
                 else if (rank < i) then
                    call compute_force(m, r_np_i, r_right, a_right)
                    a_comm = 0._xp
  
                    ! compute on left side and communicate interaction
                 else
                    call compute_force(m, r, r_i, a_comm)
                    a = a - a_comm
                 end if
  
  