This paper presents a numerical method for predicting the atmospheric boundary layer under stable, neutral, or unstable thermal stratifications. The flow field is described by the Reynolds' averaged Navier-Stokes equations complemented by the k-is an element of turbulence model. Density variations are introduced into the momentum equation using the Boussinesq approximation, and appropriate buoyancy terms are included in the k and is an element of equations. An original expression for the closure coefficient related to the buoyancy production term is proposed in order to improve the accuracy of the simulations. The resulting mathematical model has been implemented in FLUENT. The results presented in this paper include comparisons with respect to the Monin-Obukhov similarity theory, measurements, and earlier numerical solutions based on k-is an element of turbulence models available in the literature. It is shown that the proposed version of the k-is an element of model significantly improves the accuracy of the simulations for the stable atmospheric boundary layer. In neutral and unstable thermal stratifications, it is shown that the version of the k-is an element of models available in the literature also produce accurate simulations.