diff --git a/src/srcinfo.h b/src/srcinfo.h
index 9ed53d3da7efc874f514b84e458301a06fa6501f..37401921917253200fe476e0e73ef57cad15ad0f 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
 character(len=20) AUTHOR
 character(len=40) EXECDATE
 character(len=40) HOST
-parameter (VERSION='7ee703f-dirty')
+parameter (VERSION='7642d96-dirty')
 parameter (BRANCH='master')
 parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Wed Apr 13 18:49:20 CEST 2022')
+parameter (EXECDATE='Thu Apr 28 15:46:58 CEST 2022')
 parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 9ed53d3da7efc874f514b84e458301a06fa6501f..37401921917253200fe476e0e73ef57cad15ad0f 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
 character(len=20) AUTHOR
 character(len=40) EXECDATE
 character(len=40) HOST
-parameter (VERSION='7ee703f-dirty')
+parameter (VERSION='7642d96-dirty')
 parameter (BRANCH='master')
 parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Wed Apr 13 18:49:20 CEST 2022')
+parameter (EXECDATE='Thu Apr 28 15:46:58 CEST 2022')
 parameter (HOST ='spcpc606')
diff --git a/testcases/benchmark_HeLaZ_molix.m b/testcases/benchmark_HeLaZ_molix.m
index 053ab242ad725a04bb26c0cefac5802f19d83a88..18a1275a5e2434a38ff0af73c78df701b4430988 100644
--- a/testcases/benchmark_HeLaZ_molix.m
+++ b/testcases/benchmark_HeLaZ_molix.m
@@ -7,7 +7,7 @@ cd ..
 system('make');
 cd wk
 SIMDIR = '../results/benchmarks/shearless_molix_Kin_e/';
-system(['cd ',SIMDIR,';',' ./../../../bin/helaz3.04 0; cd ../../../wk'])
+system(['cd ',SIMDIR,';',' ./../../../bin/helaz3 0; cd ../../../wk'])
 
 % Run molix
 % cd ../../molix
@@ -33,7 +33,7 @@ cd ..
 system('make');
 cd wk
 SIMDIR = '../results/benchmarks/shearless_molix_Adiab_e/';
-system(['cd ',SIMDIR,';',' ./../../../bin/helaz 0; cd ../../../wk'])
+system(['cd ',SIMDIR,';',' ./../../../bin/helaz3; cd ../../../wk'])
 % Run molix
 % cd ../../molix
 % system('make');
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 9608760e51e290ebf73e463d6467b634209a0325..a0ecc04369e97d37ca06373f1f1ae90873b342da 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -21,7 +21,7 @@ FMT = '.fig';
 
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 300; TAVG_1 = 400; % Averaging times duration
+TAVG_0 = 600; TAVG_1 = 700; % Averaging times duration
 compz  = 'avg';
 % chose your field to plot in spacetime diag (uzf,szf,Gx)
 fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'phi',1,compz);
@@ -45,12 +45,12 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'kxky';
+options.PLAN      = 'RZ';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 9;%'avg';
+options.COMP      = 'avg';
 % options.TIME      = dat.Ts5D;
-options.TIME      = 800:820;
+options.TIME      = 0:1:20;
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -63,17 +63,17 @@ options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 % options.NAME      = '\phi';
-% options.NAME      = 'n_i';
+options.NAME      = 'n_i';
 % options.NAME      = 'N_i^{00}';
-options.NAME      = 'T_i';
+% options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
 options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 'avg';
-options.TIME      = [300 350 400];
-data.a = data.EPS * 2000;
+options.COMP      = 8;
+options.TIME      = [500 700 900];
+data.a = data.EPS * 1000;
 fig = photomaton(data,options);
 save_figure(data,fig)
 end
@@ -82,7 +82,7 @@ if 0
 %% 3D plot on the geometry
 options.INTERP    = 1;
 options.NAME      = 'n_i';
-options.PLANES    = 1:3:30;
+options.PLANES    = 1:3:15;
 options.TIME      = [100];
 options.PLT_MTOPO = 1;
 data.rho_o_R      = 2e-3; % Sound larmor radius over Machine size ratio
@@ -107,8 +107,8 @@ end
 if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
-options.P2J  = 1;
-options.ST   = 1;
+options.P2J  = 0;
+options.ST   = 0;
 options.NORMALIZED = 0;
 fig = show_moments_spectrum(data,options);
 save_figure(data,fig)
@@ -151,10 +151,11 @@ end
 if 0
 %% Mode evolution
 options.NORMALIZED = 1;
-options.K2PLOT = 0.01:0.01:1.0;
-options.TIME   = 20:1:60;
+options.K2PLOT = 1;
+options.TIME   = 5:1:15;
 options.NMA    = 1;
-options.NMODES = 20;
+options.NMODES = 5;
+options.iz     = 8;
 fig = mode_growth_meter(data,options);
 save_figure(data,fig)
 end
@@ -181,7 +182,7 @@ end
 
 if 0
 %% linear growth rate for 3D Zpinch
-trange = [30 40];
+trange = [5 15];
 options.keq0 = 1; % chose to plot planes at k=0 or max
 options.kxky = 1;
 options.kzkx = 0;
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index c92702190410cedcba6886d4be2a96dd9b7ae9dc..7ec3b2090bef969b96e954bbd320c1cb65fa4697 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -3,14 +3,8 @@ helazdir = '/home/ahoffman/HeLaZ/';
 % Directory of the simulation
 % if 1% Local results
 outfile ='';
-% outfile ='3D_Zpinch/dbg';
-% outfile ='3D_Zpinch/entropy_mode';
-% outfile ='3D_Zpinch/ITG_Npol_gt_1';
-% outfile ='3D_Zpinch/Ivanov_ITG';
-% outfile ='3D_Zpinch/NL_ITG';
-% outfile ='debug/test_zpar/';
+outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK';
 % outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0/';
-% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_0.8/';
-outfile ='shearless_cyclone/marconi/';
-% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8/';
-call analysis_3D
\ No newline at end of file
+% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0/';
+% outfile ='pedestal/128x128x16x4x2_L_120_LnT_40_nuDG_0.1';
+run analysis_3D
\ No newline at end of file
diff --git a/wk/gene_analysis_3D.m b/wk/gene_analysis_3D.m
index cfde12edd34a82bd438fc6ff1d3805dcf3b29827..700916aa8b9c747b1c5bd071f09dcc8730ddb589 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -1,8 +1,8 @@
 % folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
-% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
-% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.8/';
-folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
+folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
+% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
+% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
 % folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
 % folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
@@ -43,18 +43,18 @@ end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-options.INTERP    = 1;
+options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.NAME      = '\phi';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xy';
+options.PLAN      = 'xz';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 9;
-options.TIME      = 0:700;
+options.COMP      = 'avg';
+options.TIME      = 400:700;
 gene_data.a = data.EPS * 2000;
 create_film(gene_data,options,'.gif')
 end
diff --git a/wk/quick_run.m b/wk/quick_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..c78fa5aee8a924e6d44b1e8910e5c5a769311633
--- /dev/null
+++ b/wk/quick_run.m
@@ -0,0 +1,131 @@
+%% QUICK RUN SCRIPT
+% This script create a directory in /results and run a simulation directly
+% from matlab framework. It is meant to run only small problems in linear
+% for benchmark and debugging purpose since it makes matlab "busy"
+%
+RUN = 1; % To run or just to load
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+HELAZDIR = '/home/ahoffman/HeLaZ/';
+EXECNAME = 'helaz3';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU      = 0.1;   % Collision frequency
+TAU     = 1.0;    % e/i temperature ratio
+K_N     = 20;%1.9;%2.22;   % Density gradient drive
+K_T     = 20;%0.25*K_N;   % Temperature '''
+K_E     = 0.0;   % Electrostat '''
+SIGMA_E = 0.05;%0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+KIN_E   = 1;     % 1: kinetic electrons, 2: adiabatic electrons
+%% GRID PARAMETERS
+PMAXE   = 2;     % Hermite basis size of electrons
+JMAXE   = 1;     % Laguerre "
+PMAXI   = 2;     % " ions
+JMAXI   = 1;     % "
+NX      = 32;    % real space x-gridpoints
+NY      = 32;     %     ''     y-gridpoints
+LX      = 300;   % Size of the squared frequency domain
+LY      = 300;     % Size of the squared frequency domain
+NZ      = 16;     % number of perpendicular planes (parallel grid)
+NPOL    = 1;
+SG      = 0;     % Staggered z grids option
+%% GEOMETRY
+% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+GEOMETRY= 's-alpha';
+Q0      = 2.5;    % safety factor
+SHEAR   = 0.0;    % magnetic shear (Not implemented yet)
+EPS     = 0.18;    % inverse aspect ratio
+%% TIME PARMETERS
+TMAX    = 50;  % Maximal time unit
+DT      = 1e-3;   % Time step
+SPS0D   = 1;      % Sampling per time unit for 2D arrays
+SPS2D   = 0;      % Sampling per time unit for 2D arrays
+SPS3D   = 1;      % Sampling per time unit for 2D arrays
+SPS5D   = 1;    % Sampling per time unit for 5D arrays
+SPSCP   = 0;    % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+SIMID   = 'quick_run';  % Name of the simulation
+LINEARITY = 'linear';   % activate non-linearity (is cancelled if KXEQ0 = 1)
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO      = 'DG';
+GKCO    = 1; % gyrokinetic operator
+ABCO    = 1; % interspecies collisions
+INIT_ZF = 0; ZF_AMP = 0.0;
+CLOS    = 0;   % Closure model (0: =0 truncation, 1: gyrofluid 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
+%% OUTPUTS
+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;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% unused
+HD_CO   = 0.0;    % Hyper diffusivity cutoff ratio
+kmax    = NX*pi/LX;% Highest fourier mode
+MU      = 0.0; % Hyperdiffusivity coefficient
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+MU_X    = MU;     % 
+MU_Y    = MU;     % 
+MU_Z    = 0.0;     %
+MU_P    = 0.0;     %
+MU_J    = 0.0;     % 
+LAMBDAD = 0.0;
+NOISE0  = 1.0e-5; % Init noise amplitude
+BCKGD0  = 0.0;    % Init background
+GRADB   = 1.0;
+CURVB   = 1.0;
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',HELAZDIR,'bin/',EXECNAME,' 2 1 1 0; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 3 2 0; cd ../../../wk'])
+end
+
+%% Load results
+%%
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR  = [HELAZDIR,'results/',filename,'/'];
+% Load outputs from jobnummin up to jobnummax
+JOBNUMMIN = 00; JOBNUMMAX = 00; 
+data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+
+%% Short analysis
+if 1
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+trange = [0.5 1]*data.Ts3D(end);
+nplots = 1;
+lg = compute_fluxtube_growth_rate(data,trange,nplots);
+end
+
+if 0
+%% Ballooning plot
+options.time_2_plot = data.Ts3D(end);
+options.normalized  = 1;
+options.field       = 'phi';
+fig = plot_ballooning(data,options);
+end
+
+if 1
+%% linear growth rate for 3D Zpinch (kz fourier transform)
+trange = [0.5 1]*data.Ts3D(end);
+options.keq0 = 1; % chose to plot planes at k=0 or max
+options.kxky = 1;
+options.kzkx = 0;
+options.kzky = 0;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+save_figure(data,fig)
+end