A three-dimensional magnetohydrodynamic (MHD) simulation code has been newly developed. In this code, the physical variables are defined in the staggered mesh system and each of the MHD equations is split into the advection phase and the non-advection phase. The constrained interpolation profile method, which gives a highly accurate solution with less numerical diffusion, is applied to the calculation for the advection phase. In contrast, the non-advection phase is calculated using the finite difference method. The developed code is tested for the relaxation processes of the mass density, plasma pressure and magnetic field in the Large Helical Device (LHD).