The thermal denaturation (melting) of DNA molecules is a thermodynamic and conformational order-disorder transition from a double-stranded to a single-stranded state. Modeling DNA melting at the atomistic level is severely hampered due to the time and length scales of this transition. We propose a coarse-grained model where DNA strands are represented as oligomer chains of N beads using the single-site bond-fluctuation model on a cubic lattice. This approach incorporates physically relevant characteristics such as the sequence and orientation dependence of base-stacking and base-pairing interactions, as well as the semiflexibility of the chains. We perform parallel tempering Monte Carlo simulations of dilute solutions of short DNA strands in the canonical ensemble. Due to the strong short-ranged and anisotropic nature of the interactions, we employ various biased trials to improve the phase-space sampling. Feedback optimization of the temperature distribution and multihistogram reweighting techniques were used to obtain accurate estimates of the transition temperature. This procedure allows the direct calculation of thermodynamic and conformational properties across the thermally induced order-disorder transition. We explore how the interaction heterogeneity, broad stacking transition and chain stiffness may induce specific heat capacity effects, shift the location of the melting temperature and broaden the transition. Overall, the phenomenological behavior predicted is in qualitative agreement with experimental observations.