We describe a module for the Athena code that solves the gray equations of radiation hydrodynamics (RHD), based on the first two moments of the radiative transfer equation. We combine explicit Godunov methods to advance the gas and radiation variables including non-stiff source terms with a local implicit method to integrate stiff source terms. We adopt the M₁ closure relation, including leading source terms to $mathcal{O}(betatau)$ and employ the reduced speed of light approximation (RSLA) with subcycling of the radiation variables to reduce computational costs. We consider self-gravitating fragmentation and evolution of turbulent gaseous clouds, modeling the propagation and interaction of radiation from embedded star clusters that form with the surrounding gas. To model the luminosity sources, we use the star particle algorithm of Gong & Ostriker (2013) based on the particle mesh method combined with an efficient open boundary condition Poisson solver for the self-gravitational potential. Our code is dimensionally unsplit in one, two, and three space dimensions and is parallelized using MPI. The streaming and diffusion limits are well-described by the M₁ closure model, and our implementation shows excellent behavior for a problem with a concentrated radiation source containing both regimes simultaneously. Our operator-split method is ideally suited for problems with a slowly-varying radiation field and dynamical gas flows in which the effect of the RSLA is minimal. We present an analysis of the dispersion relation of RHD linear waves highlighting the conditions of applicability for the RSLA. To demonstrate the accuracy of our method, we utilize a suite of radiation and RHD tests covering a broad range of regimes, including RHD waves, shocks, and equilibria, showing second-order convergence in most cases, and a test to demonstrate the accuracy of particle orbits obtained using our method. Applying our method to the study of feedback-regulated star formation in models of giant molecular clouds, we conclude that the radiation force on dust from reprocessed radiation is an efficient mechanism for cloud disruption, which may be particularly important in super star clusters with deep gravitational potential wells.