  Dust plays a key role in the evolution of the ISM and its correct modelling
in numerical simulations is therefore fundamental. We present a new and
self-consistent model that treats grain thermal coupling with the gas,
radiation balance, and surface chemistry for molecular hydrogen. This method
can be applied to any dust distribution with an arbitrary number of grain types
without affecting the overall computational cost. In this paper we describe in
detail the physics and the algorithm behind our approach, and in order to test
the methodology, we present some examples of astrophysical interest, namely (i)
a one-zone collapse with complete gas chemistry and thermochemical processes,
(ii) a 3D model of a low-metallicity collapse of a minihalo starting from
cosmological initial conditions, and (iii) a turbulent molecular cloud with
H-C-O chemistry (277 reactions), together with self-consistent cooling and
heating solved on the fly. Although these examples employ the publicly
available code KROME, our approach can be easily integrated into any
computational framework.
