This work focuses on developing efficient and robust implementation methods for hybridizable discontinuous Galerkin (HDG) schemes for fluid and ocean dynamics. In the first part, we compare choices in weak formulations and their numerical consequences. We address details in making the leap from the mathematical formulation to the implementation, including the different spaces and mappings, discretization of the integral operators, boundary conditions, and assembly of the linear systems. We provide a flexible mapping procedure amenable to both quadrature-free and quadrature-based discretizations, and compare the accuracy of the two on different problem geometries. We verify the quadrature-free approach, demonstrating that optimal orders of convergence can be obtained, even on non-affine and curvilinear geometries. The second part of the work investigates the scalability of HDG schemes, identifying memory and time-to-solution bottlenecks. The form of the quadrature-free integral operators is exploited to develop a novel and efficient matrix-free approach to solving the global linear system that arises from HDG discretizations. Additional manipulations to improve numerical robustness are discussed. To mitigate the complexity of the implementation, we provide an automated and computationally efficient verification procedure for the HDG methodologies discussed, using a hierarchical approach to provide diagnostic information and isolate problems. Finally, challenges related to the effective visualization of high-order, discontinuous HDG-FEM data for fluid and ocean applications are illustrated and strategies are provided to address them.