To simulate and study ocean phenomena involving complex dynamics over a wider range of scales, from regional to small scales (e.g., thousands of kilometers to meters), resolving submesoscale features, nonlinear internal waves, subduction, and overturning where they occur, non-hydrostatic (NHS) ocean models are needed, at least locally. The main computational burden for NHS models arises from solving a globally coupled 3D elliptic PDE for the NHS pressure. To address this challenge, we start with a high-order hybridizable discontinuous Galerkin (HDG) (Nguyen et al. 2009) finite element NHS ocean solver (Ueckermann and Lermusiaux 2016) that is well suited for multidynamics systems. We present a new adaptive algorithm to decompose a domain into NHS and HS dynamics subdomains and solve their corresponding equations, thereby reducing the cost associated with the NHS pressure solution step. The NHS/HS subdomains are adapted based on new numerical NHS estimators, such that NHS dynamics is used only where needed. We compare and explore choices of boundary conditions imposed on the internal boundaries between subdomains of different dynamics. We evaluate and analyze the computational costs and accuracy of the adaptive NHS-HS solver using three idealized NHS dynamics test cases, (i) idealized internal waves (Vitousek and Fringer 2011), (ii)  tidally-forced oscillatory flow over seamounts and (iii)  bottom gravity currents. We then complete more realistic NHS-HS simulations of Rayleigh-Taylor instability-driven subduction events by nesting with our MSEAS realistic and operational data-assimilative HS ocean modeling system. Finally, we discuss DG-FEM-based numerical techniques to stabilize and accelerate the high-order ocean solvers by leveraging the high aspect ratio characteristic of ocean domains.