diff --git a/src/inital.F90 b/src/inital.F90
index 2dfd084..a2273ff 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -1,285 +1,350 @@
 !******************************************************************************!
 !!!!!! initialize the moments and load/build coeff tables
 !******************************************************************************!
 SUBROUTINE inital
 
   USE basic
   USE model, ONLY : CO, NON_LIN
   USE initial_par
   USE prec_const
   USE time_integration
   USE array, ONLY : Sepj,Sipj
   USE collision
   USE closure
   USE ghosts
   USE restarts
+  USE numerics, ONLY: wipe_turbulence, wipe_zonalflow
   IMPLICIT NONE
 
   CALL set_updatetlevel(1)
 
   !!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!!
   ! through loading a previou state
   IF ( RESTART ) THEN
     IF (my_id .EQ. 0) WRITE(*,*) 'Load moments'
     CALL load_moments ! get N_0
     CALL poisson ! compute phi_0=phi(N_0)
   ! through initialization
   ELSE
     ! set phi with noise and set moments to 0
     IF (INIT_NOISY_PHI) THEN
+      IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
       CALL init_phi
-
     ! set moments_00 (GC density) with noise and compute phi afterwards
     ELSE
-      IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments'
-      CALL init_moments ! init noisy N_0
+      IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy gyrocenter density'
+      CALL init_gyrodens ! init only gyrocenter density
+      ! CALL init_moments ! init all moments randomly (unadvised)
       CALL poisson ! get phi_0 = phi(N_0)
     ENDIF
   ENDIF
 
+  ! Option for wiping the ZF modes (ky==0)
+  IF ( WIPE_ZF .GT. 0 ) THEN
+    IF (my_id .EQ. 0) WRITE(*,*) '-Wiping ZF modes'
+    CALL wipe_zonalflow
+  ENDIF
+
   ! Option for wiping the turbulence and check growth of secondary inst.
-  IF ( WIPE_TURB ) THEN
+  IF ( WIPE_TURB .GT. 0 ) THEN
     IF (my_id .EQ. 0) WRITE(*,*) '-Wiping turbulence'
     CALL wipe_turbulence
   ENDIF
+
   ! Option for initializing a gaussian blob on the zonal profile
   IF ( INIT_BLOB ) THEN
     IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
-    CALL put_blob
+    CALL initialize_blob
   ENDIF
 
   IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
   CALL apply_closure_model
 
   IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
   CALL update_ghosts
 
   !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
   IF ( NON_LIN ) THEN;
     IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
     CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
   ENDIF
 
   !!!!!! Load the COSOlver collision operator coefficients !!!!!!
   IF (ABS(CO) .GT. 1) THEN
     CALL load_COSOlver_mat
     ! Compute collision
     CALL compute_TColl ! compute C_0 = C(N_0)
   ENDIF
 END SUBROUTINE inital
 !******************************************************************************!
 
 !******************************************************************************!
-!!!!!!! Initialize the moments randomly
+!!!!!!! Initialize all the moments randomly
 !******************************************************************************!
 SUBROUTINE init_moments
   USE basic
   USE grid
   USE fields
   USE prec_const
   USE utility, ONLY: checkfield
   USE initial_par
   USE model, ONLY : NON_LIN
   IMPLICIT NONE
 
   REAL(dp) :: noise
   REAL(dp) :: kx, ky, sigma, gain, ky_shift
   INTEGER, DIMENSION(12) :: iseedarr
 
   ! Seed random number generator
   iseedarr(:)=iseed
   CALL RANDOM_SEED(PUT=iseedarr+my_id)
 
     !**** Broad noise initialization *******************************************
     DO ip=ips_e,ipe_e
       DO ij=ijs_e,ije_e
 
         DO ikx=ikxs,ikxe
           DO iky=ikys,ikye
             DO iz=izs,ize
               CALL RANDOM_NUMBER(noise)
               moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
             END DO
           END DO
         END DO
 
         IF ( contains_kx0 ) THEN
           DO iky=2,Nky/2 !symmetry at kx = 0 for all z
             moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
           END DO
         ENDIF
 
       END DO
     END DO
 
     DO ip=ips_i,ipe_i
       DO ij=ijs_i,ije_i
 
         DO ikx=ikxs,ikxe
           DO iky=ikys,ikye
             DO iz=izs,ize
               CALL RANDOM_NUMBER(noise)
               moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
             END DO
           END DO
         END DO
 
         IF ( contains_kx0 ) THEN
           DO iky=2,Nky/2 !symmetry at kx = 0 for all z
             moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
           END DO
         ENDIF
 
       END DO
     END DO
 
     ! Putting to zero modes that are not in the 2/3 Orszag rule
     IF (NON_LIN) THEN
       DO ikx=ikxs,ikxe
       DO iky=ikys,ikye
       DO iz=izs,ize
         DO ip=ips_e,ipe_e
         DO ij=ijs_e,ije_e
           moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
         ENDDO
         ENDDO
         DO ip=ips_i,ipe_i
         DO ij=ijs_i,ije_i
           moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
         ENDDO
         ENDDO
       ENDDO
       ENDDO
       ENDDO
     ENDIF
 END SUBROUTINE init_moments
 !******************************************************************************!
 
