We have implemented a fully implicit time integration scheme into a version of Los Alamos National Laboratories' popular ocean general circulation model POP (Parallel Ocean Program). POP has traditionally used a semi-implicit time integration scheme where key linear terms in the equations are treated implicitly and the remaining terms are treated explicitly. In our fully implicit scheme, we use Crank-Nicholson on all terms of the equations and Jacobian free Newton-Krylov iterations to solve the resulting system. We substantially reduce the number of iterations needed by POP's semi-implicit integration scheme as a preconditioner. In idealized geometries, we can achieve timesteps 800 times larger than allowed by the semi-implicit scheme. For simulations of the North Atlantic and the global ocean, we can obtain time steps up to 25 times larger than allowed by the semi-implicit scheme. These realistic simulations use actual topography, continental geometry, initial Levitus temperature and salinity, observed annual averaged wind stress and restoring boundary conditions at the surface for temperature and salinity.