Numerical simulations of stratified flows at high Reynolds number have several important applications in engineering and physical oceanography and help us better understand the dynamics of turbulence in geophysical flows. Applications like these motivate tackling the significant computational and numerical challenges that high-Reynolds-number simulations can pose. This work presents the components of a high-order accurate flow solver designed to simulate stratified flows. The Navier-Stokes solver utilizes a Fourier pseudo-spectral method in the horizontal direction and a modal spectral element discretization in the vertical. We adopt an implicit-explicit time discretization scheme. At each time step, the scheme involves solving several one-dimensional Helmholtz problems. Static condensation and modal boundary-adapted basis functions result in an inexpensive algorithm based on solving a large number of small tridiagonal systems. Particular care is taken to ensure that the model is suitable to simulate stratified turbulent wakes at a body-based Reynolds number of the order of 10^6. The numerical model is well-suited for other stratified, highly turbulent flows developing in long, high aspect-ratio domains. A description of the key features of the proposed numerical model is presented, as well as the approach followed to test the correctness of the Navier-Stokes solver. The series of benchmark problems concludes with a turbulent stratified wake generated by a sphere in linear stratification.