diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index b58fcb545ab3a36d8f7d7d3228032ea3d3992f6b..1d2f1bb101ba5f4fe1b3283a9db73b1d8044fe60 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -238,7 +238,7 @@ save_figure
 end
 
 %%
-t0    = 50;
+t0    = 0;
 skip_ = 1; 
 DELAY = 0.01*skip_;
 FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
diff --git a/wk/marconi_scaling.m b/wk/marconi_scaling.m
index 624b85d096bf7a02a30095dbbabb7d5156224dc6..b0ea6a977b83160a26888c023cb50933d3884cd3 100644
--- a/wk/marconi_scaling.m
+++ b/wk/marconi_scaling.m
@@ -7,9 +7,9 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '12:00:00'; % allocation time hh:mm:ss
 CLUSTER.NODES = '1';        % MPI process
 CLUSTER.CPUPT = '1';        % CPU per task
-CLUSTER.NTPN  = '32';       % N tasks per node
+CLUSTER.NTPN  = '20';       % N tasks per node
 CLUSTER.PART  = 'prod';      % dbg or prod
-CLUSTER.MEM   = '10GB';     % Memory
+CLUSTER.MEM   = '16GB';     % Memory
 %% PHYSICAL PARAMETERS
 NU      = 1e-1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
@@ -21,14 +21,14 @@ NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
 N       = 1024;     % Frequency gridpoints (Nkr = N/2)
 L       = 100;     % Size of the squared frequency domain
-PMAXE   = 6;     % Highest electron Hermite polynomial degree
-JMAXE   = 4;     % Highest ''       Laguerre ''
-PMAXI   = 6;     % Highest ion      Hermite polynomial degree
-JMAXI   = 4;     % Highest ''       Laguerre ''
+PMAXE   = 1;     % Highest electron Hermite polynomial degree
+JMAXE   = 1;     % Highest ''       Laguerre ''
+PMAXI   = 1;     % Highest ion      Hermite polynomial degree
+JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS 
-TMAX    = 1;  % Maximal time unit
-DT      = 1e-2;   % Time step
-SPS0D   = 0;    % Sampling per time unit for profiler
+TMAX    = 5;  % Maximal time unit
+DT      = 5e-2;   % Time step
+SPS0D   = 1/DT;    % Sampling per time unit for profiler
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
 SPS5D   = 0;    % Sampling per time unit for 5D arrays
 RESTART = 0;      % To restart from last checkpoint
diff --git a/wk/parallel_scaling.m b/wk/parallel_scaling.m
index 4aa5c842d69bdd96508dddd58fd6341979fc2959..b1502db684cb808a733fe096f653c1e7f9ccc102 100644
--- a/wk/parallel_scaling.m
+++ b/wk/parallel_scaling.m
@@ -10,21 +10,21 @@ Results_256_21.time  = [2450, 1346, 0680, 0389,  0323,  0307];
 Results_512_21.np    = [   1,    2,    4,    8,    16,   20,   24];
 Results_512_21.time  = [3429, 1680, 0842, 0443,  0292, 0322, 0362];
 
-% Handwritten results for 512x256, P,J=3,2, Tmax = 2
-Results_512_32.np    = [   1,    2,    4,    8,    16,   20];
-Results_512_32.time  = [4450, 2267, 1136, 0595,  0363, 0323];
+% Handwritten results for 512x256, P,J=3,2, Tmax = 2, mu=0, dt 1e-2?
+Results_512_32.np    = [   1,    2,    4,    8,    16,   20,   24];
+Results_512_32.time  = [4450, 2267, 1136, 0595,  0363, 0323, 0000];
 
 % Handwritten results for 1024x512, P,J=1,1, Tmax = 5 dt = 0.05, mu = 0
-Results_1024_11.np    = [   1,    2,    4,    6,    8,   12,   16];
-Results_1024_11.time  = [1568, 1046, 0490, 0347, 0257, 0221,  219];
+Results_1024_11.np    = [   1,    2,    4,    6,    8,   12,   16,   20,   24];
+Results_1024_11.time  = [1568, 1046, 0490, 0347, 0257, 0221, 0219, 0000,  0000];
 
 % Handwritten results for 1024x512, P,J=2,2, Tmax = 2 dt = 0.05, mu = 0
 Results_1024_22.np    = [   1,    2,    4,    6,    8,   10,   12,   16,   20];
 Results_1024_22.time  = [2391, 1373, 0654, 0457, 0343, 0297, 0274, 0219, 0206];
 
 % Handwritten results for 1024x512, P,J=6,4, Tmax = 2 dt = 0.05, mu = 0
