With increasing resolution in numerical ocean models, nonhydrostatic pressure effects have to be accounted for. In sigma-coordinate mode split ocean models, this pressure may be regarded as a pressure correction. An elliptic equation must be solved for the nonhydrostatic pressure, and the gradients are used to correct the provisional hydrostatic velocity components in each time step. The focus in the present work is on the surface boundary condition for the elliptic equation. In the literature, both Dirichlet and Neumann boundary conditions are suggested and applied. To investigate the sensitivity of the numerical results to the choice of boundary condition, three numerical experiments are performed. The first and second experiments are studies of the propagation and steepening of nonlinear internal waves. The first study is on tank scale and the second experiment is on ocean scale. In the tank-scale experiment, the density and the flow fields are very robust to the choice of boundary condition. In the ocean-scale experiment, the waves produced with a Dirichlet boundary condition become more damped than the waves produced with a Neumann boundary condition. The third study involves a surface buoyant jet. It is shown that well-known characteristics of the plume front are reproduced with a Neumann boundary condition, but the rotating turbulent core of this front is lost with a Dirichlet condition. It is accordingly argued that the appropriate surface boundary condition in mode split nonhydrostatic ocean models is the Neumann condition.