A finite element algorithm is described which implements the Galerkin approximation to the Navier-Stokes equations and incorporates five predominate features. Although none of these five features is unique to this algorithm, their orchestration, as described in this paper, results in an algorithm which is not only easy to implement but also stable, accurate, and robust, as well as computationally efficient. The zero stress natural boundary condition is implemented which permits calculation of the outflow velocity distribution. A nine-node, Lagrangian, isoparametric, quadrilateral elementis used to represent the velocity while the pressure uses a four-node, Lagrangian, superparametric element coincident with the velocity element. The easily implemented, computationally efficient frontal solution technique is used to assemble the element coefficient matrices, impose the boundary conditions, and solve the resulting linear system of equations. An implicit backward Euler time integration rule provides a very stable solution method for time dependent problems. A Picard scheme with a relatively large radius of convergence is used for iteration of the non-linear equations at each time step. Results are given from the calculation of two dimensional, steady-state and time-dependent convection dominated flows of viscous incompressible fluids.