diff --git a/src/inital.F90 b/src/inital.F90
index 73bc3cb5838e3feb71f70a50edcbce6a7ef13eee..91f05c860ed68319fbfc79fac939a9a714ad6dfd 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -197,21 +197,16 @@ SUBROUTINE init_phi
       !**** Cancel previous moments initialization
       moments_e = 0._dp; moments_i = 0._dp
 
-    ELSEIF(INIT_ZF .GT. 0) THEN
-
-      !**** Zonal Flow initialization *******************************************
-      ! Every mode is zero
-      phi       = 0._dp
-      moments_e = 0._dp
-      moments_i = 0._dp
-      ! Except at ikr = mode number + 1, symmetry is already included since kr>=0
-      IF( (INIT_ZF+1 .GT. ikrs) .AND. (INIT_ZF+1 .LT. ikre) ) THEN
-        phi(INIT_ZF+1,ikz_0) = (2._dp*PI)**2*init_background/deltakr/deltakz/2._dp
+      IF(INIT_ZF .GT. 0) THEN
+
+        !**** Zonal Flow initialization *******************************************
+        ! put a mode at ikr = mode number + 1, symmetry is already included since kr>=0
+        IF( (INIT_ZF+1 .GT. ikrs) .AND. (INIT_ZF+1 .LT. ikre) ) THEN
+          phi(INIT_ZF+1,ikz_0) = ZF_AMP*(2._dp*PI)**2/deltakr/deltakz/2._dp
+          moments_i(1,1,INIT_ZF+1,ikz_0,:) = krarray(INIT_ZF+1)**2*phi(INIT_ZF+1,ikz_0)
+          moments_e(1,1,INIT_ZF+1,ikz_0,:) = 0._dp
+        ENDIF
       ENDIF
-      !**** Init ZF in density
-      moments_e(1,1,INIT_ZF+1,ikz_0,:) = -(2._dp*PI)**2*init_background/deltakr/deltakz * imagu/2._dp
-      moments_i(1,1,INIT_ZF+1,ikz_0,:) = -(2._dp*PI)**2*init_background/deltakr/deltakz * imagu/2._dp
-
     ELSE ! we compute phi from noisy moments and poisson
 
       CALL poisson
diff --git a/src/initial_par_mod.F90 b/src/initial_par_mod.F90
index 56a94133f86627138a559e8ab2b7c6c6404694be..1c619aa81ab10c1e414b85a5b053ebe4755603b9 100644
--- a/src/initial_par_mod.F90
+++ b/src/initial_par_mod.F90
@@ -9,6 +9,7 @@ MODULE initial_par
   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
   ! Initial background level
   REAL(dp), PUBLIC, PROTECTED :: init_background=0._dp
   ! Initial noise amplitude
@@ -22,9 +23,6 @@ MODULE initial_par
   REAL(dp), PUBLIC, PROTECTED ::  init_ampli_density=0.1_dp ! Oscillation amplitude
   REAL(dp), PUBLIC, PROTECTED ::  init_ampli_temp=0.1_dp
 
-  CHARACTER(len=128), PUBLIC :: selfmat_file  ! COSOlver matrix file names
-  CHARACTER(len=128), PUBLIC :: iemat_file  ! COSOlver matrix file names
-  CHARACTER(len=128), PUBLIC :: eimat_file  ! COSOlver matrix file names
   CHARACTER(len=128), PUBLIC :: mat_file    ! COSOlver matrix file names
 
   PUBLIC :: initial_outputinputs, initial_readinputs
@@ -44,9 +42,6 @@ CONTAINS
     NAMELIST /INITIAL_CON/ init_background
     NAMELIST /INITIAL_CON/ init_noiselvl
     NAMELIST /INITIAL_CON/ iseed
-    NAMELIST /INITIAL_CON/ selfmat_file
-    NAMELIST /INITIAL_CON/ iemat_file
-    NAMELIST /INITIAL_CON/ eimat_file
     NAMELIST /INITIAL_CON/ mat_file
 
     READ(lu_in,initial_con)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 83f6cabaf2ed4294fbadd58c599c067cdaee8e9b..f06cf81828b15d97eb1d8079c2920dfc22786776 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -37,6 +37,8 @@ MODULE model
   REAL(dp), PUBLIC, PROTECTED :: nu_ee, nu_ie         ! e-e, i-e coll. frequ.
   REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui   ! factor of the gammaD sum
 
+  REAL(dp), PUBLIC, PROTECTED :: INV_LIN_SYS = 1._dp  ! Invert the sign of the linear RHS (for measuring damped eigenvalues)
+
   PUBLIC :: model_readinputs, model_outputinputs
 
 CONTAINS
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index deeca00f82269be414ca222c308fde0d9d3397f7..d56a0bfeff849ff9878843f018c8d7fbb38a195c 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -162,10 +162,10 @@ SUBROUTINE moments_eq_rhs_e
 
           !! Sum of all linear terms
           moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
-              -i_kz  * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
+              - INV_LIN_SYS * i_kz  * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
               - mu*kperp2**2 * moments_e(ip,ij,ikr,ikz,updatetlevel) &
               + mu_p * Hyper_diff_p + mu_j * Hyper_diff_j &
-              + TColl
+              + INV_LIN_SYS * TColl
 
           !! Adding non linearity
           IF ( NON_LIN ) THEN
@@ -347,10 +347,10 @@ SUBROUTINE moments_eq_rhs_i
 
           !! Sum of linear terms
           moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
-              -i_kz  * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
+              -INV_LIN_SYS * i_kz  * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
               - mu*kperp2**2 * moments_i(ip,ij,ikr,ikz,updatetlevel) &
               + mu_p * Hyper_diff_p + mu_j * Hyper_diff_j &
-              + TColl
+              +INV_LIN_SYS * TColl
 
           !! Adding non linearity
           IF ( NON_LIN ) THEN