-
 !******************************************************************************!
-!!!!!!! Initialize a noisy ES potential and cancel the moments
+!!!!!!! Initialize the gyrocenter density randomly
 !******************************************************************************!
-SUBROUTINE init_phi
+SUBROUTINE init_gyrodens
   USE basic
   USE grid
   USE fields
   USE prec_const
+  USE utility, ONLY: checkfield
   USE initial_par
+  USE model, ONLY : NON_LIN
   IMPLICIT NONE
 
   REAL(dp) :: noise
   REAL(dp) :: kx, ky, sigma, gain, ky_shift
   INTEGER, DIMENSION(12) :: iseedarr
 
-  IF (INIT_NOISY_PHI) THEN
-    IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
-    ! Seed random number generator
-    iseedarr(:)=iseed
-    CALL RANDOM_SEED(PUT=iseedarr+my_id)
+  ! Seed random number generator
+  iseedarr(:)=iseed
+  CALL RANDOM_SEED(PUT=iseedarr+my_id)
 
-      !**** noise initialization *******************************************
+    !**** Broad noise initialization *******************************************
+    DO ip=ips_e,ipe_e
+      DO ij=ijs_e,ije_e
 
-      DO ikx=ikxs,ikxe
-        DO iky=ikys,ikye
-          DO iz=izs,ize
-            CALL RANDOM_NUMBER(noise)
-            phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))*AA_x(ikx)*AA_y(iky)
-          ENDDO
+        DO ikx=ikxs,ikxe
+          DO iky=ikys,ikye
+            DO iz=izs,ize
+              CALL RANDOM_NUMBER(noise)
+              IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
+                moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+              ELSE
+                moments_e(ip,ij,ikx,iky,iz,:) = 0._dp
+              ENDIF
+            END DO
+          END DO
         END DO
+
+        IF ( contains_kx0 ) THEN
+          DO iky=2,Nky/2 !symmetry at kx = 0 for all z
+            moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
+          END DO
+        ENDIF
+
       END DO
+    END DO
 
-      !symmetry at kx = 0 to keep real inverse transform
-      IF ( contains_kx0 ) THEN
-        DO iky=2,Nky/2
-          phi(ikx_0,iky,:) = phi(ikx_0,Nky+2-iky,:)
-        END DO
-        phi(ikx_0,Ny/2,:) = REAL(phi(ikx_0,Ny/2,:)) !origin must be real
-      ENDIF
+    DO ip=ips_i,ipe_i
+      DO ij=ijs_i,ije_i
 
-      !**** ensure no previous moments initialization
-      moments_e = 0._dp; moments_i = 0._dp
+        DO ikx=ikxs,ikxe
+          DO iky=ikys,ikye
+            DO iz=izs,ize
+              CALL RANDOM_NUMBER(noise)
+              IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
+                moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+              ELSE
+                moments_i(ip,ij,ikx,iky,iz,:) = 0._dp
+              ENDIF
+            END DO
+          END DO
+        END DO
 
