diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index ecf557383e32e7c152ab778a4812ae57c3cc229c..b4c3baebecd9442b548b6a8c335ef842f515a153 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -94,15 +94,21 @@ while(CONTINUE)
             PHI_  = cat(4,PHI_,PHI);
         end
         if W_NA00
+            if KIN_E
+             Ne00_ = cat(4,Ne00_,Ne00);
+            end
             Ni00_ = cat(4,Ni00_,Ni00);
-            Ne00_ = cat(4,Ne00_,Ne00);
         end
         if W_DENS
+            if KIN_E
             DENS_E_ = cat(4,DENS_E_,DENS_E);
+            end
             DENS_I_ = cat(4,DENS_I_,DENS_I);
         end
         if W_TEMP
-            TEMP_E_ = cat(4,TEMP_E_,TEMP_E);
+            if KIN_E 
+                TEMP_E_ = cat(4,TEMP_E_,TEMP_E);
+            end
             TEMP_I_ = cat(4,TEMP_I_,TEMP_I);
         end
         if W_NAPJ || W_SAPJ
@@ -110,11 +116,15 @@ while(CONTINUE)
         end
         if W_NAPJ
             Nipj_ = cat(6,Nipj_,Nipj);
-            Nepj_ = cat(6,Nepj_,Nepj);
+            if KIN_E
+                Nepj_ = cat(6,Nepj_,Nepj);
+            end
         end
         if W_SAPJ
      	  Sipj_ = cat(6,Sipj_,Sipj);
-          Sepj_ = cat(6,Sepj_,Sepj);
+          if KIN_E
+           Sepj_ = cat(6,Sepj_,Sepj);
+          end
         end
 
         % Evolution of simulation parameters
diff --git a/matlab/load_results.m b/matlab/load_results.m
index ad16253186baa5d97f83f08c35f7b7f72b5911c2..39fd519f4df8be55b49ca3078a9094d6e0fb2dc8 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -14,8 +14,8 @@ W_NAPJ    = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
 W_SAPJ    = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
 W_DENS    = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
 W_TEMP    = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
-% KIN_E     = strcmp(h5readatt(filename,'/data/input',     'KIN_E' ),'y');
-KIN_E     = 1;
+KIN_E     = strcmp(h5readatt(filename,'/data/input',     'KIN_E' ),'y');
+% KIN_E     = 1;
 
 
 if W_GAMMA
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
index 1cccc49ff54b7d090c4c3ab127f40c78b5b99927..f30b11dfbb6d6439a1f00982570a14f139f3cc44 100644
--- a/matlab/post_processing.m
+++ b/matlab/post_processing.m
@@ -52,8 +52,22 @@ dx2phi = zeros(Nx,Ny,Nz,Ns3D); % "
 
 for it = 1:numel(Ts3D)
     for iz = 1:numel(z)
-        NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
-        ne00  (:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
+        if KIN_E
+            NE_ = Ne00(:,:,iz,it); 
+            ne00  (:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
+            if(W_DENS)
+            DENS_E_ = DENS_E(:,:,iz,it);
+            dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny)));
+            Z_n_e  (:,:,iz,it) = real(fftshift(ifft2((DENS_E_.*(KY==0)),Nx,Ny)));
+            end
+            if(W_TEMP)
+            TEMP_E_ = TEMP_E(:,:,iz,it);
+            dyTe(:,:,iz,it)  = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny)));
+            temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny)));
+            Z_T_e  (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_.*(KY==0)),Nx,Ny)));
+            end
+        end
+        NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
         ni00  (:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny)));
         phi   (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny)));
         Z_phi (:,:,iz,it) = real(fftshift(ifft2((PH_.*(KY==0)),Nx,Ny)));
@@ -61,20 +75,15 @@ for it = 1:numel(Ts3D)
         dx2phi(:,:,iz,it) = real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny)));
         dyphi (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny)));
         if(W_DENS)
-        DENS_E_ = DENS_E(:,:,iz,it); DENS_I_ = DENS_I(:,:,iz,it);
+        DENS_I_ = DENS_I(:,:,iz,it);
         dyni   (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny)));
-        dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny)));
         dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny)));
-        Z_n_e  (:,:,iz,it) = real(fftshift(ifft2((DENS_E_.*(KY==0)),Nx,Ny)));
         Z_n_i  (:,:,iz,it) = real(fftshift(ifft2((DENS_I_.*(KY==0)),Nx,Ny)));
         end
         if(W_TEMP)
