We propose a new turbulence closure model based on the budget equations for the key second moments: turbulent kinetic energy (TKE), turbulent potential energy (TPE) and vertical turbulent fluxes of momentum and buoyancy (proportional to potential temperature). Besides the concept of the turbulent total energy (TTE = TKE + TPE), we take into account the non-gradient correction to the traditional buoyancy flux formulation. The proposed model permits the existence of turbulence at any gradient Richardson number, Ri. For the stationary, homogeneous regime the turbulence closure model yields universal dependencies of the flux Richardson number, turbulent Prandtl number, anisotropy of turbulence, and normalized vertical fluxes of momentum and heat on the gradient Richardson number, Ri. We also take into account an additional vertical flux of momentum and additional productions of turbulent kinetic energy, turbulent potential energy and turbulent flux of potential temperature due to large-scale internal gravity waves (IGW). Accounting for the internal gravity waves, the Ri-dependencies of the flux Richardson number, turbulent Prandtl number, anisotropy of turbulence, vertical fluxes of momentum and heat lose their universality. In particular, with increasing wave energy, the maximal value of the flux Richardson number (attained at very large Ri) decreases. In contrast to the mean wind shear which generates only the horizontal TKE, IGW generate both horizontal and vertical TKE, and thus lead to a more isotropic turbulence at very large Ri. IGW also increase the share of TPE in the turbulent total energy. A well-known effect of IGW is their direct contribution to the vertical transport of momentum. Depending on the direction (downward or upward), it either strengthens of weakens the total vertical flux of momentum. Predictions from the proposed model are consistent with available data from atmospheric and laboratory experiments, direct numerical simulations (DNS) and large-eddy simulations (LES).