We present in this paper a pressure stabilized, finite element method for the numerical approximation of three-dimensional, non-hydrostatic mesoscale ocean flows. The model considered here incorporates surface wind stress, bottom friction and Coriolis acceleration, and it is applicable to irregular bottom topographies. An implicit unconditionally stable scheme is employed for the time advancement and an anisotropic stabilization technique is used for the spatial finite element discretization. The numerical results obtained on test cases demonstrate the robustness and accuracy of the method proposed here.