diff --git a/src/io/io.cpp b/src/io/io.cpp
index 32fb8804d..c9817cc6d 100644
--- a/src/io/io.cpp
+++ b/src/io/io.cpp
@@ -15,6 +15,7 @@
 #endif  // HDF5
 #include "../grid/grid3D.h"
 #include "../io/io.h"
+#include "../utils/cuda_utilities.h"
 #include "../utils/timing_functions.h"  // provides ScopedTimer
 #ifdef MPI_CHOLLA
   #include "../mpi/mpi_routines.h"
@@ -1865,6 +1866,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     dataset_buffer_my = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real));
     dataset_buffer_mz = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real));
     dataset_buffer_E  = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real));
+  #ifdef MHD
+    std::vector<Real> dataset_buffer_magnetic_x(H.nx_real * H.ny_real);
+    std::vector<Real> dataset_buffer_magnetic_y(H.nx_real * H.ny_real);
+    std::vector<Real> dataset_buffer_magnetic_z(H.nx_real * H.ny_real);
+  #endif  // MHD
   #ifdef DE
     dataset_buffer_GE = (Real *)malloc(H.nx_real * H.ny_real * sizeof(Real));
   #endif
@@ -1875,19 +1881,38 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     // Copy the xy slices to the memory buffers
     for (j = 0; j < H.ny_real; j++) {
       for (i = 0; i < H.nx_real; i++) {
-        id     = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + zslice * H.nx * H.ny;
+        id     = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost, zslice, H.nx, H.ny);
         buf_id = j + i * H.ny_real;
+  #ifdef MHD
+        int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, j + H.n_ghost, zslice, H.nx, H.ny);
+        int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost - 1, zslice, H.nx, H.ny);
+        int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost, zslice - 1, H.nx, H.ny);
+  #endif  // MHD
   #ifdef MPI_CHOLLA
         // When there are multiple processes, check whether this slice is in
         // your domain
         if (zslice >= nz_local_start && zslice < nz_local_start + nz_local) {
-          id = (i + H.n_ghost) + (j + H.n_ghost) * H.nx + (zslice - nz_local_start + H.n_ghost) * H.nx * H.ny;
-  #endif  // MPI_CHOLLA
+          id = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost, zslice - nz_local_start + H.n_ghost, H.nx,
+                                              H.ny);
+    #ifdef MHD
+          int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, j + H.n_ghost,
+                                                      zslice - nz_local_start + H.n_ghost, H.nx, H.ny);
+          int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost - 1,
+                                                      zslice - nz_local_start + H.n_ghost, H.nx, H.ny);
+          int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, j + H.n_ghost,
+                                                      zslice - nz_local_start + H.n_ghost - 1, H.nx, H.ny);
+    #endif  // MHD
+  #endif    // MPI_CHOLLA
           dataset_buffer_d[buf_id]  = C.density[id];
           dataset_buffer_mx[buf_id] = C.momentum_x[id];
           dataset_buffer_my[buf_id] = C.momentum_y[id];
           dataset_buffer_mz[buf_id] = C.momentum_z[id];
           dataset_buffer_E[buf_id]  = C.Energy[id];
+  #ifdef MHD
+          dataset_buffer_magnetic_x[buf_id] = 0.5 * (C.magnetic_x[id] + C.magnetic_x[id_xm1]);
+          dataset_buffer_magnetic_y[buf_id] = 0.5 * (C.magnetic_y[id] + C.magnetic_y[id_ym1]);
+          dataset_buffer_magnetic_z[buf_id] = 0.5 * (C.magnetic_z[id] + C.magnetic_z[id_zm1]);
+  #endif  // MHD
   #ifdef DE
           dataset_buffer_GE[buf_id] = C.GasEnergy[id];
   #endif
@@ -1905,6 +1930,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
           dataset_buffer_my[buf_id] = 0;
           dataset_buffer_mz[buf_id] = 0;
           dataset_buffer_E[buf_id]  = 0;
+    #ifdef MHD
+          dataset_buffer_magnetic_x[buf_id] = 0;
+          dataset_buffer_magnetic_y[buf_id] = 0;
+          dataset_buffer_magnetic_z[buf_id] = 0;
+    #endif  // MHD
     #ifdef DE
           dataset_buffer_GE[buf_id] = 0;
     #endif
@@ -1924,6 +1954,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_my, "/my_xy");
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_mz, "/mz_xy");
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_E, "/E_xy");
+  #ifdef MHD
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_x.data(), "/magnetic_x_xy");
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_y.data(), "/magnetic_y_xy");
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_z.data(), "/magnetic_z_xy");
+  #endif  // MHD
   #ifdef DE
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_GE, "/GE_xy");
   #endif
