The time evolution of angular momentum and surface rotation of massive stars is strongly influenced by fossil magnetic fields via magnetic braking. We present a new module containing a simple, comprehensive implementation of such a field at the surface of a massive star within the Modules for Experiments in Stellar Astrophysics (mesa) software instrument. We test two limiting scenarios for magnetic braking: distributing the angular momentum loss throughout the star in the first case, and restricting the angular momentum loss to a surface reservoir in the second case. We perform a systematic investigation of the rotational evolution using a grid of OB star models with surface magnetic fields (M = 5 − 60 M , Ω/Ω crit = 0.2 − 1.0, B p = 1 − 20 kG). We then employ a representative grid of B-type star models (M = 5, 10, 15 M , Ω/Ω crit = 0.2, 0.5, 0.8, B p = 1, 3, 10, 30 kG) to compare to the results of a recent selfconsistent analysis of the sample of known magnetic B-type stars. We infer that magnetic massive stars arrive at the zero age main sequence with a range of rotation rates, rather than with one common value. In particular, some stars are required to have close-to-critical rotation at the ZAMS. However, magnetic braking yields surface rotation rates converging to a common low value, making it difficult to infer the initial rotation rates of evolved, slowly-rotating stars.