Summary: Numerical methods for the primitive equations (PEs) of oceanic flow are presented in this paper. First, a two-dimensional Poisson equation with a suitable boundary condition is derived to solve the surface pressure. Consequently, we derive a new formulation of the PEs in which the surface pressure Poisson equation replaces the nonlocal incompressibility constraint, which is known to be inconvenient to implement. Based on this new formulation, backward Euler and Crank-Nicolson schemes are presented. The marker and cell scheme, which gives values of physical variables on staggered mesh grid points, are chosen as spatial discretization. The convergence analysis of the fully discretized scheme is established in detail. The accuracy check for the scheme is also shown.