This paper describes a numerical method with two time levels for solving partial differential equations that represent layered models of ocean circulation. The method is designed to be used with a barotropic-baroclinic splitting that separates the fast and slow motions into subsystems that are solved by different techniques. With this method, some of the dependent variables are predicted and then corrected. After the initial prediction steps, all steps involve centered differencing and averaging about the midpoint of the time interval in question. The Coriolis and pressure terms are evaluated at the same time levels to avoid a first-order truncation error in the geostrophic balance between those terms. Compared to the three-level leapfrog method that is widely used in geophysical fluid dynamics, the present method does not admit a computational mode, and the maximum permissible time step is at least twice as large. In addition, with the two-level method it is possible to use a nearly nonoscillatory advection algorithm to solve the equations for mass and momentum. In a simple test problem, the present method gives less phase error than the leapfrog method and two other methods, and it does not give any amplitude error. The differences are especially large when comparing the present method to the leapfrog method. The two-level method also gives good results in a nonlinear test problem involving vanishing layers at the top and bottom of the fluid domain and an interface moving along sloping bottom topography.