A three-dimensional, non-hydrostatic pressure, numerical model for free surface flows is presented. By decomposing the pressure term into hydrostatic and non-hydrostatic parts, the numerical model uses an integrated time step with two fractional steps. In the first fractional step, the momentum equations are solved without the hydrostatic pressure term using Newton's method in conjunction with the generalised minimal residual (gmres) method. This combined method does not require the determination of a Jacobian matrix explicitly but simply the product of the Jacobian and a vector, thereby reducing the amount of storage required and significantly decreasing the overall * Applied Mathematics, The University of Adelaide, South Australia, Australia. C624 computational time required. By using Newton's method, the numerical model can handle implicitly almost all variables, unlike many other numerical models. Hence numerical stability is achieved effectively. In the second fractional step, the pressure-Poisson equation is solved iteratively with a preconditioned linear gmres method. It is shown that preconditioning reduces the processing time dramatically. After the new pressure field is obtained the intermediate velocities, which are calculated from the previous fractional step, are updated and then these updated velocities preserve the local mass conservation. The newly developed model is verified against analytical solutions, with good agreement.