diff --git a/matlab/compute/invert_kxky_to_kykx_gene_results.m b/matlab/compute/invert_kxky_to_kykx_gene_results.m
index 462846b7a85aa00d4f964901258455899898df68..3f320d653472aa52e6583aee1097355680323f1e 100644
--- a/matlab/compute/invert_kxky_to_kykx_gene_results.m
+++ b/matlab/compute/invert_kxky_to_kykx_gene_results.m
@@ -9,19 +9,19 @@ data.TPER_I = rotate(data.TPER_I);
 data.TPAR_I = rotate(data.TPAR_I);
 data.TEMP_I = rotate(data.TEMP_I);
 
-data.Ny = data.Nky*2-1;
-data.Nx = data.Nkx;
+data.grids.Ny = data.grids.Nky*2-1;
+data.grids.Nx = data.grids.Nkx;
 
-if numel(data.kx)>1
-    dkx = data.kx(2); 
+if numel(data.grids.kx)>1
+    dkx = data.grids.kx(2); 
 else
     dkx = 1;
 end
 
-dky = data.ky(2);
+dky = data.grids.ky(2);
 Lx = 2*pi/dkx;   Ly = 2*pi/dky;
-x = linspace(-Lx/2,Lx/2,data.Nx+1); x = x(1:end-1);
-y = linspace(-Ly/2,Ly/2,data.Ny+1); y = y(1:end-1);
-data.x = x; data.y = y; data.Lx = Lx; data.Ly = Ly;
+x = linspace(-Lx/2,Lx/2,data.grids.Nx+1); x = x(1:end-1);
+y = linspace(-Ly/2,Ly/2,data.grids.Ny+1); y = y(1:end-1);
+data.grids.x = x; data.grids.y = y; data.grids.Lx = Lx; data.grids.Ly = Ly;
 end
 
diff --git a/matlab/load/compile_results_2D.m b/matlab/load/compile_results_2D.m
new file mode 100644
index 0000000000000000000000000000000000000000..b2ec59459f7093786433df59687a77e403555c2d
--- /dev/null
+++ b/matlab/load/compile_results_2D.m
@@ -0,0 +1,37 @@
+function [field, Ts2D] = compile_results_2D(DIRECTORY,JOBNUMMIN,JOBNUMMAX,fieldname)
+CONTINUE = 1;
+JOBNUM   = JOBNUMMIN; JOBFOUND = 0;
+% field
+field    = [];
+Ts2D    = [];
+ii = 1;
+while(CONTINUE)
+    filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
+    % Check presence and jobnummax
+    if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
+        %test if it is corrupted or currently running
+        try
+            openable = ~isempty(h5read(filename,'/data/var2d/time'));
+        catch
+            openable = 0;
+        end
+        if openable
+        % load field %%
+        [ F, T, ~] = load_2D_data(filename, fieldname);
+        field  = cat(3,field,F);
+        Ts2D  = cat(1,Ts2D,T);
+        ii = ii + 1;
+        JOBFOUND = JOBFOUND + 1;
+        end
+    elseif (JOBNUM > JOBNUMMAX)
+        CONTINUE = 0;
+    end
+    JOBNUM   = JOBNUM + 1;
+end
+
+if(JOBFOUND == 0)
+    disp('no results found, please verify the paths');
+    return;
+end
+
+end
\ No newline at end of file
diff --git a/matlab/load/compile_results_5Da.m b/matlab/load/compile_results_5Da.m
new file mode 100644
index 0000000000000000000000000000000000000000..914de52aeb7af3d7515a422f3ec87e4212c2e243
--- /dev/null
+++ b/matlab/load/compile_results_5Da.m
@@ -0,0 +1,37 @@
+function [field, Ts5D] = compile_results_5Da(DIRECTORY,JOBNUMMIN,JOBNUMMAX,fieldname)
+CONTINUE = 1;
+JOBNUM   = JOBNUMMIN; JOBFOUND = 0;
+% field
+field    = [];
+Ts5D    = [];
+ii = 1;
+while(CONTINUE)
+    filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
+    % Check presence and jobnummax
+    if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
+        %test if it is corrupted or currently running
+        try
+            openable = ~isempty(h5read(filename,'/data/var5d/time'));
+        catch
+            openable = 0;
+        end
+        if openable
+        % load field %%
+        [ F, T, ~] = load_5Da_data(filename, fieldname);
+        field  = cat(7,field,F);
+        Ts5D  = cat(1,Ts5D,T);
+        ii = ii + 1;
+        JOBFOUND = JOBFOUND + 1;
+        end
+    elseif (JOBNUM > JOBNUMMAX)
+        CONTINUE = 0;
+    end
+    JOBNUM   = JOBNUM + 1;
+end
+
+if(JOBFOUND == 0)
+    disp('no results found, please verify the paths');
+    return;
+end
+
+end
\ No newline at end of file
diff --git a/matlab/load/compile_results_low_mem.m b/matlab/load/compile_results_low_mem.m
index cac486ac58f1d3c5863e3f0a75f2853d79fa4000..0c80bd79b489899c6d3289c367cbf1eb921e2066 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -1,6 +1,8 @@
 function [DATA] = compile_results_low_mem(DATA,DIRECTORY,JOBNUMMIN,JOBNUMMAX)
 CONTINUE = 1;
 JOBNUM   = JOBNUMMIN; JOBFOUND = 0;
