In this method we present a fractional step discretization of the time-dependent incompressible Navier-Stokes equations. The method is based on a projection formulation in which we first solve diffusion-convection equations to predict intermediate velocities which are then projected onto the space of divergence-free vector fields. Our treatment of the diffusion-convection step uses a specialized second-order upwind method for differencing the nonlinear convective terms that provides a robust treatment of these terms at high Reynolds number. In contrast to conventional projection-type discretizations that impose a discrete form of the divergence-free constraint, we only approximately impose the constraint; i.e., the velocity field we compute is not exactly divergence-free. The approximate projection is computed using a conventional discretization of the Laplacian and the resulting linear system is solved using conventional multigrid methods. Numerical examples are presented to validate the second-order convergence of the method for Euler, finite Reynolds number, and Stokes flow. A second example illustrating the behavior of the algorithm on an unstable shear layer is also presented.