diff --git a/matlab/plot/plot_snapshot.m b/matlab/plot/plot_snapshot.m
new file mode 100644
index 0000000000000000000000000000000000000000..d8c897279205acf5902865fbce26de1173cd4d82
--- /dev/null
+++ b/matlab/plot/plot_snapshot.m
@@ -0,0 +1,35 @@
+dir = '/home/ahoffman/HeLaZ/results/CBC/NM_F4_kT_4.5_192x64x24x6x4/';
+fname = 'check_phi.out';
+
+filename = [dir,fname];
+startRow = 2;
+
+formatSpec = '%24f%2s%f%[^\n\r]';
+
+fileID = fopen(filename,'r');
+
+DIMS = fscanf (fileID,'%d %d %d', [3,1]); % indicate the programme in bar chart “header”
+
+dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false, 'EndOfLine', '\r\n');
+
+dataArray{2} = strtrim(dataArray{2});
+
+fclose(fileID);
+
+field_c = dataArray{:, 1} + 1i * dataArray{:, 3};
+
+% Clear temporary variables
+clearvars filename startRow formatSpec fileID dataArray ans;
+
+field_c = reshape(field_c,DIMS');
+
+% Plot the snapshot
+
+field_r = ifourier_GENE(mean(field_c,3));
+
+%
+% toplot = abs(fftshift(field_c(:,:,1),2));
+toplot = real(fftshift(ifourier_GENE(field_c(:,:,1))));
+figure
+pclr = pcolor(toplot); set(pclr,'EdgeColor','none');
+colormap(bluewhitered); shading interp;
\ No newline at end of file
diff --git a/matlab/process_field.m b/matlab/process_field.m
index f574fa2fe4bc8e6bb9b24354228b0404bcc08839..ed7fe9e3683dd2184f1270f06d10be0e4b8a8e78 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -89,18 +89,18 @@ switch OPTIONS.COMP
             i = OPTIONS.COMP;
             compr = @(x) x(i,:,:);
             if REALP
-                COMPNAME = sprintf(['x=','%2.1f'],DATA.x(i));
+                COMPNAME = sprintf(['y=','%2.1f'],DATA.x(i));
             else
-                COMPNAME = sprintf(['k_x=','%2.1f'],DATA.kx(i));
+                COMPNAME = sprintf(['k_y=','%2.1f'],DATA.kx(i));
             end
             FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
         case 2
             i = OPTIONS.COMP;
             compr = @(x) x(:,i,:);
             if REALP
-                COMPNAME = sprintf(['y=','%2.1f'],DATA.y(i));
+                COMPNAME = sprintf(['x=','%2.1f'],DATA.y(i));
             else
-                COMPNAME = sprintf(['k_y=','%2.1f'],DATA.ky(i));
+                COMPNAME = sprintf(['k_x=','%2.1f'],DATA.ky(i));
             end
             FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
         case 3
diff --git a/matlab/setup.m b/matlab/setup.m
index cc4959d6f72e1656a4916fbc91513a7e6f66a578..edba5849f4ff211eacc15cef4f587541f978ac47 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -99,8 +99,8 @@ else
 end
 % temp. dens. drives
 drives_ = [];
-if abs(K_N) > 0; drives_ = [drives_,'_kN_',num2str(K_N)]; end;
-if abs(K_T) > 0; drives_ = [drives_,'_kT_',num2str(K_T)]; end;
+if abs(K_Ni) > 0; drives_ = [drives_,'_kN_',num2str(K_Ni)]; end;
+if abs(K_Ti) > 0; drives_ = [drives_,'_kT_',num2str(K_Ti)]; end;
 % collision
 coll_ = ['_nu_%1.1e_',CONAME];
 coll_   = sprintf(coll_,NU);
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index b6291344786e375ae4db577bef27c67a823cc2c1..da2009fbd175cea0ce6458e4ba91a8867e1ada9f 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -55,7 +55,6 @@ CONTAINS
       CASE DEFAULT
         ERROR STOP 'Error stop: collision model not recognized!!'
     END SELECT
-    print*, collision_kcut
 
   END SUBROUTINE collision_readinputs
 
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 74ec32c551b504574e2b95763c569f86dc0a62ad..665f18ba08fcc5ccb39437f953d76855a32b8618 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -441,12 +441,12 @@ CONTAINS
  IMPLICIT NONE
  REAL :: shift, kx_shift
  ! For periodic CHI BC or 0 dirichlet
- LOGICAL :: PERIODIC_CHI_BC = .false.
+ LOGICAL :: PERIODIC_CHI_BC = .TRUE.
  ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe))
  ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe))
 
  ! No periodic connection for extension of the domain
