diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 4185a1a1b7080a16ea79d708546cbe624aeb53a8..1b4dabc587037c4a34a6b51e48617b6bcd558bc0 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -32,8 +32,8 @@ SUBROUTINE apply_closure_model
             moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
             moments_e(ipeg_e  ,ij,ikx,iky,iz,updatetlevel) = 0._dp
           ENDDO
-          kernel_e(ijsg_e,ikx,iky)      = 0._dp
-          kernel_e(ijeg_e,ikx,iky)      = 0._dp
+          kernel_e(ijsg_e,ikx,iky,iz)      = 0._dp
+          kernel_e(ijeg_e,ikx,iky,iz)      = 0._dp
 
           DO ip = ipsg_i,ipeg_i
             moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
@@ -45,8 +45,8 @@ SUBROUTINE apply_closure_model
             moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
             moments_i(ipeg_i  ,ij,ikx,iky,iz,updatetlevel) = 0._dp
           ENDDO
-          kernel_i(ijsg_i,ikx,iky)      = 0._dp
-          kernel_i(ijeg_i,ikx,iky)      = 0._dp
+          kernel_i(ijsg_i,ikx,iky,iz)      = 0._dp
+          kernel_i(ijeg_i,ikx,iky,iz)      = 0._dp
         ENDDO
       ENDDO
     ENDDO
@@ -66,8 +66,8 @@ SUBROUTINE apply_closure_model
             moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
             moments_e(ipeg_e  ,ij,ikx,iky,iz,updatetlevel) = 0._dp
           ENDDO
-          kernel_e(ijsg_e,ikx,iky)      = 0._dp
-          kernel_e(ijeg_e,ikx,iky)      = 0._dp
+          kernel_e(ijsg_e,ikx,iky,iz)      = 0._dp
+          kernel_e(ijeg_e,ikx,iky,iz)      = 0._dp
 
           DO ip = ipsg_i,ipeg_i
             moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
@@ -79,8 +79,8 @@ SUBROUTINE apply_closure_model
             moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
             moments_i(ipeg_i  ,ij,ikx,iky,iz,updatetlevel) = 0._dp
           ENDDO
-          kernel_i(ijsg_i,ikx,iky)      = 0._dp
-          kernel_i(ijeg_i,ikx,iky)      = 0._dp
+          kernel_i(ijsg_i,ikx,iky,iz)      = 0._dp
+          kernel_i(ijeg_i,ikx,iky,iz)      = 0._dp
         ENDDO
       ENDDO
     ENDDO
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index b6dca64698d3a81f74a91b5b42f8487668a93fda..5fe8fa392ecc6e9210f249007618d689b5e3fc02 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -43,7 +43,7 @@ CONTAINS
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
       ! Get adiabatic moment
-      TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*be_2) * qe_taue * Kernel_e(ij_,ikx_,iky_)*phi(ikx_,iky_,iz_)
+      TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*be_2) * qe_taue * Kernel_e(ij_,ikx_,iky_,iz_)*phi(ikx_,iky_,iz_)
         !** build required fluid moments **
         n_     = 0._dp
         upar_  = 0._dp; uperp_ = 0._dp
@@ -51,9 +51,9 @@ CONTAINS
         DO in_ = 1,jmaxe+1
           n_dp = REAL(in_-1,dp)
           ! Store the kernels for sparing readings
