diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 6d2dfa7652705810e0d94cb24c7910a60b5c902c..015406b7e6d8c1a9e9fba748f88544ae924a915b 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -15,7 +15,7 @@
 % title('gyroLES results for H&P fig. 2c')
 
 %% Plot linear results
-path = '/home/ahoffman/gene/linear_zpinch_results/';
+% path = '/home/ahoffman/gene/linear_zpinch_results/';
 % fname ='tmp.txt';
 % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.047_32x16.txt';
 % fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.047_32x16.txt';
@@ -29,7 +29,7 @@ path = '/home/ahoffman/gene/linear_zpinch_results/';
 % fname ='GENE_LIN_Kn_1.8_KT_0.4_nuLDDK_0.0047_64x32.txt';
 % fname ='GENE_LIN_Kn_1.8_KT_0.45_nu_0_32x16.txt';
 % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.0235_64x32.txt';
-fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
 % fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.235_64x32.txt';
 % fname ='GENE_LIN_Kn_1.7_KT_0.425_nuSG_0.235_64x32.txt';
 % fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSGDK_0.047_32x16.txt';
@@ -40,6 +40,8 @@ fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
 % fname ='GENE_LIN_Kn_2.0_KT_0.5_nuSGDK_0.0235_32x16.txt';
 % fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
 % fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
+path = '/home/ahoffman/gene/linear_CBC_results/';
+fname = 'CBC_100_ky_1e-1_to_1e1.txt';
 data_ = load([path,fname]);
 
 figure
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index 49d5924876e609537fa9cce8f2d068734d9540f1..435573a688a5d16f2b99149aeba7ced167d3724a 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -1,17 +1,27 @@
 function [FIG] = plot_ballooning(data,options)
     FIG.fig = figure; iplot = 1;
-    [~,it] = min(abs(data.Ts3D - options.time_2_plot));
+    [~,it0] = min(abs(data.Ts3D - options.time_2_plot(1)));
+    [~,it1] = min(abs(data.Ts3D - options.time_2_plot(end)));
     [~,ikyarray] = min(abs(data.ky - options.kymodes));
-    phi_real=(real(data.PHI(:,:,:,it)));
-    phi_imag=(imag(data.PHI(:,:,:,it)));
+    phi_real=mean(real(data.PHI(:,:,:,it0:it1)),4);
+    phi_imag=mean(imag(data.PHI(:,:,:,it0:it1)),4);
     % Apply baollooning tranform
     for iky=ikyarray
         dims = size(phi_real);
-        phib_real = zeros(  dims(2)*dims(3)  ,1);
+
+        if options.sheared
+            idx = -Nkx:1:Nkx;
+            ikxlim = dims(2);
+        else
+            idx = 0;
+            ikxlim = 1;
+        end
+        
+        phib_real = zeros(  ikxlim*dims(3)  ,1);
         phib_imag= phib_real;
         b_angle = phib_real;
-
-        for ikx =1: dims(2)
+        
+        for ikx =1:ikxlim
             start_ =  (ikx -1)*dims(3) +1;
             end_ =  ikx*dims(3);
             phib_real(start_:end_)  = phi_real(iky,ikx,:);
@@ -20,8 +30,7 @@ function [FIG] = plot_ballooning(data,options)
 
         % Define ballooning angle
         Nkx = numel(data.kx)-1; coordz = data.z;
-        idx = -Nkx:1:Nkx;
-        for ikx =1: dims(2)
+        for ikx =1: ikxlim
             for iz=1:dims(3)
                 ii = dims(3)*(ikx -1) + iz;
                 b_angle(ii) =coordz(iz) + 2*pi*idx(ikx);
@@ -47,8 +56,8 @@ function [FIG] = plot_ballooning(data,options)
         legend('real','imag','norm')
         xlabel('$\chi / \pi$')
         ylabel('$\phi_B (\chi)$');
