A procedure for the numerical solution of nonorthogonal boundary-fitted coordinate systems, i.e., a coordinate line coincides with the boundary, is presented. This method generates curvilinear coordinates as the solution of two elliptic partial differential equations with Dirichlet boundary conditions, one coordinate being specified to be constant on each of the boundaries, and a distribution of the other specified along the boundaries. No restrictions are placed on the irregularity of the boundaries, which are even allowed to be time dependent, such as might occur in problems involving the computation of the location of the free surface, flooding boundaries, and streambank erosion problems. In addition, fields containing multiple bodies or branches can be handled as easily as simple geometries. Regardless of the shape and number of bodies and regardless of the spacing of the curvilinear coordinate lines, all numerical computations, both to generate the coordinate system and to subsequently solve the system of partial differential equations of interest, e.g., the vertically integrated hydrodynamic equations, are done on a rectangular grid with square mesh. Since the boundary-fitted coordinate system has coordinate lines coincident with all boundaries, all boundary conditions may be expressed at grid points, and normal derivatives on the bodies may be represented using only finite differences between grid points on coordinate lines. No interpolation is needed even though the coordinate system is not orthogonal at the boundary. Several example sketches of coordinate systems for numerical modeling problems in the estuarine, riverine, and reservoir environments are presented. In addition, an actual computer generated plot of a boundary-fitted coordinate system for a region representative of the shape of Charleston Harbor is also presented. Several features of boundary-fitted coordinate systems are especially suited to the numerical modeling of hydraulic problems. Some of these are: (A.) The boundary geometry is completely arbitrary and is specified entirely by input. (B.) Complicated configurations, such as channels with branches and islands, can be treated as easily as simple configurations. (C.) Moving boundaries can be treated naturally without interpolation being required in the handling of boundary conditions or the expression of partial derivatives. (D.) Grid points can be concentrated in regions of rapid flow changes. (E.) General codes can be written that are applicable to different locations with different configurations since the code generated to approximate the solution of a given set of partial differential equations is independent of the physical geometry of the problem.