On the decadal to centennial timescale, changes in climate are controlled strongly by changes in ocean circulation. However, because of limitations inherent to the time integration schemes used in present-day ocean models, state-of-the-art climate change simulations resolve the oceans only very coarsely. With an aim to enable long-term simulations of ocean circulation at the high resolutions required for a better representation of global ocean dynamics, we have implemented fully-implicit time integration schemes in a version of the popular ocean general circulation model POP (Parallel Ocean Program), employing Jacobian free Newton-Krylov techniques. Here, we describe the numerical principles underlying iPOP in some detail and present a few computational results. While there are many advantages to this approach, including a consistent and uniform treatment of the terms in the governing equations, the primary advantage lies in the ability to take time steps that are of relevance to the physical phenomenon that is being studied. The time step is not limited (for stability reasons) by the fastest modes of the system.