  In circumstellar disks around T Tauri stars, visible and near-infrared
stellar irradiation is intercepted by dust at the disk's optical surface and
reprocessed into thermal infrared; this subsequently undergoes radiative
diffusion through the optically thick bulk of the disk. The gas component --
overwhelmingly dominant by mass, but contributing little to the opacity -- is
heated primarily by gas-grain collisions. In hydrodynamical simulations,
however, typical models for this heating process (local isothermality,
$\beta$-cooling, two-temperature radiation hydrodynamics) incorporate
simplifying assumptions that limit their ranges of validity. To build on these
methods, we develop a ``three-temperature" numerical scheme, which
self-consistently models energy exchange between gas, dust, and radiation, as a
part of the PLUTO radiation-hydrodynamics code. With a range of test problems
in 0D, 1D, 2D, and 3D, we demonstrate the efficacy of our method, and make the
case for its applicability to a wide range of problems in disk physics,
including hydrodynamic instabilities and disk-planet interaction.
