A numerical approach for the modeling of thermal effects in turbulent mixing of large flows with uniform eddy viscosities and eddy diffusivities is presented. The distinctive features of this approach are the simultaneous compliance of mass and momentum equations, and accurate predictions of the surface waves and pressure distribution in flow fields with arbitrary bottom geometry. The numerical approach employs centered three‐point finite difference expressions for the first and second spatial derivatives with unequal mesh size and a forward time integration. Numerical stability analysis of the governing equations is utilized to select optimum space and time increments for numerical integration.