-Results_1024_64.np    = [   1,    2,    4,    6,    8,   10,   12,   16];
-Results_1024_64.time  =[38957,22675,10886, 7424, 5768, 0000, 3904, 2947];
+Results_1024_64.np    = [   1,    2,    4,    6,    8,   10,   12,   16,   20,  24];
+Results_1024_64.time  =[38957,22675,10886, 7424, 5768, 0000, 3904, 2947, 2417,2075];
 
 %
 fig = figure;
@@ -50,7 +50,7 @@ title('Strong scaling')
 grid on  
     
 
-FIGNAME = [SIMDIR,'strong_scaling.png'];
+FIGNAME = '/home/ahoffman/HeLaZ/results/strong_scaling_old.pdf';
 saveas(fig,FIGNAME);
 disp(['Figure saved @ : ',FIGNAME])
 
@@ -88,7 +88,7 @@ title('Weak scaling')
 grid on  
     
 
-FIGNAME = [SIMDIR,'weak_scaling.png'];
+FIGNAME = '/home/ahoffman/HeLaZ/results/weak_scaling.pdf';
 saveas(fig,FIGNAME);
 disp(['Figure saved @ : ',FIGNAME])
 end
\ No newline at end of file
diff --git a/wk/parallel_scaling_new.m b/wk/parallel_scaling_new.m
new file mode 100644
index 0000000000000000000000000000000000000000..47def3cdc8082dc168676efb1e3ee193c095d2b8
--- /dev/null
+++ b/wk/parallel_scaling_new.m
@@ -0,0 +1,51 @@
+default_plots_options
+
+%% Strong scaling measurement
+
+% Handwritten results for 512x256, P,J=2,1, Tmax = 5
+Results_512_21.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
+Results_512_21.time  = [0000, 0000, 0000, 0000,  0000, 0000, 0000, 0000];
+
+% Handwritten results for 512x256, P,J=3,2, Tmax = 2, mu=0, dt 1e-2?
+Results_512_32.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
+Results_512_32.time  = [0000, 0000, 0000, 0000,  0000, 0000, 0000, 0000];
+
+% Handwritten results for 1024x512, P,J=1,1, Tmax = 5 dt = 0.05, mu = 0
+Results_1024_11.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
+Results_1024_11.time  = [0000, 0000, 0000, 0000,  0000, 0000, 0000, 0000];
+
+% Handwritten results for 1024x512, P,J=2,2, Tmax = 2 dt = 0.05, mu = 0
+Results_1024_22.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
+Results_1024_22.time  = [0000, 0000, 0000, 0000,  0000, 0000, 0000, 0000];
+
+% Handwritten results for 1024x512, P,J=6,4, Tmax = 2 dt = 0.05, mu = 0
+Results_1024_64.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
+Results_1024_64.time  = [0000, 0000, 0000, 0000,  0000, 0000, 0000, 0000];
+
+%
+fig = figure;
+
+plot(1:24,1:24,'-k','DisplayName','Ideal')
+hold on
+% res = Results_256_21;
+% plot(res.np,res.time(1)./(res.time),'o--','DisplayName','$256\times128$, $P,J=2,1$');
+res = Results_512_21;
+plot(res.np,res.time(1)./(res.time),'v-','DisplayName','$512\times256$, $P,J=2,1$');
+res = Results_512_32;
+plot(res.np,res.time(1)./(res.time),'>-','DisplayName','$512\times256$, $P,J=3,2$');
+res = Results_1024_11;
+plot(res.np,res.time(1)./(res.time),'o-','DisplayName','$1024\times512$, $P,J=1,1$');
+res = Results_1024_22;
+plot(res.np,res.time(1)./(res.time),'s-','DisplayName','$1024\times512$, $P,J=2,2$');xlim([1,max(res.np)]);
+res = Results_1024_64;
+plot(res.np,res.time(1)./(res.time),'d-','DisplayName','$1024\times512$, $P,J=6,4$');xlim([1,max(res.np)]);
+xlabel('$N_p$'); ylabel('speedup')
+xlim([1,24]); ylim([1,24])
+legend('show')
+title('Strong scaling')
+grid on  
+    
+
+FIGNAME = '/home/ahoffman/HeLaZ/results/strong_scaling_new.pdf';
+saveas(fig,FIGNAME);
+disp(['Figure saved @ : ',FIGNAME])
diff --git a/wk/parallel_scaling_old.m b/wk/parallel_scaling_old.m
new file mode 100644
index 0000000000000000000000000000000000000000..b1502db684cb808a733fe096f653c1e7f9ccc102
--- /dev/null
+++ b/wk/parallel_scaling_old.m
@@ -0,0 +1,94 @@
+default_plots_options
+
+%% Strong scaling measurement
+
+% Handwritten results for 256x128, P,J=2,1, Tmax = 20
+Results_256_21.np    = [   1,    2,    4,    8,    10,    12];
+Results_256_21.time  = [2450, 1346, 0680, 0389,  0323,  0307];
+
+% Handwritten results for 512x256, P,J=2,1, Tmax = 5
+Results_512_21.np    = [   1,    2,    4,    8,    16,   20,   24];
+Results_512_21.time  = [3429, 1680, 0842, 0443,  0292, 0322, 0362];
+
+% Handwritten results for 512x256, P,J=3,2, Tmax = 2, mu=0, dt 1e-2?
+Results_512_32.np    = [   1,    2,    4,    8,    16,   20,   24];
+Results_512_32.time  = [4450, 2267, 1136, 0595,  0363, 0323, 0000];
+
+% Handwritten results for 1024x512, P,J=1,1, Tmax = 5 dt = 0.05, mu = 0
+Results_1024_11.np    = [   1,    2,    4,    6,    8,   12,   16,   20,   24];
+Results_1024_11.time  = [1568, 1046, 0490, 0347, 0257, 0221, 0219, 0000,  0000];
+
+% Handwritten results for 1024x512, P,J=2,2, Tmax = 2 dt = 0.05, mu = 0
+Results_1024_22.np    = [   1,    2,    4,    6,    8,   10,   12,   16,   20];
+Results_1024_22.time  = [2391, 1373, 0654, 0457, 0343, 0297, 0274, 0219, 0206];
+
+% Handwritten results for 1024x512, P,J=6,4, Tmax = 2 dt = 0.05, mu = 0
+Results_1024_64.np    = [   1,    2,    4,    6,    8,   10,   12,   16,   20,  24];
+Results_1024_64.time  =[38957,22675,10886, 7424, 5768, 0000, 3904, 2947, 2417,2075];
+
+%
+fig = figure;
+
+plot(1:24,1:24,'-k','DisplayName','Ideal')
+hold on
+% res = Results_256_21;
+% plot(res.np,res.time(1)./(res.time),'o--','DisplayName','$256\times128$, $P,J=2,1$');
+res = Results_512_21;
+plot(res.np,res.time(1)./(res.time),'v-','DisplayName','$512\times256$, $P,J=2,1$');
+res = Results_512_32;
+plot(res.np,res.time(1)./(res.time),'>-','DisplayName','$512\times256$, $P,J=3,2$');
+res = Results_1024_11;
+plot(res.np,res.time(1)./(res.time),'o-','DisplayName','$1024\times512$, $P,J=1,1$');
+res = Results_1024_22;
+plot(res.np,res.time(1)./(res.time),'s-','DisplayName','$1024\times512$, $P,J=2,2$');xlim([1,max(res.np)]);
+res = Results_1024_64;
+plot(res.np,res.time(1)./(res.time),'d-','DisplayName','$1024\times512$, $P,J=6,4$');xlim([1,max(res.np)]);
+xlabel('$N_p$'); ylabel('speedup')
+xlim([1,24]); ylim([1,24])
+legend('show')
+title('Strong scaling')
+grid on  
+    
+
+FIGNAME = '/home/ahoffman/HeLaZ/results/strong_scaling_old.pdf';
+saveas(fig,FIGNAME);
+disp(['Figure saved @ : ',FIGNAME])
+
+if 0
+%% Weak scaling
+% Handwritten results for P,J=2,1, Tmax = 5, dt = 0.01, Nz = Nr
+Results_1_64.np    = [   1,    2,    4,    8];
+Results_1_64.Nr    = [  64,   90,  128,  180];
+Results_1_64.time  = [0064, 0074, 0082, 0101];
+
+% Handwritten results for P,J=2,1, Tmax = 5, dt = 0.01, Nz = 128
+Results_1_128.np    = [   1,    2,    4,    8,   16];
+Results_1_128.Nr    = [  32,   64,  128,  256,  512];
+Results_1_128.time  = [0032, 0037, 0043, 0049, 0070];
+
+% Handwritten results for Tmax = 5, dt = 0.05, mu = 0, etab =0, Pi=Ji=Pe=Je=1
+Results_1_128.np    = [   1,    2,    4,    6,   16];
+Results_1_128.N     = [ 256,  360,  512,  720, 1024];
+Results_1_128.time  = [0059, 0072, 0000, 0153, 0070];
+
+%
+fig = figure;
+
+plot(Results_1_64.np,Results_1_64.time,'ob','DisplayName','$256\times128$');
+hold on
+plot(Results_1_64.np,Results_1_64.time(1)*ones(numel(Results_1_64.np)),'--b','DisplayName','Ideal')
+
+plot(Results_1_128.np,Results_1_128.time,'or','DisplayName','$256\times128$');
+plot(Results_1_128.np,Results_1_128.time(1)*ones(numel(Results_1_128.np)),'--r','DisplayName','Ideal')
+
+xlim([1,max(res.np)]);
+xlabel('$N_p$'); ylabel('CPU time [s]') 
+legend('show')
+title('Weak scaling')
+grid on  
+    
+
+FIGNAME = '/home/ahoffman/HeLaZ/results/weak_scaling.pdf';
+saveas(fig,FIGNAME);
+disp(['Figure saved @ : ',FIGNAME])
+end
\ No newline at end of file
diff --git a/wk/test_parallel.m b/wk/test_parallel.m
index 2eb69bbc1a29e3ffce633b592f9c95b791abb637..70eadc474a6c2e72b7ddc8c2c115866cacdc9832 100644
--- a/wk/test_parallel.m
+++ b/wk/test_parallel.m
@@ -12,22 +12,22 @@ ETAT    = 0.0;    % Temperature gradient
 MU      = 5e-4;   % Hyper diffusivity coefficient
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 512;     % Frequency gridpoints (Nkr = N/2)
-L       = 100;     % Size of the squared frequency domain
+N       = 128;     % Frequency gridpoints (Nkr = N/2)
+L       = 33;     % Size of the squared frequency domain
 PMAXE   = 2;     % Highest electron Hermite polynomial degree
 JMAXE   = 1;     % Highest ''       Laguerre ''
 PMAXI   = 2;     % Highest ion      Hermite polynomial degree
 JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS 
-TMAX    = 100;  % Maximal time unit
-DT      = 1e-2;   % Time step
-SPS0D   = 1;    % Sampling per time unit for profiler
-SPS2D   = 1;      % Sampling per time unit for 2D arrays
-SPS5D   = 0.1;    % Sampling per time unit for 5D arrays
-RESTART = 1;      % To restart from last checkpoint
+TMAX    = 40;  % Maximal time unit
+DT      = 2e-2;   % Time step
+SPS0D   = 1/DT;    % Sampling per time unit for profiler
+SPS2D   = 2;      % Sampling per time unit for 2D arrays
+SPS5D   = 2;    % Sampling per time unit for 5D arrays
+RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 1;
 %% OPTIONS
-SIMID   = 'test_parallel';  % Name of the simulation
+SIMID   = 'test_sapj_opt';  % Name of the simulation
 CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -40,6 +40,6 @@ KPAR    = 0.0;    % Parellel wave vector component
 LAMBDAD = 0.0; 
 NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 CANCEL_ODD_P = 0;% Cancels the odd polynomials degree
-MARCONI = 0;
+LOAD_MARCONI = 0;
 
 setup