In this work, we consider the primitive equations of an ocean circulation model for the southern pacific, which consists of the time-dependent Navier-Stokes equations in the β-plane coupled with the temperature transport equation. Specifically, the full three-dimensional equations are adopted and formulated as a monolithic system of nonstationary, nonlinear, coupled partial differential equations. The El Nino phenomenon is simulated by the action of given wind stresses on the ocean surface. We present an approximation scheme with Crank-Nicolson finite differences in time, and in space we take inf-sup stable Galerkin finite elements for the Navier-Stokes part and bilinear elements for the temperature. We solve the resulting, nonlinear monolithic discrete system by Newton's method. Numerical experiments with realistic geometry and material data are conducted, which show the...