Extragalactic nebular emission has long been a workhorse probe of the processes driving galaxy evolution, but the richness of JWST spectroscopy has shifted the bottleneck from data acquisition to physical interpretation and modelling. In this context, we present a major update to the Monte Carlo radiative transfer code COLT to facilitate self-consistent modelling of nebular line and continuum emission from simulated galaxies. We introduce a new thermal equilibrium solver that iteratively couples to the existing ionization solver and radiation field to compute effective gas temperatures by accurately balancing photoionization heating, radiative and dielectronic recombination, collisional ionization, charge exchange, metal and primordial line cooling, free-free emission, and Compton scattering. To prevent over-cooling where non-equilibrium hydrodynamics dominate, we introduce a Courant-limited cooling prescription tied to each cell's sound-crossing time, preserving temperatures in the diffuse halo while allowing physically motivated cooling in the interstellar medium (ISM). Applied to an isolated local galaxy simulation, the equilibrium solver reshapes the ISM phase space by reducing spuriously excessive lukewarm (T=10³−10⁴K) gas and better resolving warm ionized and cold neutral phases, while leaving the CGM largely intact. We further implement a level population solver based on modern atomic data, enabling accurate cooling and emissivities for a large library of UV to infrared metal lines, together with newly implemented primordial nebular continuum emission from free-free, free-bound, and two-photon processes. Finally, by applying COLT to the high-redshift THESAN-ZOOM simulations, we reproduce observed emission-line ratios, establishing COLT as a robust framework for forward modelling nebular emission across cosmic time.