We examine a nonlinear elimination method for the free-surface ocean equations based on barotropic-baroclinic decomposition. The two dimensional scalar continuity equation is treated implicitly with a preconditioned Jacobian-free Newton-Krylov method (JFNK) and the remaining three dimensional equations are subcycled explicitly within the JFNK residual evaluation. In this approach, the footprint of the underlying Krylov vector is greatly reduced over that required by fully coupled implicit methods. The method is second order accurate and scales algorithmically with timesteps much larger than explicit methods. Moreover, the heterogeneous nature of the algorithm lends itself readily to emerging architectures.