@@ -1957,6 +1992,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     dataset_buffer_my = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real));
     dataset_buffer_mz = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real));
     dataset_buffer_E  = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real));
+  #ifdef MHD
+    dataset_buffer_magnetic_x.resize(H.nx_real * H.nz_real);
+    dataset_buffer_magnetic_y.resize(H.nx_real * H.nz_real);
+    dataset_buffer_magnetic_z.resize(H.nx_real * H.nz_real);
+  #endif  // MHD
   #ifdef DE
     dataset_buffer_GE = (Real *)malloc(H.nx_real * H.nz_real * sizeof(Real));
   #endif
@@ -1967,19 +2007,38 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     // Copy the xz slices to the memory buffers
     for (k = 0; k < H.nz_real; k++) {
       for (i = 0; i < H.nx_real; i++) {
-        id     = (i + H.n_ghost) + yslice * H.nx + (k + H.n_ghost) * H.nx * H.ny;
+        id     = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice, k + H.n_ghost, H.nx, H.ny);
         buf_id = k + i * H.nz_real;
+  #ifdef MHD
+        int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, yslice, k + H.n_ghost, H.nx, H.ny);
+        int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - 1, k + H.n_ghost, H.nx, H.ny);
+        int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice, k + H.n_ghost - 1, H.nx, H.ny);
+  #endif  // MHD
   #ifdef MPI_CHOLLA
         // When there are multiple processes, check whether this slice is in
         // your domain
         if (yslice >= ny_local_start && yslice < ny_local_start + ny_local) {
-          id = (i + H.n_ghost) + (yslice - ny_local_start + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny;
-  #endif  // MPI_CHOLLA
+          id = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - ny_local_start + H.n_ghost, k + H.n_ghost, H.nx,
+                                              H.ny);
+    #ifdef MHD
+          int id_xm1 = cuda_utilities::compute1DIndex(i + H.n_ghost - 1, yslice - ny_local_start + H.n_ghost,
+                                                      k + H.n_ghost, H.nx, H.ny);
+          int id_ym1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - ny_local_start + H.n_ghost - 1,
+                                                      k + H.n_ghost, H.nx, H.ny);
+          int id_zm1 = cuda_utilities::compute1DIndex(i + H.n_ghost, yslice - ny_local_start + H.n_ghost,
+                                                      k + H.n_ghost - 1, H.nx, H.ny);
+    #endif  // MHD
+  #endif    // MPI_CHOLLA
           dataset_buffer_d[buf_id]  = C.density[id];
           dataset_buffer_mx[buf_id] = C.momentum_x[id];
           dataset_buffer_my[buf_id] = C.momentum_y[id];
           dataset_buffer_mz[buf_id] = C.momentum_z[id];
           dataset_buffer_E[buf_id]  = C.Energy[id];
+  #ifdef MHD
+          dataset_buffer_magnetic_x[buf_id] = 0.5 * (C.magnetic_x[id] + C.magnetic_x[id_xm1]);
+          dataset_buffer_magnetic_y[buf_id] = 0.5 * (C.magnetic_y[id] + C.magnetic_y[id_ym1]);
+          dataset_buffer_magnetic_z[buf_id] = 0.5 * (C.magnetic_z[id] + C.magnetic_z[id_zm1]);
+  #endif  // MHD
   #ifdef DE
           dataset_buffer_GE[buf_id] = C.GasEnergy[id];
   #endif
@@ -1997,6 +2056,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
           dataset_buffer_my[buf_id] = 0;
           dataset_buffer_mz[buf_id] = 0;
           dataset_buffer_E[buf_id]  = 0;
+    #ifdef MHD
+          dataset_buffer_magnetic_x[buf_id] = 0;
+          dataset_buffer_magnetic_y[buf_id] = 0;
+          dataset_buffer_magnetic_z[buf_id] = 0;
+    #endif  // MHD
     #ifdef DE
           dataset_buffer_GE[buf_id] = 0;
     #endif
@@ -2016,6 +2080,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_my, "/my_xz");
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_mz, "/mz_xz");
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_E, "/E_xz");
+  #ifdef MHD
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_x.data(), "/magnetic_x_xz");
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_y.data(), "/magnetic_y_xz");
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_z.data(), "/magnetic_z_xz");
+  #endif  // MHD
   #ifdef DE
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_GE, "/GE_xz");
   #endif
@@ -2050,6 +2119,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     dataset_buffer_my = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real));
     dataset_buffer_mz = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real));
     dataset_buffer_E  = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real));
