We present starburst models for far-infrared/sub-millimeter/millimeter (FIR/sub-mm/mm) line emission of molecular and atomic gas in an evolving starburst region, which is treated as an ensemble of non-interacting hot bubbles which drive spherical shells of swept-up gas into a surrounding uniform gas medium. These bubbles and shells are driven by stellar winds and supernovae within massive star clusters formed during an instantaneous starburst. The underlying stellar radiation from the evolving clusters affects the properties and structure of photodissociation regions (PDRs) in the shells, and hence the spectral energy distributions (SEDs) of the molecular and atomic line emission from these swept-up shells and the associated parent giant molecular clouds (GMCs) contains a signature of the stage of evolution of the starburst. The physical and chemical properties of the shells and their structure are computed using a a simple well known similarity solution for the shell expansion, a stellar population synthesis code, and a time-dependent PDR chemistry model. The SEDs for several molecular and atomic lines ($^{12}$CO and its isotope $^{13}$CO, HCN, HCO$^+$, C, O, and C$^+$) are computed using a non-local thermodynamic equilibrium (non-LTE) line radiative transfer model. By comparing our models with the available observed data of nearby infrared bright galaxies, especially M 82, we constrain the models and in the case of M 82, we provide estimates for the ages (5 - 6 Myr, 10 Myr) of recent starburst activity. We also derive a total H$_2$ gas mass of $\sim$ 2 - 3.4 $\times$ 10$^8$ M$_{\odot}$ for the observed regions of the central 1 kpc starburst disk of M 82.