We present a thorough theoretical study of the title reaction using a wave packet based statistical model, which is ideally suited for reactions dominated by long-lived complexes. Both state-to-state integral cross sections and thermal rate constants were calculated for the three isotopic (H-2, D-2, and HD) reactions using a recent ab initio potential energy surface. Product (CH or CD) vibrational state distributions were found to decrease monotonically with the vibrational quantum number. The product rotational degree of freedom is highly excited up to the highest energetically available states, in good agreement with experimental results. For the C + HD reaction, the CD + H channel was found to be more populated than the CH + D channel, consistent with experimental observations. The thermal rate constants depend weakly on temperature and agree very well with recent experimental measurements at the room temperature.