-        TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it);
-        dyTe(:,:,iz,it)  = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny)));
+        TEMP_I_ = TEMP_I(:,:,iz,it);
         dyTi(:,:,iz,it)  = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny)));
-        temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny)));
         temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny)));
-        Z_T_e  (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_.*(KY==0)),Nx,Ny)));
         Z_T_i  (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_.*(KY==0)),Nx,Ny)));
         end
     end
@@ -119,7 +128,6 @@ phi_kp_t = zeros(numel(kp_ip),Nz,Ns3D);
 %
 for it = 1:numel(Ts3D) % Loop over 2D aX_XYays
     for iz = 1:numel(z)
-    NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
     phi_maxx_maxy(iz,it)   =  max( max(squeeze(phi(:,:,iz,it))));
     phi_avgx_maxy(iz,it)   =  max(mean(squeeze(phi(:,:,iz,it))));
     phi_maxx_avgy(iz,it)   = mean( max(squeeze(phi(:,:,iz,it))));
@@ -147,7 +155,9 @@ end
 %
 for it = 1:numel(Ts5D) % Loop over 5D aX_XYays
     [~, it2D] = min(abs(Ts3D-Ts5D(it)));
-    Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky;
+    if KIN_E
+        Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky;
+    end
     Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky;
 end
 
diff --git a/matlab/setup.m b/matlab/setup.m
index a9a1f15ba5d1f0e8f5468264b2ea9bf42c07db47..aa8683ce69ffce482afbf7f1c38030f292c3da9b 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -13,7 +13,7 @@ GRID.Nz    = NZ;    % z resolution
 GRID.q0    = Q0;    % q factor
 GRID.shear = SHEAR; % shear
 GRID.eps   = EPS;   % inverse aspect ratio
-
+if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end;
 % Model parameters
 MODEL.CO      = CO;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 MODEL.CLOS    = CLOS;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 84259b0f77cca3641911abaf13592dc1cf25f820..58092ab4055d785053046bb7c0635d86ecbe3dd6 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -23,6 +23,7 @@ fprintf(fid,['  Nz     = ', num2str(GRID.Nz),'\n']);
 fprintf(fid,['  q0     = ', num2str(GRID.q0),'\n']);
 fprintf(fid,['  shear  = ', num2str(GRID.shear),'\n']);
 fprintf(fid,['  eps    = ', num2str(GRID.eps),'\n']);
