diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 1fcc187f8d6070ff4eb27493bd01190d8c9317ff..30ab59ffbb14c46218425227c3b91a1f3176213e 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -48,12 +48,12 @@ for i = 1:2
                 MODESTR = '$\max_{k_y}$';
             otherwise
                 plt_ = @(x,ik) movmean(abs(squeeze(x(OPTIONS.ik,ik,FRAMES))),Nma);
-                MODESTR = ['$k_y=$',num2str(DATA.ky(OPTIONS.ik))];
+                MODESTR = ['$k_y=$',num2str(DATA.grids.ky(OPTIONS.ik))];
         end
         kstr = 'k_x';
         % Number max of modes to plot is kx>0 (1/2) of the non filtered modes (2/3)
-        Nmax = ceil(DATA.Nkx*1/3);
-        k = DATA.kx;
+        Nmax = ceil(DATA.grids.Nkx*1/3);
+        k = DATA.grids.kx;
     elseif MODES_SELECTOR == 1
         switch OPTIONS.ik
             case 'sum'
@@ -64,12 +64,12 @@ for i = 1:2
                 MODESTR = '$\max_{k_x}$';
             otherwise
                 plt_ = @(x,ik) movmean(abs(squeeze(x(ik,OPTIONS.ik,FRAMES))),Nma);
-                MODESTR = ['$k_x=$',num2str(DATA.kx(OPTIONS.ik))];
+                MODESTR = ['$k_x=$',num2str(DATA.grids.kx(OPTIONS.ik))];
         end
         kstr = 'k_y';
         % Number max of modes to plot is ky>0 (1/1) of the non filtered modes (2/3)
-        Nmax = ceil(DATA.Nky*2/3);
-        k = DATA.ky;
+        Nmax = ceil(DATA.grids.Nky*2/3);
+        k = DATA.grids.ky;
     end 
     if NORMALIZED
         plt = @(x,ik) plt_(x,ik)./max(plt_(x,ik));
