Configurational thermodynamics of C in body-centered cubic (bcc) or tetragonal (bct) Fe is investigated combining several computational techniques. Pairwise CC interaction energies in bcc Fe are calculated by density functional theory (DFT) and embedded atom method (EAM) potential respectively. The interaction between C atom and homogeneous strain is calculated assuming C acts as force dipoles in linear elastic medium of Fe. The CC and C–strain interactions are input into Monte Carlo (MC) simulations to find equilibrium C configuration on octahedral interstitial sublattices (OISs) in bcc/bct Fe and corresponding thermodynamic properties. In bcc Fe, DFT-MC and EAM-MC both give a single-phase region with C distributed with disorder on all the three OISs (α) at high temperature, and a two-phase region with ferrite and an ordered compound (α‴) at low temperature. The compound is Fe16C1 according to DFT or Fe16C2 according to EAM inputs, both having two OISs occupied. When a homogeneous tetragonal lattice strain is applied, the disordered phase exhibits a preferential sublattice occupation (Zener ordered α′), which is primarily caused by C–strain interaction and is mitigated by CC interactions. The ordered compound may also have two (α‴) or one (α″) OIS occupied. C clustering in bcc/bct Fe follows a conditional spinodal mechanism, namely long-range order being a prerequisite of spinodal decomposition. This is verified by the difference in solution thermodynamics of α/α′ (ordering-type) and α‴/α″ (clustering-type), as well as kinetic Ising model simulations which reveal a temporal sequence of short-range ordering, long-range ordering, and eventually C clustering.