+DATA.CODENAME = 'GYAC';
+DATA.folder = DIRECTORY;
 DATA.params_evol.TJOB_SE  = []; % Start and end times of jobs
 DATA.params_evol.NU_EVOL  = []; % evolution of parameter nu between jobs
 DATA.params_evol.CO_EVOL  = []; % evolution of CO
@@ -51,8 +53,8 @@ while(CONTINUE)
                     [ GGAMMA_R, Ts0D, ~] = load_0D_data(filename, 'gflux_xi');
                     PGAMMA_R             = load_0D_data(filename, 'pflux_xi');
                 end
-                GGAMMA_ = cat(1,GGAMMA_,GGAMMA_R); clear GGAMMA_R
-                PGAMMA_ = cat(1,PGAMMA_,PGAMMA_R); clear PGAMMA_R
+                GGAMMA_ = cat(2,GGAMMA_,GGAMMA_R); clear GGAMMA_R
+                PGAMMA_ = cat(2,PGAMMA_,PGAMMA_R); clear PGAMMA_R
             end
 
             if W_HF
@@ -62,7 +64,7 @@ while(CONTINUE)
                     % old version
                     [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_xi');
                 end
-                HFLUX_ = cat(1,HFLUX_,HFLUX_X); clear HFLUX_X
+                HFLUX_ = cat(2,HFLUX_,HFLUX_X); clear HFLUX_X
             end
 
             Ts0D_   = cat(1,Ts0D_,Ts0D);
diff --git a/matlab/load/load_0D_data.m b/matlab/load/load_0D_data.m
index a4d20bd59a2f8d139eda8e338f57bc6a311db639..c9d75ee1d21e08fdeb79f72e4e4ac2e50cdb1a44 100644
--- a/matlab/load/load_0D_data.m
+++ b/matlab/load/load_0D_data.m
@@ -5,6 +5,6 @@ function [ data, time, dt ] = load_0D_data( filename, variablename )
     Na    = h5readatt(filename,'/data/input/model','Na');
     cstart= h5readatt(filename,'/data/input/basic','start_iframe0d'); 
     data  = h5read(filename,['/data/var0d/',variablename]);
-%     data  = reshape(data,[Na, numel(time)]);
+    data  = reshape(data,[Na, numel(time)]);
 end
 
diff --git a/matlab/load/load_2D_data.m b/matlab/load/load_2D_data.m
index b013a63d4be724587abf43b7913cba5421387647..67da60c9dd1931e73cde73c7520ecb6370a0c153 100644
--- a/matlab/load/load_2D_data.m
+++ b/matlab/load/load_2D_data.m
@@ -1,15 +1,28 @@
 function [ data, time, dt ] = load_2D_data( filename, variablename )
-%LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ
+%LOAD_3D_DATA load a 3D variable stored in a hdf5 result file from HeLaZ
     time     = h5read(filename,'/data/var2d/time');
-    kx       = h5read(filename,'/data/grid/coordkx');
-    ky       = h5read(filename,'/data/grid/coordky');
-    dt    = h5readatt(filename,'/data/input','dt');
-    cstart= h5readatt(filename,'/data/input','start_iframe2d'); 
+    dt    = h5readatt(filename,'/data/input/basic','dt');
+    cstart= h5readatt(filename,'/data/input/basic','start_iframe3d'); 
+    
+    % Find array size by loading the first output
+    tmp   = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+1,'%06d')]);
+    try % check if it is complex or real
+        sz  = size(tmp.real);
+        cmpx = 1;
+    catch
+        sz  = size(tmp);
+        cmpx = 0;
+    end
+
+    data     = zeros(sz);
     
