In this paper we describe the numerical code for a three-dimensional ground- water flow model. The model is developed for the case of variably saturated porous media, applicable to both the unsaturated (soil) zone and the saturated (groundwater) zone. The governing equation is nonlinear, and is linearized using either Picard or Newton iteration. The large sparse systems of linear equations generated by the finite element discretization are solved using efficient preconditioned conjugate gradient-like methods. Tetrahedral elements and linear basis functions are used for the discretization in space, and a weighted finite difference formula is used for the discretization in time. The code handles: temporally and spatially variable boundary conditions, including seepage faces and evaporation/precipitation inputs; heterogeneous material properties and hydraulic characteristics, including saturated conductivities, porosities, and storage coefficients; and various expressions to describe the moisture content-pressure head and relative conductivity-pressure head relationships. The model solves for nodal pressure heads, and uses these values to compute the water saturations and velocities over the flow domain. The water saturation and velocity values can be used ad input in the LEA3D and NONLEA3D transport codes, which are describede in a companion paper.