-      !**** Zonal Flow initialization *******************************************
-      ! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0
-      IF(INIT_ZF .GT. 0) THEN
-      IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
-        IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
-          DO iz = izs,ize
-            phi(INIT_ZF+1,iky_0,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
-            moments_i(1,1,INIT_ZF+1,iky_0,iz,:) = kxarray(INIT_ZF+1)**2*phi(INIT_ZF+1,iky_0,iz)* COS((iz-1)/Nz*2._dp*PI)
-            moments_e(1,1,INIT_ZF+1,iky_0,iz,:) = 0._dp
-          ENDDO
+        IF ( contains_kx0 ) THEN
+          DO iky=2,Nky/2 !symmetry at kx = 0 for all z
+            moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
+          END DO
         ENDIF
-      ENDIF
-    ELSE ! we compute phi from noisy moments and poisson
-      CALL poisson
-    ENDIF
 
-END SUBROUTINE init_phi
+      END DO
+    END DO
+
+    ! Putting to zero modes that are not in the 2/3 Orszag rule
+    IF (NON_LIN) THEN
+      DO ikx=ikxs,ikxe
+      DO iky=ikys,ikye
+      DO iz=izs,ize
+        DO ip=ips_e,ipe_e
+        DO ij=ijs_e,ije_e
+          moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+        ENDDO
+        ENDDO
+        DO ip=ips_i,ipe_i
+        DO ij=ijs_i,ije_i
+          moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+        ENDDO
+        ENDDO
+      ENDDO
+      ENDDO
+      ENDDO
+    ENDIF
+END SUBROUTINE init_gyrodens
 !******************************************************************************!
 
 !******************************************************************************!
-!!!!!!! Remove all ky!=0 modes to conserve only zonal modes in a restart
+!!!!!!! Initialize a noisy ES potential and cancel the moments
 !******************************************************************************!
-SUBROUTINE wipe_turbulence
-  USE fields
+SUBROUTINE init_phi
+  USE basic
   USE grid
+  USE fields
+  USE prec_const
+  USE initial_par
   IMPLICIT NONE
-  DO ikx=ikxs,ikxe
-  DO iky=ikys,ikye
-  DO iz=izs,ize
-    DO ip=ips_e,ipe_e
-    DO ij=ijs_e,ije_e
-      IF( iky .NE. iky_0) THEN
-        moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
-      ELSE
-        moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
-      ENDIF
-    ENDDO
-    ENDDO
-    DO ip=ips_i,ipe_i
-    DO ij=ijs_i,ije_i
-      IF( iky .NE. iky_0) THEN
-        moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
-      ELSE
-        moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
+
+  REAL(dp) :: noise
+  REAL(dp) :: kx, ky, sigma, gain, ky_shift
+  INTEGER, DIMENSION(12) :: iseedarr
+
+  ! Seed random number generator
+  iseedarr(:)=iseed
+  CALL RANDOM_SEED(PUT=iseedarr+my_id)
+
+    !**** noise initialization *******************************************
+
+    DO ikx=ikxs,ikxe
+      DO iky=ikys,ikye
+        DO iz=izs,ize
+          CALL RANDOM_NUMBER(noise)
+          phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
+        ENDDO
+      END DO
+    END DO
+
+    !symmetry at kx = 0 to keep real inverse transform
+    IF ( contains_kx0 ) THEN
+      DO iky=2,Nky/2
+        phi(ikx_0,iky,:) = phi(ikx_0,Nky+2-iky,:)
+      END DO
+      phi(ikx_0,Ny/2,:) = REAL(phi(ikx_0,Ny/2,:)) !origin must be real
+    ENDIF
+
+    !**** ensure no previous moments initialization
+    moments_e = 0._dp; moments_i = 0._dp
+
+    !**** Zonal Flow initialization *******************************************
+    ! put a mode at ikx = mode number + 1, symmetry is already included since kx>=0
+    IF(INIT_ZF .GT. 0) THEN
+      IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
+      IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
+        DO iz = izs,ize
+          phi(INIT_ZF+1,iky_0,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
+          moments_i(1,1,INIT_ZF+1,iky_0,iz,:) = kxarray(INIT_ZF+1)**2*phi(INIT_ZF+1,iky_0,iz)* COS((iz-1)/Nz*2._dp*PI)
+          moments_e(1,1,INIT_ZF+1,iky_0,iz,:) = 0._dp
+        ENDDO
       ENDIF
-    ENDDO
-    ENDDO
-  ENDDO
-  ENDDO
-  ENDDO
-END SUBROUTINE
+    ENDIF
+
+END SUBROUTINE init_phi
+!******************************************************************************!
+
 !******************************************************************************!
 !******************************************************************************!
 !!!!!!! Initialize an ionic Gaussian blob on top of the preexisting modes
 !******************************************************************************!
-SUBROUTINE put_blob
+SUBROUTINE initialize_blob
   USE fields
   USE grid
   USE model, ONLY: sigmai2_taui_o2
   IMPLICIT NONE
   REAL(dp) ::kx, ky, sigma, gain
   sigma = 0.5_dp
   gain  = 5e2_dp
 
   DO ikx=ikxs,ikxe
     kx = kxarray(ikx)
   DO iky=ikys,ikye
     ky = kyarray(iky)
   DO iz=izs,ize
     DO ip=ips_i,ipe_i
     DO ij=ijs_i,ije_i
       IF( (iky .NE. iky_0) .AND. (ip .EQ. 1) .AND. (ij .EQ. 1)) THEN
         moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :) &
         + gain*sigma/SQRT2 * exp(-(kx**2+ky**2)*sigma**2/4._dp) &
           * AA_x(ikx)*AA_y(iky)!&
           ! * exp(sigmai2_taui_o2*(kx**2+ky**2))
       ENDIF
     ENDDO
     ENDDO
   ENDDO
   ENDDO
   ENDDO
-END SUBROUTINE put_blob
+END SUBROUTINE initialize_blob
 !******************************************************************************!
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index fefa911..de90b69 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -1,74 +1,77 @@
 MODULE initial_par
   !   Module for initial parameters
 
   USE prec_const
   IMPLICIT NONE
   PRIVATE
 
   ! Initialization through a noisy phi
   LOGICAL,  PUBLIC, PROTECTED :: INIT_NOISY_PHI = .false.
   ! Initialization through a zonal flow phi
   INTEGER,  PUBLIC, PROTECTED :: INIT_ZF    = 0
   REAL(DP), PUBLIC, PROTECTED :: ZF_AMP     = 1E+3_dp
-  ! Wipe turbulence in the restart
-  LOGICAL,  PUBLIC, PROTECTED :: WIPE_TURB = .false.
+  ! Wipe zonal flow in the restart (=1) or at each step (=2)
+  INTEGER,  PUBLIC, PROTECTED :: WIPE_ZF   = 0
+  ! Wipe turbulence in the restart (=1) or at each step (=2)
+  INTEGER,  PUBLIC, PROTECTED :: WIPE_TURB = 0
   ! Init a Gaussian blob density in the middle
   LOGICAL,  PUBLIC, PROTECTED :: INIT_BLOB = .false.
   ! Initial background level
   REAL(dp), PUBLIC, PROTECTED :: init_background=0._dp
   ! Initial noise amplitude
   REAL(dp), PUBLIC, PROTECTED :: init_noiselvl=1E-6_dp
   ! Initialization for random number generator
   INTEGER,  PUBLIC, PROTECTED :: iseed=42
 
   CHARACTER(len=128), PUBLIC :: mat_file    ! COSOlver matrix file names
 
   PUBLIC :: initial_outputinputs, initial_readinputs
 
 CONTAINS
 
 
   SUBROUTINE initial_readinputs
     ! Read the input parameters
 
     USE basic, ONLY : lu_in, RESTART
     USE prec_const
     IMPLICIT NONE
 
     NAMELIST /INITIAL_CON/ INIT_NOISY_PHI
     NAMELIST /INITIAL_CON/ INIT_ZF
+    NAMELIST /INITIAL_CON/ WIPE_ZF
     NAMELIST /INITIAL_CON/ WIPE_TURB
     NAMELIST /INITIAL_CON/ INIT_BLOB
     NAMELIST /INITIAL_CON/ init_background
     NAMELIST /INITIAL_CON/ init_noiselvl
     NAMELIST /INITIAL_CON/ iseed
     NAMELIST /INITIAL_CON/ mat_file
 
     READ(lu_in,initial_con)
     !WRITE(*,initial_con)
 
   END SUBROUTINE initial_readinputs
 
 
   SUBROUTINE initial_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), "INIT_NOISY_PHI", INIT_NOISY_PHI)
 
     CALL attach(fidres, TRIM(str), "INIT_ZF", INIT_ZF)
 
     CALL attach(fidres, TRIM(str), "init_background", init_background)
 
     CALL attach(fidres, TRIM(str), "init_noiselvl", init_noiselvl)
 
     CALL attach(fidres, TRIM(str), "iseed", iseed)
 
   END SUBROUTINE initial_outputinputs
 
 END MODULE initial_par
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 1289550..e8f1b48 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -1,256 +1,331 @@
 MODULE numerics
     USE basic
     USE prec_const
     USE grid
     USE utility
     USE coeff
     implicit none
 
     PUBLIC :: compute_derivatives, build_dnjs_table, evaluate_kernels, compute_lin_coeff