-    data     = zeros(numel(kx),numel(ky),numel(time));
     for it = 1:numel(time)
         tmp         = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+it,'%06d')]);
-        data(:,:,it) = tmp.real + 1i * tmp.imaginary;
+        if cmpx
+            data(:,:,it) = tmp.real + 1i * tmp.imaginary;
+        else
+            data(:,:,it) = tmp;
+        end
     end
 
 end
diff --git a/matlab/load/load_5Da_data.m b/matlab/load/load_5Da_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..2077b119d61e69e5eda85b105bd5465bf23241e7
--- /dev/null
+++ b/matlab/load/load_5Da_data.m
@@ -0,0 +1,36 @@
+function [ data, time, dt ] = load_5Da_data( filename, variablename )
+%LOAD_3D_DATA load a 3D variable stored in a hdf5 result file from HeLaZ
+    time  = h5read(filename,'/data/var5d/time');
+    dt    = h5readatt(filename,'/data/input/basic','dt');
+    cstart= h5readatt(filename,'/data/input/basic','start_iframe3d'); 
+    Na    = h5readatt(filename,'/data/input/model','Na'); 
+    Np    = h5readatt(filename,'/data/input/grid', 'Np'); 
+    Nj    = h5readatt(filename,'/data/input/grid', 'Nj'); 
+    Nky   = h5readatt(filename,'/data/input/grid', 'Nky'); 
+    Nkx   = h5readatt(filename,'/data/input/grid', 'Nkx'); 
+    Nz    = h5readatt(filename,'/data/input/grid', 'Nz'); 
+
+    
+    % Find array size by loading the first output
+    tmp   = h5read(filename,['/data/var5d/',variablename,'/', num2str(cstart+1,'%06d')]);
+    try % check if it is complex or real
+        sz  = size(tmp.real);
+        cmpx = 1;
+    catch
+        sz  = size(tmp);
+        cmpx = 0;
+    end
+    % add time dimension
+    data     = zeros([Na,Np,Nj,Nky,Nkx,Nz,numel(time)]);
+    
+    for it = 1:numel(time)
+        tmp         = h5read(filename,['/data/var5d/',variablename,'/', num2str(cstart+it,'%06d')]);
+        if cmpx
+            data(:,:,:,:,:,:,it) = reshape(tmp.real,sz) + 1i * reshape(tmp.imaginary,sz);
+        else
+            data(:,:,:,:,:,:,it) = reshape(tmp,sz);
+        end
+    end
+
+end
+
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index 8c877d2e5c235a922efaf812d17ec03ce216e05a..196e899617f3ed7af330823bf805d1bb3f7e4c1a 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -5,35 +5,35 @@ DATA.namelist = namelist;
 DATA.folder   = folder;
 %% Grid
 coofile = 'coord.dat.h5';
-DATA.vp  = h5read([folder,coofile],'/coord/vp');
-DATA.Nvp = numel(DATA.vp);
+DATA.grids.vp  = h5read([folder,coofile],'/coord/vp');
+DATA.grids.Nvp = numel(DATA.grids.vp);
 
-DATA.mu  = h5read([folder,coofile],'/coord/mu');
-DATA.Nmu = numel(DATA.mu);
+DATA.grids.mu  = h5read([folder,coofile],'/coord/mu');
+DATA.grids.Nmu = numel(DATA.grids.mu);
 
-DATA.kx  = h5read([folder,coofile],'/coord/kx');
-DATA.Nkx = numel(DATA.kx);
-DATA.Nx  = DATA.Nkx;
+DATA.grids.kx  = h5read([folder,coofile],'/coord/kx');
+DATA.grids.Nkx = numel(DATA.grids.kx);
+DATA.grids.Nx  = DATA.grids.Nkx;
 
-DATA.ky  = h5read([folder,coofile],'/coord/ky');
-DATA.Nky = numel(DATA.ky);
-DATA.Ny  = DATA.Nky*2-1;
+DATA.grids.ky  = h5read([folder,coofile],'/coord/ky');
+DATA.grids.Nky = numel(DATA.grids.ky);
+DATA.grids.Ny  = DATA.grids.Nky*2-1;
 
-DATA.z   = h5read([folder,coofile],'/coord/z');
-DATA.Nz  = numel(DATA.z);
+DATA.grids.z   = h5read([folder,coofile],'/coord/z');
+DATA.grids.Nz  = numel(DATA.grids.z);
 