-          Knp0 =  Kernel_e(in_,ikx_,iky_)
-          Knp1 =  Kernel_e(in_+1,ikx_,iky_)
-          Knm1 =  Kernel_e(in_-1,ikx_,iky_)
+          Knp0 =  Kernel_e(in_  ,ikx_,iky_,iz_)
+          Knp1 =  Kernel_e(in_+1,ikx_,iky_,iz_)
+          Knm1 =  Kernel_e(in_-1,ikx_,iky_,iz_)
           ! Nonadiabatic moments (only different from moments when p=0)
           nadiab_moment_0j   = moments_e(1,in_  ,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
           ! Density
@@ -67,10 +67,10 @@ CONTAINS
         ENDDO
       T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
       ! Add energy restoring term
-      TColl_ = TColl_ + T_* 4._dp *  j_dp          * Kernel_e(ij_  ,ikx_,iky_)
-      TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_)
-      TColl_ = TColl_ - T_* 2._dp *  j_dp          * Kernel_e(ij_-1,ikx_,iky_)
-      TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_) - Kernel_e(ij_-1,ikx_,iky_))
+      TColl_ = TColl_ + T_* 4._dp *  j_dp          * Kernel_e(ij_  ,ikx_,iky_,iz_)
+      TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_,iz_)
+      TColl_ = TColl_ - T_* 2._dp *  j_dp          * Kernel_e(ij_-1,ikx_,iky_,iz_)
+      TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_,iz_) - Kernel_e(ij_-1,ikx_,iky_,iz_))
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -79,9 +79,9 @@ CONTAINS
       upar_  = 0._dp
       DO in_ = 1,jmaxe+1
         ! Parallel velocity
-         upar_  = upar_  + Kernel_e(in_,ikx_,iky_) * moments_e(2,in_,ikx_,iky_,iz_,updatetlevel)
+         upar_  = upar_  + Kernel_e(in_,ikx_,iky_,iz_) * moments_e(2,in_,ikx_,iky_,iz_,updatetlevel)
       ENDDO
-      TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_)
+      TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,iz_)
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -93,9 +93,9 @@ CONTAINS
       DO in_ = 1,jmaxe+1
         n_dp = REAL(in_-1,dp)
         ! Store the kernels for sparing readings
-        Knp0 =  Kernel_e(in_  ,ikx_,iky_)
-        Knp1 =  Kernel_e(in_+1,ikx_,iky_)
-        Knm1 =  Kernel_e(in_-1,ikx_,iky_)
+        Knp0 =  Kernel_e(in_  ,ikx_,iky_,iz_)
+        Knp1 =  Kernel_e(in_+1,ikx_,iky_,iz_)
+        Knm1 =  Kernel_e(in_-1,ikx_,iky_,iz_)
         ! Nonadiabatic moments (only different from moments when p=0)
         nadiab_moment_0j = moments_e(1,in_,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
         ! Density
@@ -106,7 +106,7 @@ CONTAINS
         Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
       ENDDO
       T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
-      TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_)
+      TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_,iz_)
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
     ENDIF
@@ -149,7 +149,7 @@ CONTAINS
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
       ! Get adiabatic moment
-      TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bi_2) * qi_taui * Kernel_i(ij_,ikx_,iky_)*phi(ikx_,iky_,iz_)
+      TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bi_2) * qi_taui * Kernel_i(ij_,ikx_,iky_,iz_)*phi(ikx_,iky_,iz_)
         !** build required fluid moments **
         n_     = 0._dp
         upar_  = 0._dp; uperp_ = 0._dp
@@ -157,9 +157,9 @@ CONTAINS
         DO in_ = 1,jmaxi+1
           n_dp = REAL(in_-1,dp)
           ! Store the kernels for sparing readings
-          Knp0 =  Kernel_i(in_,ikx_,iky_)
-          Knp1 =  Kernel_i(in_+1,ikx_,iky_)
-          Knm1 =  Kernel_i(in_-1,ikx_,iky_)
+          Knp0 =  Kernel_i(in_  ,ikx_,iky_,iz_)
+          Knp1 =  Kernel_i(in_+1,ikx_,iky_,iz_)
+          Knm1 =  Kernel_i(in_-1,ikx_,iky_,iz_)
           ! Nonadiabatic moments (only different from moments when p=0)
           nadiab_moment_0j   = moments_i(1,in_  ,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
           ! Density
@@ -173,10 +173,10 @@ CONTAINS
         ENDDO
         T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
       ! Add energy restoring term
-      TColl_ = TColl_ + T_* 4._dp *  j_dp          * Kernel_i(ij_  ,ikx_,iky_)
-      TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_)
-      TColl_ = TColl_ - T_* 2._dp *  j_dp          * Kernel_i(ij_-1,ikx_,iky_)
-      TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_) - Kernel_i(ij_-1,ikx_,iky_))
+      TColl_ = TColl_ + T_* 4._dp *  j_dp          * Kernel_i(ij_  ,ikx_,iky_,iz_)
+      TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_,iz_)
+      TColl_ = TColl_ - T_* 2._dp *  j_dp          * Kernel_i(ij_-1,ikx_,iky_,iz_)
+      TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_,iz_) - Kernel_i(ij_-1,ikx_,iky_,iz_))
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -185,9 +185,9 @@ CONTAINS
       upar_  = 0._dp
       DO in_ = 1,jmaxi+1
         ! Parallel velocity
