A two-part presentation is given in which a combination of domain splitting, coordinate transformation, and an implicit algorithm is suggested as a versatile three-dimensional program for accurate simulation of geometrically complex environmental flow problems. In part 1, a numerical method for solving the time-dependent laminar Navier-Stokes equations posed in general curvilinear coordinates is presented. The issues of equation form, consistency, and accuracy are considered as they relate to numerical simulation of incompressible flows on curvilinear grid systems. The method is applied to flows in square and polar cavities and evaluated in terms of corresponding benchmark physical and numerical experiments.