-DATA.paramshort = [num2str(DATA.Nkx),'x',num2str(DATA.Nky),'x',num2str(DATA.Nz),...
-                    'x',num2str(DATA.Nvp),'x',num2str(DATA.Nmu)];
-if numel(DATA.kx)>1
-    dkx = DATA.kx(2); 
+DATA.params_string = [num2str(DATA.grids.Nkx),'x',num2str(DATA.grids.Nky),'x',num2str(DATA.grids.Nz),...
+                    'x',num2str(DATA.grids.Nvp),'x',num2str(DATA.grids.Nmu)];
+if numel(DATA.grids.kx)>1
+    dkx = DATA.grids.kx(2); 
 else
     dkx = 1;
 end
-dky = DATA.ky(2);
+dky = DATA.grids.ky(2);
 Lx = 2*pi/dkx;   Ly = 2*pi/dky;
-x = linspace(-Lx/2,Lx/2,DATA.Nx+1); x = x(1:end-1);
-y = linspace(-Ly/2,Ly/2,DATA.Ny+1); y = y(1:end-1);
-DATA.x = x; DATA.y = y; DATA.Lx = Lx; DATA.Ly = Ly;
+x = linspace(-Lx/2,Lx/2,DATA.grids.Nx+1); x = x(1:end-1);
+y = linspace(-Ly/2,Ly/2,DATA.grids.Ny+1); y = y(1:end-1);
+DATA.grids.x = x; DATA.grids.y = y; DATA.grids.Lx = Lx; DATA.grids.Ly = Ly;
 %% Transport
 nrgfile           = 'nrg.dat.h5';
 % nrgfile           = 'nrg_1.h5';
@@ -45,10 +45,10 @@ DATA.HFLUX_X   = h5read([folder,nrgfile],'/nrgions/Q_es');
 phifile   = 'field.dat.h5';
 % phifile   = 'field_1.h5';
 DATA.Ts3D = h5read([folder,phifile],'/field/time');
-DATA.DENS_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
-DATA.TPER_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
-DATA.TPAR_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
-DATA.PHI    = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
+DATA.DENS_I = zeros(DATA.grids.Nkx,DATA.grids.Nky,DATA.grids.Nz,numel(DATA.Ts3D));
+DATA.TPER_I = zeros(DATA.grids.Nkx,DATA.grids.Nky,DATA.grids.Nz,numel(DATA.Ts3D));
+DATA.TPAR_I = zeros(DATA.grids.Nkx,DATA.grids.Nky,DATA.grids.Nz,numel(DATA.Ts3D));
+DATA.PHI    = zeros(DATA.grids.Nkx,DATA.grids.Nky,DATA.grids.Nz,numel(DATA.Ts3D));
 
 momfile = 'mom_ions.dat.h5';
 % momfile = 'mom_ions_1.h5';
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index 02291c8e2e11b50772b0763ae61e5f833c4ecb39..d4d71fccd12d133015d88c0d4da74ccce9fb0964 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -6,98 +6,83 @@ Ja = DATA.grids.Jarray;
 Time_ = DATA.Ts3D;
 FIGURE.fig = figure; FIGURE.FIGNAME = ['mom_spectrum_',DATA.params_string];
 set(gcf, 'Position',  [100 50 1000 400])
+if OPTIONS.ST == 0
+    OPTIONS.LOGSCALE = 0;
+end
 if OPTIONS.LOGSCALE
-%     compress = @(x,ia) log(sum(abs(squeeze(x(ia,:,:,:))),3));
-    compress = @(x,ia) log(sum(abs(squeeze(x(:,:,:,:))),3));
+    logname = 'log';
+    compress = @(x,ia) squeeze(log(sum(abs(x(ia,:,:,:,:)),4)));
+%     compress = @(x,ia) log(sum(abs(squeeze(x(:,:,:,:))),3));
 else
-%     compress = @(x,ia) sum(abs(squeeze(x(ia,:,:,:))),3);
-    compress = @(x,ia) sum(abs(squeeze(x(:,:,:,:))),3);
+    logname = '';
+    compress = @(x,ia) squeeze(sum(abs(x(ia,:,:,:,:)),4));
+%     compress = @(x,ia) sum(abs(squeeze(x(:,:,:,:))),3);
 end
 for ia = 1:DATA.inputs.Na
 %     Napjz  = sum(abs(squeeze(DATA.Napjz(ia,:,:,:,:))),3);
-    Napjz  =compress(DATA.Napjz);
+    Napjz  =compress(DATA.Napjz,ia);
     subplot(double(DATA.inputs.Na),1,double(ia))
