diff --git a/sources/integrals.F90 b/sources/integrals.F90
index e3d2227..47db729 100644
--- a/sources/integrals.F90
+++ b/sources/integrals.F90
@@ -684,9 +684,9 @@ module integrals
 ! the integral of |Bx|
 !
 #if NDIMS == 3
-        inarr(11) = inarr(11) + sum(abs(pdata%u(ibx,nb:ne,nb:ne,nb:ne))) * dvol
+        inarr(11) = inarr(11) + sum(abs(pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * dvol
 #else /* NDIMS == 3 */
-        inarr(11) = inarr(11) + sum(abs(pdata%u(ibx,nb:ne,nb:ne,  :  ))) * dvol
+        inarr(11) = inarr(11) + sum(abs(pdata%q(ibx,nb:ne,nb:ne,  :  ))) * dvol
 #endif /* NDIMS == 3 */
 
 ! magnetic energy
@@ -701,10 +701,10 @@ module integrals
 !
 #if NDIMS == 3
           inarr(12) = inarr(12)                                                &
-                    + 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nb,nb:ne))) * dxz
+                    + 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,nb,nb:ne))) * dxz
 #else /* NDIMS == 3 */
           inarr(12) = inarr(12)                                                &
-                    + 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nb,  :  ))) * dxz
+                    + 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,nb,  :  ))) * dxz
 #endif /* NDIMS == 3 */
 
 ! advection of Bx along Y
@@ -765,12 +765,12 @@ module integrals
 !
 #if NDIMS == 3
           inarr(22) = inarr(22)                                                &
-                    + 2.5d-01 * sum((pdata%u(ibx,nb:ne,nb,nb:ne)**2            &
-                                   + pdata%u(ibz,nb:ne,nb,nb:ne)**2)) * dxz
+                    + 2.5d-01 * sum((pdata%q(ibx,nb:ne,nb,nb:ne)**2            &
+                                   + pdata%q(ibz,nb:ne,nb,nb:ne)**2)) * dxz
 #else /* NDIMS == 3 */
           inarr(22) = inarr(22)                                                &
-                    + 2.5d-01 * sum((pdata%u(ibx,nb:ne,nb,  :  )**2            &
-                                   + pdata%u(ibz,nb:ne,nb,  :  )**2)) * dxz
+                    + 2.5d-01 * sum((pdata%q(ibx,nb:ne,nb,  :  )**2            &
+                                   + pdata%q(ibz,nb:ne,nb,  :  )**2)) * dxz
 #endif /* NDIMS == 3 */
 
 ! advection of magnetic field along Y
@@ -861,12 +861,12 @@ module integrals
 !
 #if NDIMS == 3
           inarr(22) = inarr(22)                                                &
-                    + 2.5d-01 * sum((pdata%u(ibx,nb:ne,ne,nb:ne)**2            &
-                                   + pdata%u(ibz,nb:ne,ne,nb:ne)**2)) * dxz
+                    + 2.5d-01 * sum((pdata%q(ibx,nb:ne,ne,nb:ne)**2            &
+                                   + pdata%q(ibz,nb:ne,ne,nb:ne)**2)) * dxz
 #else /* NDIMS == 3 */
           inarr(22) = inarr(22)                                                &
-                    + 2.5d-01 * sum((pdata%u(ibx,nb:ne,ne,  :  )**2            &
-                                   + pdata%u(ibz,nb:ne,ne,  :  )**2)) * dxz
+                    + 2.5d-01 * sum((pdata%q(ibx,nb:ne,ne,  :  )**2            &
+                                   + pdata%q(ibz,nb:ne,ne,  :  )**2)) * dxz
 #endif /* NDIMS == 3 */
 
 ! advection of magnetic field along Y