@@ -120,7 +120,7 @@ for i = 1:2
 %     subplot(2+d,3,3+3*(i-1))
     FIGURE.axes(3+3*(i-1)) = subplot(2+d,3,3+3*(i-1),'parent',FIGURE.fig);
         plot(k(MODES),gamma,...
-                'DisplayName',['(',num2str(DATA.Pmaxi-1),',',num2str(DATA.Jmaxi-1),')']); hold on;
+                'DisplayName',['(',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']); hold on;
         for i_ = 1:numel(mod2plot)
             plot(k(MODES(mod2plot(i_))),gamma(mod2plot(i_)),'x','color',clr_(i_,:));
         end
@@ -132,10 +132,10 @@ for i = 1:2
 end
 
 if d
-    [~,ikx] = min(abs(DATA.kx-OPTIONS.fftz.kx));
-    [~,iky] = min(abs(DATA.ky-OPTIONS.fftz.ky));
+    [~,ikx] = min(abs(DATA.grids.kx-OPTIONS.fftz.kx));
+    [~,iky] = min(abs(DATA.grids.ky-OPTIONS.fftz.ky));
     sz_=size(DATA.PHI);nkz = sz_(3)/2;
-    k = [(0:nkz/2), (-nkz/2+1):-1]/DATA.Npol;
+    k = [(0:nkz/2), (-nkz/2+1):-1]/DATA.grids.Npol;
     % Spectral treatment of z-dimension
     Y = fft(DATA.PHI(iky,ikx,:,:),[],3);
     phi = squeeze(Y(1,1,2:2:end,:)); 
@@ -158,7 +158,7 @@ if d
         set(gca,'YDir','normal')
     %     xlim([t(1) t(end)]); %ylim([1e-5 1])
         xlabel('$k_\parallel$'); ylabel('$t c_s /\rho_s$');
-        title(['$k_x=$',num2str(DATA.kx(ikx)),', $k_y=$',num2str(DATA.ky(iky))]);  
+        title(['$k_x=$',num2str(DATA.grids.kx(ikx)),', $k_y=$',num2str(DATA.grids.ky(iky))]);  
 
     subplot(3,3,8)
         for i_ = 1:numel(mod2plot)
@@ -174,7 +174,7 @@ if d
 
     subplot(3,3,9)
         plot(k(MODES),gamma,...
-                'DisplayName',['(',num2str(DATA.Pmaxi-1),',',num2str(DATA.Jmaxi-1),')']); hold on;
+                'DisplayName',['(',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']); hold on;
         for i_ = 1:numel(mod2plot)
             plot(k(MODES(mod2plot(i_))),gamma(mod2plot(i_)),'x','color',clr_(i_,:));
         end
@@ -185,5 +185,5 @@ if d
         title('Growth rates')   
     
 end
-
+suptitle(DATA.paramshort)
 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 ea3275d81288033d6b047794774194aee2d8fc93..cac486ac58f1d3c5863e3f0a75f2853d79fa4000 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -43,14 +43,25 @@ while(CONTINUE)
             BETA      = h5readatt(filename,'/data/input/model','beta');
 
             if W_GAMMA
-                [ GGAMMA_R, Ts0D, ~] = load_0D_data(filename, 'gflux_x');
-                PGAMMA_R            = load_0D_data(filename, 'pflux_x');
+                try
+                    [ GGAMMA_R, Ts0D, ~] = load_0D_data(filename, 'gflux_x');
+                    PGAMMA_R             = load_0D_data(filename, 'pflux_x');
+                catch
+                    % old version
+                    [ 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
             end
 
             if W_HF
-                [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_x');
+                try
+                    [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_x');
+                catch
+                    % old version
+                    [ HFLUX_X, Ts0D, ~] = load_0D_data(filename, 'hflux_xi');
+                end
                 HFLUX_ = cat(1,HFLUX_,HFLUX_X); clear HFLUX_X
             end
 
diff --git a/matlab/load/load_0D_data.m b/matlab/load/load_0D_data.m
index ee55c6fa0f56ae6657c31f525b57b649ffa720c2..a4d20bd59a2f8d139eda8e338f57bc6a311db639 100644
--- a/matlab/load/load_0D_data.m
+++ b/matlab/load/load_0D_data.m
@@ -1,9 +1,10 @@
 function [ data, time, dt ] = load_0D_data( filename, variablename )
 %LOAD_0D_DATA load a 0D variable stored in a hdf5 result file from HeLaZ
-    time     = h5read(filename,'/data/var0d/time');
+    time  = h5read(filename,'/data/var0d/time');
     dt    = h5readatt(filename,'/data/input/basic','dt');
+    Na    = h5readatt(filename,'/data/input/model','Na');
     cstart= h5readatt(filename,'/data/input/basic','start_iframe0d'); 
-    data     = h5read(filename,['/data/var0d/',variablename]);
-
+    data  = h5read(filename,['/data/var0d/',variablename]);
+%     data  = reshape(data,[Na, numel(time)]);
 end
 
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 0234d379eabedebbec1aee428fb04d1ef101e966..58305dae1f943add87b63787eae14dbfcd4ca686 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -26,11 +26,12 @@ clr_ = lines(20);
 mvm = @(x) movmean(x,OPTIONS.NMVA);
     FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.params_string]; %set(gcf, 'Position',  [500, 1000, 1000, 600])
     FIGURE.ax1 = subplot(2,1,1,'parent',FIGURE.fig);
-        plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
-            'color',clr_((DATA.grids.Np-1)/2,:),...
+    for ia = 1:DATA.inputs.Na
+        plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI(ia,:)*SCALE),'--',...
+            'color',clr_((DATA.grids.Np-1)/2+(ia-1),:),...
             'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
-        plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'-',...
-            'color',clr_((DATA.grids.Np-1)/2,:),...
+        plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X(ia,:)*SCALE),'-',...
+            'color',clr_((DATA.grids.Np-1)/2+(ia-1),:),...
             'DisplayName',['$Q_x$ ',DATA.paramshort]); hold on;
         ylabel('Transport')  
         if(~isnan(Qx_infty_avg))
@@ -44,6 +45,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
             catch
             end
         end
+    end
         xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
     grid on; set(gca,'xticklabel',[]); 
     title({DATA.param_title,...
@@ -82,4 +84,6 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
         subplot(311)
         plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
     end
+    suptitle(DATA.paramshort)
+
 end
\ No newline at end of file
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index 1f93b35d707c9c5e696cb7f99549f952ca7208d8..02291c8e2e11b50772b0763ae61e5f833c4ecb39 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -6,9 +6,16 @@ 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.LOGSCALE
+%     compress = @(x,ia) log(sum(abs(squeeze(x(ia,:,:,:))),3));
+    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);
+end
 for ia = 1:DATA.inputs.Na
-    Napjz  = sum(abs(squeeze(DATA.Napjz(ia,:,:,:,:))),3);
+%     Napjz  = sum(abs(squeeze(DATA.Napjz(ia,:,:,:,:))),3);
+    Napjz  =compress(DATA.Napjz);
     subplot(double(DATA.inputs.Na),1,double(ia))
     if OPTIONS.P2J
         plotname = ['$\langle\sum_k |N_',species_name{ia},'^{pj}|\rangle_{p+2j=const}$'];
@@ -92,5 +99,7 @@ for ia = 1:DATA.inputs.Na
     title(plotname)
 
 end
+suptitle(DATA.paramshort)
+
 end