diff --git a/sources/integrals.F90 b/sources/integrals.F90
index b57a5e7..2eb063b 100644
--- a/sources/integrals.F90
+++ b/sources/integrals.F90
@@ -666,7 +666,7 @@ module integrals
 ! current density Z component
 !
             tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,nb ,nb:ne)             &
-                                 -  pdata%u(iby,nbl:nel,nbl,nb:ne)) / dh(1)    &
+                                 -  pdata%u(iby,nbl:nel,nb ,nb:ne)) / dh(1)    &
                                  - (pdata%u(ibx,nb :ne ,nb ,nb:ne)             &
                                   - pdata%u(ibx,nb :ne ,nbl,nb:ne)) / dh(2)
 
@@ -676,7 +676,7 @@ module integrals
 ! current density Z component
 !
             tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,nb ,  :  )             &
-                                 -  pdata%u(iby,nbl:nel,nbl,  :  )) / dh(1)    &
+                                 -  pdata%u(iby,nbl:nel,nb ,  :  )) / dh(1)    &
                                  - (pdata%u(ibx,nb :ne ,nb ,  :  )             &
                                   - pdata%u(ibx,nb :ne ,nbl,  :  )) / dh(2)
 
@@ -735,7 +735,7 @@ module integrals
 #if NDIMS == 3
 ! current density Z component
 !
-            tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,neu,nb:ne)             &
+            tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,ne ,nb:ne)             &
                                  -  pdata%u(iby,nbl:nel,ne ,nb:ne)) / dh(1)    &
                                  - (pdata%u(ibx,nb :ne ,neu,nb:ne)             &
                                   - pdata%u(ibx,nb :ne ,ne ,nb:ne)) / dh(2)
@@ -745,7 +745,7 @@ module integrals
 #else /* NDIMS == 3 */
 ! current density Z component
 !
-            tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,neu,  :  )             &
+            tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,ne ,  :  )             &
                                  -  pdata%u(iby,nbl:nel,ne ,  :  )) / dh(1)    &
                                  - (pdata%u(ibx,nb :ne ,neu,  :  )             &
                                   - pdata%u(ibx,nb :ne ,ne ,  :  )) / dh(2)
@@ -785,7 +785,7 @@ module integrals
 ! current density Y component
 !
               tmp(:,:,1) = 0.5d+00 * (pdata%u(ibz,nbu:neu,nb:ne,nb )           &
-                                   -  pdata%u(ibz,nbl:nel,nb:ne,nbl)) / dh(1)  &
+                                   -  pdata%u(ibz,nbl:nel,nb:ne,nb )) / dh(1)  &
                                    - (pdata%u(ibx,nb :ne ,nb:ne,nb )           &
                                     - pdata%u(ibx,nb :ne ,nb:ne,nbl)) / dh(3)
 
@@ -819,7 +819,7 @@ module integrals
 
 ! current density Y component
 !
-              tmp(:,:,1) = 0.5d+00 * (pdata%u(ibz,nbu:neu,nb:ne,neu)           &
+              tmp(:,:,1) = 0.5d+00 * (pdata%u(ibz,nbu:neu,nb:ne,ne )           &
                                    -  pdata%u(ibz,nbl:nel,nb:ne,ne )) / dh(1)  &
                                    - (pdata%u(ibx,nb :ne ,nb:ne,neu)           &
                                     - pdata%u(ibx,nb :ne ,nb:ne,ne )) / dh(3)