In climate modeling, the three-dimensional incompressible Navier-Stokes equations formulated in a rotational frame of reference are a very good representation of the dynamics of ocean circulation. As such, it is paramount to further advance numerical schemes to approximate their solution. In this presentation we describe how the First Order System Least Squares Finite Element Method (FOSLS) can be used to solve the model and the nature of the discrete linear system arising from its application, which are optimally solved by application of Multigrid methods.