Starting from an atomic interaction model of an alloy with a clustering tendency, a mesoscopic thermodynamic simulation method is developed. The alloy is represented by a mesoscopic rigid lattice of cells characterized by the number of solute atoms they contain and a discrete cellular Monte Carlo (CMC) algorithm is designed to perform simulations at this scale. The central quantity is a mesoscopic free energy, whose components (local free energy and stiffness parameter) are extracted from the equilibrium properties of single-cell and two-cell systems. Particular attention has been paid to properly incorporate, at mesoscale, the correlations between neighboring cells with identical concentrations. Finally, the ability of the CMC method to compute thermodynamic properties such as phase diagrams, interface profiles, and equilibrium fluctuation spectra is illustrated on a model alloy as well as on an Ising model with interactions up to the second-nearest neighbors, previously developed for body-centered-cubic Fe-Cu alloys. The very good agreement between the CMC calculation results for different cell sizes and direct computations at the atomic scale indicate that the proposed scale change is thermodynamically consistent.