-    if OPTIONS.P2J
-        plotname = ['$\langle\sum_k |N_',species_name{ia},'^{pj}|\rangle_{p+2j=const}$'];
-        [JJ,PP] = meshgrid(Ja,Pa);
-        P2J = PP + 2*JJ;
-        % form an axis of velocity ordered moments
-        p2j = unique(P2J); yname = '$p+2j$';
-        % weights to average
-        weights = zeros(size(p2j));
-        % space time of moments amplitude wrt to p+2j degree
-        Na_ST = zeros(numel(p2j),numel(Time_));
-        % fill the st diagramm by averaging every p+2j=n moments together
-        for ip = 1:numel(Pa)
-            for ij = 1:numel(Ja)
-                [~,ip2j] = min(abs(p2j-(Pa(ip)+2*Ja(ij))));
-                Na_ST(ip2j,:) = Na_ST(ip2j,:) + transpose(squeeze(Napjz(ip,ij,:)));
-                weights(ip2j) = weights(ip2j) + 1;
-            end
-        end
-        % doing the average
-        for ip2j = 1:numel(p2j)
-            Na_ST(ip2j,:) = Na_ST(ip2j,:)/weights(ip2j);
-        end
-        ticks_labels = {};
-    else % We just order our moments w.r.t. to the convention ..
-         % (0,0) (1,0) (2,0) (0,1) (3,0) (1,1) (4,0) (2,1) (0,2) etc.
-        plotname = '$\langle\sum_k |N_i^{pj}|^2\rangle_z$';
-        Nmoments = numel(Napjz(:,:,1)); % total number of moments
-        HL_deg   = zeros(Nmoments,2);   % list of degrees, first column Hermite, second Laguerre
-        im  = 2;
-        deg = 1; % the zero degree is always here first place
-        ticks_labels = cell(10,1);
-        ticks_labels{1} = '(0,0)';
-        while(im<=Nmoments)
-        FOUND = 1;
-            while(FOUND) % As we find a pair matching the degree we retry
-            FOUND = 0;
-            for ij = 1:DATA.grids.Nj
-                for ip = 1:DATA.grids.Np
-                    if((ip-1)+2*(ij-1) == deg)
-                        % Check if the pair is already added
-                        check_ = HL_deg == [DATA.grids.Parray(ip) DATA.grids.Jarray(ij)];
-                        check_ = sum(check_(:,1) .* check_(:,2));
-                        if ~check_ % if not add it
-                            HL_deg(im,1) = DATA.grids.Parray(ip);
-                            HL_deg(im,2) = DATA.grids.Jarray(ij);
-                            ticks_labels{im} = ['(',num2str(DATA.grids.Parray(ip)),',',num2str(DATA.grids.Jarray(ij)),')'];
-                            im = im + 1; FOUND = 1;
-                        end
-                    end
+    plotname = [logname,'$\langle\sum_k |N_',species_name{ia},'^{pj}|\rangle_{p+2j=const}$'];
+    %We order the moments (0,0) (1,0) (2,0) (0,1) (3,0) (1,1) (4,0) (2,1) (0,2) etc.
+    Nmoments = numel(Napjz(:,:,1)); % total number of moments
+    HL_deg   = zeros(Nmoments,2);   % list of degrees, first column Hermite, second Laguerre
+    im  = 2;
+    deg = 1; % the zero degree is always here first place
+    ticks_labels = cell(10,1);
+    ticks_labels{1} = '(0,0)';
+    while(im<=Nmoments)
+    FOUND = 1;
+        while(FOUND) % As we find a pair matching the degree we retry
+        FOUND = 0;
+        for ij = 1:DATA.grids.Nj
+            for ip = 1:DATA.grids.Np
+                if((ip-1)+2*(ij-1) == deg)
+                    % Check if the pair is already added
+                    check_ = HL_deg == [DATA.grids.Parray(ip) DATA.grids.Jarray(ij)];
+                    check_ = sum(check_(:,1) .* check_(:,2));
+                    if ~check_ % if not add it
+                        HL_deg(im,1) = DATA.grids.Parray(ip);
+                        HL_deg(im,2) = DATA.grids.Jarray(ij);
+                        ticks_labels{im} = ['(',num2str(DATA.grids.Parray(ip)),',',num2str(DATA.grids.Jarray(ij)),')'];
+                        im = im + 1; FOUND = 1;
                     end
                 end
+                end
             end
