The ocean general circulation models are essential numerical tools for studying marine systems, managing industrial offshore activities, and forecasting climate change. Density-driven flows play a key role in many ocean phenomena where nonhydrostatic effects are significant, such as internal wave formation, deep-water convection, and slope bottom currents (Galán del Sastre and Bermejo, 2016). In these flows, differences in water mass density can result in a sharp gradient in the temperature and salinity functions. Without sufficient spatial resolution, finite element ocean models may exhibit spurious oscillations due to the Gibbs phenomenon for such steep gradients, potentially leading to numerical instability and unphysical results (Frerichs and John, 2021). The standard finite element method typically requires a very fine mesh to achieve stable solutions, leading to prohibitively high computational costs. The high-order finite element method has well-established accuracy and efficiency advantages over traditional computational techniques (Fehn et al., 2017). We propose a semi-Lagrangian high-order finite element method for calculating the numerical solution of a nonhydrostatic ocean model. The multi-resolution weighted essentially non-oscillatory (WENO) is applied to get high order accuracy in smooth regions and to avoid spurious oscillations near sharp gradients (Zhu and Shu, 2018). We have tested the new numerical scheme on two-dimensional convection-diffusion problems and in flows where the nonhydrostatic effects are important. The results indicate that this scheme stabilizes the solution by significantly reducing spurious oscillations without requiring mesh refinement.