+    PUBLIC :: wipe_turbulence, wipe_zonalflow
 
 CONTAINS
 
 ! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
 SUBROUTINE compute_derivatives
 
 END SUBROUTINE compute_derivatives
 
 !******************************************************************************!
 !!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin
 !******************************************************************************!
 SUBROUTINE build_dnjs_table
   USE array, Only : dnjs
   USE coeff
   IMPLICIT NONE
 
   INTEGER :: in, ij, is, J
   INTEGER :: n_, j_, s_
 
   J = max(jmaxe,jmaxi)
 
   DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry
     n_ = in - 1
     DO ij = in,J+1
       j_ = ij - 1
       DO is = ij,J+1
         s_ = is - 1
 
         dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0))
         ! By symmetry
         dnjs(in,is,ij) = dnjs(in,ij,is)
         dnjs(ij,in,is) = dnjs(in,ij,is)
         dnjs(ij,is,in) = dnjs(in,ij,is)
         dnjs(is,ij,in) = dnjs(in,ij,is)
         dnjs(is,in,ij) = dnjs(in,ij,is)
       ENDDO
     ENDDO
   ENDDO
 END SUBROUTINE build_dnjs_table
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Evaluate the kernels once for all
 !******************************************************************************!
 SUBROUTINE evaluate_kernels
   USE basic
   USE array, Only : kernel_e, kernel_i, gxy, gyy
   USE grid
   USE model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2
   IMPLICIT NONE
 
   REAL(dp)    :: factj, j_dp, j_int
   REAL(dp)    :: be_2, bi_2, alphaD
   REAL(dp)    :: kx, ky, kperp2
 
   !!!!! Electron kernels !!!!!
   !Precompute species dependant factors
   factj = 1.0 ! Start of the recursive factorial
   DO ij = 1, jmaxe+1
     j_int = jarray_e(ij)
     j_dp = REAL(j_int,dp) ! REAL of degree
 
     ! Recursive factorial
     IF (j_dp .GT. 0) THEN
       factj = factj * j_dp
     ELSE
       factj = 1._dp
     ENDIF
 
     DO ikx = ikxs,ikxe
       kx     = kxarray(ikx)
       DO iky = ikys,ikye
         ky    = kyarray(iky)
         DO iz = izs,ize
           kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
           be_2  =  kperp2 * sigmae2_taue_o2
           kernel_e(ij, ikx, iky,iz) = be_2**j_int * exp(-be_2)/factj
         ENDDO
       ENDDO
     ENDDO
   ENDDO
   ! Kernels closure
   DO ikx = ikxs,ikxe
     kx     = kxarray(ikx)
     DO iky = ikys,ikye
       ky    = kyarray(iky)
       DO iz = izs,ize
         kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
         be_2  = kperp2 * sigmae2_taue_o2
         ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1)
         kernel_e(ijeg_e,ikx,iky,iz) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky,iz)
         ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
         IF ( be_2 .NE. 0 ) THEN
           kernel_e(ijsg_e,ikx,iky,iz) = (real(ijsg_e,dp))/be_2*kernel_e(ijs_e,ikx,iky,iz)
         ELSE
           kernel_e(ijsg_e,ikx,iky,iz) = 0._dp
         ENDIF
       ENDDO
     ENDDO
   ENDDO
 
   !!!!! Ion kernels !!!!!
   factj = 1.0 ! Start of the recursive factorial
   DO ij = 1, jmaxi+1
     j_int = jarray_e(ij)
     j_dp = REAL(j_int,dp) ! REAL of degree
 
     ! Recursive factorial
     IF (j_dp .GT. 0) THEN
       factj = factj * j_dp
     ELSE
       factj = 1._dp
     ENDIF
 
     DO ikx = ikxs,ikxe
       kx     = kxarray(ikx)
       DO iky = ikys,ikye
         ky    = kyarray(iky)
         DO iz = izs,ize
           kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
           bi_2  = kperp2 * sigmai2_taui_o2
           kernel_i(ij, ikx, iky,iz) = bi_2**j_int * exp(-bi_2)/factj
         ENDDO
       ENDDO
     ENDDO
   ENDDO
   ! Kernels closure
   DO ikx = ikxs,ikxe
     kx     = kxarray(ikx)
     DO iky = ikys,ikye
       ky    = kyarray(iky)
       DO iz = izs,ize
         kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
         bi_2  =  kperp2 * sigmai2_taui_o2
         ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj
         kernel_i(ijeg_i,ikx,iky,iz) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky,iz)
         ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
         IF ( bi_2 .NE. 0 ) THEN
           kernel_i(ijsg_i,ikx,iky,iz) = (real(ijsg_i,dp))/bi_2*kernel_i(ijs_i,ikx,iky,iz)
         ELSE
           kernel_i(ijsg_i,ikx,iky,iz) = 0._dp
         ENDIF
       ENDDO
     ENDDO
   ENDDO
 END SUBROUTINE evaluate_kernels
 !******************************************************************************!
 
 SUBROUTINE compute_lin_coeff
   USE array
   USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, eta_T, eta_n
   USE prec_const
   USE grid,  ONLY: parray_e, parray_i, jarray_e, jarray_i, &
                    ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i
   IMPLICIT NONE
   INTEGER     :: p_int, j_int ! polynom. degrees
   REAL(dp)    :: p_dp, j_dp
   REAL(dp)    :: kx, ky, z
   !! Electrons linear coefficients for moment RHS !!!!!!!!!!
   DO ip = ips_e, ipe_e
     p_int= parray_e(ip)   ! Hermite degree
     p_dp = REAL(p_int,dp) ! REAL of Hermite degree
     DO ij = ijs_e, ije_e
       j_int= jarray_e(ij)   ! Laguerre degree
       j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
       xnepj(ip,ij) = taue_qe * 2._dp * (p_dp + j_dp + 1._dp)
       ynepp1j  (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp+1._dp)
       ynepm1j  (ip,ij) = -SQRT(tau_e)/sigma_e * (j_dp+1)*SQRT(p_dp)
       ynepp1jm1(ip,ij) = +SQRT(tau_e)/sigma_e *     j_dp*SQRT(p_dp+1._dp)
       ynepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e *     j_dp*SQRT(p_dp)
       zNepm1j  (ip,ij) = -SQRT(tau_e)/sigma_e * (2._dp*j_dp+1_dp)*SQRT(p_dp)
       zNepm1jp1(ip,ij) = +SQRT(tau_e)/sigma_e *       (j_dp+1_dp)*SQRT(p_dp)
       zNepm1jm1(ip,ij) = +SQRT(tau_e)/sigma_e *              j_dp*SQRT(p_dp)
     ENDDO
   ENDDO
   DO ip = ips_e, ipe_e
     p_int= parray_e(ip)   ! Hermite degree
     p_dp = REAL(p_int,dp) ! REAL of Hermite degree
     xnepp1j(ip) = sqrtTaue_qe * SQRT(p_dp + 1_dp)
     xnepm1j(ip) = sqrtTaue_qe * SQRT(p_dp)
     xnepp2j(ip) = taue_qe * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
     xnepm2j(ip) = taue_qe * SQRT(p_dp * (p_dp - 1._dp))
   ENDDO
   DO ij = ijs_e, ije_e
     j_int= jarray_e(ij)   ! Laguerre degree
     j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
     xnepjp1(ij) = -taue_qe * (j_dp + 1._dp)
     xnepjm1(ij) = -taue_qe * j_dp
   ENDDO
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !! Ions linear coefficients for moment RHS !!!!!!!!!!
   DO ip = ips_i, ipe_i
     p_int= parray_i(ip)   ! Hermite degree
     p_dp = REAL(p_int,dp) ! REAL of Hermite degree
     DO ij = ijs_i, ije_i
       j_int= jarray_i(ij)   ! Laguerre degree
       j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
       xnipj(ip,ij) = taui_qi * 2._dp * (p_dp + j_dp + 1._dp)
       ynipp1j  (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp+1._dp)
       ynipm1j  (ip,ij) = -SQRT(tau_i)/sigma_i * (j_dp+1)*SQRT(p_dp)
       ynipp1jm1(ip,ij) = +SQRT(tau_i)/sigma_i *     j_dp*SQRT(p_dp+1._dp)
       ynipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i *     j_dp*SQRT(p_dp)
       zNipm1j  (ip,ij) = -SQRT(tau_i)/sigma_i * (2._dp*j_dp+1_dp)*SQRT(p_dp)
       zNipm1jp1(ip,ij) = +SQRT(tau_i)/sigma_i *       (j_dp+1_dp)*SQRT(p_dp)
       zNipm1jm1(ip,ij) = +SQRT(tau_i)/sigma_i *              j_dp*SQRT(p_dp)
     ENDDO
   ENDDO
   DO ip = ips_i, ipe_i
     p_int= parray_i(ip)   ! Hermite degree
     p_dp = REAL(p_int,dp) ! REAL of Hermite degree
     xnipp1j(ip) = sqrtTaui_qi * SQRT(p_dp + 1._dp)
     xnipm1j(ip) = sqrtTaui_qi * SQRT(p_dp)
     xnipp2j(ip) = taui_qi * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
     xnipm2j(ip) = taui_qi * SQRT(p_dp * (p_dp - 1._dp))
   ENDDO
   DO ij = ijs_i, ije_i
     j_int= jarray_i(ij)   ! Laguerre degree
     j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
     xnipjp1(ij) = -taui_qi * (j_dp + 1._dp)
     xnipjm1(ij) = -taui_qi * j_dp
   ENDDO
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !! ES linear coefficients for moment RHS !!!!!!!!!!
   DO ip = ips_i, ipe_i
     p_int= parray_i(ip)   ! Hermite degree
     DO ij = ijs_i, ije_i
       j_int= jarray_i(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(ip,ij)    =+eta_n + 2.*j_dp*eta_T
         xphijp1(ip,ij)  =-eta_T*(j_dp+1._dp)
         xphijm1(ip,ij)  =-eta_T* j_dp
       ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
         xphij(ip,ij)    =+eta_T/SQRT2
         xphijp1(ip,ij)  = 0._dp; xphijm1(ip,ij)  = 0._dp;
       ELSE
         xphij(ip,ij)    = 0._dp; xphijp1(ip,ij)  = 0._dp
         xphijm1(ip,ij)  = 0._dp;
       ENDIF
     ENDDO
   ENDDO
 
 END SUBROUTINE compute_lin_coeff
 
+!******************************************************************************!
+!!!!!!! Remove all ky!=0 modes to conserve only zonal modes in a restart
+!******************************************************************************!
+SUBROUTINE wipe_turbulence
+  USE fields
+  USE grid
+  IMPLICIT NONE
+  DO ikx=ikxs,ikxe
+  DO iky=ikys,ikye
+  DO iz=izs,ize
+    DO ip=ips_e,ipe_e
+    DO ij=ijs_e,ije_e
+      IF( iky .NE. iky_0) THEN
+        moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
+      ELSE
+        moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
+      ENDIF
+    ENDDO
+    ENDDO
+    DO ip=ips_i,ipe_i
+    DO ij=ijs_i,ije_i
+      IF( iky .NE. iky_0) THEN
+        moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
+      ELSE
+        moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
+      ENDIF
+    ENDDO
+    ENDDO
+    IF( iky .NE. iky_0) THEN
+      phi(ikx,iky,iz) = 0e-3_dp*phi(ikx,iky,iz)
+    ELSE
+      phi(ikx,iky,iz) = 1e+0_dp*phi(ikx,iky,iz)
+    ENDIF
+  ENDDO
+  ENDDO
+  ENDDO
+END SUBROUTINE
+!******************************************************************************!
+!!!!!!! Remove all ky==0 modes to conserve only non zonal modes in a restart
+!******************************************************************************!
+SUBROUTINE wipe_zonalflow
+  USE fields
+  USE grid
+  IMPLICIT NONE
+  DO ikx=ikxs,ikxe
+  DO iky=ikys,ikye
+  DO iz=izs,ize
+    DO ip=ips_e,ipe_e
+    DO ij=ijs_e,ije_e
+      IF( iky .EQ. iky_0) THEN
+        moments_e( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_e( ip,ij,ikx,iky,iz, :)
+      ELSE
+        moments_e( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_e( ip,ij,ikx,iky,iz, :)
+      ENDIF
+    ENDDO
+    ENDDO
+    DO ip=ips_i,ipe_i
+    DO ij=ijs_i,ije_i
+      IF( iky .EQ. iky_0) THEN
+        moments_i( ip,ij,ikx,iky,iz, :) = 0e-3_dp*moments_i( ip,ij,ikx,iky,iz, :)
+      ELSE
+        moments_i( ip,ij,ikx,iky,iz, :) = 1e+0_dp*moments_i( ip,ij,ikx,iky,iz, :)
+      ENDIF
+    ENDDO
+    ENDDO
+    IF( iky .EQ. iky_0) THEN
+      phi(ikx,iky,iz) = 0e-3_dp*phi(ikx,iky,iz)
+    ELSE
+      phi(ikx,iky,iz) = 1e+0_dp*phi(ikx,iky,iz)
+    ENDIF
+  ENDDO
+  ENDDO
+  ENDDO
+END SUBROUTINE
 
 END MODULE numerics
diff --git a/src/stepon.F90 b/src/stepon.F90
index 882e127..b7a3678 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -1,142 +1,147 @@
 SUBROUTINE stepon
   !   Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
   USE advance_field_routine, ONLY: advance_time_level, advance_field, advance_moments
   USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj
   USE basic
   USE closure
   USE collision, ONLY : compute_TColl
   USE fields, ONLY: moments_e, moments_i, phi
+  USE initial_par, ONLY: WIPE_ZF, WIPE_TURB
   USE ghosts
   USE grid
   USE model
   use prec_const
   USE time_integration
+  USE numerics, ONLY: wipe_zonalflow, wipe_turbulence
   USE utility, ONLY: checkfield
 
   IMPLICIT NONE
 
   INTEGER :: num_step
   LOGICAL :: mlend
 
    DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
    !----- BEFORE: All fields are updated for step = n
       ! Compute right hand side from current fields
       ! N_rhs(N_n,phi_n, S_n, Tcoll_n)
       CALL moments_eq_rhs_e
       CALL moments_eq_rhs_i
       ! ---- step n -> n+1 transition
       ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
       CALL advance_time_level
       ! Update moments with the hierarchy RHS (step by step)
       ! N_n+1 = N_n + N_rhs(n)
       CALL advance_moments
       ! Closure enforcement of N_n+1
       CALL apply_closure_model
       ! Exchanges the ghosts values of N_n+1
       CALL update_ghosts
       ! Update collision C_n+1 = C(N_n+1)
       CALL compute_TColl
       ! Update electrostatic potential phi_n = phi(N_n+1)
       CALL poisson
       ! Update nonlinear term S_n -> S_n+1(phi_n+1,N_n+1)
-      IF ( NON_LIN ) CALL compute_Sapj
-
+      IF ( NON_LIN )         CALL compute_Sapj
+      ! Cancel zonal modes artificially
+      IF ( WIPE_ZF .EQ. 2)   CALL wipe_zonalflow
+      ! Cancel non zonal modes artificially
+      IF ( WIPE_TURB .EQ. 2) CALL wipe_turbulence
       !-  Check before next step
       CALL checkfield_all()
       IF( nlend ) EXIT ! exit do loop
 
       CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
    !----- AFTER: All fields are updated for step = n+1
    END DO
 
    CONTAINS
 
       SUBROUTINE checkfield_all ! Check all the fields for inf or nan
         ! Execution time start
         CALL cpu_time(t0_checkfield)
 
         IF(NON_LIN) CALL anti_aliasing   ! ensure 0 mode for 2/3 rule
         IF(NON_LIN) CALL enforce_symetry ! Enforcing symmetry on kx = 0
 
         mlend=.FALSE.
         IF(.NOT.nlend) THEN
            mlend=mlend .or. checkfield(phi,' phi')
            DO ip=ips_e,ipe_e
              DO ij=ijs_e,ije_e
               mlend=mlend .or. checkfield(moments_e(ip,ij,:,:,:,updatetlevel),' moments_e')
              ENDDO
            ENDDO
            DO ip=ips_i,ipe_i
              DO ij=ijs_i,ije_i
               mlend=mlend .or. checkfield(moments_i(ip,ij,:,:,:,updatetlevel),' moments_i')
              ENDDO
            ENDDO
            CALL MPI_ALLREDUCE(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
         ENDIF
 
         ! Execution time end
         CALL cpu_time(t1_checkfield)
         tc_checkfield = tc_checkfield + (t1_checkfield - t0_checkfield)
       END SUBROUTINE checkfield_all
 
       SUBROUTINE anti_aliasing
         DO ip=ips_e,ipe_e
           DO ij=ijs_e,ije_e
             DO ikx=ikxs,ikxe
               DO iky=ikys,ikye
                 DO iz=izs,ize
                   moments_e( ip,ij,ikx,iky,iz,:) = AA_x(ikx)* AA_y(iky) * moments_e( ip,ij,ikx,iky,iz,:)
                 END DO
               END DO
             END DO
           END DO
         END DO
         DO ip=ips_i,ipe_i
           DO ij=ijs_i,ije_i
             DO ikx=ikxs,ikxe
               DO iky=ikys,ikye
                 DO iz=izs,ize
                   moments_i( ip,ij,ikx,iky,iz,:) = AA_x(ikx)* AA_y(iky) * moments_i( ip,ij,ikx,iky,iz,:)
                 END DO
               END DO
             END DO
           END DO
         END DO
       END SUBROUTINE anti_aliasing
 
       SUBROUTINE enforce_symetry ! Force X(k) = X(N-k)* complex conjugate symmetry
         IF ( contains_kx0 ) THEN
           ! Electron moments
           DO ip=ips_e,ipe_e
             DO ij=ijs_e,ije_e
               DO iz=izs,ize
                 DO iky=2,Nky/2 !symmetry at kx = 0
                   moments_e( ip,ij,ikx_0,iky,iz, :) = CONJG(moments_e( ip,ij,ikx_0,Nky+2-iky,iz, :))
                 END DO
                 ! must be real at origin
                 moments_e(ip,ij, ikx_0,iky_0,iz, :) = REAL(moments_e(ip,ij, ikx_0,iky_0,iz, :))
               END DO
             END DO
           END DO
           ! Ion moments
           DO ip=ips_i,ipe_i
             DO ij=ijs_i,ije_i
               DO iz=izs,ize
                 DO iky=2,Nky/2 !symmetry at kx = 0
                   moments_i( ip,ij,ikx_0,iky,iz, :) = CONJG(moments_i( ip,ij,ikx_0,Nky+2-iky,iz, :))
                 END DO
                 ! must be real at origin and top right
                 moments_i(ip,ij, ikx_0,iky_0,iz, :) = REAL(moments_i(ip,ij, ikx_0,iky_0,iz, :))
               END DO
             END DO
           END DO
           ! Phi
           DO iky=2,Nky/2 !symmetry at kx = 0
             phi(ikx_0,iky,:) = phi(ikx_0,Nky+2-iky,:)
           END DO
           ! must be real at origin
           phi(ikx_0,iky_0,:) = REAL(phi(ikx_0,iky_0,:))
         ENDIF
       END SUBROUTINE enforce_symetry
 
 END SUBROUTINE stepon