  We show a numerical scheme to solve the moment equations of the radiative
transfer, i.e., M1 model which follows the evolution of the energy density, $ E
$, and the energy flux, $ \mbox{\boldmath$F$} $. In our scheme we reconstruct
the intensity from $ E $ and $ \mbox{\boldmath$F$} $ so that it is consistent
with the closure relation, relation, $ \chi = (3 + 4 f ^2)/(5 + 2 \sqrt{4 - 3 f
^2}) $.
  Here the symbols, $ \chi $, $ f = |\mbox{\boldmath$F$}|/(cE) $, and $ c $,
denote the Eddington factor, the reduced flux, and the speed of light,
respectively. We evaluate the numerical flux across the cell surface from the
kinetically reconstructed intensity. It is an explicit function of $ E $ and $
\mbox{\boldmath$F$} $ in the neighboring cells across the surface considered.
We include absorption and reemission within a numerical cell in the evaluation
of the numerical flux. The numerical flux approaches to the diffusion
approximation when the numerical cell itself is optically thick. Our numerical
flux gives a stable solution even when some regions computed are very optically
thick. We show the advantages of the numerical flux with examples. They include
flash of beamed photons and irradiated protoplanetary disks.
