The use of σ-coordinate ocean models has historically been considered a disadvantage for large-scale climate studies. The main reason resides in the non-alignment of the vertical coordinate isosurfaces with either geopotential surfaces or isopycnals making it harder to accurately compute the horizontal pressure gradient, advection, and isoneutral tracer diffusion. Moreover, this class of model requires a vertical mixing parameterization robust to large changes in the vertical resolution between shallow and deep areas. In this paper we show that, with some adjustments of the tracer advection, the surface boundary layer parameterization and the vertical grid, a σ-coordinate model can achieve an accurate representation of the oceanic interior and mixed-layer dynamics. To do so, a new way of handling the temporal discretization of the rotated biharmonic operator is used to achieve tracer variance dissipation in an adiabatic and computationally efficient way. Furthermore, a redesign of the K-Profile surface layer Parameterization (KPP) to improve the regularity of the solution and the overall numerical efficiency of the scheme is introduced. To validate the new algorithmic developments, we perform a set of coarse-resolution realistic basin-scale Pacific simulations. Besides improving the conservation of water mass properties, the use of an isoneutral tracer hyperdiffusion is shown to have a negative feedback on the circulation error growth rate, thus significantly reducing the sensitivity of the model solution to the degree of topographic smoothing. The overall validation of our simulations, focusing on the key characteristics of the circulation in the Pacific Ocean, provides some evidence of the efficacy of a terrain-following coordinate for large scale applications.