  We explore the performance and advantages/disadvantages of using
unconditionally stable explicit super time-stepping (STS) algorithms versus
implicit schemes with Krylov solvers for integrating parabolic operators in
thermodynamic MHD models of the solar corona. Specifically, we compare the
second-order Runge-Kutta Legendre (RKL2) STS method with the implicit backward
Euler scheme computed using the preconditioned conjugate gradient (PCG) solver
with both a point-Jacobi and a non-overlapping domain decomposition ILU0
preconditioner. The algorithms are used to integrate anisotropic Spitzer
thermal conduction and artificial kinematic viscosity at time-steps much larger
than classic explicit stability criteria allow. A key component of the
comparison is the use of an established MHD model (MAS) to compute a real-world
simulation on a large HPC cluster. Special attention is placed on the parallel
scaling of the algorithms. It is shown that, for a specific problem and model,
the RKL2 method is comparable or surpasses the implicit method with PCG solvers
in performance and scaling, but suffers from some accuracy limitations. These
limitations, and the applicability of RKL methods are briefly discussed.