-            % No pair found anymore, increase the degree
-            deg = deg + 1;
-        end
-
-        % form an axis of velocity ordered moments
-        p2j = 1:Nmoments; yname = '$(P,J)$';
-        % space time of moments amplitude wrt to p+2j degree
-        Na_ST = zeros(Nmoments,numel(Time_));
-        for i_ = 1:Nmoments
-           Na_ST(i_,:) = Napjz((HL_deg(i_,1)/DATA.grids.deltap)+1,HL_deg(i_,2)+1,:); 
         end
-    end    
+        % No pair found anymore, increase the degree
+        deg = deg + 1;
+    end
+    % form an axis of velocity ordered moments
+    p2j = 1:Nmoments; yname = '$(P,J)$';
+    % space time of moments amplitude wrt to p+2j degree
+    Na_ST = zeros(Nmoments,numel(Time_));
+    for i_ = 1:Nmoments
+       Na_ST(i_,:) = Napjz((HL_deg(i_,1)/DATA.grids.deltap)+1,HL_deg(i_,2)+1,:); 
+    end  
     % plots
     % scaling
     if OPTIONS.NORMALIZED
-    plt = @(x,i) x(i,:)./max(x(i,:));
+    plt = @(x,i) squeeze(x(i,:)./max(x(i,:)));
     else
-    plt = @(x,i) x;
-    end
-
-    imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j))); 
-    set(gca,'YDir','normal')        
-    xlabel('$t$'); ylabel(yname);
-    if ~isempty(ticks_labels)
-        yticks(p2j);
-        yticklabels(ticks_labels)
+    plt = @(x,i) squeeze(x(i,:));
     end
+    if OPTIONS.ST
+        imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j))); 
+        set(gca,'YDir','normal')        
+        xlabel('$t$'); ylabel(yname);
+        if ~isempty(ticks_labels)
+            yticks(p2j);
+            yticklabels(ticks_labels)
+        end
+    else
+        colors_ = jet(numel(p2j));
+        for i = 1:numel(p2j)
+           semilogy(Time_,plt(Na_ST,i),...
+               'DisplayName',ticks_labels{i},...
+               'color',colors_(i,:)); hold on;
+        end
     title(plotname)
-
 end
 suptitle(DATA.paramshort)
 
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 5056f61e4e709e2f65f094e9dd41484c57c8bfb0..43a3713a1c24cdb63d077c38b9d0e56dc72311fb 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -24,7 +24,8 @@ for i = 1:numel(data.outfilenames)
     STETC = [ STETC; h5read(outfilename,'/profiler/Tc_step')];
     TS0TC = [ TS0TC; h5read(outfilename,'/profiler/time')];
 end
-
+CPUTI = CPUTI(end);
+DTSIM = mean(DTSIM);
 N_T          = 12;
 
 missing_Tc   = STETC - RHSTC - ADVTC - GHOTC - CLOTC ...
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index b3aa355a68f1fb929cdf35527a539bf252465c4b..a968d3db896445b50f2370645eb044d329366f29 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -47,14 +47,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
 %Paper 2
 % folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
 % folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
-% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
+folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
 % folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_00/';
 % folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_01/';
 
 % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/';
 % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_01/';
 % folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
-folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
+% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
 % folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
 % folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
 % folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
@@ -69,10 +69,12 @@ folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
 gene_data = load_gene_data(folder);
 gene_data.FIGDIR = folder;
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
-gene_data.Pmaxi = gene_data.Nvp-1;
-gene_data.Jmaxi = gene_data.Nmu-1;
+gene_data.grids.Np = gene_data.grids.Nvp-1;
+gene_data.grids.Nj = gene_data.grids.Nmu-1;
 gene_data.CODENAME = 'GENE';
-
+gene_data.inputs = gene_data.grids;
+gene_data.inputs.Na = 1;
+gene_data.paramshort = gene_data.params_string;
 %% Dashboard (Compilation of main plots of the sim)
 dashboard(gene_data);
 
@@ -86,7 +88,7 @@ options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag
 options.INTERP   = 1;
 options.NCUT     = 4;              % Number of cuts for averaging and error estimation
 options.RESOLUTION = 256;
-fig = plot_radial_transport_and_spacetime(gene_data,options,'GENE');
+fig = plot_radial_transport_and_spacetime(gene_data,options);
 save_figure(gene_data,fig,'.png')
 end
 
