The surface pressure reformulation of the Bryan-Cox-Semtner global ocean model, which we developed recently, has been further improved by eliminating the rigid-lid approximation in favor of an implicit free-surface method for the barotropic mode. This retains the advantages of the surface pressure formulation and in addition (1) substantially improves the computational efficiency of the barotropic mode; (2) computes surface height directly, allowing comparison with and assimilation of alimetry data (with a rigid-lid case the surface height is accurately computed only in the limit of a steady state); (3) improves accuracy in the computation of long-wavelength Rossby waves and includes the effects of long-wavelength surface gravity waves; (4) provides robustness by selectively damping computational modes; (5) simplifies satisfying global energetic balances required for energetic consistency; and (6) alleviates ''checkerboarding'' in the surface height or pressure fields. A variable thickness top-layer model is developed to account for the free surface, and a generalized form of implicit time discretization of the barotropic equations is introduced. This discretization is analyzed in detail from the point of view of accuracy, stability, and damping properties, and certain cases are selected for their advantageous properties. The performance of the method is demonstrated on the Connection Machine, CM-200, and the results of numerical simulations are compared using the stream function, rigid-lid surface pressure, and free-surface formulations.