- IF(Nexc .GT. 1) PERIODIC_CHI_BC = .false.
+ IF(Nexc .GT. 1) PERIODIC_CHI_BC = .TRUE.
 
  !! No shear case (simple id mapping)
  !3            | 1    2    3    4    5    6 |  ky = 3 dky
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
index 7d6bae2e193a3ab1e14a47737f448d89b9b7ea1b..225819204a168c0ac7fec165ca8620f299c0e015 100644
--- a/wk/CBC_kT_PJ_scan.m
+++ b/wk/CBC_kT_PJ_scan.m
@@ -3,117 +3,164 @@ default_plots_options
 HELAZDIR = '/home/ahoffman/HeLaZ/';
 EXECNAME = 'helaz3';
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% KT_a = [9:2:17];
-KT_a = 7;
-g_max= KT_a*0;
-g_avg= KT_a*0;
-g_std= KT_a*0;
-k_max= KT_a*0;
+KT_a = [3:0.5:5];
+P_a_6   = [6 6 6 6 6 6 6];
+J_a_6   = [0 1 2 3 4 5 6];
+P_a_10  = 10*ones(1,6);
+J_a_10  = 5:10;
+% P_a_10b = 10*ones(1,7);
+% J_a_10b = 10:17;
+P_a_20  = 20*ones(1,11);
+J_a_20  = 10:20;
+
+P_a = [P_a_20]; J_a = [J_a_20];
+% KT_a = 5.0; P_a = 20; J_a = 20;
+
+g_max= zeros(numel(P_a),numel(KT_a));
+g_avg= g_max*0;
+g_std= g_max*0;
+k_max= g_max*0;
 
 CO    = 'DG'; GKCO = 0;
-NU    = 0.01;
-DT    = 1e-2;
-TMAX  = 25;
-ky_   = 0.3;
-SIMID = 'linear_CBC_kT_scan_ky_0.3';  % Name of the simulation
-RUN   = 1;
-figure
-P = 4;
-% for P = [2 4 6]
-J = P/2;
+NU    = 0.0;
+DT    = 7e-3;
+TMAX  = 40;
+ky_   = 0.15;
+SIMID = 'linear_CBC_kT_threshold';  % Name of the simulation
+RUN   = 0;
 