-        title(['ky=',sprintf('%1.1f',data.ky(iky)),...
-               ',t=',sprintf('%1.1f',data.Ts3D(it))]);
+        title(['$k_y=',sprintf('%1.1f',data.ky(iky)),...
+               ',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
         iplot = iplot + 1;
     end
 end
diff --git a/matlab/plot/save_figure.m b/matlab/plot/save_figure.m
index 1aa1dd64952fd181fced9e760a7ddfc761a139d5..d1ce56e2c59f6dfa9b6df4e6d62fac4d4dafd2b3 100644
--- a/matlab/plot/save_figure.m
+++ b/matlab/plot/save_figure.m
@@ -1,10 +1,10 @@
 %% Auxiliary script to save figure using a dir (FIGDIR), name (FIGNAME) 
 %  and parameters
 function save_figure(data,FIGURE)
-if ~exist([data.localdir,'/fig'], 'dir')
-   mkdir([data.localdir,'/fig'])
+if ~exist([data.FIGDIR,'/fig'], 'dir')
+   mkdir([data.FIGDIR,'/fig'])
 end
-saveas(FIGURE.fig,[data.localdir,'/fig/', FIGURE.FIGNAME,'.fig']);
-saveas(FIGURE.fig,[data.localdir, FIGURE.FIGNAME,'.png']);
-disp(['Figure saved @ : ',[data.localdir, FIGURE.FIGNAME,'.png']])
+saveas(FIGURE.fig,[data.FIGDIR,'/fig/', FIGURE.FIGNAME,'.fig']);
+saveas(FIGURE.fig,[data.FIGDIR, FIGURE.FIGNAME,'.png']);
+disp(['Figure saved @ : ',[data.FIGDIR, FIGURE.FIGNAME,'.png']])
 end
\ No newline at end of file
diff --git a/matlab/setup.m b/matlab/setup.m
index 5233fad4be26fae7a2537e84e84f66cc3c755764..41d31c2112651cc91b02239a586e584f91504fd2 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -25,7 +25,7 @@ MODEL.KIN_E   = KIN_E;
 if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
 MODEL.mu_x    = MU_X;
 MODEL.mu_y    = MU_Y;
-MODEL.mu_z    = 0;
+MODEL.mu_z    = MU_Z;
 MODEL.mu_p    = MU_P;
 MODEL.mu_j    = MU_J;
 MODEL.nu      = NU; % hyper diffusive coefficient nu for HW
diff --git a/wk/analysis_3D.m b/wk/analysis_HeLaZ.m
similarity index 85%
rename from wk/analysis_3D.m
rename to wk/analysis_HeLaZ.m
index e9ab6b1776f0368871a18153f1704c897b397ae6..1d0382518b1aac8bd7a56a7129af36367271a619 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_HeLaZ.m
@@ -10,16 +10,16 @@ system(['mkdir -p ',MISCDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 08;
+JOBNUMMIN = 00; JOBNUMMAX = 20;
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
-
+data.FIGDIR = LOCALDIR;
 
 %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 default_plots_options
 disp('Plots')
 FMT = '.fig';
 
-if 1
+if 0
 %% Space time diagramm (fig 11 Ivanov 2020)
 options.TAVG_0   = 0.8*data.Ts3D(end); 
 options.TAVG_1   = data.Ts3D(end); % Averaging times duration
@@ -40,7 +40,7 @@ end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-options.INTERP    = 0;
+options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.NAME      = '\phi';
 % options.NAME      = 'N_i^{00}';
@@ -48,12 +48,12 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xy';
+options.PLAN      = 'yz';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 9;
+options.COMP      = 'avg';
 % options.TIME      = dat.Ts5D;
-options.TIME      = 920:1:1250;
+options.TIME      = 700:1:960;
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -62,35 +62,33 @@ end
 if 0
 %% 2D snapshots
 % Options
-options.INTERP    = 1;
+options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
-% options.NAME      = '\phi';
+options.NAME      = '\phi';
 % options.NAME      = 'n_i';
-options.NAME      = 'N_i^{00}';
+% options.NAME      = 'N_i^{00}';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'xy';
+options.PLAN      = 'yz';
 % options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
-options.COMP      = 'avg';
-options.TIME      = [900 923 927 990];
+options.COMP      = 1;
+options.TIME      =[800 900 950];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 save_figure(data,fig)
 end
 
-if 0
-%% 3D plot on the geometry
-options.INTERP    = 1;
-options.NAME      = 'n_i';
-options.PLANES    = [1:2:12];
-options.TIME      = [60];
-options.PLT_MTOPO = 1;
-data.rho_o_R      = 2e-3; % Sound larmor radius over Machine size ratio
-fig = show_geometry(data,options);
-save_figure(data,fig)
+if 1
+%% Ballooning plot
+options.time_2_plot = [800 900];
+options.kymodes     = [0.5];
+options.normalized  = 1;
+options.sheared     = 0;
+options.field       = 'phi';
+fig = plot_ballooning(data,options);
 end
 
 if 0
@@ -99,9 +97,9 @@ if 0
 % options.XPERP     = linspace( 0,6,64);
 options.SPAR      = gene_data.vp';
 options.XPERP     = gene_data.mu';
-options.iz        = 13;
-options.T         = 20;
-options.PLT_FCT   = 'pcolor';
+options.iz        = 'avg';
+options.T         = 1000;
+options.PLT_FCT   = 'contour';
 options.ONED      = 0;
 options.non_adiab = 1;
 options.SPECIE    = 'i';
@@ -164,12 +162,12 @@ if 0
 %% Mode evolution
 options.NORMALIZED = 1;
 options.K2PLOT = 1;
-options.TIME   = 2:20;
+options.TIME   = [0.5 1]*data.Ts3D(end);
 options.NMA    = 1;
 options.NMODES = 15;
 options.iz     = 9;
 fig = mode_growth_meter(data,options);
-save_figure(data,fig)
+save_figure(gbms_dat,fig)
 end
 
 if 0
@@ -202,3 +200,15 @@ options.kzky = 1;
 [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
 save_figure(data,fig)
 end
+
+if 0
+%% 3D plot on the geometry
+options.INTERP    = 1;
+options.NAME      = 'n_i';
+options.PLANES    = [1:2:12];
+options.TIME      = [60];
+options.PLT_MTOPO = 1;
+data.rho_o_R      = 2e-3; % Sound larmor radius over Machine size ratio
+fig = show_geometry(data,options);
+save_figure(data,fig)
+end
\ No newline at end of file
diff --git a/wk/analysis_gbms.m b/wk/analysis_gbms.m
new file mode 100644
index 0000000000000000000000000000000000000000..f3a35be5dd2010050b51271abc41a00e9b589989
--- /dev/null
+++ b/wk/analysis_gbms.m
@@ -0,0 +1,81 @@
+% cd /home/ahoffman/Documents/gbms/scan_test
+% addpath(genpath('/home/ahoffman/Documents/gbms/matlab_scripts'));
+% res = gbms_get_scandir('/home/ahoffman/Documents/gbms/scan_test/scan_test/');
+% figure; plot(res.paramscan,res.growth_rate)
+
+%%
+resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/shearless_linear_cyclone/';
+% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/linear_cyclone/';
+% resdir = '/home/ahoffman/molix/';
+outfile = [resdir,'field.dat.h5'];
+
+gbms_dat.Ts3D    = h5read(outfile,'/data/var2d/time');
+gbms_dat.Nt      = unique(numel(gbms_dat.Ts3D));
+gbms_dat.kx      = unique(h5read(outfile,'/data/var2d/phi/coordkx'));
+gbms_dat.ky      = unique(h5read(outfile,'/data/var2d/phi/coordky'));
+gbms_dat.z       = unique(h5read(outfile,'/data/var2d/phi/coordz'));
+gbms_dat.Nx = numel(gbms_dat.kx); gbms_dat.Nkx = numel(gbms_dat.kx); 
+gbms_dat.Ny = numel(gbms_dat.ky); gbms_dat.Nky = numel(gbms_dat.ky); 
+gbms_dat.Nz = numel(gbms_dat.z);
+
+dky = min(gbms_dat.ky(gbms_dat.ky>0)); Ly = 2*pi/dky;
+gbms_dat.y  = linspace(-Ly/2,Ly/2,gbms_dat.Ny+1); gbms_dat.y = gbms_dat.y(1:end-1);
+gbms_dat.x  = 0;
+gbms_dat.PHI = zeros(gbms_dat.Ny,gbms_dat.Nx,gbms_dat.Nz,gbms_dat.Nt);
+gbms_dat.param_title = 'GBMS';
+for it = 1:gbms_dat.Nt
+    
+    tmp = h5read(outfile,['/data/var2d/phi/',sprintf('%.6d',it-1)]);
+    gbms_dat.PHI(:,:,:,it) = permute(tmp.real + 1i * tmp.imaginary,[2 1 3]);
+    
+end
+
+gbms_dat.localdir = resdir;
+
+%%
+if 0
+%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Options
+options.INTERP    = 1;
+options.POLARPLOT = 0;
+options.NAME      = '\phi';
+options.PLAN      = 'yz';
+options.COMP      = 1;
+options.TIME      = 0:200;
+gbms_dat.EPS          = 0.1
+gbms_dat.a = gbms_dat.EPS * 2000;
+create_film(gbms_dat,options,'.gif')
+end
+
+if 0
+%% 2D snapshots
+% Options
+options.INTERP    = 0;
+options.POLARPLOT = 0;
+options.AXISEQUAL = 1;
+options.NAME      = '\phi';
+options.PLAN      = 'yz';
+options.COMP      = 1;
+options.TIME      = 100;
+gbms_dat.EPS = 1e-3;
+gbms_dat.a = gbms_dat.EPS * 2000;
+fig = photomaton(gbms_dat,options);
+save_figure(gbms_dat,fig)
+end
+
+if 0
+%% linear growth rate for 3D fluxtube
+trange = [10 200];
+nplots = 1;
+lg = compute_fluxtube_growth_rate(gbms_dat,trange,nplots);
+end
+
+if 1
+%% Ballooning plot
+options.time_2_plot = data.Ts3D(end);
+options.kymodes     = [0.5];
+options.normalized  = 1;
+options.sheared     = 0;
+options.field       = 'phi';
+fig = plot_ballooning(gbms_dat,options);
+end
\ No newline at end of file
diff --git a/wk/gene_analysis_3D.m b/wk/analysis_gene.m
similarity index 84%
rename from wk/gene_analysis_3D.m
rename to wk/analysis_gene.m
index 31652e6e02748c7eeb11b462fed77d82d5620c3d..0fe05469c77ecd4444c4c8922bc0ddaa687fc4d8 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/analysis_gene.m
@@ -1,7 +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/linear_s_alpha_CBC_100/';
+folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.2/';
+% folder = '/misc/gene_results/shearless_cyclone/linear_s_alpha_CBC_100/';
+% folder = '/misc/gene_results/shearless_cyclone/allmodes_CBC_100/';
 % 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/';
@@ -16,8 +17,9 @@ options.TAVG_1   = gene_data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
 options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
 options.INTERP   = 1;
+gene_data.FIGDIR = folder;;
 fig = plot_radial_transport_and_spacetime(gene_data,options);
-% save_figure(data,fig)
+save_figure(gene_data,fig)
 end
 
 if 0
@@ -34,8 +36,8 @@ options.NAME      = '\phi';
 options.PLAN      = 'kxky';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
-options.COMP      = 'avg';
-options.TIME      = [0 500];
+options.COMP      = 9;
+options.TIME      = [0 50 100 200 300];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig)
@@ -55,7 +57,7 @@ options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = 000:700;
+options.TIME      = 000:300;
 gene_data.a = data.EPS * 2000;
 create_film(gene_data,options,'.gif')
 end
@@ -88,11 +90,11 @@ end
 
 if 0
 %% Show f_i(vpar,mu)
-options.times   = 200:600;
+options.times   = 200:300;
 options.specie  = 'i';
-options.PLT_FCT = 'pcolor';
+options.PLT_FCT = 'contour';
 options.folder  = folder;
-options.iz      = 1;
+options.iz      = 9;
 options.FIELD   = '<f_>';
 options.ONED    = 0;
 % options.FIELD   = 'Q_es';
@@ -101,12 +103,12 @@ end
 
 if 0
 %% Mode evolution
-options.NORMALIZED = 1;
+options.NORMALIZED = 0;
 options.K2PLOT = 1;
-options.TIME   = 50:150;
+options.TIME   = 100:200;
 options.NMA    = 1;
-options.NMODES = 5;
+options.NMODES = 15;
 options.iz     = 9;
 fig = mode_growth_meter(gene_data,options);
-save_figure(data,fig)
-end
\ No newline at end of file
+save_figure(gene_data,fig)
+end
diff --git a/wk/analysis_header_2DZP.m b/wk/header_2DZP_results.m
similarity index 99%
rename from wk/analysis_header_2DZP.m
rename to wk/header_2DZP_results.m
index cd4f356ac323fb751a737043e92399b2fcbc9e41..9007b24197b022fb10e398ae757b68b4cce65a58 100644
--- a/wk/analysis_header_2DZP.m
+++ b/wk/header_2DZP_results.m
@@ -168,4 +168,4 @@ outfile ='';
 % MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
 % end
 
-analysis_3D
\ No newline at end of file
+analysis_HeLaZ
\ No newline at end of file
diff --git a/wk/analysis_header.m b/wk/header_3D_results.m
similarity index 89%
rename from wk/analysis_header.m
rename to wk/header_3D_results.m
index 3ff76bdcb862c14adbc219320007d520a20714a2..113d91144a967916ca693fcd1b9683d5cbaee6b2 100644
--- a/wk/analysis_header.m
+++ b/wk/header_3D_results.m
@@ -7,11 +7,11 @@ outfile ='';
 outfile ='';
 % outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
 % outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
-outfile ='shearless_cyclone/linear_CBC_120';
+outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0';
 % outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
 % outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1';
 % outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK';
 % outfile ='shearless_cyclone/128x128x16x8x4_L_120_CTC_1.0/';
 % 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
+run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index f5a8c35d3367a5850a2ab1112a06574d467f67cd..8d437c867bd14636a22241ce8e5604aa2b1e381d 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -13,7 +13,7 @@ EXECNAME = 'helaz3';
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.0;   % Collision frequency
+NU      = 0.2;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 K_N     = 2.22;   % Density gradient drive
 K_T     = 6.96;   % Temperature '''
@@ -21,15 +21,15 @@ K_E     = 0.0;   % Electrostat '''
 SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 KIN_E   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
 %% GRID PARAMETERS
-PMAXE   = 4;     % Hermite basis size of electrons
-JMAXE   = 2;     % Laguerre "
-PMAXI   = 4;     % " ions
-JMAXI   = 2;     % "
+PMAXE   = 8;     % Hermite basis size of electrons
+JMAXE   = 4;     % Laguerre "
+PMAXI   = 8;     % " ions
+JMAXI   = 4;     % "
 NX      = 1;    % real space x-gridpoints
-NY      = 24;     %     ''     y-gridpoints
+NY      = 64;     %     ''     y-gridpoints
 LX      = 100;   % Size of the squared frequency domain
-LY      = 100;     % Size of the squared frequency domain
-NZ      = 32;     % number of perpendicular planes (parallel grid)
+LY      = 60;     % Size of the squared frequency domain
+NZ      = 16;     % number of perpendicular planes (parallel grid)
 NPOL    = 1;
 SG      = 0;     % Staggered z grids option
 %% GEOMETRY
@@ -39,8 +39,8 @@ Q0      = 1.4;    % safety factor
 SHEAR   = 0.0;    % magnetic shear (Not implemented yet)
 EPS     = 0.18;    % inverse aspect ratio
 %% TIME PARMETERS
-TMAX    = 10;  % Maximal time unit
-DT      = 2e-2;   % Time step
+TMAX    = 30;  % Maximal time unit
+DT      = 2*5e-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
@@ -48,12 +48,12 @@ 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
+SIMID   = 'dbg';  % 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
+GKCO    = 0; % gyrokinetic operator
 ABCO    = 1; % interspecies collisions
 INIT_ZF = 0; ZF_AMP = 0.0;
 CLOS    = 1;   % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))s
@@ -89,10 +89,10 @@ setup
 system(['rm fort*.90']);
 % Run linear simulation
 if RUN
-    system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 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'])
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
 end
 
 %% Load results
@@ -111,16 +111,17 @@ nplots = 1;
 lg = compute_fluxtube_growth_rate(data,trange,nplots);
 end
 
-if 0
+if 1
 %% Ballooning plot
 options.time_2_plot = data.Ts3D(end);
-options.kymodes     = [0.4 0.5 0.6];
+options.kymodes     = [0.5];
 options.normalized  = 1;
+options.sheared     = 0;
 options.field       = 'phi';
 fig = plot_ballooning(data,options);
 end
 
-if 1
+if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
 options.P2J        = 1;
@@ -133,8 +134,8 @@ options.JOBNUM     = 0;
 options.TIME       = [0 50];
 options.specie     = 'i';
 options.compz      = 'avg';
-% fig = show_moments_spectrum(data,options);
-fig = show_napjz(data,options);
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
 save_figure(data,fig)
 end
 
@@ -148,3 +149,15 @@ options.kzky = 0;
 [lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
 save_figure(data,fig)
 end
+
+if 1
+%% Mode evolution
+options.NORMALIZED = 1;
+options.K2PLOT = 1;
+options.TIME   = [0.8 1]*data.Ts3D(end);
+options.NMA    = 1;
+options.NMODES = 5;
+options.iz     = 9;
+fig = mode_growth_meter(data,options);
+save_figure(gbms_dat,fig)
+end