A parallel finite-difference algorithm for the solution of diffusive, three-dimensional (3D) transient electromagnetic field simulations is presented. The purpose of the scheme is the simulation of both electric fields and the time derivative of magnetic fields generated by galvanic sources (grounded wires) over arbitrarily complicated distributions of conductivity and magnetic permeability. Using a staggered grid and a modified DuFort-Frankel method, the scheme steps Maxwell's equations in time. Electric field initialization is done by a conjugate-gradient solution of a 3D Poisson problem, as is common in 3D resistivity modeling. Instead of calculating the initial magnetic field directly, its time derivative and curl are employed in order to advance the electric field in time. A divergence-free condition is enforced for both the magnetic-field time derivative and the total conduction-current density, providing accurate results at late times. In order to simulate large realistic earth models, the algorithm has been designed to run on parallel computer platforms. The upward continuation boundary condition for a stable solution in the infinitely resistive air layer involves a two-dimensional parallel fast Fourier transform. Example simulations are compared with analytical, integral-equation and spectral Lanczos decomposition solutions and demonstrate the accuracy of the scheme.