We report a high-precision numerical estimation of the critical exponent $alpha$ of the specific heat of the random-field Ising model in four dimensions. Our result $alpha = 0.12(1)$ indicates a diverging specific-heat behavior and is consistent with the estimation coming from the modified hyperscaling relation using our estimate of $theta$ via the anomalous dimensions $eta$ and $bar{eta}$. Our analysis benefited form a high-statistics zero-temperature numerical simulation of the model for two distributions of the random fields, namely a Gaussian and Poissonian distribution, as well as recent advances in finite-size scaling and reweighting methods for disordered systems. An original estimate of the critical slowing down exponent $z$ of the maximum-flow algorithm used is also provided.