diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 2645c259f08c9089c23f589adb4d40144dfcaac7..0445ae1bee607b7306c2752f9cf0499495be1f3d 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -43,7 +43,8 @@ MODULE grid
   REAL(dp), DIMENSION(:),   ALLOCATABLE, PUBLIC :: zarray_full
   ! local z weights for computing simpson rule
   INTEGER,  DIMENSION(:),   ALLOCATABLE, PUBLIC :: zweights_SR
-  REAL(dp), PUBLIC, PROTECTED  ::  deltax,  deltay, deltaz, inv_deltaz, diff_dz_coeff
+  REAL(dp), PUBLIC, PROTECTED  ::  deltax,  deltay, deltaz, inv_deltaz
+  REAL(dp), PUBLIC, PROTECTED  ::  diff_kx_coeff, diff_ky_coeff, diff_dz_coeff
   INTEGER,  PUBLIC, PROTECTED  ::  ixs,  ixe,  iys,  iye,  izs,  ize
   INTEGER,  PUBLIC, PROTECTED  ::  izgs, izge ! ghosts
   LOGICAL,  PUBLIC, PROTECTED  ::  SG = .true.! shifted grid flag
@@ -292,6 +293,7 @@ CONTAINS
       deltaky = 2._dp*PI/Ly
       ky_max  = Nky*deltaky
       ky_min  = deltaky
+      diff_ky_coeff= (1._dp/ky_max)**4
     ENDIF
     ! Build the full grids on process 0 to diagnose it without comm
     ALLOCATE(kyarray_full(1:Nky))
@@ -370,9 +372,10 @@ CONTAINS
       kxarray_full(1) = 0._dp
       local_kxmax     = 0._dp
     ELSE ! Build apprpopriate grid
-      deltakx     = 2._dp*PI/Lx
-      kx_max      = (Ny/2)*deltakx
-      kx_min      = deltakx
+      deltakx      = 2._dp*PI/Lx
+      kx_max       = (Nkx/2)*deltakx
+      kx_min       = deltakx
+      diff_kx_coeff= (1._dp/kx_max)**4
       ! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
       local_kxmax = 0._dp
       DO ikx = ikxs,ikxe
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 803ca3d23aa794eb2d767a89c3e825148987d320..3a493643875aa862251991ed6eed83655f2d8625 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -110,8 +110,9 @@ SUBROUTINE moments_eq_rhs_e
                 ! Drives (density + temperature gradients)
                 - i_ky * Tphi &
                 ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
-                ! - (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
-                - (mu_x*kx**2 + mu_y*ky**2)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
+                - mu_x*diff_kx_coeff*kx**4*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
+                - mu_y*diff_ky_coeff*ky**4*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
+                ! - (mu_x*kx**2 + mu_y*ky**2)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
                 ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"  see Pueschel 2010 (eq 25)
                 + mu_z * diff_dz_coeff * ddz4_Nepj(ip,ij,iky,ikx,iz) &
                 ! Collision term
@@ -241,8 +242,8 @@ SUBROUTINE moments_eq_rhs_i
                   ! Drives (density + temperature gradients)
                   - i_ky * Tphi &
                   ! Numerical hyperdiffusion (totally artificial, for stability purpose)
-                  ! - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
-                  - (mu_x*kx**2 + mu_y*ky**2)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
+                  - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
+                  ! - (mu_x*kx**2 + mu_y*ky**2)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
                   ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
                   + mu_z * diff_dz_coeff * ddz4_Nipj(ip,ij,iky,ikx,iz) &
                   ! Collision term