This paper reports on three-dimensional, time-dependent numerical simulations of dynamics and thermodynamics in the diffuse interstellar medium (ISM). Our models are local, account for sheared galactic rotation, magnetic fields, realistic cooling, and thermal conduction, and resolve scales ~ 1-200 pc. This combination permits the study of quasi-steady-state turbulence in a cloudy medium representing the warm/cold atomic ISM. Turbulence is driven by the magnetorotational instability (MRI); our models are the first to study the saturated state of MRI under strongly inhomogeneous conditions, with cloud/intercloud density and temperature contrasts ~ 100. For volume-averaged densities n=0.25 - 4 particles/cc, the mean saturated-state velocity dispersion ranges from 8-1 km/s, with a scaling log(delta v) proportional to -0.77 log(n). The MRI is therefore likely quite important in driving turbulence in low-density regions of the ISM, both away from the midplane in the inner Galaxy (as observed at high latitudes), and throughout the far outer Galaxy (where the mean density drops and the disk flares). The MRI may even be key to suppressing star formation at large radii in spiral galaxies, where the pressure can be high enough that without MRI-driven turbulence, a gravitationally-unstable cold layer would form. As expected, we find that turbulence affects the thermal structure of the ISM. In all our simulations, the fraction of thermally-unstable gas increases as the MRI develops, and in the saturated state is largest in high velocity dispersion models. The mass fractions of warm-stable and unstable gas are typically comparable, in agreement with observations. While inclusion of resistive dissipation of magnetic fields could enhance the amount of thermally-unstable gas compared to current models, our present results indicate that even high levels of turbulence cannot wipe out the signature of thermal instability, and that a shift to a ``phase continuum'' description is probably unwarranted. Instead, we find that temperature and density PDFs are broadened (and include extreme departures from equilibrium), but retain the bimodal character of the classical two-phase description. Our presentation also includes results on the distribution of clump masses (the mass spectrum peaks at ~ 100 solar masses), comparisons of saturated-state MRI scalings with single-phase simulation results (we find the mean magnetic field strength is independent of the mean density), and examples of synthetic HI line profile maps (showing that physical clumps are not easily distinguished in velocity components, and vice versa).

For the associated movie in MPEG format (27 MB), click here.