  We present a new magnetohydrodynamic (MHD) simulation code with the aim of
providing accurate numerical solutions to astrophysical phenomena where
discontinuities, shock waves, and turbulence are inherently important. The code
implements the HLLD approximate Riemann solver, the
fifth-order-monotonicity-preserving interpolation (MP5) scheme, and the
hyperbolic divergence cleaning method for a magnetic field. This choice of
schemes significantly improved numerical accuracy and stability, and saved
computational costs in multidimensional problems. Numerical tests of one- and
two-dimensional problems showed the advantages of using the high-order scheme
by comparing with results from a standard second-order TVD MUSCL scheme. The
present code enabled us to explore long-term evolution of a three-dimensional
accretion disk around a black hole, in which compressible MHD turbulence caused
continuous mass accretion via nonlinear growth of the magneto-rotational
instability (MRI). Numerical tests with various computational cell sizes
exhibited a convergent picture of the early nonlinear growth of the MRI in a
global model, and indicated that the MP5 scheme has more than twice the
resolution of the MUSCL scheme in practical applications.
