diff --git a/.gitignore b/.gitignore
index 963a466aad33ec20d3a50f2a45305702af86d82d..c03a314529fbe9cbacdbedcf7bb6102da2ceac79 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,7 +34,7 @@ iCa/
 *.out
 src/srcinfo.h
 src/srcinfo/srcinfo.h
-fort*.90
+*/fort*.90
 local/
-*.sh
+*/*.sh
 Gallery/
diff --git a/batch_script.sh b/batch_script.sh
new file mode 100644
index 0000000000000000000000000000000000000000..34ffd2e77898150641c4307868a4771bb5694f8b
--- /dev/null
+++ b/batch_script.sh
@@ -0,0 +1,13 @@
+#!/bin/bash
+#SBATCH --job-name=HeLaZ
+#SBATCH --time=24:00:00
+#SBATCH --nodes=1
+#SBATCH --cpus-per-task=1
+#SBATCH --ntasks-per-node=48
+#SBATCH --mem=64GB
+#SBATCH --error=err.txt
+#SBATCH --output=out.txt
+#SBATCH --account=FUA36_TSVVT422
+#SBATCH --partition=skl_fua_dbg
+module load autoload hdf5 fftw
+srun --cpu-bind=cores ./../../../bin/helaz3 2 12 2
diff --git a/fort.90 b/fort.90
new file mode 100644
index 0000000000000000000000000000000000000000..97a735157d7d992f2ef536b8391a3e604e47e284
--- /dev/null
+++ b/fort.90
@@ -0,0 +1,80 @@
+&BASIC
+  nrun   = 100000000
+  dt     = 0.01
+  tmax   = 200
+  maxruntime = 356400
+/
+&GRID
+  pmaxe  = 4
+  jmaxe  = 2
+  pmaxi  = 4
+  jmaxi  = 2
+  Nx     = 128
+  Lx     = 100
+  Ny     = 128
+  Ly     = 120
+  Nz     = 16
+  q0     = 1.4
+  shear  = 0
+  eps    = 0.18
+  SG     = .t.
+/
+&OUTPUT_PAR
+  nsave_0d = 200
+  nsave_1d = -1
+  nsave_2d = 200
+  nsave_3d = 200
+  nsave_5d = 2000
+  write_doubleprecision = .f.
+  write_gamma = .t.
+  write_hf    = .t.
+  write_phi   = .t.
+  write_Na00  = .t.
+  write_Napj  = .t.
+  write_Sapj  = .f.
+  write_dens  = .t.
+  write_temp  = .t.
+  job2load    = -1
+/
+&MODEL_PAR
+  ! Collisionality
+  CLOS    = 0
+  NL_CLOS = 0
+  LINEARITY = 'nonlinear'
+  KIN_E   = .f.
+  mu_x    = 0.0
+  mu_y    = 0.0
+  mu_z    = 0.2
+  mu_p    = 0
+  mu_j    = 0
+  nu      = 0.01
+  tau_e   = 1
+  tau_i   = 1
+  sigma_e = 0.023338
+  sigma_i = 1
+  q_e     = -1
+  q_i     = 1
+  K_n     = 2.22
+  K_T     = 6.96
+  K_E     = 0
+  GradB     = 1
+  CurvB     = 1
+  lambdaD = 0
+/
+&COLLISION_PAR
+  collision_model = 'DG'
+  gyrokin_CO      = .true.
+  interspecies    = .true.
+!  mat_file        = '../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+  mat_file        = '../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+  INIT_OPT    = 'ppj'
+  ACT_ON_MODES       = 'donothing'
+  init_background  = 0
+  init_noiselvl = 0.00005
+  iseed         = 42
+/
+&TIME_INTEGRATION_PAR
+  numerical_scheme = 'RK4'
+/
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index bafdc82eaea46b86504498e9337bea8c10d31cbe..841f80d2f6ec0a858bd14c0cba09c375d0c0cd9f 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -46,7 +46,8 @@ MODULE array
   REAL(dp), DIMENSION(:),   ALLOCATABLE :: xnipp1j, xnipm1j, xnipp2j, xnipm2j, xnipjp1, xnipjm1
   REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ynipp1j, ynipm1j,   ynipp1jm1, ynipm1jm1 ! mirror lin coeff for non adiab mom
   REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zNipm1j, zNipm1jp1, zNipm1jm1            ! mirror lin coeff for adiab mom
-  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_e, xphijp1_e, xphijm1_e
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_i, xphijp1_i, xphijm1_i
   ! Geometrical operators
   ! Curvature
   REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky  ! dimensions: kx, ky, z, odd/even p
diff --git a/src/memory.F90 b/src/memory.F90
index 72d8331fa510c12ff98de194aecbe9c848d22851..5c94ee7e4e4b9255b38625153c74557b5b593c4f 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -86,9 +86,14 @@ SUBROUTINE memory
   CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1)
 
   ! elect. pot. linear terms
-  CALL allocate_array( xphij,   ips_i,ipe_i, ijs_i,ije_i)
-  CALL allocate_array( xphijp1, ips_i,ipe_i, ijs_i,ije_i)
-  CALL allocate_array( xphijm1, ips_i,ipe_i, ijs_i,ije_i)
+  IF (KIN_E) THEN
+    CALL allocate_array( xphij_e,   ips_e,ipe_e, ijs_e,ije_e)
+    CALL allocate_array( xphijp1_e, ips_e,ipe_e, ijs_e,ije_e)
+    CALL allocate_array( xphijm1_e, ips_e,ipe_e, ijs_e,ije_e)
+  ENDIF
+  CALL allocate_array( xphij_i,   ips_i,ipe_i, ijs_i,ije_i)
+  CALL allocate_array( xphijp1_i, ips_i,ipe_i, ijs_i,ije_i)
+  CALL allocate_array( xphijm1_i, ips_i,ipe_i, ijs_i,ije_i)
 
   ! Curvature and geometry
   CALL allocate_array( Ckxky,   ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 87380472cbea3fd93d0fed2c4ea574d2c5b68c7f..0934b5336b234a9977646d663bd04f2eefb80536 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -99,9 +99,9 @@ SUBROUTINE moments_eq_rhs_e
           Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1
           !! Electrical potential term
           IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
-            Tphi = (xphij  (ip,ij)*kernel_e(ij  ,ikx,iky,iz,eo) &
-                  + xphijp1(ip,ij)*kernel_e(ij+1,ikx,iky,iz,eo) &
-                  + xphijm1(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
+            Tphi = (xphij_i  (ip,ij)*kernel_e(ij  ,ikx,iky,iz,eo) &
+                  + xphijp1_i(ip,ij)*kernel_e(ij+1,ikx,iky,iz,eo) &
+                  + xphijm1_i(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
           ELSE
             Tphi = 0._dp
           ENDIF
@@ -227,9 +227,9 @@ SUBROUTINE moments_eq_rhs_i
 
           !! Electrical potential term
           IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
-            Tphi = (xphij  (ip,ij)*kernel_i(ij  ,ikx,iky,iz,eo) &
-                  + xphijp1(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) &
-                  + xphijm1(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
+            Tphi = (xphij_i  (ip,ij)*kernel_i(ij  ,ikx,iky,iz,eo) &
+                  + xphijp1_i(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) &
+                  + xphijm1_i(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
           ELSE
             Tphi = 0._dp
           ENDIF
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index d8e2b178aa48249d68fbeefc814b4577d6239d4b..ba8d0653e120c44d9ba44c48e120751cf89b3616 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -227,6 +227,27 @@ SUBROUTINE compute_lin_coeff
   ENDDO
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !! ES linear coefficients for moment RHS !!!!!!!!!!
+  IF (KIN_E) THEN
+    DO ip = ips_e, ipe_e
+      p_int= parray_e(ip)   ! Hermite degree
+      DO ij = ijs_e, ije_e
+        j_int= jarray_e(ij)   ! REALof Laguerre degree
+        j_dp = REAL(j_int,dp) ! REALof Laguerre degree
+        !! Electrostatic potential pj terms
+        IF (p_int .EQ. 0) THEN ! kronecker p0
+          xphij_e(ip,ij)    =+K_n + 2.*j_dp*K_T
+          xphijp1_e(ip,ij)  =-K_T*(j_dp+1._dp)
+          xphijm1_e(ip,ij)  =-K_T* j_dp
+        ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
+          xphij_e(ip,ij)    =+K_T/SQRT2
+          xphijp1_e(ip,ij)  = 0._dp; xphijm1_e(ip,ij)  = 0._dp;
+        ELSE
+          xphij_e(ip,ij)    = 0._dp; xphijp1_e(ip,ij)  = 0._dp
+          xphijm1_e(ip,ij)  = 0._dp;
+        ENDIF
+      ENDDO
+    ENDDO
+  ENDIF
   DO ip = ips_i, ipe_i
     p_int= parray_i(ip)   ! Hermite degree
     DO ij = ijs_i, ije_i
@@ -234,19 +255,18 @@ SUBROUTINE compute_lin_coeff
       j_dp = REAL(j_int,dp) ! REALof Laguerre degree
       !! Electrostatic potential pj terms
       IF (p_int .EQ. 0) THEN ! kronecker p0
-        xphij(ip,ij)    =+K_n + 2.*j_dp*K_T
-        xphijp1(ip,ij)  =-K_T*(j_dp+1._dp)
-        xphijm1(ip,ij)  =-K_T* j_dp
+        xphij_i(ip,ij)    =+K_n + 2.*j_dp*K_T
+        xphijp1_i(ip,ij)  =-K_T*(j_dp+1._dp)
+        xphijm1_i(ip,ij)  =-K_T* j_dp
       ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
-        xphij(ip,ij)    =+K_T/SQRT2
-        xphijp1(ip,ij)  = 0._dp; xphijm1(ip,ij)  = 0._dp;
+        xphij_i(ip,ij)    =+K_T/SQRT2
+        xphijp1_i(ip,ij)  = 0._dp; xphijm1_i(ip,ij)  = 0._dp;
       ELSE
-        xphij(ip,ij)    = 0._dp; xphijp1(ip,ij)  = 0._dp
-        xphijm1(ip,ij)  = 0._dp;
+        xphij_i(ip,ij)    = 0._dp; xphijp1_i(ip,ij)  = 0._dp
+        xphijm1_i(ip,ij)  = 0._dp;
       ENDIF
     ENDDO
   ENDDO
-
 END SUBROUTINE compute_lin_coeff
 
 !******************************************************************************!
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index cc69a2b1845375eb87bffede814f6669edf55335..613cf37f8175cc1b1a118e7fb867b39ee48f11ec 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -72,10 +72,10 @@ CONTAINS
           ! Read state of system from checkpoint file
           IF (KIN_E) THEN
           WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
-          CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3/))
+          CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3,5/))
           ENDIF
           WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
-          CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3/))
+          CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3,5/))
 
           CALL closef(fidrst)
 
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index a6f008afc17c076f8e24309f841c00d4113cd20f..c7a82f47017566a668dc9dd95f2a252823f6b314 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -11,7 +11,7 @@ system(['mkdir -p ',MISCDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 10; 
+JOBNUMMIN = 00; JOBNUMMAX = 0; 
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 
 
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index a20170243f45717dbe2bd3467f97fa96fff3bef9..657b52a8a128af26748710fbcdae93caf5735bf5 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -3,6 +3,6 @@
 outfile ='';
 outfile ='linear_shearless_cyclone/dbg';
 % outfile ='debug/test_zpar/';
-% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0_colless/';
+outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0/';
 % outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_0.8/';
 analysis_3D
\ No newline at end of file