diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 3dca4e309433da22a41cbeae65d4477232297d87..8e921396915b92da861e5c717496e5362c283026 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -62,6 +62,7 @@ MODULE array
   REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_i
   ! Poisson operator (ikx,iky,iz)
   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
+  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_pol_ion
 
   ! Gyrocenter density for electron and ions (ikx,iky,iz)
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ne00
diff --git a/src/memory.F90 b/src/memory.F90
index c9509278094ed98800c478c2e1469bf47a0728db..afe8aebc7c988bd596db2c0fa62f1df176c3d91d 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -17,6 +17,7 @@ SUBROUTINE memory
   CALL allocate_array(        phi_ZF, ikxs,ikxe, izs,ize)
   CALL allocate_array(        phi_EM, ikys,ikye, izs,ize)
   CALL allocate_array(inv_poisson_op, ikys,ikye, ikxs,ikxe, izs,ize)
+  CALL allocate_array(   inv_pol_ion, ikys,ikye, ikxs,ikxe, izs,ize)
 
   !Electrons arrays
   IF(KIN_E) THEN
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 77d0d35749e853075968b23cdf9ac5d5f74ff357..7b2e743584f4e168bb84a40072d8ea22427345b0 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -94,7 +94,7 @@ END SUBROUTINE evaluate_kernels
 !******************************************************************************!
 SUBROUTINE evaluate_poisson_op
   USE basic
-  USE array, Only : kernel_e, kernel_i, inv_poisson_op
+  USE array, Only : kernel_e, kernel_i, inv_poisson_op, inv_pol_ion
   USE grid
   USE model, ONLY : tau_e, tau_i, q_e, q_i, KIN_E
   IMPLICIT NONE
@@ -125,9 +125,10 @@ SUBROUTINE evaluate_poisson_op
       END DO
     ! Adiabatic model
     ELSE
-      pol_e = 1._dp - qe2_taue
+      pol_e = qe2_taue - 1._dp
     ENDIF
     inv_poisson_op(iky, ikx, iz) =  1._dp/(qe2_taue + qi2_taui - pol_i - pol_e)
+    inv_pol_ion   (iky, ikx, iz) =  1._dp/(qi2_taui - pol_i)
   ENDIF
   END DO zloop
   END DO kyloop
diff --git a/src/poisson.F90 b/src/poisson.F90
index 75e2b887a670919ba094a6c80df19f307761d65f..49f1db404180e9eb9582e85b546087e72e14f498 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -15,9 +15,9 @@ SUBROUTINE poisson
   IMPLICIT NONE
 
   INTEGER     :: ini,ine, i_, root_bcast
-  COMPLEX(dp) :: fa_phi, intf_   ! current flux averaged phi
+  COMPLEX(dp) :: fsa_phi, intf_   ! current flux averaged phi
   INTEGER     :: count !! mpi integer to broadcast the electric potential at the end
-  COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e     ! charge density q_a n_a
+  COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e, integrant  ! charge density q_a n_a and aux var
 
   ! Execution time start
   CALL cpu_time(t0_poisson)
@@ -27,36 +27,35 @@ SUBROUTINE poisson
 
     kxloop: DO ikx = ikxs,ikxe
       kyloop: DO iky = ikys,ikye
-        !!!! Compute ion gyro charge density
+        !!!! Compute ion particle charge density q_i n_i
         rho_i = 0._dp
         DO ini=1,jmaxi+1
-          DO iz = izs,ize
-          rho_i(iz) = rho_i(iz) &
-           +q_i*kernel_i(ini,iky,ikx,iz,0)*moments_i(ip0_i,ini,iky,ikx,iz,updatetlevel)
-          ENDDO
+          rho_i(izs:ize) = rho_i(izs:ize) &
+           +q_i*kernel_i(ini,iky,ikx,izs:ize,0)*moments_i(ip0_i,ini,iky,ikx,izs:ize,updatetlevel)
         END DO
 
-        !!!! Compute electron gyro charge density
+        !!!! Compute electron particle charge density q_e n_e
         rho_e = 0._dp
         IF (KIN_E) THEN ! Kinetic electrons
           DO ine=1,jmaxe+1
-            DO iz = izs,ize
-            rho_e(iz) = rho_e(iz) &
-             +q_e*kernel_e(ine,iky,ikx,iz,0)*moments_e(ip0_e,ine,iky,ikx,iz,updatetlevel)
-            ENDDO
+            rho_e(izs:ize) = rho_e(izs:ize) &
+             +q_e*kernel_e(ine,iky,ikx,izs:ize,0)*moments_e(ip0_e,ine,iky,ikx,izs:ize,updatetlevel)
           END DO
         ELSE  ! Adiabatic electrons
-          ! Adiabatic charge density (linked to flux averaged phi)
-          fa_phi = 0._dp
-          IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (average)
-            DO ini=1,jmaxi+1 ! some over FLR contributions
-                rho_e(izs:ize) = Jacobian(izs:ize,0)*moments_i(ip0_i,ini,iky,ikx,izs:ize,updatetlevel)&
-                           *kernel_i(ini,iky,ikx,izs:ize,0)*(inv_poisson_op(iky,ikx,izs:ize)-1._dp)
-                call simpson_rule_z(rho_e(izs:ize),intf_) ! integrate over z
-                fa_phi = fa_phi + intf_
-            ENDDO
-            rho_e(izs:ize) = fa_phi*iInt_Jacobian !Normalize by 1/int(Jxyz)dz
+          ! Adiabatic charge density (linked to flux surface averaged phi)
+          fsa_phi = 0._dp
+          ! We compute the flux surface average solving a flux surface averaged
+          ! Poisson equation, i.e.
+          !
+          ! [qi^2(1-sum_j K_j^2)/tau_i] <phi>_psi = <q_i n_i >_psi
+          !       inv_pol_ion^-1         fsa_phi  = simpson(Jacobian rho_i ) * iInt_Jacobian
+          IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (y-average)
+            ! Prepare integrant for z-average
+            integrant(izs:ize) = Jacobian(izs:ize,0)*rho_i(izs:ize)*inv_pol_ion(iky,ikx,izs:ize)
+            call simpson_rule_z(integrant(izs:ize),intf_) ! get the flux averaged phi
+            fsa_phi = intf_
           ENDIF
+          rho_e(izs:ize) = fsa_phi * iInt_Jacobian !Normalize by 1/int(Jxyz)dz
         ENDIF
 
         !!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!