diff --git a/wk/benchmark and scan scripts/Ajay_scan_CH4_lin_ITG.m b/wk/benchmark and scan scripts/old/Ajay_scan_CH4_lin_ITG.m
similarity index 100%
rename from wk/benchmark and scan scripts/Ajay_scan_CH4_lin_ITG.m
rename to wk/benchmark and scan scripts/old/Ajay_scan_CH4_lin_ITG.m
diff --git a/wk/benchmark and scan scripts/CBC_P_J_scan.m b/wk/benchmark and scan scripts/old/CBC_P_J_scan.m
similarity index 100%
rename from wk/benchmark and scan scripts/CBC_P_J_scan.m
rename to wk/benchmark and scan scripts/old/CBC_P_J_scan.m
diff --git a/wk/benchmark and scan scripts/CBC_hypcoll_PJ_scan.m b/wk/benchmark and scan scripts/old/CBC_hypcoll_PJ_scan.m
similarity index 100%
rename from wk/benchmark and scan scripts/CBC_hypcoll_PJ_scan.m
rename to wk/benchmark and scan scripts/old/CBC_hypcoll_PJ_scan.m
diff --git a/wk/benchmark and scan scripts/CBC_kT_PJ_scan.m b/wk/benchmark and scan scripts/old/CBC_kT_PJ_scan.m
similarity index 100%
rename from wk/benchmark and scan scripts/CBC_kT_PJ_scan.m
rename to wk/benchmark and scan scripts/old/CBC_kT_PJ_scan.m
diff --git a/wk/benchmark and scan scripts/CBC_kT_nu_scan.m b/wk/benchmark and scan scripts/old/CBC_kT_nu_scan.m
similarity index 100%
rename from wk/benchmark and scan scripts/CBC_kT_nu_scan.m
rename to wk/benchmark and scan scripts/old/CBC_kT_nu_scan.m
diff --git a/wk/benchmark and scan scripts/CBC_nu_PJ_scan.m b/wk/benchmark and scan scripts/old/CBC_nu_PJ_scan.m
similarity index 100%
rename from wk/benchmark and scan scripts/CBC_nu_PJ_scan.m
rename to wk/benchmark and scan scripts/old/CBC_nu_PJ_scan.m
diff --git a/wk/benchmark and scan scripts/lin_3D_Zpinch.m b/wk/benchmark and scan scripts/old/lin_3D_Zpinch.m
similarity index 100%
rename from wk/benchmark and scan scripts/lin_3D_Zpinch.m
rename to wk/benchmark and scan scripts/old/lin_3D_Zpinch.m
diff --git a/wk/benchmark and scan scripts/lin_ETPY.m b/wk/benchmark and scan scripts/old/lin_ETPY.m
similarity index 100%
rename from wk/benchmark and scan scripts/lin_ETPY.m
rename to wk/benchmark and scan scripts/old/lin_ETPY.m
diff --git a/wk/benchmark and scan scripts/lin_ITG.m b/wk/benchmark and scan scripts/old/lin_ITG.m
similarity index 100%
rename from wk/benchmark and scan scripts/lin_ITG.m
rename to wk/benchmark and scan scripts/old/lin_ITG.m
diff --git a/wk/benchmark and scan scripts/lin_KBM.m b/wk/benchmark and scan scripts/old/lin_KBM.m
similarity index 100%
rename from wk/benchmark and scan scripts/lin_KBM.m
rename to wk/benchmark and scan scripts/old/lin_KBM.m
diff --git a/wk/benchmark and scan scripts/lin_MTM.m b/wk/benchmark and scan scripts/old/lin_MTM.m
similarity index 100%
rename from wk/benchmark and scan scripts/lin_MTM.m
rename to wk/benchmark and scan scripts/old/lin_MTM.m
diff --git a/wk/benchmark and scan scripts/lin_RHT.m b/wk/benchmark and scan scripts/old/lin_RHT.m
similarity index 100%
rename from wk/benchmark and scan scripts/lin_RHT.m
rename to wk/benchmark and scan scripts/old/lin_RHT.m
diff --git a/wk/benchmark and scan scripts/lin_TEM.m b/wk/benchmark and scan scripts/old/lin_TEM.m
similarity index 100%
rename from wk/benchmark and scan scripts/lin_TEM.m
rename to wk/benchmark and scan scripts/old/lin_TEM.m
diff --git a/wk/dashboard.m b/wk/dashboard.m
index 6c4f65b69c028cad173114e510fddc2d97595a24..7937ef7697747be45643c452e5f694f558542d07 100644
--- a/wk/dashboard.m
+++ b/wk/dashboard.m
@@ -22,7 +22,7 @@ options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag
 options.INTERP   = 1;
 options.NCUT     = 4;              % Number of cuts for averaging and error estimation
 options.RESOLUTION = 256;
