  The random field Ising model with Gaussian disorder is studied using a new
Monte Carlo algorithm. The algorithm combines the advantanges of the replica
exchange method and the two-replica cluster method and is much more efficient
than the Metropolis algorithm for some disorder realizations. Three-dimensional
sytems of size $24^3$ are studied. Each realization of disorder is simulated at
a value of temperature and uniform field that is adjusted to the phase
transition region for that disorder realization. Energy and magnetization
distributions show large variations from one realization of disorder to
another. For some realizations of disorder there are three well separated peaks
in the magnetization distribution and two well separated peaks in the energy
distribution suggesting a first-order transition.
