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_star=5-60$ M$_odot$, $Omega/Omega_{rm crit} =0.2-1.0$, $B_{rm p} =1-20$ kG). We then employ a representative grid of B-type star models ($M_star=5, 10, 15$ M$_odot$, $Omega/Omega_{rm crit} =0.2 , 0.5, 0.8$, $B_{rm p} = 1, 3 ,10, 30$ kG) to compare to the results of a recent self-consistent 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.