A mayor problem that arises in the computation of stellar atmosphere models is the self consistent determination of the temperature distribution via the constraint of energy conservation. The energy balance includes the gains due to the absorption of radiation: int χ(ν) J(ν) dν, and the losses due to emission: int χ(ν) S(ν) dν . It is well known that, within each one of the two above integrals, the part corresponding to spectral ranges whose opacity χ(ν) is huge can overcome by many orders of magnitude the part that corresponds to the remaining frequencies. On the other hand, at those frequencies where χ(ν) is very large, the mean intensity J(ν) of the radiation field shall be equal, up to many significant digits, to the source function S(ν) and consequently to the Planck function B(ν,T). Then their net share to the energy balance shall be null, albeit separately their contributions to the gain and loss integrals are the most important numerically. Thus, the spectral range whose physical contribution to the overall balance is null will dominate numerically both sides of the energy balance equation, and consequently the errors on the determination of J(ν) and S(ν) at these frequencies will falsify the balance. It is possible to circumvent the numerical problem brought about by the foregoing circumstances by solving the radiative transfer equation for the variable I(n,ν) - S(ν), instead of the customary intensity I(n,ν). We present here a novel iterative algorithm, based on iteration factors already employed by us with success, which makes it possible a fast correction of the temperature by computing directly the difference between the radiative losses and gains at each step of the iterations.