  We devise a physical model of formation and distribution of molecular gas
clouds in galaxies. We use the model to predict the intensities of rotational
transition lines of carbon monoxide (CO) and the molecular hydrogen (H$_{\rm
2}$) abundance. Using the outputs of Illustris-TNG cosmological simulations, we
populate molecular gas clouds of unresolved sizes in individual simulated
galaxies, where the effect of the interstellar radiation field with dust
attenuation is also taken into account. We then use the publicly available code
DESPOTIC to compute the CO line luminosities and H$_{\rm 2}$ densities without
assuming the CO-to-H$_{\rm 2}$ conversion factor ($\alpha_{\rm CO}$). Our
method allows us to study the spatial and kinematic structures traced by
CO(1-0) and higher transition lines. We compare the CO luminosities and H$_{\rm
2}$ masses with recent observations of galaxies at low and high redshifts. Our
model reproduces well the observed CO-luminosity function and the estimated
H$_{\rm 2}$ mass in the local Universe. About ten per cent of molecules in the
Universe reside in dwarf galaxies with stellar masses lower than $10^9~{\rm
M_\odot}$, but the galaxies are generally `CO-dark' and have typically high
$\alpha_{\rm CO}$. Our model predicts generally lower CO line luminosities than
observations at redshifts $z\gtrsim 1$--$2$. We argue that the difference can
be explained by the highly turbulent structure suggested for the high-redshift
star-forming galaxies.
