Accurate numerical modeling of nonhydrostatic dynamics is essential to properly describe the physical phenomena in stratified ocean flows, including internal waves, Kelvin-Helmholtz instabilities and flows encountering steep topography. However, resolving these dynamics can often be computationally expensive due to the complex mathematical formulations involved. Nonhydrostatic models consist of the full incompressible Navier-Stokes equations with the Boussinesq approximation, convection-diffusion equations for salinity and temperature, and an equation of state for the density. The high-order finite element method for solving these equations has well-established accuracy and efficiency advantages over traditional computational techniques [1,2]. It is well known that the spatial and temporal discretization of convection terms may lead to spurious oscillations in the numerical solution when the model resolution is not fine enough [3]. This problematic behavior can produce unphysical values such as negative density. In some cases, they generate unstable solutions. Moreover, these oscillations grow as the order of the basis functions increases. We introduce a high-order finite element method in combination with a second order BDF semi-Lagragian scheme with a stabilization technique in a nonhydrostatic ocean model. Mesh adaptation strategies are used to identify elements that can exhibit spurious oscillations, applying to them the so-called quasi-monotone scheme. The numerical results show that for convection-dominated convection-diffusion problems, the spurious oscillations are eliminated. We also carry out several numerical experiments in two-dimensions bi-phase ocean flows where the nonhydrostatic effects are important. The results illustrate that the method is an effective numerical tool to stabilize the numerical solution and to increase the accuracy.