A new scalable and efficient implementation of the mesoscopic distinct element method for massively parallel numerical simulations of carbon nanotube systems is introduced. Carbon nanotubes are represented as chains of rigid bodies, linked by elastic bonds and dispersive van der Waals (vdW) forces. The enhanced vector model formalism of the elastic bond between rigid bodies, developed recently, is employed here to capture the elastic deformation of nanotubes. Dispersive interactions between the neighboring nanotubes are described with the coarse-grained vdW potential. Time integration is performed using a velocity Verlet integration scheme with tunable damping in order to describe the energy dissipation to the implicit degrees of freedom. Due to the scalable message passing interface (MPI) parallelization, enabled by rigid particle dynamics module (PE) of the waLBerla multiphysics framework, our method is capable of modeling large assemblies of carbon nanotubes. This advance enables us to move closer to the length and time scales required to extract representative mechanics of carbon nanotube materials. The promising scalability of the new implementation is probed in two examples of self-assembly of ultrathin carbon nanotube films and carbon nanotube buckypapers, where formation of hierarchical networks of carbon nanotube bundles, storing both elastic and vdW adhesion energies is observed. The relaxation of one cubic micrometer of buckypaper illustrates the code scalability.