+fprintf(fid,['  SG     = ',           GRID.SG,'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&OUTPUT_PAR\n');
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index eff2e01673390d6acd0e33b22cbfbb27d21e645c..5515dca97f0dad23d250e825bfa81bf437c754be 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -4,8 +4,8 @@ outfile ='';
 %% Directory of the simulation
 if 1% Local results
 outfile ='';
-outfile ='';
-outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_lin';
+outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_adiabe';
+% outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK_lin';
 % outfile ='fluxtube_salphaB_s0/64x64x16_5x3_L_200_q0_2.7_e_0.18_kN_2.22_kT_6_nu_1e-01_DGGK';
 % outfile ='simulation_A/1024x256_3x2_L_120_kN_1.6667_nu_1e-01_DGGK';
 % outfile ='Linear_Device/64x64x20_5x2_Lx_20_Ly_150_q0_25_kN_0.24_kT_0.03_nu_1e-02_DGGK';
@@ -15,16 +15,15 @@ outfile ='fluxtube_salphaB_s0/50x100x20_5x3_L_300_q0_2.7_e_0.18_kN_2.22_kT_6_nu_
     CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
     system(CMD);
 else% Marconi results
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x300_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/500x500_L_120_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x200_L_200_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x300_L_120_P_8_J_4_eta_0.6_nu_1e-01_PAGK_mu_0e+00/out.txt';
 BASIC.RESDIR      = ['../',outfile(46:end-8),'/'];
 BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
 end
 
 %% Load the results
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20; 
+JOBNUMMIN = 00; JOBNUMMAX = 00; 
 % JOBNUMMIN = 07; JOBNUMMAX = 20; % For CO damping sim A 
 compile_results %Compile the results from first output found to JOBNUMMAX if existing
 
@@ -46,7 +45,7 @@ end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-t0    =000; iz = 2; ix = 1; iy = 1;
+t0    =000; iz = 1; ix = 1; iy = 1;
 skip_ =1; FPS = 30; DELAY = 1/FPS;
 [~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
 [~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
@@ -64,12 +63,12 @@ FIELD = phi;          NAME = 'phi'; FIELDNAME = '\phi';
 % FIELD = Gamma_x;      NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
 
 % Sliced
-plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; 
+% plt = @(x) real(x(ix, :, :,:)); X = Y_YZ; Y = Z_YZ; XNAME = 'y'; YNAME = 'z'; 
 % plt = @(x) real(x( :,iy, :,:)); X = X_XZ; Y = Z_XZ; XNAME = 'x'; YNAME = 'z'; 
-% plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; 
+plt = @(x) real(x( :, :,iz,:)); X = X_XY; Y = Y_XY; XNAME = 'x'; YNAME = 'y'; 
 
 % K-space
-% FIELD = PHI;          NAME = 'PHI'; FIELDNAME = '\tilde \phi';
+% FIELD = PHI.*(KY~=0);          NAME = 'NZPHI'; FIELDNAME = '\tilde \phi_{k_y\neq0}';
 % FIELD = Ne00;         NAME = 'Ne00'; FIELDNAME = 'N_e^{00}';
 % FIELD = Ni00;         NAME = 'Ni00'; FIELDNAME = 'N_i^{00}';
 % plt = @(x) fftshift((abs(x( :, :,1,:))),2); X = fftshift(KX,2); Y = fftshift(KY,2); XNAME = 'k_x'; YNAME = 'k_y'; 
@@ -145,13 +144,14 @@ if 0
 
 % Chose the field to plot
 % FIELD = Ni00;   FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-% FIELD = Ne00;   FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+FIELD = Ne00;   FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
 % FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
 % FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
 
 % Chose when to plot it
-tf = [0 15 27 28 30];
+% tf = [0 15 27 28 30];
+tf = 200:200:1200;
 % tf = 8000;
 
 % Planar plot: choose a plane to plot at x0/y0/z0 coordinates
@@ -178,7 +178,7 @@ fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position',
     subplot(1,numel(tf),i_)
         DATA = plt_2(squeeze(plt(FIELD,it)));
         pclr = pcolor(plt_x(X),plt_y(Y),plt_z(DATA)); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1])
-%         caxis([0 1]*5e3);
+        caxis([0 1]*1e3);
 %         caxis([-1 1]*5);
         xlabel(latexize(XNAME)); ylabel(latexize(YNAME));
         if(i_ > 1); set(gca,'ytick',[]); end;
@@ -217,7 +217,8 @@ end
 
 if 0
 %% zonal vs nonzonal energies for phi(t)
-it0 = 01; itend = Ns3D;
+t0 = 200;  [~, it0] = min(abs(Ts3D-t0)); 
+itend = Ns3D;
 trange = it0:itend;
 pltx = @(x) x;%-x(1);
 plty = @(x) x./max(squeeze(x));
@@ -237,9 +238,10 @@ subplot(121)
     xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
 
 subplot(122)
-    plot(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)));
+%     plot(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)));
+    plot3(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)),Ts3D(trange));
     title('Phase space'); legend(CONAME)
-    xlabel('$E_Z$'); ylabel('$E_{NZ}$'); grid on;% xlim([0 500]);
+    xlabel('$E_Z$'); ylabel('$E_{NZ}$'); zlabel('time'); grid on;% xlim([0 500]);
 end
 
 if 0
@@ -281,7 +283,7 @@ FIELD = Ne00.*conj(Ne00);   FNAME = 'Ne002'; FIELDLTX = '|N_e^{00}|^2'
 % FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
 
 % Chose when to plot it
-tf = 1500:200:2500;
+tf = 200:200:1200;
 % tf = 8000;
 
 % Sliced
@@ -308,7 +310,7 @@ end
 if 0
 %% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
 % tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
-tstart = 1000; tend = 4000;
+tstart = 0000; tend = 1000;
 % tstart = 10000; tend = 12000;
 % Chose the field to plot
 % FIELD = Ni00;   FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
diff --git a/wk/benchmark_HeLaZ_molix.m b/wk/benchmark_HeLaZ_molix.m
index a7394a7eee3c29452b6e46c88f4ad29a383df702..bb6da7ab40c14a323db5dcb13abca47c5e88e306 100644
--- a/wk/benchmark_HeLaZ_molix.m
+++ b/wk/benchmark_HeLaZ_molix.m
@@ -17,7 +17,8 @@ system(['cd ',SIMDIR,';',' ./../../../bin/helaz 0; cd ../../../wk'])
 % Compare the results with molix at a given time
 time_2_plot = 5.0;
 [Y_,X_]=molix_plot_phi([SIMDIR,'molix_phi.h5'],time_2_plot);