-i=1;
-for K_T = KT_a
-    
-%Set Up parameters
-for j = 1
-    CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
-    TAU     = 1.0;    % e/i temperature ratio
-    K_N = 2.22; K_Ne = K_N;
-    K_Te     = K_T;            % Temperature '''
-    SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-    KIN_E   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
-    BETA    = 0e-1;     % electron plasma beta 
-    PMAXE   = P; JMAXE   = J;
-    PMAXI   = P; JMAXI   = J;
-    NX      = 12;    % real space x-gridpoints
-    NY      = 2;     %     ''     y-gridpoints
-    LX      = 2*pi/0.1;   % Size of the squared frequency domain
-    LY      = 2*pi/ky_;
-    NZ      = 16;    % number of perpendicular planes (parallel grid)
-    NPOL    = 1; SG = 0;
-    GEOMETRY= 's-alpha';
-    Q0      = 1.4;    % safety factor
-    SHEAR   = 0.8;    % magnetic shear (Not implemented yet)
-    EPS     = 0.18;    % inverse aspect ratio
-    SPS0D = 1; SPS2D = 0; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
-    JOB2LOAD= -1;
-    LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
-    ABCO    = 1; % interspecies collisions
-    INIT_ZF = 0; ZF_AMP = 0.0;
-    CLOS    = 0;   % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
-    NL_CLOS = 0;   % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
-    KERN    = 0;   % Kernel model (0 : GK)
-    INIT_OPT= 'phi';   % Start simulation with a noisy mom00/phi/allmom
-    W_DOUBLE = 1;
-    W_GAMMA  = 1; W_HF     = 1;
-    W_PHI    = 1; W_NA00   = 1;
-    W_DENS   = 1; W_TEMP   = 1;
-    W_NAPJ   = 1; W_SAPJ   = 0;
-    HD_CO   = 0.0;    % Hyper diffusivity cutoff ratio
-    MU      = 0.0; % Hyperdiffusivity coefficient
-    INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
-    MU_X    = MU;     %
-    MU_Y    = MU;  N_HD    = 4;
-    MU_Z    = 2.0; MU_P    = 0.0;     %
-    MU_J    = 0.0; LAMBDAD = 0.0;
-    NOISE0  = 0.0e-5; % Init noise amplitude
-    BCKGD0  = 1.0;    % Init background
-GRADB   = 1.0;CURVB   = 1.0;
-end
-%%-------------------------------------------------------------------------
-% RUN
-setup
-if RUN
-    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
-end
 
-% Load results
-filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR  = [HELAZDIR,'results/',filename,'/'];
-data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+for i = 1:numel(P_a)
+P = P_a(i); J = J_a(i);
+    j=1;
+    %Set Up parameters
+    for K_Ti = KT_a
+        CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
+        TAU     = 1.0;    % e/i temperature ratio
+        K_Ni = 2.22; K_Ne = K_Ni;
+        K_Te     = K_Ti;            % Temperature '''
+        SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+        KIN_E   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
+        BETA    = 0e-1;     % electron plasma beta 
+        PMAXE   = P; JMAXE   = J;
+        PMAXI   = P; JMAXI   = J;
+        NX      = 8;    % real space x-gridpoints
+        NY      = 6;     %     ''     y-gridpoints
+        LX      = 2*pi/0.15;   % Size of the squared frequency domain
+        LY      = 2*pi/ky_;
+        NZ      = 24;    % number of perpendicular planes (parallel grid)
+        NPOL    = 1; SG = 0; NEXC = 1;
+        GEOMETRY= 's-alpha';
+        Q0      = 1.4;    % safety factor
+        SHEAR   = 0.8;    % magnetic shear (Not implemented yet)
+        EPS     = 0.18;    % inverse aspect ratio
+        SPS0D = 1; SPS2D = 0; SPS3D = 5;SPS5D= 1/5; SPSCP = 0;
+        JOB2LOAD= -1;
+        LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
+        ABCO    = 1; % interspecies collisions
+        INIT_ZF = 0; ZF_AMP = 0.0;
+        CLOS    = 0;   % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+        NL_CLOS = 0;   % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+        KERN    = 0;   % Kernel model (0 : GK)
+        INIT_OPT= 'phi';   % Start simulation with a noisy mom00/phi/allmom
+        W_DOUBLE = 1;
+        W_GAMMA  = 1; W_HF     = 1;
+        W_PHI    = 1; W_NA00   = 1;
+        W_DENS   = 1; W_TEMP   = 1;
+        W_NAPJ   = 1; W_SAPJ   = 0;
+        HD_CO   = 0.0;    % Hyper diffusivity cutoff ratio
+        MU      = 0.0; % Hyperdiffusivity coefficient
+        INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+        MU_X    = MU;     %
+        MU_Y    = MU;  N_HD    = 4;
+        MU_Z    = 1.0; MU_P    = 0.0;     %
+        MU_J    = 0.0; LAMBDAD = 0.0;
+        NOISE0  = 1.0e-4; % Init noise amplitude
+        BCKGD0  = 0.0;    % Init background
+        GRADB   = 1.0;CURVB   = 1.0;
+
+        %%-------------------------------------------------------------------------
+        % RUN
+        setup
+        if RUN
+            system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 2 2 1 0; cd ../../../wk'])
+        end
+
+        % Load results
+        filename = [SIMID,'/',PARAMS,'/'];
+        LOCALDIR  = [HELAZDIR,'results/',filename,'/'];
+        data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+
+        %linear growth rate (adapted for 2D zpinch and fluxtube)
+        trange = [0.5 1]*data.Ts3D(end);
+        nplots = 0;
+        lg = compute_fluxtube_growth_rate(data,trange,nplots);
+        [gmax,     kmax] = max(lg.g_ky(:,end));
+        [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+        msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+        msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
+
+        g_max(i,j) = gmax;
+        k_max(i,j) = kmax;
+        
+        [g_avg(i,j), ik_] = max(lg.avg_g);
+        g_std(i,j) = max(lg.std_g(ik_));
 
-%linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 0;
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
-[gmax,     kmax] = max(lg.g_ky(:,end));
-[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
-msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
-msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
-    
-    g_max(i) = gmax;
-    k_max(i) = kmax;
-    
-    g_avg(i) = lg.avg_g;
-    g_std(i) = lg.std_g;
-    
-    i = i + 1;
+        j = j + 1;
+        if 0
+        %% Verify gamma time trace
+        figure
+        for ik_ = 1:numel(lg.ky)
+            plot(lg.trange(2:end),lg.g_ky(ik_,:)','DisplayName',['$k_y=',num2str(lg.ky(ik_)),'$']); hold on;
+        end
+        xlabel('$t$'); ylabel('$\gamma$');
+        title(data.param_title); legend('show');
+        drawnow
+        end
+    end
 end
-%%
 
-% plot(KT_a,max(g_max,0));
-y_ = g_avg; 
-e_ = g_std;
+if 1
+%% PLOTS
+ERR_WEIGHT = 1/3; %weight of the error to compute marginal stability
+%% Superposed 1D plots
+sz_ = size(g_max);
+figure
+for i = 1:sz_(1)
+    y_ = g_avg(i,:); 
+    e_ = g_std(i,:);
 
-y_ = y_.*(y_-e_>0);
-e_ = e_ .* (y_>0);
-errorbar(KT_a,y_,e_,...
-    'LineWidth',1.2,...
-    'DisplayName',['(',num2str(P),',',num2str(J),')']); 
-hold on;
-title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
-legend('show'); xlabel('$K_T$'); ylabel('$\gamma$');
-drawnow
+    y_ = y_.*(y_-e_*ERR_WEIGHT>0);
+    e_ = e_ .* (y_>0);
+        errorbar(KT_a,y_,e_,...
+            'LineWidth',1.2,...
+            'DisplayName',['(',num2str(P_a(i)),',',num2str(J_a(i)),')']); 
+%     plot(KT_a,y_,...
+%         'LineWidth',1.2,...
+%         'DisplayName',['(',num2str(P_a(i)),',',num2str(J_a(i)),')']); 
+    hold on;
 
+end
+title('Linear CBC $K_T$ threshold');
+legend('show'); xlabel('$K_T$'); ylabel('$\max_{k_y}(\gamma_k)$');
+drawnow
+%% Color map
+[NP__, KT__] = meshgrid(P_a+2*J_a, KT_a);
+% GG_ = g_avg;
+GG_ = g_avg .* (g_avg-g_std*ERR_WEIGHT > 0);
+figure;
+% pclr = pcolor(KT__,NP__,g_max'); set(pclr,'EdgeColor','none');
+pclr = imagesc(KT_a,1:numel(P_a),GG_);
+LABELS = [];
+for i_ = 1:numel(P_a)
+    LABELS = [LABELS; '(',sprintf('%2.0f',P_a(i_)),',',sprintf('%2.0f',J_a(i_)),')'];
+end
+yticks(1:numel(P_a));
+yticklabels(LABELS);
+xlabel('$\kappa_T$'); ylabel('$(P,J)$');
+title('Linear ITG threshold in CBC');
+colormap(bluewhitered);
+%%
+%%
+end
 
diff --git a/wk/Dimits_fig3.m b/wk/Dimits_fig3.m
index 9528aa88453dcae5da7172da19c92ff73dd27d56..21ad776c7a93746f156f11310845501680f5642f 100644
--- a/wk/Dimits_fig3.m
+++ b/wk/Dimits_fig3.m
@@ -37,6 +37,7 @@
 	     4.5 1.1e+0 4.0e-1;...%192x96x24x13x7 kymin=0.05
 	     4.5 9.6e-1 1.5e-1;...%128x64x16x13x2 kymin=0.05
          4.5 7.9e-1 1.8e-1;...%128x64x16x13x7 kymin=0.05
+         4.5 1.2e+0 5.4e-1;...%192x64x24x6x4  kymin=0.05 ! Lx is too small... (weird oscillations)
 	    ];
 	%-------------- GENE ---------------
 	kT_Qi_GENE = ...
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index cc348d028b9958709ef7b9edc0c7880f048f5c76..370c61dd64bb1df607e3ec3109728a6ced6a39d7 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -23,8 +23,8 @@ FMT = '.fig';
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
 % data.scale = 1;%/(data.Nx*data.Ny)^2;
-options.TAVG_0   = 400;%0.4*data.Ts3D(end);
-options.TAVG_1   = 600;%0.9*data.Ts3D(end); % Averaging times duration
+options.TAVG_0   = 350;%0.4*data.Ts3D(end);
+options.TAVG_1   = 1000;%0.9*data.Ts3D(end); % Averaging times duration
 options.NCUT     = 4;              % Number of cuts for averaging and error estimation
 options.NMVA     = 1;              % Moving average for time traces
 % options.ST_FIELD = '\Gamma_x';   % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -45,8 +45,8 @@ if 0
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
-% options.NAME      = '\phi';
-options.NAME      = 'N_i^{00}';
+options.NAME      = '\phi';
+% options.NAME      = 'N_i^{00}';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
@@ -57,7 +57,7 @@ options.PLAN      = 'xy';
 options.COMP      = 'avg';
 % options.TIME      = data.Ts5D(end-30:end);
 % options.TIME      =  data.Ts3D;
-options.TIME      = [00:1:800];
+options.TIME      = [1:0.2:500];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -66,13 +66,13 @@ end
 if 1
 %% 2D snapshots
 % Options
-options.INTERP    = 1;
+options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 0;
-options.NAME      = '\phi';
+% options.NAME      = '\phi';
 % options.NAME      = '\psi';
 % options.NAME      = 'n_e';
-% options.NAME      = 'N_i^{00}';
+options.NAME      = 'N_i^{00}';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
@@ -80,7 +80,7 @@ options.PLAN      = 'kxky';
 % options.NAME      'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [100 200 500];
+options.TIME      = [1000 1100 1200];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 % save_figure(data,fig)
@@ -102,13 +102,13 @@ end
 
 if 0
 %% Kinetic distribution function sqrt(<f_a^2>xy) (GENE vsp)
-% options.SPAR      = linspace(-3,3,32)+(6/127/2);
-% options.XPERP     = linspace( 0,6,32);
-options.SPAR      = gene_data.vp';
-options.XPERP     = gene_data.mu';
+options.SPAR      = linspace(-3,3,32)+(6/127/2);
+options.XPERP     = linspace( 0,6,32);
+% options.SPAR      = gene_data.vp';
+% options.XPERP     = gene_data.mu';
 options.iz        = 'avg';
 options.T         = [250 600];
-options.PLT_FCT   = 'contour';
+options.PLT_FCT   = 'pcolor';
 options.ONED      = 0;
 options.non_adiab = 0;
 options.SPECIE    = 'i';
@@ -123,7 +123,7 @@ if 0
 options.P2J        = 0;
 options.ST         = 1;
 options.PLOT_TYPE  = 'space-time';
-options.NORMALIZED = 1;
+options.NORMALIZED = 0;
 options.JOBNUM     = 0;
 options.TIME       = [1000];
 options.specie     = 'i';
@@ -170,10 +170,10 @@ end
 if 0
 %% Mode evolution
 options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME   = [00:800];
+options.K2PLOT = [0.1 0.2 0.3 0.4];
+options.TIME   = [00:1200];
 options.NMA    = 1;
-options.NMODES = 1;
+options.NMODES = 5;
 options.iz     = 'avg';
 fig = mode_growth_meter(data,options);
 save_figure(data,fig,'.png')
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index 9007b24197b022fb10e398ae757b68b4cce65a58..1d41457a568551b88b8fbceb0207a2a6bfc57f5a 100644
--- a/wk/header_2DZP_results.m
+++ b/wk/header_2DZP_results.m
@@ -1,4 +1,5 @@
 %% Directory of the simulation
+helazdir = '/home/ahoffman/HeLaZ/';
 % if 1% Local results
 outfile ='';
 outfile ='';
@@ -168,4 +169,14 @@ outfile ='';
 % MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
 % end
 
-analysis_HeLaZ
\ No newline at end of file
+%% nu scan
+% outfile = 'Zpinch_rerun/kN_2.2_coll_scan_128x48x5x3';
+% outfile = 'Zpinch_rerun/Ultra_HD_312x196x5x3';
+% outfile = 'Zpinch_rerun/UHD_nu_001_LDGK';
+outfile = 'Zpinch_rerun/UHD_nu_01_LDGK';
+% outfile = 'Zpinch_rerun/UHD_nu_1_LDGK';
+
+%%
+JOBNUMMIN = 01; JOBNUMMAX = 10;
+
+run analysis_HeLaZ
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index a7f051d6af91614932f89399922d25c73c967371..294359d0d28089b0c891acdb3de936bb284feb5e 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -57,10 +57,17 @@ helazdir = '/home/ahoffman/HeLaZ/';
 
 % outfile = 'CBC/kT_scan_128x64x16x5x3';
 % outfile = 'CBC/kT_scan_192x96x16x3x2';
+% outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
+% outfile = 'dbg/nexc_dbg';
+outfile = 'CBC/NM_F4_kT_4.5_192x64x24x6x4';
 
-outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
+% outfile = 'CBC_Ke_EM/192x96x24x5x3';
+% outfile = 'CBC_Ke_EM/96x48x16x5x3';
+% outfile = 'CBC_Ke_EM/minimal_res';
+%% KBM
+% outfile = 'NL_KBM/192x64x24x5x3';
 %% Linear CBC
 % outfile = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
 
-JOBNUMMIN = 00; JOBNUMMAX = 20;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
 run analysis_HeLaZ
diff --git a/wk/lin_TEM.m b/wk/lin_TEM.m
index 6396527a7c4e0bb3988d3439697900690b8d9cfe..35eab68d66563e90babc3acd4b018c931dd50830 100644
--- a/wk/lin_TEM.m
+++ b/wk/lin_TEM.m
@@ -26,19 +26,19 @@ K_Ti    = 6.96;            % ion Temperature '''
 SIGMA_E = 0.05196152422706632;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 % SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 KIN_E   = 1;         % 1: kinetic electrons, 0: adiabatic electrons
-BETA    = 0.0;     % electron plasma beta
+BETA    = 0.001;     % electron plasma beta
 %% GRID PARAMETERS
-P = 4;
-J = P/2;
+P = 3;
+J = 2;
 PMAXE   = P;     % Hermite basis size of electrons
 JMAXE   = J;     % Laguerre "
 PMAXI   = P;     % " ions
 JMAXI   = J;     % "
-NX      = 11;    % real space x-gridpoints
-NY      = 2;     %     ''     y-gridpoints
-LX      = 2*pi/0.1;   % Size of the squared frequency domain
-LY      = 2*pi/0.9;     % Size of the squared frequency domain
-NZ      = 32;    % number of perpendicular planes (parallel grid)
+NX      = 32;    % real space x-gridpoints
+NY      = 16;     %     ''     y-gridpoints
+LX      = 64;%2*pi/0.1;   % Size of the squared frequency domain
+LY      = 200;%2*pi/0.9;     % Size of the squared frequency domain
+NZ      = 16;    % number of perpendicular planes (parallel grid)
 NPOL    = 1;
 SG      = 0;     % Staggered z grids option
 %% GEOMETRY
@@ -47,7 +47,7 @@ GEOMETRY= 's-alpha';
 % GEOMETRY= 'circular';
 Q0      = 1.4;    % safety factor
 SHEAR   = 0.8;    % magnetic shear
-NEXC    = 1;      % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+NEXC    = 4;      % To extend Lx if needed (Lx = Nexc/(kymin*shear))
 EPS     = 0.18;   % inverse aspect ratio
 %% TIME PARMETERS
 TMAX    = 7;  % Maximal time unit