  The three-dimensional bimodal random-field Ising model is investigated using
the N-fold version of the Wang-Landau algorithm. The essential energy subspaces
are determined by the recently developed critical minimum energy subspace
technique, and two implementations of this scheme are utilized. The random
fields are obtained from a bimodal discrete $(\pm\Delta)$ distribution, and we
study the model for various values of the disorder strength $\Delta$,
$\Delta=0.5, 1, 1.5$ and 2, on cubic lattices with linear sizes $L=4-24$. We
extract information for the probability distributions of the specific heat
peaks over samples of random fields. This permits us to obtain the phase
diagram and present the finite-size behavior of the specific heat. The question
of saturation of the specific heat is re-examined and it is shown that the open
problem of universality for the random-field Ising model is strongly influenced
by the lack of self-averaging of the model. This property appears to be
substantially depended on the disorder strength.
