This paper deals with the development of efficient numerical solvers for the primitive equations of the ocean in turbulent regime. We derive the numerical approximation of a reduced model by the Smagorinsky turbulence model that includes stabilization of the pressure discretization by a penalty technique. We perform the numerical analysis of this discretization (stability, convergence, error estimates), obtaining error estimates of at most first order in natural norms, due to the penalty structure of the Smagorinsky eddy viscosity. We finally perform some numerical tests for the primitive and Navier–Stokes equations, that confirm the theoretical convergence expectations.