In this paper, we present a new linear system solver for use in a fully-implicit ocean model. The new solver allows to perform bifurcation analysis of relatively high-resolution primitive-equation ocean-climate models. It is based on a block-ILU approach and takes special advantage of the mathematical structure of the governing equations. In implicit models Jacobian matrices have to be constructed. Analytical construction is hard for complicated but more realistic representations of mixing. This is overcome by evaluating the Jacobian in part numerically. The performance of the new implicit ocean model is demonstrated using (i) a high-resolution model of the wind-forced double-gyre flow problem in a (relatively small) midlatitude spherical b...