This work proposes a thermophysical model for realistic surface layers on airless small bodies (RSTPM) for the use of interpreting their multiepoch thermal light curves (e.g., WISE/NEOWISE). RSTPM considers the real orbital cycle, rotation cycle, rough surface, temperature-dependent thermal parameters, as well as contributions of sunlight reflection to observations. It is thus able to produce a precise temperature distribution and thermal emission of airless small bodies regarding the variations on orbital timescales. Details of the physics, mathematics, and numerical algorithms of RSTPM are presented. When used to interpret multiepoch thermal light curves by WISE/NEOWISE, RSTPM can give constraints on the spin orientation and surface physical properties, such as the mean thermal inertia or the mean size of dust grains, the roughness fraction, and the albedo via a radiometric procedure. As an application example, we apply this model to the main-belt object (24) Themis, the largest object of the Themis family, which is believed to be the source region of many main-belt comets. We find multiepoch (2010, 2014-2018) observations of Themis by WISE/NEOWISE, yielding 18 thermal light curves. By fitting these data with RSTPM, the best-fit spin orientation of Themis is derived to be (λ = 137°, β = 59°) in ecliptic coordinates, and the mean radius of dust grains on the surface is estimated to be b = 140-114+500 (6∼ 640) μm, indicating that the surface thermal inertia varies from ∼3 Jm-2 s-0.5 K-1 to ∼60 Jm-2 s -0.5 K-1 due to seasonal temperature variation. A more detailed analysis found that the thermal light curves of Themis show a weak feature that depends on the rotation phase, which is indicative of heterogeneous thermal properties or imperfections of the light-curve inversion shape model.