  In this paper, we propose, analyze, and test a new fully discrete, efficient,
decoupled, stable, and practically second-order time-stepping algorithm for
computing MHD ensemble flow averages under uncertainties in the initial
conditions and forcing. For each viscosity and magnetic diffusivity pair, the
algorithm picks the largest possible parameter $\theta\in[0,1]$ to avoid the
instability that arises due to the presence of some explicit viscous terms. At
each time step, the algorithm shares the same system matrix with all $J$
realizations but with different right-hand-side vectors. That saves assembling
time and computer memory, allows the reuse of the same preconditioner, and can
take the advantage of block linear solvers. For the proposed algorithm, we
prove stability and convergence rigorously. To illustrate the predicted
convergence rates of our analysis, numerical experiments with manufactured
solutions are given on a unit square domain. Finally, we test the scheme on a
benchmark channel flow over a step problem and it performs well.