-         upar_  = upar_  + Kernel_i(in_,ikx_,iky_) * moments_i(2,in_,ikx_,iky_,iz_,updatetlevel)
+         upar_  = upar_  + Kernel_i(in_,ikx_,iky_,iz_) * moments_i(2,in_,ikx_,iky_,iz_,updatetlevel)
       ENDDO
-      TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_)
+      TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,iz_)
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -199,9 +199,9 @@ CONTAINS
       DO in_ = 1,jmaxi+1
         n_dp = REAL(in_-1,dp)
         ! Store the kernels for sparing readings
-        Knp0 =  Kernel_i(in_  ,ikx_,iky_)
-        Knp1 =  Kernel_i(in_+1,ikx_,iky_)
-        Knm1 =  Kernel_i(in_-1,ikx_,iky_)
+        Knp0 =  Kernel_i(in_  ,ikx_,iky_,iz_)
+        Knp1 =  Kernel_i(in_+1,ikx_,iky_,iz_)
+        Knm1 =  Kernel_i(in_-1,ikx_,iky_,iz_)
         ! Nonadiabatic moments (only different from moments when p=0)
         nadiab_moment_0j = moments_i(1,in_,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
         ! Density
@@ -212,7 +212,7 @@ CONTAINS
         Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
       ENDDO
       T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
-      TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_)
+      TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_,iz_)
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
     ENDIF
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index e66080a71e5cc09d09bb39307d3a55f680718839..107e334979fb64a4770e25ccc79e7019dd8cc5e0 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -54,7 +54,7 @@ zloop: DO iz = izs,ize
           kyloope: DO iky = ikys,ikye ! Loop over ky
             kx     = kxarray(ikx)
             ky     = kyarray(iky)
-            kerneln = kernel_e(in, ikx, iky)
+            kerneln = kernel_e(in, ikx, iky, iz)
 
             ! First convolution terms
             Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
@@ -145,7 +145,7 @@ zloop: DO iz = izs,ize
           kyloopi: DO iky = ikys,ikye ! Loop over ky
             kx      = kxarray(ikx)
             ky      = kyarray(iky)
-            kerneln = kernel_i(in, ikx, iky)
+            kerneln = kernel_i(in, ikx, iky, iz)
 
             ! First convolution terms
             Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
diff --git a/src/poisson.F90 b/src/poisson.F90
index 6ec3a2fd38ee42899ade6480f7c1365f1a69f15f..00ad9de7ad31981141db5a8fcefc4431590d5275 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -37,7 +37,7 @@ SUBROUTINE poisson
           sum_kernel2_e    = 0._dp
           ! loop over n only if the max polynomial degree
           DO ine=1,jmaxe+1 ! ine = n+1
-            Kne = kernel_e(ine, ikx, iky)
+            Kne = kernel_e(ine,ikx,iky,iz)
             sum_kernel_mom_e  = sum_kernel_mom_e  + Kne * moments_e(1, ine, ikx, iky, iz, updatetlevel)
             sum_kernel2_e     = sum_kernel2_e     + Kne**2 ! ... sum recursively ...
           END DO
@@ -47,7 +47,7 @@ SUBROUTINE poisson
             sum_kernel2_i    = 0._dp
             ! loop over n only if the max polynomial degree
             DO ini=1,jmaxi+1
