In this paper, we present a new approach to treating many-body molecular dynamics on coupled electronic surfaces. The method is based on a semiclassical limit of the quantum Liouville equation. The formal result is a set of coupled classical-like partial differential equations for generalized distribution functions which describe both the nuclear probability densities on the coupled surfaces and the coherences between the electronic states. The Hamiltonian dynamics underlying the evolution of these distributions is augmented by nonclassical source and sink terms, which allow the flow of probability between the coupled surfaces and the corresponding formation and decay of electronic coherences. The formal results are shown analytically to reproduce the well-known Rabi and Landau-Zener results in appropriate limits. In addition, a direct numerical solution of the phase space partial differential equations is performed, and the results compared with exact quantum solutions for a model curve-crossing problem, yielding excellent agreement. Future trajectory-based implementation of the method in molecular dynamics simulations is also discussed.