A new software framework using a well-established high-order spectral element discretization is presented for solving the compressible Navier–Stokes equations for purposes of research in atmospheric dynamics in bounded and unbounded limited-area domains, with a view toward capturing spatiotemporal intermittency that may be particularly challenging to attain using low order schemes. A review of the discretization is provided, emphasizing properties such as the matrix product formalism and other design considerations that will facilitate its effective use on emerging exascale platforms, and a new geometry-independent, element boundary exchange method is described to maintain continuity. A variety of test problems are presented that demonstrate accuracy of the implementation primarily in wave-dominated or transitional flow regimes; conservation properties are also demonstrated. A strong scaling CPU study in a three-dimensional domain without using threading shows an average parallel efficiency of ≳ 99% up to 2×104 MPI tasks that is not affected negatively by expansion polynomial order. On-node performance is also examined and reveals that, while the primary numerical operations achieve their theoretical arithmetic intensity, the application performance is largely limited by available memory bandwidth.