-              Kni = kernel_i(ini, ikx, iky)
+              Kni = kernel_i(ini,ikx,iky,iz)
               sum_kernel_mom_i  = sum_kernel_mom_i  + Kni * moments_i(1, ini, ikx, iky, iz, updatetlevel)
               sum_kernel2_i     = sum_kernel2_i     + Kni**2 ! ... sum recursively ...
             END DO
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 129a4c1755ac8cf5aed20d95305d4f8a4ffcc00e..b3e5dd94975b5f62f6aba9a658dfca0d587272f2 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -34,7 +34,7 @@ SUBROUTINE compute_radial_ion_transport
                       imagu * ky_ * moments_i(1,1,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
                   DO ij = ijs_i, ije_i
                       pflux_local = pflux_local - &
-                          imagu * ky_ * kernel_i(ij,ikx,iky) * moments_i(1,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
+                          imagu * ky_ * kernel_i(ij,ikx,iky,iz) * moments_i(1,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
                   ENDDO
               ENDDO
           ENDDO
@@ -81,14 +81,14 @@ SUBROUTINE compute_density
             ! electron density
             dens_e(ikx,iky,iz) = 0._dp
             DO ij = ijs_e, ije_e
-                dens_e(ikx,iky,iz) = dens_e(ikx,iky,iz) + kernel_e(ij,ikx,iky) * &
-                 (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky)*phi(ikx,iky,iz))
+                dens_e(ikx,iky,iz) = dens_e(ikx,iky,iz) + kernel_e(ij,ikx,iky,iz) * &
+                 (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz)*phi(ikx,iky,iz))
             ENDDO
             ! ion density
             dens_i(ikx,iky,iz) = 0._dp
             DO ij = ijs_i, ije_i
-                dens_i(ikx,iky,iz) = dens_i(ikx,iky,iz) + kernel_i(ij,ikx,iky) * &
-                (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky)*phi(ikx,iky,iz))
+                dens_i(ikx,iky,iz) = dens_i(ikx,iky,iz) + kernel_i(ij,ikx,iky,iz) * &
+                (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iz))
             ENDDO
           ENDDO
         ENDDO
@@ -118,9 +118,9 @@ SUBROUTINE compute_temperature
             DO ij = ijs_e, ije_e
               j_dp = REAL(ij-1,dp)
               temp_e(ikx,iky,iz) = temp_e(ikx,iky,iz) + &
-                2._dp/3._dp * (2._dp*j_dp*kernel_e(ij,ikx,iky) - (j_dp+1)*kernel_e(ij+1,ikx,iky) - j_dp*kernel_e(ij-1,ikx,iky))&
-                 * (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky)*phi(ikx,iky,iz)) + &
-                SQRT2/3._dp * kernel_e(ij,ikx,iky) * moments_e(3,ij,ikx,iky,iz,updatetlevel)
+                2._dp/3._dp * (2._dp*j_dp*kernel_e(ij,ikx,iky,iz) - (j_dp+1)*kernel_e(ij+1,ikx,iky,iz) - j_dp*kernel_e(ij-1,ikx,iky,iz))&
+                 * (moments_e(1,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz)*phi(ikx,iky,iz)) &
+                + SQRT2/3._dp * kernel_e(ij,ikx,iky,iz) * moments_e(3,ij,ikx,iky,iz,updatetlevel)
             ENDDO
 
             ! ion temperature
@@ -128,9 +128,9 @@ SUBROUTINE compute_temperature
             DO ij = ijs_i, ije_i
               j_dp = REAL(ij-1,dp)
               temp_i(ikx,iky,iz) = temp_i(ikx,iky,iz) + &
-                2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky) - (j_dp+1)*kernel_i(ij+1,ikx,iky) - j_dp*kernel_i(ij-1,ikx,iky))&
-                 * (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky)*phi(ikx,iky,iz)) + &
-                SQRT2/3._dp * kernel_i(ij,ikx,iky) * moments_i(3,ij,ikx,iky,iz,updatetlevel)
+                2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky,iz) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz) - j_dp*kernel_i(ij-1,ikx,iky,iz))&
+                 * (moments_i(1,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz)*phi(ikx,iky,iz)) &
+                + SQRT2/3._dp * kernel_i(ij,ikx,iky,iz) * moments_i(3,ij,ikx,iky,iz,updatetlevel)
             ENDDO
           ENDDO
         ENDDO