diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 56f1217c77c1a6c268e2d7acddf6362ee4d15b35..f8785a89748999b4592f211238ef45a106faa8f1 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -266,7 +266,7 @@ END SUBROUTINE fft1D_plans
     !   module variable (convolution theorem)
     SUBROUTINE poisson_bracket_and_sum( ky_array, kx_array, inv_Ny, inv_Nx, AA_y, AA_x,&
                                         local_nky, total_nkx, F_, G_,&
-                                        ExB, ExB_NL_factor,sky_ExB,sum_real_)
+                                        ExB_NL_CORRECTION, ExB_NL_factor,sky_ExB,sum_real_)
         IMPLICIT NONE
         INTEGER,                                  INTENT(IN) :: local_nky,total_nkx
         REAL(xp),                                 INTENT(IN) :: inv_Nx, inv_Ny
@@ -277,7 +277,7 @@ END SUBROUTINE fft1D_plans
                                                   INTENT(IN) :: F_, G_
         COMPLEX(xp), DIMENSION(total_nkx,local_nky), &
                                                   INTENT(IN) :: ExB_NL_factor
-        LOGICAL, INTENT(IN) :: ExB
+        LOGICAL, INTENT(IN) :: ExB_NL_CORRECTION
         REAL(xp),DIMENSION(local_nky),            INTENT(IN) :: sky_ExB
         real(c_xp_r), pointer,                 INTENT(INOUT) :: sum_real_(:,:)
         ! local variables
@@ -297,7 +297,7 @@ END SUBROUTINE fft1D_plans
                         ikxG(ikx,iky) = imagu*kxs*G_(iky,ikx)
                 ENDDO
         ENDDO
-        IF(ExB) THEN 
+        IF(ExB_NL_CORRECTION) THEN 
             ! Apply the ExB shear correction factor exp(ixkySJdT)
             CALL apply_ExB_NL_factor(ikxF,ExB_NL_factor)
             CALL apply_ExB_NL_factor(ikyG,ExB_NL_factor)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index f468a0cd83091b5b34f508c060c12223b5b1e2cc..4447969b5cd7a93a6dd60371d3bbd8dab0fb0ca6 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -35,6 +35,8 @@ MODULE model
   LOGICAL,  PUBLIC, PROTECTED :: RM_LD_T_EQ = .false.
   ! Flag to force the reality condition symmetry for the kx at ky=0
   LOGICAL,  PUBLIC, PROTECTED :: FORCE_SYMMETRY = .false.
+  ! Add or remove the ExB nonlinear correction (Mcmillan 2019)
+  LOGICAL,  PUBLIC, PROTECTED :: ExB_NL_CORRECTION = .true.
 
   ! Module's routines
   PUBLIC :: model_readinputs, model_outputinputs
@@ -51,7 +53,8 @@ CONTAINS
     NAMELIST /MODEL/ KERN, LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, &
                          Na, ADIAB_E, ADIAB_I, tau_i, &
                          mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, &
-                         nu, k_gB, k_cB, lambdaD, beta, ExBrate, ikxZF, ZFamp
+                         nu, k_gB, k_cB, lambdaD, beta, ExBrate, ExB_NL_CORRECTION,&
+                         ikxZF, ZFamp
 
     READ(lu_in,model)
 
@@ -78,6 +81,9 @@ CONTAINS
       EM = .FALSE.
     ENDIF
 
+    IF(abs(ExBrate) .LT. epsilon(ExBrate))&
+      ExB_NL_CORRECTION = .false.
+
   END SUBROUTINE model_readinputs
 
   SUBROUTINE model_outputinputs(fid)
@@ -106,6 +112,7 @@ CONTAINS
     CALL attach(fid, TRIM(str),   "lambdaD", lambdaD)
     CALL attach(fid, TRIM(str),    "MHD_PD",  MHD_PD)
     CALL attach(fid, TRIM(str),      "beta",    beta)
+    CALL attach(fid, TRIM(str),   "ExBrate", ExBrate)
     CALL attach(fid, TRIM(str),   "ADIAB_E", ADIAB_E)
     CALL attach(fid, TRIM(str),   "ADIAB_I", ADIAB_I)
     CALL attach(fid, TRIM(str),     "tau_i",   tau_i)
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 2440078e749ff0ea70fbdc10176ac687f67e2eda..38dd61bd32eead22cb775ac6e36ae021aacd7732 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -9,12 +9,12 @@ MODULE nonlinear
                          kyarray, AA_y, local_nky, inv_Ny,&
                          total_nkx,kxarray, AA_x, inv_Nx,&
                          local_nz,ngz,zarray,nzgrid, deltakx, iky0, contains_kx0, contains_ky0
-  USE model,       ONLY : LINEARITY, EM, ikxZF, ZFamp
+  USE model,       ONLY : LINEARITY, EM, ikxZF, ZFamp, ExB_NL_CORRECTION
   USE closure,     ONLY : evolve_mom, nmaxarray
   USE prec_const,  ONLY : xp
   USE species,     ONLY : sqrt_tau_o_sigma
   USE time_integration, ONLY : updatetlevel
-  USE ExB_shear_flow,   ONLY : ExB_NL_factor, inv_ExB_NL_factor, ExB, sky_ExB
+  USE ExB_shear_flow,   ONLY : ExB_NL_factor, inv_ExB_NL_factor, sky_ExB
   use, intrinsic :: iso_c_binding
 
   IMPLICIT NONE
@@ -100,7 +100,7 @@ SUBROUTINE compute_nonlinear
               ! this function adds its result to bracket_sum_r
                 CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
                                               local_nky,total_nkx,F_cmpx,G_cmpx,&
-                                              ExB, ExB_NL_factor, sky_ExB, bracket_sum_r)
+                                              ExB_NL_CORRECTION, ExB_NL_factor, sky_ExB, bracket_sum_r)
   !-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
               IF(EM) THEN
                 ! First convolution terms
@@ -116,12 +116,11 @@ SUBROUTINE compute_nonlinear
                 ! this function adds its result to bracket_sum_r
                 CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
                                               local_nky,total_nkx,F_cmpx,G_cmpx,&
-                                              ExB, ExB_NL_factor, sky_ExB, bracket_sum_r)
+                                              ExB_NL_CORRECTION, ExB_NL_factor, sky_ExB, bracket_sum_r)
               ENDIF
             ENDDO n
             ! Apply the ExB shearing rate factor before going back to k-space
-            IF (ExB) THEN
-              ! print*, SUM(bracket_sum_r)
+            IF (ExB_NL_CORRECTION) THEN
               CALL apply_inv_ExB_NL_factor(bracket_sum_r,inv_ExB_NL_factor)
             ENDIF
             ! Put the real nonlinear product back into k-space