diff --git a/Makefile b/Makefile
index d104f738ba9f4dac3c2fa764d9496e568c895798..6a4c8eecf38798e780ad771287d138d02d9941cb 100644
--- a/Makefile
+++ b/Makefile
@@ -23,7 +23,7 @@ dirs:
 	mkdir -p $(BINDIR)
 	mkdir -p $(OBJDIR)
 	mkdir -p $(MODDIR)
-	mkdir -p $(BACKUPDIR)
+	mkdir -p $(CHCKPTDIR)
 
 src/srcinfo.h:
 	( cd src/srcinfo; $(MAKE);)
@@ -45,24 +45,24 @@ cleanbin:
 $(OBJDIR)/diagnose.o : src/srcinfo.h
 
 FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \
-$(OBJDIR)/coeff_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/convolution_mod.o \
+$(OBJDIR)/coeff_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
 $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \
 $(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o \
 $(OBJDIR)/model_mod.o $(OBJDIR)/mkl_dfti.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o \
 $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o \
-$(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o \
+$(OBJDIR)/grid_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o \
 $(OBJDIR)/utility_mod.o
 
  $(EXEC): $(FOBJ)
 	$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
 
- $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@
 
  $(OBJDIR)/array_mod.o : src/array_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/array_mod.F90 -o $@
 
- $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o  $(OBJDIR)/fourier_grid_mod.o
+ $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o  $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
 
  $(OBJDIR)/basic_mod.o : src/basic_mod.F90 $(OBJDIR)/prec_const_mod.o
@@ -71,16 +71,16 @@ $(OBJDIR)/utility_mod.o
  $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o
 		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@
 
- $(OBJDIR)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/convolution_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
 		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/compute_Sapj.F90 -o $@
 
  $(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/control.F90 -o $@
 
- $(OBJDIR)/convolution_mod.o : src/convolution_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/mkl_dfti.o
-	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/convolution_mod.F90 -o $@
+ $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/mkl_dfti.o $(OBJDIR)/grid_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@
 
- $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnose.F90 -o $@
 
  $(OBJDIR)/diagnostics_par_mod.o : src/diagnostics_par_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
@@ -92,10 +92,10 @@ $(OBJDIR)/utility_mod.o
  $(OBJDIR)/fields_mod.o : src/fields_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@
 
- $(OBJDIR)/fourier_grid_mod.o : src/fourier_grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
-		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_grid_mod.F90 -o $@
+ $(OBJDIR)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@
 
- $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@
 
  $(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
@@ -104,7 +104,7 @@ $(OBJDIR)/utility_mod.o
  $(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@
 
- $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/fourier_grid_mod.o
+ $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
 
  $(OBJDIR)/mkl_dfti.o : $(MKLROOT)/include/mkl_dfti.f90
@@ -113,25 +113,25 @@ $(OBJDIR)/utility_mod.o
  $(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
 
- $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs.F90 -o $@
 
- $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/poisson.F90 -o $@
 
  $(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppexit.F90 -o $@
 
- $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppinit.F90 -o $@
 
  $(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_mod.F90 -o $@
 
- $(OBJDIR)/readinputs.o : src/readinputs.F90  $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/readinputs.o : src/readinputs.F90  $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
 
- $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
+ $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
 
  $(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
@@ -140,5 +140,5 @@ $(OBJDIR)/utility_mod.o
  $(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@
 
- $(OBJDIR)/utility_mod.o : src/utility_mod.F90  $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/fourier_grid_mod.o
+ $(OBJDIR)/utility_mod.o : src/utility_mod.F90  $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
diff --git a/local/dirs.inc b/local/dirs.inc
index cafb3079daeddfe957ff6c88e4020487c5b4ba4f..6ff72f9689ab58aa7d192eb00508a67ee6ea58e2 100644
--- a/local/dirs.inc
+++ b/local/dirs.inc
@@ -1,10 +1,12 @@
 #Local code, binaries, pputils library
-PREFIX   = $(HOME)/Documents/HeLaZ#write the path to the root of your MONOLIZ distro here
+PREFIX   = $(HOME)/Documents/HeLaZ#write the path to the root of your HeLaZ distrib here
 SRCDIR   = $(PREFIX)/src
 BINDIR   = $(PREFIX)/bin
 OBJDIR   = $(PREFIX)/obj
 LIBDIR   = $(PREFIX)/lib
 MODDIR   = $(PREFIX)/mod
 FMDIR    = $(PREFIX)/FM
-FFTWDIR  = /usr/local/fftw-3.3.5
-BACKUPDIR  = $(PREFIX)/backup
+FFTW3DIR = /usr/local/fftw3
+CHCKPTDIR  = $(PREFIX)/checkpoint
+
+# Naming ideas : HeLaZ, MoNoLiT (Moment Non Linear Torroidal)
diff --git a/local/make.inc b/local/make.inc
index 5a62151b39a1333468a861eee36d59c233dcaf26..b1943c29340c8b72027fb998df7f7cc2e5dc1e74 100644
--- a/local/make.inc
+++ b/local/make.inc
@@ -129,6 +129,12 @@ ifdef MKL_LIB
 
 endif
 
+ifdef FFTW3DIR
+      LIBS  += -lfftw3
+      LDIRS += -L/usr/local/fftw3/lib
+      IDIRS += -I/usr/local/fftw3/include
+endif
+
 #
 # Add mumps libraries if preprocessor set
 #
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index 709686d01648e355e7cd85b5f38fb74a3cd2da3b..c615d87be7eb2313904dcaf883ff55beb6d9ad4d 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -23,17 +23,17 @@ CONTAINS
     USE basic
     USE time_integration
     USE array
-    USE fourier_grid
+    USE grid
     use prec_const
     IMPLICIT NONE
 
     COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f
     COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f_rhs
-    INTEGER :: ikr, ikz, istage
+    INTEGER :: istage
 
-    SELECT CASE (updatetlevel)   
+    SELECT CASE (updatetlevel)
     CASE(1)
-      DO ikz=ikzs,ikze  
+      DO ikz=ikzs,ikze
         DO ikr=ikrs,ikre
           DO istage=1,ntimelevel
             f(ikr,ikz,1) = f(ikr,ikz,1) + dt*b_E(istage)*f_rhs(ikr,ikz,istage)
@@ -41,7 +41,7 @@ CONTAINS
         END DO
       END DO
     CASE DEFAULT
-      DO ikz=ikzs,ikze  
+      DO ikz=ikzs,ikze
         DO ikr=ikrs,ikre
           f(ikr,ikz,updatetlevel) = f(ikr,ikz,1);
           DO istage=1,updatetlevel-1
diff --git a/src/auxval.F90 b/src/auxval.F90
index 412809bdde03363898d4759d16d6f12ea5a28b09..85ef9e9b833e10aaa2592ed654fcbd0b153f50c0 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -2,12 +2,12 @@ subroutine auxval
   !   Set auxiliary values, at beginning of simulation
 
   USE basic
-  USE fourier_grid
+  USE grid
   USE array
   ! use mumps_bsplines, only: putrow_mumps => putrow, updtmat_mumps => updtmat, factor_mumps => factor, bsolve_mumps => bsolve ! from BSPLINES
   use prec_const
   IMPLICIT NONE
-  
+
   INTEGER :: irows,irowe, irow, icol
   WRITE(*,*) '=== Set auxiliary values ==='
 
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 859e81ac6c112dd9a07b1f756dc6a5c7bf330e76..735a3b457efa4df3d55a8d211bc8664f631f8aa4 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -4,22 +4,22 @@ MODULE basic
   use prec_const
   IMPLICIT none
 
-  LOGICAL  :: RESTART=.true.    ! To restart the simulation or look for saved state
-  INTEGER  :: nrun=1            ! Number of time steps to run
-  real(dp) :: tmax=100000.0     ! Maximum simulation time
-  real(dp) :: dt=1.0            ! Time step
-  real(dp) :: time=0            ! Current simulation time (Init from restart file)
-
-  INTEGER :: jobnum=0                  ! Job number
-  INTEGER :: step=0                    ! Calculation step of this run
-  INTEGER :: cstep=0                   ! Current step number (Init from restart file)
-  LOGICAL :: nlend=.FALSE.             ! Signal end of run
-
-  INTEGER :: ierr                      ! flag for MPI error
-  INTEGER :: iframe1d                  ! counting the number of times 1d datasets are outputed (for diagnose)
-  INTEGER :: iframe2d                  ! counting the number of times 2d datasets are outputed (for diagnose)
-  INTEGER :: iframe3d                  ! counting the number of times 3d datasets are outputed (for diagnose)
-  INTEGER :: iframe5d                  ! counting the number of times 5d datasets are outputed (for diagnose)
+  INTEGER  :: nrun   = 1           ! Number of time steps to run
+  real(dp) :: tmax   = 100000.0    ! Maximum simulation time
+  real(dp) :: dt     = 1.0         ! Time step
+  real(dp) :: time   = 0           ! Current simulation time (Init from restart file)
+
+  INTEGER :: jobnum  = 0           ! Job number
+  INTEGER :: step    = 0           ! Calculation step of this run
+  INTEGER :: cstep   = 0           ! Current step number (Init from restart file)
+  LOGICAL :: RESTART = .FALSE.     ! Signal end of run
+  LOGICAL :: nlend   = .FALSE.     ! Signal end of run
+
+  INTEGER :: ierr                  ! flag for MPI error
+  INTEGER :: iframe1d              ! counting the number of times 1d datasets are outputed (for diagnose)
+  INTEGER :: iframe2d              ! counting the number of times 2d datasets are outputed (for diagnose)
+  INTEGER :: iframe3d              ! counting the number of times 3d datasets are outputed (for diagnose)
+  INTEGER :: iframe5d              ! counting the number of times 5d datasets are outputed (for diagnose)
 
   !  List of logical file units
   INTEGER :: lu_in   = 90              ! File duplicated from STDIN
@@ -43,10 +43,9 @@ CONTAINS
     use prec_const
     IMPLICIT NONE
 
-    NAMELIST /BASIC/  nrun, dt, tmax
+    NAMELIST /BASIC/  nrun, dt, tmax, RESTART
 
     READ(lu_in,basic)
-    !WRITE(*,basic)
 
   END SUBROUTINE basic_data
   !================================================================================
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index ca2787e34e7db2be358941af43ff0716a0fafed7..f6c2fd2f2a8e2c4e3c3823ae0a39a1a50468c0a0 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -2,9 +2,9 @@ SUBROUTINE compute_Sapj(Sepj, Sipj)
 
   USE array, ONLY : dnjs
   USE basic
-  USE convolution
+  USE fourier
   USE fields!, ONLY : phi, moments_e, moments_i
-  USE fourier_grid
+  USE grid
   USE model
   USE prec_const
   USE time_integration!, ONLY : updatetlevel
@@ -13,7 +13,7 @@ SUBROUTINE compute_Sapj(Sepj, Sipj)
   COMPLEX(dp), DIMENSION(Nkr,Nkz) :: F_, G_, CONV
   COMPLEX(dp), DIMENSION(ips_e:ipe_e, ijs_e:ije_e, Nkr, Nkz) :: Sepj
   COMPLEX(dp), DIMENSION(ips_i:ipe_i, ijs_i:ije_i, Nkr, Nkz) :: Sipj
-  INTEGER :: ip, ij, in, is, ikr, ikz
+  INTEGER :: in, is
   REAL(dp):: kr, kz, kernel, be_2, bi_2, factn
   REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
 
@@ -32,7 +32,7 @@ SUBROUTINE compute_Sapj(Sepj, Sipj)
             kz     = kzarray(ikz)
             be_2   = sigmae2_taue_o2 * (kr**2 + kz**2)
             kernel = be_2**(in-1)/factn * EXP(-be_2)
-            
+
             ! First convolution term
             F_(ikr,ikz) = (kz - kr)  * phi(ikr,ikz) * kernel
             ! Second convolution term
@@ -105,5 +105,5 @@ SUBROUTINE compute_Sapj(Sepj, Sipj)
     ENDDO jloopi
   ENDDO ploopi
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- 
+
 END SUBROUTINE compute_Sapj
diff --git a/src/convolution_mod.F90 b/src/convolution_mod.F90
index b50e917363162b5cb019fc29c6f1431dbea6e596..7ca4cb8073e38ba2d0c2b1886efe59d271dd74d0 100644
--- a/src/convolution_mod.F90
+++ b/src/convolution_mod.F90
@@ -9,7 +9,7 @@ MODULE convolution
   SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D )
 
       USE prec_const
-      USE fourier_grid, ONLY : Nkr, Nkz, Pad
+      USE grid, ONLY : Nkr, Nkz, Pad
       USE MKL_DFTI
 
       IMPLICIT NONE
@@ -53,7 +53,7 @@ MODULE convolution
 
       !CALL backward_FFT( F_2D )
       !CALL backward_FFT( G_2D )
-      
+
       !WRITE(*,*) 'C =  F/(Mkr *Mkz) x G/(Mkr *Mkz)..'
       DO ix=1,Mkr*Mkz
         C_1D(ix) = F_1D(ix)/Mkr/Mkz * G_1D(ix)/Mkr/Mkz
@@ -81,7 +81,7 @@ MODULE convolution
   SUBROUTINE convolve_2D_F2R( F_2D, G_2D, C_2D )
 
     USE prec_const
-    USE fourier_grid, ONLY : Nkr, Nkz, Pad
+    USE grid, ONLY : Nkr, Nkz, Pad
     USE MKL_DFTI
 
     IMPLICIT NONE
@@ -125,12 +125,12 @@ MODULE convolution
   END SUBROUTINE convolve_2D_F2R
 
 
-  !! Compute a forward FFT to go back to Fourier space 
+  !! Compute a forward FFT to go back to Fourier space
   !  - used to save computation after the sum of convolution_2D_F2R in compute_Sapj
   SUBROUTINE forward_FFT( C_2D )
 
     USE prec_const
-    USE fourier_grid, ONLY : Nkr, Nkz, Pad
+    USE grid, ONLY : Nkr, Nkz, Pad
     USE MKL_DFTI
 
     IMPLICIT NONE
@@ -177,7 +177,7 @@ MODULE convolution
   SUBROUTINE backward_FFT( C_2D )
 
     USE prec_const
-    USE fourier_grid, ONLY : Nkr, Nkz, Pad
+    USE grid, ONLY : Nkr, Nkz, Pad
     USE MKL_DFTI
 
     IMPLICIT NONE
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index a532329168b8a6a4a62f065d50d5dff5bdd1c581..47e9dacf784463330e3a5bc86bbc64c67e33a1a0 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -2,9 +2,9 @@ SUBROUTINE diagnose(kstep)
   !   Diagnostics, writing simulation state to disk
 
   USE basic
-  USE fourier_grid
+  USE grid
   USE diagnostics_par
-  USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach
+  USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf
   USE model
   USE initial_par
   USE fields
@@ -23,8 +23,7 @@ SUBROUTINE diagnose(kstep)
   !_____________________________________________________________________________
   !                   1.   Initial diagnostics
 
-  IF (kstep .EQ. 0) THEN
-
+  IF ((kstep .EQ. 0)) THEN
      ! jobnum is either 0 from initialization or some integer values read from chkrst(0)
      IF(jobnum .LE. 99) THEN
          WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
@@ -52,11 +51,11 @@ SUBROUTINE diagnose(kstep)
 
      ! Initialize counter of number of saves for each category
      IF (cstep==0) THEN
-        iframe2d=0
-     END IF
+         iframe2d=0
+     ENDIF
      CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
      IF (cstep==0) THEN
-      iframe5d=0
+         iframe5d=0
      END IF
      CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
 
@@ -70,15 +69,15 @@ SUBROUTINE diagnose(kstep)
      CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
 
      IF (write_Ni00) THEN
-      CALL creatg(fidres, "/data/var2d/Ni00", "Ni00")
-      CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
-      CALL putarr(fidres, "/data/var2d/Ni00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+       CALL creatg(fidres, "/data/var2d/Ni00", "Ni00")
+       CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+       CALL putarr(fidres, "/data/var2d/Ni00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      IF (write_phi) THEN
-      CALL creatg(fidres, "/data/var2d/phi", "phi")
-      CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
-      CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+       CALL creatg(fidres, "/data/var2d/phi", "phi")
+       CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+       CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      !  var5d group (moments)
@@ -86,31 +85,31 @@ SUBROUTINE diagnose(kstep)
      CALL creatd(fidres, rank, dims,  "/data/var5d/time",     "Time t*c_s/R")
      CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
      IF (write_moments) THEN
-        CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
-        CALL putarr(fidres,  "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e),       "p_e",ionode=0)
-        CALL putarr(fidres,  "/data/var5d/moments_e/coordj", jarray_e(ijs_e:ije_e),       "j_e",ionode=0)
-        CALL putarr(fidres, "/data/var5d/moments_e/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
-        CALL putarr(fidres, "/data/var5d/moments_e/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
-
-        CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
-        CALL putarr(fidres,  "/data/var5d/moments_i/coordp", parray_i(ips_i:ipe_i),       "p_i",ionode=0)
-        CALL putarr(fidres,  "/data/var5d/moments_i/coordj", jarray_i(ijs_i:ije_i),       "j_i",ionode=0)
-        CALL putarr(fidres, "/data/var5d/moments_i/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
-        CALL putarr(fidres, "/data/var5d/moments_i/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+       CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
+       CALL putarr(fidres,  "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e),       "p_e",ionode=0)
+       CALL putarr(fidres,  "/data/var5d/moments_e/coordj", jarray_e(ijs_e:ije_e),       "j_e",ionode=0)
+       CALL putarr(fidres, "/data/var5d/moments_e/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+       CALL putarr(fidres, "/data/var5d/moments_e/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+
+       CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
+       CALL putarr(fidres,  "/data/var5d/moments_i/coordp", parray_i(ips_i:ipe_i),       "p_i",ionode=0)
+       CALL putarr(fidres,  "/data/var5d/moments_i/coordj", jarray_i(ijs_i:ije_i),       "j_i",ionode=0)
+       CALL putarr(fidres, "/data/var5d/moments_i/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+       CALL putarr(fidres, "/data/var5d/moments_i/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      IF (write_non_lin) THEN
-        CALL creatg(fidres, "/data/var5d/Sepj", "Sepj")
-        CALL putarr(fidres,  "/data/var5d/Sepj/coordp", parray_e(ips_e:ipe_e),       "p_e",ionode=0)
-        CALL putarr(fidres,  "/data/var5d/Sepj/coordj", jarray_e(ijs_e:ije_e),       "j_e",ionode=0)
-        CALL putarr(fidres, "/data/var5d/Sepj/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
-        CALL putarr(fidres, "/data/var5d/Sepj/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
-
-        CALL creatg(fidres, "/data/var5d/Sipj", "Sipj")
-        CALL putarr(fidres,  "/data/var5d/Sipj/coordp", parray_i(ips_i:ipe_i),       "p_i",ionode=0)
-        CALL putarr(fidres,  "/data/var5d/Sipj/coordj", jarray_i(ijs_i:ije_i),       "j_i",ionode=0)
-        CALL putarr(fidres, "/data/var5d/Sipj/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
-        CALL putarr(fidres, "/data/var5d/Sipj/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+       CALL creatg(fidres, "/data/var5d/Sepj", "Sepj")
+       CALL putarr(fidres,  "/data/var5d/Sepj/coordp", parray_e(ips_e:ipe_e),       "p_e",ionode=0)
+       CALL putarr(fidres,  "/data/var5d/Sepj/coordj", jarray_e(ijs_e:ije_e),       "j_e",ionode=0)
+       CALL putarr(fidres, "/data/var5d/Sepj/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+       CALL putarr(fidres, "/data/var5d/Sepj/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
+
+       CALL creatg(fidres, "/data/var5d/Sipj", "Sipj")
+       CALL putarr(fidres,  "/data/var5d/Sipj/coordp", parray_i(ips_i:ipe_i),       "p_i",ionode=0)
+       CALL putarr(fidres,  "/data/var5d/Sipj/coordj", jarray_i(ijs_i:ije_i),       "j_i",ionode=0)
+       CALL putarr(fidres, "/data/var5d/Sipj/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
+       CALL putarr(fidres, "/data/var5d/Sipj/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      !  Add input namelist variables as attributes of /data/input, defined in srcinfo.h
@@ -119,11 +118,8 @@ SUBROUTINE diagnose(kstep)
      WRITE(*,*)  'AUTHOR=', AUTHOR
      WRITE(*,*)    'HOST=', HOST
 
-     IF(jobnum .LE. 99) THEN
-       WRITE(str,'(a,i2.2)') "/data/input.",jobnum
-     ELSE
-       WRITE(str,'(a,i3.2)') "/data/input.",jobnum
-     END IF
+     WRITE(str,'(a,i2.2)') "/data/input"
+
      rank=0
      CALL creatd(fidres, rank,dims,TRIM(str),'Input parameters')
      CALL attach(fidres, TRIM(str),     "version",  VERSION) !defined in srcinfo.h
@@ -132,12 +128,12 @@ SUBROUTINE diagnose(kstep)
      CALL attach(fidres, TRIM(str),    "execdate", EXECDATE) !defined in srcinfo.h
      CALL attach(fidres, TRIM(str),        "host",     HOST) !defined in srcinfo.h
      CALL attach(fidres, TRIM(str),  "start_time",     time)
-     CALL attach(fidres, TRIM(str), "start_cstep",    cstep)
+     CALL attach(fidres, TRIM(str), "start_cstep",    cstep-1)
      CALL attach(fidres, TRIM(str),          "dt",       dt)
      CALL attach(fidres, TRIM(str),        "tmax",     tmax)
      CALL attach(fidres, TRIM(str),        "nrun",     nrun)
 
-     CALL fourier_grid_outputinputs(fidres, str)
+     CALL grid_outputinputs(fidres, str)
 
      CALL output_par_outputinputs(fidres, str)
 
@@ -160,8 +156,17 @@ SUBROUTINE diagnose(kstep)
 
      CALL putfile(fidres, TRIM(str), TRIM(fname),ionode=0)
 
-  END IF
+   ELSEIF((kstep .EQ. 0)) THEN
+
+      IF(jobnum .LE. 99) THEN
+         WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
+      ELSE
+         WRITE(resfile,'(a,a1,i3.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
+      END IF
 
+      CALL openf(resfile,fidres, 'D');
+
+   ENDIF
 
   !_____________________________________________________________________________
   !                   2.   Periodic diagnostics
@@ -182,6 +187,9 @@ SUBROUTINE diagnose(kstep)
      IF (nsave_2d .NE. 0) THEN
         IF (MOD(cstep, nsave_2d) == 0) THEN
            CALL diagnose_2d
+           IF (MOD(cstep, nsave_2d*10) == 0) THEN
+            WRITE(*,"(F4.0,A,F4.0)") time,"/",tmax
+           ENDIF
         END IF
      END IF
 
@@ -193,9 +201,9 @@ SUBROUTINE diagnose(kstep)
      END IF
 
      !                       2.5   Backups
-     !
-     !
-     !
+     IF (MOD(cstep, nsave_cp) == 0) THEN
+       CALL checkpoint_save
+     ENDIF
 
   !_____________________________________________________________________________
   !                   3.   Final diagnostics
@@ -254,7 +262,7 @@ CONTAINS
 
   SUBROUTINE write_field2d(field, text)
     USE futils, ONLY: attach, putarr
-    USE fourier_grid, ONLY: ikrs,ikre, ikzs,ikze
+    USE grid, ONLY: ikrs,ikre, ikzs,ikze
     USE prec_const
     IMPLICIT NONE
 
@@ -278,12 +286,11 @@ SUBROUTINE diagnose_5d
    USE futils, ONLY: append, getatt, attach, putarrnd
    USE fields
    USE array!, ONLY: Sepj, Sipj
-   USE fourier_grid, ONLY: ips_e,ipe_e, ips_i, ipe_i, &
-                           ijs_e,ije_e, ijs_i, ije_i
+   USE grid, ONLY: ips_e,ipe_e, ips_i, ipe_i, &
+                   ijs_e,ije_e, ijs_i, ije_i
    USE time_integration
    USE diagnostics_par
    USE prec_const
-   USE model, ONLY: NON_LIN
    IMPLICIT NONE
 
    CALL append(fidres,  "/data/var5d/time",           time,ionode=0)
@@ -297,7 +304,7 @@ SUBROUTINE diagnose_5d
       CALL write_field5d_i(moments_i(:,:,:,:,updatetlevel), 'moments_i')
    END IF
 
-   IF (write_non_lin .AND. NON_LIN) THEN
+   IF (write_non_lin) THEN
       CALL write_field5d_e(Sepj, 'Sepj')
       CALL write_field5d_i(Sipj, 'Sipj')
    END IF
@@ -306,7 +313,7 @@ SUBROUTINE diagnose_5d
 
    SUBROUTINE write_field5d_e(field, text)
      USE futils, ONLY: attach, putarr
-     USE fourier_grid, ONLY: ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze
+     USE grid,   ONLY: ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze
      USE prec_const
      IMPLICIT NONE
 
@@ -324,7 +331,7 @@ SUBROUTINE diagnose_5d
 
    SUBROUTINE write_field5d_i(field, text)
       USE futils, ONLY: attach, putarr
-      USE fourier_grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze
+      USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze
       USE prec_const
       IMPLICIT NONE
 
@@ -341,3 +348,35 @@ SUBROUTINE diagnose_5d
     END SUBROUTINE write_field5d_i
 
 END SUBROUTINE diagnose_5d
+
+SUBROUTINE checkpoint_save
+  USE basic
+  USE grid
+  USE diagnostics_par
+  USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf
+  USE model
+  USE initial_par
+  USE fields
+  USE time_integration
+  IMPLICIT NONE
+
+  WRITE(rstfile,'(a,a3)') TRIM(rstfile0),'.h5'
+
+  CALL creatf(rstfile, fidrst, real_prec='d')
+  CALL creatg(fidrst, '/Basic', 'Basic data')
+  CALL attach(fidrst, '/Basic', 'cstep', cstep)
+  CALL attach(fidrst, '/Basic', 'time', time)
+  CALL attach(fidrst, '/Basic', 'jobnum', jobnum)
+  CALL attach(fidrst, '/Basic', 'dt', dt)
+  CALL attach(fidrst, '/Basic', 'iframe2d', iframe2d)
+  CALL attach(fidrst, '/Basic', 'iframe5d', iframe5d)
+
+  ! Write state of system to restart file
+  CALL putarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,&
+                                                    ikrs:ikre,ikzs:ikze,1),ionode=0)
+  CALL putarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,&
+                                                    ikrs:ikre,ikzs:ikze,1),ionode=0)
+  CALL closef(fidrst)
+  WRITE(*,'(3x,a)') "Checkpoint file "//TRIM(rstfile)//" saved!"
+
+END SUBROUTINE checkpoint_save
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index c0b86ad4a36b7f9cf2c339cb7ba798dfcf1bf9d4..a9005ea86262e395fdfc0b29823d84d19c40f7cb 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -12,15 +12,15 @@ MODULE diagnostics_par
   LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision=.FALSE.
 
   INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d
-
+  INTEGER, PUBLIC, PROTECTED :: nsave_cp = 1e3
 
   !  HDF5 file
   CHARACTER(len=64), PUBLIC :: resfile0 = "results"   ! Head of main result file name
   CHARACTER(len=64), PUBLIC :: resfile                ! Main result file
-  INTEGER, PUBLIC :: fidres                           ! FID for resfile
+  INTEGER, PUBLIC           :: fidres                 ! FID for resfile
   CHARACTER(len=64), PUBLIC :: rstfile0 = "restart"   ! Head of restart file name
   CHARACTER(len=64), PUBLIC :: rstfile                ! Full restart file
-  INTEGER, PUBLIC :: fidrst                           ! FID for restart file
+  INTEGER, PUBLIC           :: fidrst                 ! FID for restart file
 
   PUBLIC :: output_par_readinputs, output_par_outputinputs
 
diff --git a/src/fourier_grid_mod.F90 b/src/fourier_grid_mod.F90
deleted file mode 100644
index 3cbcbc76e023947bb5b0bada6e0022da00c50874..0000000000000000000000000000000000000000
--- a/src/fourier_grid_mod.F90
+++ /dev/null
@@ -1,163 +0,0 @@
-MODULE fourier_grid
-  ! Grid module for spatial discretization
-
-  USE prec_const
-  IMPLICIT NONE
-  PRIVATE
-
-  !   GRID Namelist
-  INTEGER,  PUBLIC, PROTECTED :: pmaxe = 2      ! The maximal electron Hermite-moment computed
-  INTEGER,  PUBLIC, PROTECTED :: jmaxe = 2      ! The maximal electron Laguerre-moment computed
-  INTEGER,  PUBLIC, PROTECTED :: pmaxi = 2      ! The maximal ion Hermite-moment computed
-  INTEGER,  PUBLIC, PROTECTED :: jmaxi = 2      ! The maximal ion Laguerre-moment computed
-  INTEGER,  PUBLIC, PROTECTED :: maxj  = 2      ! The maximal ion Laguerre-moment computed
-  INTEGER,  PUBLIC, PROTECTED :: nkr   = 10     ! Number of total internal grid points in kr
-  REAL(dp), PUBLIC, PROTECTED :: krmin = 0._dp  ! kr coordinate for left boundary
-  REAL(dp), PUBLIC, PROTECTED :: krmax = 1._dp  ! kr coordinate for right boundary
-  INTEGER,  PUBLIC, PROTECTED :: nkz   = 10     ! Number of total internal grid points in kz
-  REAL(dp), PUBLIC, PROTECTED :: kzmin = 0._dp  ! kz coordinate for left boundary
-  REAL(dp), PUBLIC, PROTECTED :: kzmax = 1._dp  ! kz coordinate for right boundary
-  INTEGER,  PUBLIC, PROTECTED :: Pad   = 2      ! Zero padding parameter for convolution
-
-  ! Indices of s -> start , e-> END
-  INTEGER, PUBLIC, PROTECTED ::  ips_e,ipe_e, ijs_e,ije_e
-  INTEGER, PUBLIC, PROTECTED ::  ips_i,ipe_i, ijs_i,ije_i
-  INTEGER, PUBLIC, PROTECTED ::  ikrs, ikre, ikzs, ikze
-
-  ! Toroidal direction
-  REAL(dp), PUBLIC, PROTECTED ::  deltakr
-  REAL(dp), PUBLIC, PROTECTED ::  deltakz
-
-  ! Grids containing position in fourier space
-  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray
-  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray
-
-  ! Grid containing the polynomials degrees
-  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e
-  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i
-  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e
-  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i
-
-  ! Public Functions
-  PUBLIC :: set_krgrid, set_kzgrid, set_pj
-  PUBLIC :: fourier_grid_readinputs, fourier_grid_outputinputs
-  PUBLIC :: bare, bari
-
-CONTAINS
-
-  SUBROUTINE set_krgrid
-    USE prec_const
-    IMPLICIT NONE
-    INTEGER :: ikr
-    ! Start and END indices of grid
-    ikrs = 1
-    ikre = nkr
-    ! Grid spacings, precompute some inverses
-    IF ( nkr .GT. 1 ) THEN ! To avoid case with 0 intervals
-      deltakr = (krmax - krmin) / REAL(nkr-1,dp)
-    ELSE
-      deltakr = 0;
-    ENDIF
-    ! Discretized kr positions
-    ALLOCATE(krarray(ikrs:ikre))
-    DO ikr = ikrs,ikre
-      krarray(ikr) = krmin + REAL(ikr-1,dp) * deltakr
-    END DO
-  END SUBROUTINE set_krgrid
-
-  SUBROUTINE set_kzgrid
-    USE prec_const
-    IMPLICIT NONE
-    INTEGER :: ikz
-    ! Start and END indices of grid
-    ikzs = 1
-    ikze = nkz
-    ! Grid spacings, precompute some inverses
-    IF ( nkz .GT. 1 ) THEN
-      deltakz = (kzmax - kzmin) / REAL(nkz-1,dp)
-    ELSE
-      deltakz = 0;
-    ENDIF
-    ! Discretized kz positions
-    ALLOCATE(kzarray(ikzs:ikze))
-    DO ikz = ikzs,ikze
-       kzarray(ikz) = kzmin + REAL(ikz-1,dp) * deltakz
-    END DO
-  END SUBROUTINE set_kzgrid
-
-  SUBROUTINE set_pj
-    USE prec_const
-    IMPLICIT NONE
-    INTEGER :: ip, ij
-    ips_e = 1; ipe_e = pmaxe + 1
-    ips_i = 1; ipe_i = pmaxi + 1
-    ALLOCATE(parray_e(ips_e:ipe_e))
-    ALLOCATE(parray_i(ips_i:ipe_i))
-    DO ip = ips_e,ipe_e; parray_e(ip) = ip-1; END DO
-    DO ip = ips_i,ipe_i; parray_i(ip) = ip-1; END DO
-
-    ijs_e = 1; ije_e = jmaxe + 1
-    ijs_i = 1; ije_i = jmaxi + 1
-    ALLOCATE(jarray_e(ijs_e:ije_e))
-    ALLOCATE(jarray_i(ijs_i:ije_i))
-    DO ij = ijs_e,ije_e; jarray_e(ij) = ij-1; END DO
-    DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO
-    maxj  = MAX(jmaxi, jmaxe)
-  END SUBROUTINE set_pj
-
-  SUBROUTINE fourier_grid_readinputs
-    ! Read the input parameters
-
-    USE basic, ONLY : lu_in
-
-    USE prec_const
-    IMPLICIT NONE
-
-    NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
-                    nkr, krmin, krmax, nkz, kzmin, kzmax, &
-                    Pad
-
-    READ(lu_in,grid)
-    !WRITE(*,grid)
-
-  END SUBROUTINE fourier_grid_readinputs
-
-
-  SUBROUTINE fourier_grid_outputinputs(fidres, str)
-    ! Write the input parameters to the results_xx.h5 file
-
-    USE futils, ONLY: attach
-
-    USE prec_const
-    IMPLICIT NONE
-
-    INTEGER, INTENT(in) :: fidres
-    CHARACTER(len=256), INTENT(in) :: str
-    CALL attach(fidres, TRIM(str), "pmaxe", pmaxe)
-    CALL attach(fidres, TRIM(str), "jmaxe", jmaxe)
-    CALL attach(fidres, TRIM(str), "pmaxi", pmaxi)
-    CALL attach(fidres, TRIM(str), "jmaxi", jmaxi)
-    CALL attach(fidres, TRIM(str),   "nkr",   nkr)
-    CALL attach(fidres, TRIM(str), "krmin", krmin)
-    CALL attach(fidres, TRIM(str), "krmax", krmax)
-    CALL attach(fidres, TRIM(str),   "nkz",   nkz)
-    CALL attach(fidres, TRIM(str), "kzmin", kzmin)
-    CALL attach(fidres, TRIM(str), "kzmax", kzmax)
-    CALL attach(fidres, TRIM(str), "Pad", Pad)
-  END SUBROUTINE fourier_grid_outputinputs
-
-
-  FUNCTION bare(p_,j_)
-    IMPLICIT NONE
-    INTEGER :: bare, p_, j_
-    bare = (jmaxe+1)*p_ + j_ + 1
-  END FUNCTION
-
-  FUNCTION bari(p_,j_)
-    IMPLICIT NONE
-    INTEGER :: bari, p_, j_
-    bari = (jmaxi+1)*p_ + j_ + 1
-  END FUNCTION
-
-
-END MODULE fourier_grid
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..7bb0bc2bd9e2dbb2404811c93015837048f09616
--- /dev/null
+++ b/src/fourier_mod.F90
@@ -0,0 +1,129 @@
+MODULE fourier
+
+  USE prec_const
+  USE grid
+  use, intrinsic :: iso_c_binding
+  implicit none
+
+  INCLUDE 'fftw3.f03'
+
+  PRIVATE
+  PUBLIC :: fft_r2cc
+  PUBLIC :: ifft_cc2r
+  PUBLIC :: fft2_r2cc
+  PUBLIC :: ifft2_cc2r
+  PUBLIC :: convolve_2D_F2F
+  PUBLIC :: set_descriptors, free_descriptors ! Temporary to switch easily with MKL DFTI
+
+  CONTAINS
+
+  SUBROUTINE fft_r2cc(fx_in, Fk_out)
+
+    IMPLICIT NONE
+    REAL(dp),    DIMENSION(Nr),  INTENT(IN) :: fx_in
+    COMPLEX(dp), DIMENSION(Nkr), INTENT(OUT):: Fk_out
+    integer*8 plan
+
+    call dfftw_plan_dft_r2c_1d(plan,Nr,fx_in,Fk_out,FFTW_FORWARD,FFTW_ESTIMATE)
+    call dfftw_execute_dft_r2c(plan, fx_in, Fk_out)
+    call dfftw_destroy_plan(plan)
+
+  END SUBROUTINE fft_r2cc
+
+
+
+  SUBROUTINE ifft_cc2r(Fk_in, fx_out)
+
+    IMPLICIT NONE
+    COMPLEX(dp), DIMENSION(Nkr),  INTENT(IN) :: Fk_in
+    REAL(dp),    DIMENSION(Nr),   INTENT(OUT):: fx_out
+    integer*8 plan
+
+    call dfftw_plan_dft_c2r_1d(plan,Nr,Fk_in,fx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
+    call dfftw_execute_dft_c2r(plan, Fk_in, fx_out)
+    call dfftw_destroy_plan(plan)
+    fx_out = fx_out/Nr
+
+  END SUBROUTINE ifft_cc2r
+
+
+
+  SUBROUTINE fft2_r2cc( ffx_in, FFk_out )
+
+      IMPLICIT NONE
+      REAL(dp),    DIMENSION(Nr,Nz), INTENT(IN)   :: ffx_in
+      COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT):: FFk_out
+      integer*8 plan
+
+      !!! 2D Forward FFT ________________________!
+      call dfftw_plan_dft_r2c_2d(plan,Nr,Nz,ffx_in,FFk_out,FFTW_FORWARD,FFTW_ESTIMATE)
+      call dfftw_execute_dft_r2c(plan,ffx_in,FFk_out)
+      call dfftw_destroy_plan(plan)
+
+  END SUBROUTINE fft2_r2cc
+
+
+
+  SUBROUTINE ifft2_cc2r( FFk_in, ffx_out )
+
+      IMPLICIT NONE
+      COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN) :: FFk_in
+      REAL(dp),    DIMENSION(Nr,Nz), INTENT(OUT)  :: ffx_out
+      COMPLEX(dp), DIMENSION(Nkr,Nkz) :: tmp_c
+      integer*8 plan
+
+      tmp_c = FFk_in
+      !!! 2D Backward FFT ________________________!
+      call dfftw_plan_dft_c2r_2d(plan,Nr,Nz,tmp_c,ffx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
+      call dfftw_execute_dft_c2r(plan,tmp_c,ffx_out)
+      call dfftw_destroy_plan(plan)
+      ffx_out = ffx_out/Nr/Nz
+
+  END SUBROUTINE ifft2_cc2r
+
+
+
+  !!! Convolution 2D Fourier to Fourier
+  !   - Compute the convolution using the convolution theorem and MKL
+  SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D )
+
+    USE grid, ONLY : AA_r, AA_z
+    IMPLICIT NONE
+    COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN)  :: F_2D, G_2D  ! input fields
+    COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT) :: C_2D  ! output convolutioned field
+    REAL(dp), DIMENSION(Nr,Nz) :: ff, gg ! iFFT of input fields
+    REAL(dp), DIMENSION(Nr,Nz) :: ffgg  ! will receive the product of f*g in Real
+    INTEGER :: ikr, ikz
+    REAL    :: a_r
+
+    ! 2D inverse Fourier transform
+    CALL ifft2_cc2r(F_2D,ff);
+    CALL ifft2_cc2r(G_2D,gg);
+
+    ! Product in physical space
+    ffgg = ff * gg;
+
+    ! 2D Fourier tranform
+    CALL fft2_r2cc(ffgg,C_2D);
+
+    ! Anti aliasing (2/3 rule)
+    DO ikr = 1,Nkr
+      a_r = AA_r(ikr)
+      DO ikz = 1,Nkz
+         C_2D(ikr,ikz) = C_2D(ikr,ikz) * a_r * AA_z(ikz)
+      ENDDO
+    ENDDO
+
+END SUBROUTINE convolve_2D_F2F
+
+
+! Empty set/free routines to switch easily with MKL DFTI
+SUBROUTINE set_descriptors
+
+END SUBROUTINE set_descriptors
+
+SUBROUTINE free_descriptors
+
+END SUBROUTINE free_descriptors
+
+END MODULE fourier
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..c239c8df72abf0266ad3dd94454b5dfcb16ae517
--- /dev/null
+++ b/src/grid_mod.F90
@@ -0,0 +1,240 @@
+MODULE grid
+  ! Grid module for spatial discretization
+  USE prec_const
+  IMPLICIT NONE
+  PRIVATE
+
+  !   GRID Namelist
+  INTEGER,  PUBLIC, PROTECTED :: pmaxe = 1      ! The maximal electron Hermite-moment computed
+  INTEGER,  PUBLIC, PROTECTED :: jmaxe = 1      ! The maximal electron Laguerre-moment computed
+  INTEGER,  PUBLIC, PROTECTED :: pmaxi = 1      ! The maximal ion Hermite-moment computed
+  INTEGER,  PUBLIC, PROTECTED :: jmaxi = 1      ! The maximal ion Laguerre-moment computed
+  INTEGER,  PUBLIC, PROTECTED :: maxj  = 1      ! The maximal ion Laguerre-moment computed
+  INTEGER,  PUBLIC, PROTECTED :: Nr    = 16     ! Number of total internal grid points in r
+  REAL(dp), PUBLIC, PROTECTED :: Lr    = 1._dp  ! horizontal length of the spatial box
+  INTEGER,  PUBLIC, PROTECTED :: Nz    = 16     ! Number of total internal grid points in z
+  REAL(dp), PUBLIC, PROTECTED :: Lz    = 1._dp  ! vertical length of the spatial box
+  INTEGER,  PUBLIC, PROTECTED :: Nkr   = 8      ! Number of total internal grid points in kr
+  REAL(dp), PUBLIC, PROTECTED :: Lkr   = 1._dp  ! horizontal length of the fourier box
+  INTEGER,  PUBLIC, PROTECTED :: Nkz   = 16     ! Number of total internal grid points in kz
+  REAL(dp), PUBLIC, PROTECTED :: Lkz   = 1._dp  ! vertical length of the fourier box
+
+  ! For Orszag filter
+  REAL(dp), PUBLIC, PROTECTED :: two_third_krmax
+  REAL(dp), PUBLIC, PROTECTED :: two_third_kzmax
+
+  ! 1D Antialiasing arrays (2/3 rule)
+  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_r
+  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_z
+
+  ! Grids containing position in physical space
+  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: rarray
+  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray
+  REAL(dp), PUBLIC, PROTECTED ::  deltar,  deltaz
+  INTEGER,  PUBLIC, PROTECTED  ::  irs,  ire,  izs,  ize
+  INTEGER,  PUBLIC :: ir,iz ! counters
+
+  ! Grids containing position in fourier space
+  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray
+  REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray
+  REAL(dp), PUBLIC, PROTECTED ::  deltakr, deltakz
+  INTEGER,  PUBLIC, PROTECTED ::  ikrs, ikre, ikzs, ikze, a
+  INTEGER,  PUBLIC, PROTECTED :: ikr_0, ikz_0 ! Indices of k-grid origin
+  INTEGER,  PUBLIC :: ikr, ikz, ip, ij ! counters
+
+  ! Grid containing the polynomials degrees
+  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e
+  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i
+  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e
+  INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i
+  INTEGER, PUBLIC, PROTECTED ::  ips_e,ipe_e, ijs_e,ije_e
+  INTEGER, PUBLIC, PROTECTED ::  ips_i,ipe_i, ijs_i,ije_i
+
+  ! Public Functions
+  PUBLIC :: set_pj
+  PUBLIC :: set_rgrid,  set_zgrid
+  PUBLIC :: set_krgrid, set_kzgrid
+  PUBLIC :: grid_readinputs, grid_outputinputs
+  PUBLIC :: bare, bari
+CONTAINS
+
+  SUBROUTINE set_pj
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: ip, ij
+    ips_e = 1; ipe_e = pmaxe + 1
+    ips_i = 1; ipe_i = pmaxi + 1
+    ALLOCATE(parray_e(ips_e:ipe_e))
+    ALLOCATE(parray_i(ips_i:ipe_i))
+    DO ip = ips_e,ipe_e; parray_e(ip) = ip-1; END DO
+    DO ip = ips_i,ipe_i; parray_i(ip) = ip-1; END DO
+
+    ijs_e = 1; ije_e = jmaxe + 1
+    ijs_i = 1; ije_i = jmaxi + 1
+    ALLOCATE(jarray_e(ijs_e:ije_e))
+    ALLOCATE(jarray_i(ijs_i:ije_i))
+    DO ij = ijs_e,ije_e; jarray_e(ij) = ij-1; END DO
+    DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO
+    maxj  = MAX(jmaxi, jmaxe)
+  END SUBROUTINE set_pj
+
+  SUBROUTINE set_rgrid
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: ir
+    ! Start and END indices of grid
+    irs = 1
+    ire = Nr
+    ! Grid spacing
+    IF ( Nr .GT. 1 ) THEN ! To avoid case with 0 intervals
+      deltar = Lr / REAL(Nr-1,dp)
+    ELSE
+      deltar = 0;
+    ENDIF
+    ! Discretized r positions ordered as dx*(0 1 2 -3 -2 -1)
+    ALLOCATE(rarray(irs:ire))
+    DO ir = irs,ire
+      rarray(ir) = deltar*(MODULO(ir-1,Nr/2)-Nr/2*FLOOR(2.*real(ir-1)/real(Nr)))
+    END DO
+    rarray(Nr/2+1) = -rarray(Nr/2+1)
+
+  END SUBROUTINE set_rgrid
+
+  SUBROUTINE set_zgrid
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: iz
+    ! Start and END indices of grid
+    izs = 1
+    ize = Nr
+    ! Grid spacing
+    IF ( Nz .GT. 1 ) THEN ! To avoid case with 0 intervals
+      deltaz = Lz / REAL(Nz-1,dp)
+    ELSE
+      deltaz = 0;
+    ENDIF
+    ! Discretized r positions ordered as dx*(0 1 2 -3 -2 -1)
+    ALLOCATE(zarray(irs:ire))
+    DO iz = izs,ize
+      zarray(iz) = deltaz*(MODULO(iz-1,Nz/2)-Nz/2*FLOOR(2.*real(iz-1)/real(Nz)))
+    END DO
+    zarray(Nz/2+1) = -zarray(Nz/2+1)
+
+  END SUBROUTINE set_zgrid
+
+  SUBROUTINE set_krgrid
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: ikr
+
+    Nkr = Nr/2+1 ! Defined only on positive kr since fields are real
+    ! Start and END indices of grid
+    ikrs = 1
+    ikre = Nkr
+    ! Grid spacings
+    deltakr = 2._dp*PI/Nr/deltar
+
+    ! Discretized kr positions ordered as dk*(0 1 2)
+    ALLOCATE(krarray(ikrs:ikre))
+    DO ikr = ikrs,ikre
+      krarray(ikr) = REAL(ikr-1,dp) * deltakr
+      if (krarray(ikr) .EQ. 0) THEN
+        ikr_0 = ikr
+      ENDIF
+    END DO
+
+    ! Orszag 2/3 filter
+    two_third_krmax = 2._dp/3._dp*deltakr*Nkr
+    ALLOCATE(AA_r(ikrs:ikre))
+    DO ikr = ikrs,ikre
+      IF ( (krarray(ikr) .GT. -two_third_krmax) .AND. (krarray(ikr) .LT. two_third_krmax) ) THEN
+        AA_r(ikr) = 1._dp;
+      ELSE
+        AA_r(ikr) = 0._dp;
+      ENDIF
+    END DO
+  END SUBROUTINE set_krgrid
+
+  SUBROUTINE set_kzgrid
+    USE prec_const
+    IMPLICIT NONE
+
+    Nkz = Nz;
+    ! Start and END indices of grid
+    ikzs = 1
+    ikze = Nkz
+    ! Grid spacings
+    deltakz = 2._dp*PI/Nkz/deltaz
+
+    ! Discretized kz positions ordered as dk*(0 1 2 -3 -2 -1)
+    ALLOCATE(kzarray(ikzs:ikze))
+    DO ikz = ikzs,ikze
+      kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz)))
+      if (kzarray(ikz) .EQ. 0) THEN
+        ikz_0 = ikz
+      ENDIF
+    END DO
+    kzarray(Nz/2+1) = -kzarray(Nz/2+1)
+
+    ! Orszag 2/3 filter
+    two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2);
+    ALLOCATE(AA_z(ikzs:ikze))
+    DO ikz = ikzs,ikze
+      IF ( (kzarray(ikz) .GT. -two_third_kzmax) .AND. (kzarray(ikz) .LT. two_third_kzmax) ) THEN
+        AA_z(ikz) = 1._dp;
+      ELSE
+        AA_z(ikz) = 0._dp;
+      ENDIF
+    END DO
+  END SUBROUTINE set_kzgrid
+
+  SUBROUTINE grid_readinputs
+    ! Read the input parameters
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: lu_in   = 90              ! File duplicated from STDIN
+
+    NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
+                    Nr,  Lr,  Nz,  Lz
+    READ(lu_in,grid)
+
+  END SUBROUTINE grid_readinputs
+
+
+  SUBROUTINE grid_outputinputs(fidres, str)
+    ! Write the input parameters to the results_xx.h5 file
+
+    USE futils, ONLY: attach
+
+    USE prec_const
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in) :: fidres
+    CHARACTER(len=256), INTENT(in) :: str
+    CALL attach(fidres, TRIM(str), "pmaxe", pmaxe)
+    CALL attach(fidres, TRIM(str), "jmaxe", jmaxe)
+    CALL attach(fidres, TRIM(str), "pmaxi", pmaxi)
+    CALL attach(fidres, TRIM(str), "jmaxi", jmaxi)
+    CALL attach(fidres, TRIM(str),   "nr",   nr)
+    CALL attach(fidres, TRIM(str),   "Lr",   Lr)
+    CALL attach(fidres, TRIM(str),   "nz",   nz)
+    CALL attach(fidres, TRIM(str),   "Lz",   Lz)
+    CALL attach(fidres, TRIM(str),   "nkr",   nkr)
+    CALL attach(fidres, TRIM(str),   "Lkr",   Lkr)
+    CALL attach(fidres, TRIM(str),   "nkz",   nkz)
+    CALL attach(fidres, TRIM(str),   "Lkz",   Lkz)
+  END SUBROUTINE grid_outputinputs
+
+  FUNCTION bare(p_,j_)
+    IMPLICIT NONE
+    INTEGER :: bare, p_, j_
+    bare = (jmaxe+1)*p_ + j_ + 1
+  END FUNCTION
+
+  FUNCTION bari(p_,j_)
+    IMPLICIT NONE
+    INTEGER :: bari, p_, j_
+    bari = (jmaxi+1)*p_ + j_ + 1
+  END FUNCTION
+
+END MODULE grid
diff --git a/src/inital.F90 b/src/inital.F90
index d2e839c03809ccf29f2de3a04ec3beaaa5c42783..69312a55cbf39fd49bbcf3bdd7e00c2b2ffc626e 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -36,14 +36,13 @@ END SUBROUTINE inital
 !******************************************************************************!
 SUBROUTINE init_moments
   USE basic
-  USE fourier_grid
+  USE grid
   USE fields
   USE initial_par
   USE time_integration
   USE prec_const
   IMPLICIT NONE
 
-  INTEGER :: ip,ij, ikr,ikz
   INTEGER, DIMENSION(12) :: iseedarr
   REAL(dp) :: noise
 
@@ -116,7 +115,7 @@ END SUBROUTINE
 SUBROUTINE build_dnjs_table
   USE basic
   USE array, Only : dnjs
-  USE fourier_grid, Only : jmaxe, jmaxi
+  USE grid, Only : jmaxe, jmaxi
   USE coeff
   IMPLICIT NONE
 
@@ -150,14 +149,14 @@ END SUBROUTINE
 !******************************************************************************!
 SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,6
   USE basic
-  USE fourier_grid
+  USE grid
   USE initial_par
   USE array
   use model
   USE futils, ONLY : openf, getarr, closef
   IMPLICIT NONE
 
-  INTEGER :: ip,ij, ip2,ij2
+  INTEGER :: ip2,ij2
   INTEGER :: fid1, fid2, fid3, fid4
 
   !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
diff --git a/src/memory.F90 b/src/memory.F90
index 0e1e380b059d5e810757bb823c656e436d99198a..58b8b395811f931e83a600a9598b2263de8b88f6 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -4,7 +4,7 @@ SUBROUTINE memory
   USE array
   USE basic
   USE fields
-  USE fourier_grid
+  USE grid
   USE time_integration
   USE model, ONLY: CO, NON_LIN
 
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 6aa516e42cf85d7a1570251a433ccbc37a1d4eb0..0b85811c9caa211046a9e4aacc3ac040c2b726ea 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -4,12 +4,12 @@ SUBROUTINE moments_eq_rhs
   USE time_integration
   USE array
   USE fields
-  USE fourier_grid
+  USE grid
   USE model
   USE prec_const
   IMPLICIT NONE
 
-  INTEGER     :: ip,ij, ikr,ikz, ip2, ij2 ! loops indices
+  INTEGER     :: ip2, ij2 ! loops indices
   REAL(dp)    :: ip_dp, ij_dp
   REAL(dp)    :: kr, kz, kperp2
   REAL(dp)    :: taue_qe_etaB, taui_qi_etaB
diff --git a/src/poisson.F90 b/src/poisson.F90
index ee353c0c5bada4c218d79c7f58d86dff16f28937..abfa72693d9189ce1638edd1d2955055c5418f47 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -4,13 +4,13 @@ SUBROUTINE poisson
   USE time_integration, ONLY: updatetlevel
   USE array
   USE fields
-  USE fourier_grid
+  USE grid
   use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD
 
   USE prec_const
   IMPLICIT NONE
 
-  INTEGER     :: ikr,ikz, ini,ine, i0j
+  INTEGER     :: ini,ine, i0j
   REAL(dp)    :: ini_dp, ine_dp
   REAL(dp)    :: kr, kz, kperp2
   REAL(dp)    :: Kne, Kni ! sub kernel factor for recursive build
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index 879df9c14414d2076925f45e2aac1a6808f506f1..44bfcf8ecd377b0fd47aeb7541604fbc47e187bd 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -1,7 +1,7 @@
 SUBROUTINE readinputs
   ! Additional data specific for a new run
 
-  USE fourier_grid,     ONLY: fourier_grid_readinputs
+  USE grid,             ONLY: grid_readinputs
   USE diagnostics_par,  ONLY: output_par_readinputs
   USE model,            ONLY: model_readinputs
   USE initial_par,      ONLY: initial_readinputs
@@ -11,11 +11,11 @@ SUBROUTINE readinputs
   IMPLICIT NONE
 
   !WRITE(*,'(a/)') '=== Define additional data for a new run ==='
-  
+
   ! The input file must be in the same order as the following read routines commands (eg spacegrid then output then model)
 
   ! Load grid data from input file
-  CALL fourier_grid_readinputs
+  CALL grid_readinputs
 
   ! Load diagnostic options from input file
   CALL output_par_readinputs
@@ -30,4 +30,3 @@ SUBROUTINE readinputs
   CALL time_integration_readinputs
 
 END SUBROUTINE readinputs
-
diff --git a/src/stepon.F90 b/src/stepon.F90
index 678b81e0bac51bee3f974fdaa448b91929a1797f..47d0014479ec51095cee8ec4f6608a6cb6746b71 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -5,14 +5,14 @@ SUBROUTINE stepon
   USE time_integration
   USE fields, ONLY: moments_e, moments_i, phi
   USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj
-  USE fourier_grid
+  USE grid
   USE advance_field_routine, ONLY: advance_time_level, advance_field
   USE model
   USE utility, ONLY: checkfield
   use prec_const
   IMPLICIT NONE
 
-  INTEGER :: num_step, ip,ij, ikr, ikz
+  INTEGER :: num_step
 
 
    DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index 24445b854fb373b7f67c224c6e38eb67ce6fb7e2..16dd278c21772b01d4ae0288ee24c68e44538ce4 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -53,7 +53,7 @@ CONTAINS
 
   FUNCTION checkfield(field,str) RESULT(mlend)
 
-    USE fourier_grid
+    USE grid
 
     use prec_const
     IMPLICIT NONE