  In this paper, we propose and analyze a temporally second-order accurate,
fully discrete finite element method for the magnetohydrodynamic (MHD)
equations. A modified Crank--Nicolson method is used to discretize the model
and appropriate semi-implicit treatments are applied to the fluid convection
term and two coupling terms. These semi-implicit approximations result in a
linear system with variable coefficients for which the unique solvability can
be proved theoretically. In addition, we use a decoupling projection method of
the Van Kan type \cite{vankan1986} in the Stokes solver, which computes the
intermediate velocity field based on the gradient of the pressure from the
previous time level, and enforces the incompressibility constraint via the
Helmholtz decomposition of the intermediate velocity field. The energy
stability of the scheme is theoretically proved, in which the decoupled Stokes
solver needs to be analyzed in details. Optimal-order convergence of
$\mathcal{O} (\tau^2+h^{r+1})$ in the discrete $L^\infty(0,T;L^2)$ norm is
proved for the proposed decoupled projection finite element scheme, where
$\tau$ and $h$ are the time stepsize and spatial mesh size, respectively, and
$r$ is the degree of the finite elements. Existing error estimates of
second-order projection methods of the Van Kan type \cite{vankan1986} were only
established in the discrete $L^2(0,T;L^2)$ norm for the Navier--Stokes
equations. Numerical examples are provided to illustrate the theoretical
results.
