As part of a research effort by the U.S. Environmental Protection Agency to study the dispersion of pollutants in aquatic systems, a numerical model has been developed that is capable of realistically describing the hydrodynamics in lakes, embayments, nearshore marine coastal areas, and riverine and thermal outfall plumes. The model is time-dependent, three-dimensional, and variable density. Both rigid-lid and free-surface flows can be determined. The main assumptions used in the development of the model include hydrostatic pressure variation, Boussinesq approximation, and eddy coefficients to account for turbulence. A new solution procedure, which is a modification of the simplified marker and cell method, is used. The procedure permits selected terms in the equations to be treated implicitly in time. The report provides the documentation for the computer program for the numerical model.