The present study introduces a new boundary layer parameterization for weather and forecasting models. It is implemented here as a boundary layer module in Weather Research and Forecasting (WRF) model. The main novelty in the new scheme is that it includes prognostic equations for the heat flux and temperature variance, being the first WRF boundary layer scheme with that feature. This is specially aimed at improving the representation of nocturnal stable boundary layer and of its turbulence regimes, weakly and very stable. The effort is supported by previous studies that found that the two regimes and the transitions between them are better represented by simplified numerical schemes that represent the interactions between the surface and the air adjacent to it when the heat flux and temperature variance are solved prognostically. The results show that the two regimes are adequately simulated by the new scheme. Such an evaluation is presented in terms of the relationship between the turbulence velocity scale and mean wind speed, of the dependence of the potential temperature gradient near the surface and the mean wind speed, and by the relationship between flux and gradient Richardson numbers. In the new scheme, the relationship between thermal structure and the mean and turbulent flows arises naturally from the heat flux prognostic equation, not being arbitrarily imposed by an empirical stability function.