-PLOT = plot_radial_transport_and_spacetime(DATA,options,'GENE');
+PLOT = plot_radial_transport_and_spacetime(DATA,options);
 % Put on full fig
 axcp = copyobj(PLOT.ax1,full_fig); 
 set(axcp,'Position',get(axes(1),'position'));delete(axes(1));
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 6f3841cca6ad5f997056f69ce15c059a3f5fc1f6..af8a6e41f1a863aead95f9067732c1288ca0119e 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -18,14 +18,25 @@ PARTITION = '/home/ahoffman/gyacomo/';
 % resdir = 'paper_2_GYAC23/precision_study/test_3x2x128x64x24_sp_muz_2.0';
 % resdir = 'paper_2_GYAC23/precision_study/3x2x128x64x24_sp_clos_1';
 
-%% 
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0';
-% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_full_NL';
+%%
 % resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_muxy_0';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_SG';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_full_NL';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/7x4x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/9x5x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/11x6x128x64x24_dp';
 
+% resdir = 'paper_2_GYAC23/collisionless/CBC/5x3x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/7x4x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/11x6x128x64x24_dp';
+% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x192x96x32_dp';
 %% testcases
-resdir = 'testcases/zpinch_example';
+% resdir = 'testcases/zpinch_example';
 % resdir = 'testcases/cyclone_example';
+resdir = 'testcases/DLRA_zpinch/base_case';
+% resdir = 'testcases/DLRA_zpinch/nsv_filter_2';
+% resdir = 'testcases/DLRA_zpinch/nsv_filter_6';
 
  %%
 J0 = 00; J1 = 10;
@@ -35,7 +46,7 @@ DATADIR = [PARTITION,resdir,'/'];
 data    = {};
 data    = compile_results_low_mem(data,DATADIR,J0,J1);
 
-if 0
+if 1
 %% Plot transport and phi radial profile
 [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
 
@@ -49,10 +60,11 @@ options.INTERP   = 0;
 options.RESOLUTION = 256;
 fig = plot_radial_transport_and_spacetime(data,options);
 end
+
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-options.INTERP    = 0;
+options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.NAME      = '\phi';
 % options.NAME      = '\omega_z';
@@ -67,28 +79,45 @@ options.PLAN      = 'xy';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 % options.TIME      = data.Ts5D(end-30:end);
-% options.TIME      =  data.Ts3D;
-options.TIME      = [0:1500];
+options.TIME      =  data.Ts3D;
+% options.TIME      = [0:1500];
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 options.RESOLUTION = 256;
 create_film(data,options,'.gif')
 end
 
+if 0
+%% fields snapshots
+% Options
+[data.Na00, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Na00');
+
+data.Ni00 = reshape(squeeze(data.Na00(1,:,:,:,:)),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+options.INTERP    = 0;
+options.POLARPLOT = 0;
+options.AXISEQUAL = 0;
+options.NORMALIZE = 0;
+options.NAME      = 'N_i^{00}';
+% options.NAME      = '\phi';
+options.PLAN      = 'xy';
+options.COMP      = 'avg';
+options.TIME      = [800 900 1000];
+options.RESOLUTION = 256;
+fig = photomaton(data,options);
+% save_figure(data,fig)
+end
 if 0
 %% Performance profiler
 profiler(data)
 end
 
-if 1
+if 0
 %% Hermite-Laguerre spectrum
 [data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
-% data.Nipjz = log(data.Nipjz);
-% options.TIME = 'avg';
-options.P2J        = 0;
-options.ST         = 1;
+% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
+options.ST         = 0;
 options.NORMALIZED = 0;
-options.TIME       = [500:800];
+options.LOGSCALE   = 1;
 fig = show_moments_spectrum(data,options);
 % fig = show_napjz(data,options);
 % save_figure(data,fig,'.png');
@@ -109,4 +138,19 @@ options.ik     = 1; % sum, max or index
 options.fftz.flag = 0;
 fig = mode_growth_meter(data,options);
 % save_figure(data,fig,'.png')
+end
+
+
+if 0
+%% Study singular values
+[data.SV_ky_pj, data.Ts2D] = compile_results_2D(DATADIR,J0,J1,'sv_ky_pj');
+nSV = data.grids.Np * data.grids.Nj;
+colors_ = jet(nSV);
+figure
+for i = 1:nSV
+    sv = squeeze(data.SV_ky_pj(i,:));
+    semilogy(data.Ts2D,sv,...
+        'color',colors_(i,:),'DisplayName',['SV ',num2str(i)]);hold on
+end
+legend('show');
 end
\ No newline at end of file