+  #ifdef MHD
+    dataset_buffer_magnetic_x.resize(H.ny_real * H.nz_real);
+    dataset_buffer_magnetic_y.resize(H.ny_real * H.nz_real);
+    dataset_buffer_magnetic_z.resize(H.ny_real * H.nz_real);
+  #endif  // MHD
   #ifdef DE
     dataset_buffer_GE = (Real *)malloc(H.ny_real * H.nz_real * sizeof(Real));
   #endif
@@ -2060,19 +2134,37 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     // Copy the yz slices to the memory buffers
     for (k = 0; k < H.nz_real; k++) {
       for (j = 0; j < H.ny_real; j++) {
-        id     = xslice + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny;
+        id     = cuda_utilities::compute1DIndex(xslice, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny);
         buf_id = k + j * H.nz_real;
+  #ifdef MHD
+        int id_xm1 = cuda_utilities::compute1DIndex(xslice - 1, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny);
+        int id_ym1 = cuda_utilities::compute1DIndex(xslice, j + H.n_ghost - 1, k + H.n_ghost, H.nx, H.ny);
+        int id_zm1 = cuda_utilities::compute1DIndex(xslice, j + H.n_ghost, k + H.n_ghost - 1, H.nx, H.ny);
+  #endif  // MHD
   #ifdef MPI_CHOLLA
         // When there are multiple processes, check whether this slice is in
         // your domain
         if (xslice >= nx_local_start && xslice < nx_local_start + nx_local) {
-          id = (xslice - nx_local_start) + (j + H.n_ghost) * H.nx + (k + H.n_ghost) * H.nx * H.ny;
-  #endif  // MPI_CHOLLA
+          id = cuda_utilities::compute1DIndex(xslice - nx_local_start, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny);
+    #ifdef MHD
+          int id_xm1 =
+              cuda_utilities::compute1DIndex(xslice - nx_local_start - 1, j + H.n_ghost, k + H.n_ghost, H.nx, H.ny);
+          int id_ym1 =
+              cuda_utilities::compute1DIndex(xslice - nx_local_start, j + H.n_ghost - 1, k + H.n_ghost, H.nx, H.ny);
+          int id_zm1 =
+              cuda_utilities::compute1DIndex(xslice - nx_local_start, j + H.n_ghost, k + H.n_ghost - 1, H.nx, H.ny);
+    #endif  // MHD
+  #endif    // MPI_CHOLLA
           dataset_buffer_d[buf_id]  = C.density[id];
           dataset_buffer_mx[buf_id] = C.momentum_x[id];
           dataset_buffer_my[buf_id] = C.momentum_y[id];
           dataset_buffer_mz[buf_id] = C.momentum_z[id];
           dataset_buffer_E[buf_id]  = C.Energy[id];
+  #ifdef MHD
+          dataset_buffer_magnetic_x[buf_id] = 0.5 * (C.magnetic_x[id] + C.magnetic_x[id_xm1]);
+          dataset_buffer_magnetic_y[buf_id] = 0.5 * (C.magnetic_y[id] + C.magnetic_y[id_ym1]);
+          dataset_buffer_magnetic_z[buf_id] = 0.5 * (C.magnetic_z[id] + C.magnetic_z[id_zm1]);
+  #endif  // MHD
   #ifdef DE
           dataset_buffer_GE[buf_id] = C.GasEnergy[id];
   #endif
@@ -2090,6 +2182,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
           dataset_buffer_my[buf_id] = 0;
           dataset_buffer_mz[buf_id] = 0;
           dataset_buffer_E[buf_id]  = 0;
+    #ifdef MHD
+          dataset_buffer_magnetic_x[buf_id] = 0;
+          dataset_buffer_magnetic_y[buf_id] = 0;
+          dataset_buffer_magnetic_z[buf_id] = 0;
+    #endif  // MHD
     #ifdef DE
           dataset_buffer_GE[buf_id] = 0;
     #endif
@@ -2109,6 +2206,11 @@ void Grid3D::Write_Slices_HDF5(hid_t file_id)
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_my, "/my_yz");
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_mz, "/mz_yz");
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_E, "/E_yz");
+  #ifdef MHD
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_x.data(), "/magnetic_x_yz");
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_y.data(), "/magnetic_y_yz");
+    status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_magnetic_z.data(), "/magnetic_z_yz");
+  #endif  // MHD
   #ifdef DE
     status = Write_HDF5_Dataset(file_id, dataspace_id, dataset_buffer_GE, "/GE_yz");
   #endif