In this work, a spectral model is derived to investigate numerically unstably stratified homogeneous turbulence (USHT) at large Reynolds numbers. The modeling relies on an earlier work for passive scalar dynamics [Briard et al., J. Fluid Mech. 799, 159 (2016)] and can handle both shear and mean scalar gradients. The extension of this model to the case of active scalar dynamics is the main theoretical contribution of this paper. This spectral modeling is then applied at large Reynolds numbers to analyze the scaling of the kinetic energy, scalar variance, and scalar flux spectra and to study as well the temporal evolution of the mixing parameter, the Froude number, and some anisotropy indicators in USHT. A theoretical prediction for the exponential growth rate of the kinetic energy, associated with our model equations, is derived and assessed numerically. Throughout the validation part, results are compared with an analogous approach, restricted to axisymmetric turbulence, which is more accurate in term of anisotropy description, but also much more costly in terms of computational resources [Burlot et al., J. Fluid Mech. 765, 17 (2015)]. It is notably shown that our model can qualitatively recover all the features of the USHT dynamics, with good quantitative agreement on some specific aspects. In addition, some remarks are proposed to point out the similarities and differences between the physics of USHT, shear flows, and passive scalar dynamics with a mean gradient, the two latter configurations having been addressed previously with the same closure. Moreover, it is shown that the anisotropic part of the pressure spectrum in USHT scales in k^{-11/3} in the inertial range, similarly to the one in shear flows. Finally, at large Schmidt numbers, a different spectral range is found for the scalar flux: It first scales in k^{-3} around the Kolmogorov scale and then further in k^{-1} in the viscous-convective range.