The random field Ising model with Gaussian disorder is studied using a different Monte Carlo algorithm. The algorithm combines the advantages 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 systems of size 243 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.