-[ PHI, Ts3D, dt3D] = load_3D_data([SIMDIR,'outputs_00.h5'], 'phi');
+filename = [SIMDIR,'/outputs_00.h5'];
+[ PHI, Ts3D, dt3D] = load_3D_data(filename, 'phi');
 [Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
 
 plot_phi_ballooning; hold on
diff --git a/wk/benchmark_linear_1D_entropy_mode.m b/wk/benchmark_linear_1D_entropy_mode.m
index ae422643468cb28627c60acb7fbd2b8064b6c275..86528abce82fb8cbd6474d09e5d158b0d48e3761 100644
--- a/wk/benchmark_linear_1D_entropy_mode.m
+++ b/wk/benchmark_linear_1D_entropy_mode.m
@@ -5,8 +5,12 @@ SIMDIR = '../results/benchmarks/1D_linear_entropy_mode/';
 cd ..
 system('make');
 cd wk
-% Run HeLaZ
-system(['cd ',SIMDIR,';',' ./../../../bin/helaz 0; cd ../../../wk'])
+% Run HeLaZ in sequential
+system(['cd ',SIMDIR,';',' ./../../../bin/helaz 0; cd ../../../wk']);
+% Run HeLaZ in parallel (discrepencies can occur at marginal growth rate)
+% since the random seed will not be the same)
+% system(['cd ',SIMDIR,';',' mpirun -np 6 ./../../../bin/helaz 1 6 0; cd ../../../wk'])
+% system(['cd ',SIMDIR,';',' mpirun -np 6 ./../../../bin/helaz 2 3 0; cd ../../../wk'])
 %compute growth rate
 %%
 filename = [SIMDIR,'/outputs_00.h5'];
diff --git a/wk/local_run.m b/wk/local_run.m
index 0ed013d84270b3b9311354d6fc521ad08bddf5d8..10fefe6e5f12d1fbbecaa1ec4e92e4cc0b29d18e 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -6,11 +6,11 @@ CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %% PHYSICAL PARAMETERS
 NU      = 0.1;   % Collision frequency
 K_N     = 2.22;      % Density gradient drive
-K_T     = 6.0;       % Temperature '''
+K_T     = 8.0;       % Temperature '''
 K_E     = 0.00;    % Electrostat gradient
 SIGMA_E = 0.05196;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 NU_HYP  = 0.0;
-KIN_E   = 1;         % Kinetic (1) or adiabatic (2) electron model
+KIN_E   = 0;         % Kinetic (1) or adiabatic (2) electron model
 %% GRID PARAMETERS
 NX      = 50;     % Spatial radial resolution ( = 2x radial modes)
 LX      = 300;    % Radial window size
@@ -23,11 +23,12 @@ J       = 2;
 Q0      = 2.7;       % safety factor
 SHEAR   = 0.0;       % magnetic shear
 EPS     = 0.18;      % inverse aspect ratio
-GRADB   = 1.0;   % Magnetic  gradient
-CURVB   = 1.0;   % Magnetic  curvature
+GRADB   = 1.0;       % Magnetic  gradient
+CURVB   = 1.0;       % Magnetic  curvature
+SG      = 1;         % Staggered z grids option
 %% TIME PARAMETERS
-TMAX    = 10;  % Maximal time unit
-DT      = 5e-3;   % Time step
+TMAX    = 50;  % Maximal time unit
+DT      = 2e-2;   % Time step
 SPS0D   = 1;      % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS3D   = 5;      % Sampling per time unit for 3D arrays
@@ -39,13 +40,13 @@ JOB2LOAD= -1;
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK)
 CO      = 1;
 CLOS    = 0;   % Closure model (0: =0 truncation)
-NL_CLOS = 0;   % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+NL_CLOS = -1;   % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
 SIMID   = 'fluxtube_salphaB_s0';  % Name of the simulation
 % SIMID   = 'simulation_A';  % Name of the simulation
 % SIMID   = ['v3.0_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
-NON_LIN = 0;   % activate non-linearity (is cancelled if KXEQ0 = 1)
+NON_LIN = 1;   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % INIT options
-INIT_PHI= 1;   % Start simulation with a noisy phi (0= noisy moments 00)
+INIT_PHI= 0;   % Start simulation with a noisy phi (0= noisy moments 00)
 INIT_ZF   = 0; ZF_AMP = 0.0;
 INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
 %% OUTPUTS