diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 1d2f1bb101ba5f4fe1b3283a9db73b1d8044fe60..dc5198904d5581f68e0b5ac458cb358db1a61010 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,12 +1,12 @@
 %% Load results
-if LOAD_MARCONI
+if LOAD_MARCONI==1
     hostfolder = ['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS];
     localfolder= [BASIC.RESDIR,'..'];
     system(['scp -r ahoffman@login.marconi.cineca.it:',hostfolder,' ',localfolder])
 end
-JOBNUM = 0; load_results;
+% JOBNUM = 0; load_results;
 % JOBNUM = 1; load_results;
-% compile_results
+compile_results
 
 %% Retrieving max polynomial degree and sampling info
 Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
@@ -106,6 +106,7 @@ end
 
 
 %% Compute growth rate
+if NON_LIN == 0
 disp('- growth rate')
 tend   = Ts2D(end); tstart   = 0.6*tend; 
 g_          = zeros(Nkr,Nkz);
@@ -117,7 +118,7 @@ for ikr = 1:Nkr
 end
 % gkr0kz_Ni00 = max(real(g_(:,:)),[],1);
 gkr0kz_Ni00 = real(g_(ikr0KH,:));
-
+end
 %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 default_plots_options
 disp('Plots')
@@ -203,33 +204,35 @@ suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB), sprin
 save_figure
 end
 
-if 1
+if 0
 %% Photomaton : real space
-% FIELD = ni00; FNAME = 'ni';
-FIELD = phi; FNAME = 'phi';
-tf = 60;  [~,it1] = min(abs(Ts2D-tf));
-tf = 65;  [~,it2] = min(abs(Ts2D-tf)); 
-tf = 200; [~,it3] = min(abs(Ts2D-tf));
-tf = 400; [~,it4] = min(abs(Ts2D-tf));
+FIELD = ni00; FNAME = 'ni';
+% FIELD = ne00; FNAME = 'ne';
+% FIELD = phi; FNAME = 'phi';
+tf = 19;  [~,it1] = min(abs(Ts2D-tf));
+tf = 20;  [~,it2] = min(abs(Ts2D-tf)); 
+tf = 21; [~,it3] = min(abs(Ts2D-tf));
+tf = 22; [~,it4] = min(abs(Ts2D-tf));
 fig = figure; FIGNAME = [FNAME,'_snaps']; set(gcf, 'Position',  [100, 100, 1500, 400])
-    subplot(141); plt = @(x) ((x));
+plt = @(x) x;%./max(max(x));
+    subplot(141)
         DATA = plt(FIELD(:,:,it1));
-        pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
         xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]); 
         title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
-    subplot(142); plt = @(x) ((x));
+    subplot(142)
         DATA = plt(FIELD(:,:,it2));
-        pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
         xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); 
         title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
-    subplot(143); plt = @(x) ((x));
+    subplot(143)
         DATA = plt(FIELD(:,:,it3));
-        pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
         xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]); 
         title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
-    subplot(144); plt = @(x) ((x));
+    subplot(144)
         DATA = plt(FIELD(:,:,it4));
-        pclr = pcolor((RR),(ZZ),DATA./max(max(DATA))); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        pclr = pcolor((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
         xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]); 
         title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
 % suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
@@ -363,7 +366,7 @@ end
 %%
 if 0
 %% Show frame in kspace
-tf = 200; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
+tf = 100; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
 fig = figure; FIGNAME = ['krkz_frame',sprintf('t=%.0f',Ts2D(it2))];set(gcf, 'Position',  [100, 100, 700, 600])
     subplot(221); plt = @(x) fftshift((abs(x)),2);
         pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar;
diff --git a/wk/linear_study.m b/wk/linear_study.m
index b84475b511ad811c5466de0c08110e5ee22f485d..60e8ce23d5cb5f390198e0fe5243da512ec77120 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -11,17 +11,17 @@ ETAB    = 0.5;
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
 MU      = 0e-4;   % Hyper diffusivity coefficient
-LAMBDAD = 0.0; 
+LAMBDAD = 0.0;
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
 N       = 50;     % Frequency gridpoints (Nkr = N/2)
 L       = 100;     % Size of the squared frequency domain
 PMAXE   = 12;
-JMAXE   = 6; 
+JMAXE   = 6;
 PMAXI   = 12;
 JMAXI   = 6;
 KREQ0   = 1;      % put kr = 0
-%% TIME PARAMETERS 
+%% TIME PARAMETERS
 TMAX    = 300;  % Maximal time unit
 DT      = 1e-3;   % Time step
 SPS0D   = 0.5;      % Sampling per time unit for 2D arrays
@@ -36,12 +36,10 @@ CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Do
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % unused
-KR0KH   = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
 NO_E    = 0;  % Remove electrons dynamic
 % DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
 JOBNUM = 00;
-CANCEL_ODD_P = 0;% Cancels the odd polynomials degree
-KPAR    = 0.0 * (1-CANCEL_ODD_P);    % Parellel wave vector component
+KPAR    = 0.0;    % Parellel wave vector component
 
 %% PARAMETER SCANS
 if 1
@@ -70,7 +68,7 @@ for i = 1:Nparam
     system('./../bin/helaz');
     % Load and process results
     load_results
-    tend   = Ts2D(end); tstart   = 0.4*tend; 
+    tend   = Ts2D(end); tstart   = 0.4*tend;
     for ikz = 1:N
         gamma_Ni(i,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(1,ikz,:))),tstart,tend);
         Ni00_ST(i,ikz,1:numel(Ts2D)) = squeeze((Ni00(1,ikz,:)));
@@ -83,7 +81,7 @@ end
 % %% Plot
 % fig = figure; FIGNAME = 'space_time_Ni00';
 % i = 1;
-% 
+%
 % [YY, XX] = meshgrid(Ts2D,kz);
 % plt = @(x) squeeze(log(abs(x)))
 % pclr = pcolor(XX,YY,plt(Ni00_ST(i,:,:)));set(pclr, 'edgecolor','none');
@@ -149,7 +147,7 @@ for i = 1:Nparam
     run
     % Load and process results
     load_results
-    tend   = Ts2D(end); tstart   = 0.6*tend; 
+    tend   = Ts2D(end); tstart   = 0.6*tend;
     for ikz = 1:N
         gamma_Ni(i,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(1,ikz,:))),tstart,tend);
     end
@@ -177,4 +175,4 @@ title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',...
 legend('show')
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
 end
-end
\ No newline at end of file
+end
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 0c709cbe32818b159abebd5d4fc45a736c8b9512..4a2ae9df68cd32a7b06d353f2738715857497c4b 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -7,31 +7,32 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
 CLUSTER.NODES = '1';        % MPI process
 CLUSTER.CPUPT = '1';        % CPU per task
-CLUSTER.NTPN  = '16';       % N tasks per node
+CLUSTER.NTPN  = '20';       % N tasks per node
 CLUSTER.PART  = 'prod';      % dbg or prod
 CLUSTER.MEM   = '32GB';     % Memory
 %% PHYSICAL PARAMETERS
 NU      = 1e-1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.5;    % Magnetic gradient
+ETAB    = 0.4;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 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
-PMAXE   = 4;     % Highest electron Hermite polynomial degree
-JMAXE   = 2;     % Highest ''       Laguerre ''
-PMAXI   = 4;     % Highest ion      Hermite polynomial degree
-JMAXI   = 2;     % Highest ''       Laguerre ''
-%% TIME PARAMETERS 
-TMAX    = 250;  % Maximal time unit
+L       = 150;     % 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    = 200;  % Maximal time unit
 DT      = 5e-3;   % Time step
-SPS0D   = 1/DT/4;    % Sampling per time unit for profiler
+SPS0D   = 10;    % 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
+SPS5D   = 1/10;    % Sampling per time unit for 5D arrays
+SPSCP   = 1/10;    % Sampling per time unit for checkpoints
+RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
 SIMID   = 'Marconi';  % Name of the simulation
@@ -44,9 +45,8 @@ NO_E    = 0;  % Remove electrons dynamic
 % DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
 KREQ0   = 0;      % put kr = 0
 KPAR    = 0.0;    % Parellel wave vector component
-LAMBDAD = 0.0; 
+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
 
 %% Run following scripts
 setup
@@ -56,4 +56,4 @@ write_sbash
 system(['scp {fort.90,batch_script.sh,setup_and_run.sh}',...
     ' ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk']);
 LOAD_MARCONI = 1;
-disp('done');
\ No newline at end of file
+disp('done');
diff --git a/wk/marconi_scaling.m b/wk/marconi_scaling.m
index 1aecc4a58613d4d1e06ef878ed9a2bb99c9a9dba..ec5aa14864f524b8884c99c606ac1b46fd82c1a7 100644
--- a/wk/marconi_scaling.m
+++ b/wk/marconi_scaling.m
@@ -4,11 +4,11 @@ addpath(genpath('../matlab')) % ... add
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% CLUSTER PARAMETERS
-CLUSTER.TIME  = '01:30:00'; % allocation time hh:mm:ss
-CLUSTER.NODES = '1';        % MPI process
-CLUSTER.CPUPT = '1';        % CPU per task
-CLUSTER.NTPN  = '4';       % N tasks per node
-CLUSTER.PART  = 'prod';      % dbg or prod
+CLUSTER.TIME  = '01:00:00'; % allocation time hh:mm:ss
+CLUSTER.NODES = '01';       % MPI process
+CLUSTER.CPUPT = '01';       % CPU per task
+CLUSTER.NTPN  = '01';       % N tasks per node (openMP)
+CLUSTER.PART  = 'dbg';      % dbg or prod
 CLUSTER.MEM   = '16GB';     % Memory
 %% PHYSICAL PARAMETERS
 NU      = 1e-1;   % Collision frequency
@@ -19,22 +19,23 @@ ETAT    = 0.0;    % Temperature gradient
 MU      = 0e-4;   % Hyper diffusivity coefficient
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 1024;     % Frequency gridpoints (Nkr = N/2)
+N       = 512;     % Frequency gridpoints (Nkr = N/2)
 L       = 100;     % 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    = 10;  % Maximal time unit
+%% TIME PARAMETERS
+TMAX    = 5;  % Maximal time unit
 DT      = 5e-2;   % Time step
 SPS0D   = 1/DT;    % Sampling per time unit for profiler
-SPS2D   = 1;      % Sampling per time unit for 2D arrays
-SPS5D   = 0;    % Sampling per time unit for 5D arrays
+SPS2D   = 0;      % Sampling per time unit for 2D arrays
+SPS5D   = 0;      % Sampling per time unit for 5D arrays
+SPSCP   = 0;      % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
-SIMID   = ['Scaling_np',num2str(CLUSTER.NTPN)];  % Name of the simulation
+SIMID   = ['Scaling__np',num2str(CLUSTER.NTPN)];  % Name of the simulation
 CO      = -2;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -44,13 +45,12 @@ NO_E    = 0;  % Remove electrons dynamic
 % DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
 KREQ0   = 0;      % put kr = 0
 KPAR    = 0.0;    % Parellel wave vector component
-LAMBDAD = 0.0; 
+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
 
 %% Run following scripts
 setup
 
 write_sbash
 
-MARCONI = 1;
\ No newline at end of file
+MARCONI = 1;
diff --git a/wk/parallel_scaling_new.m b/wk/parallel_scaling_new.m
index 4526a5dbc9439d7b59ecbcb23a32dd7dec01c425..c9608265e7ed5bdc343f27a279464439d2bce4cb 100644
--- a/wk/parallel_scaling_new.m
+++ b/wk/parallel_scaling_new.m
@@ -7,33 +7,34 @@ if 0
 i_ = 1;
 for np = NPS
     SIM_NAME = sprintf('Scaling_np%02d',np);
-    hostfile = ['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/',SIM_NAME,'/',BASIC.PARAMS,'/outputs_00.h5'];
-    localfile= ['../results/',SIM_NAME,'/',BASIC.PARAMS,'/.'];
-    system(['scp -p ahoffman@login.marconi.cineca.it:',hostfile,' ',localfile]);
+    hostfile = ['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/',SIM_NAME,'/',BASIC.PARAMS];
+    localfile= ['../results/',SIM_NAME,'/'];
+    system(['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',localfile]);
     filename = ['../results/',SIM_NAME,'/',BASIC.PARAMS,'/outputs_00.h5'];
     TIMES(i_)   = h5readatt(filename,'/data/input','cpu_time');
     i_ = i_ + 1;
 end
+TIMES
 end
-%% Strong scaling measurement
+%% Strong scaling measurement (no diagnostics)
 
 % Handwritten results for 512x256, P,J=2,1, Tmax = 10, mu=0, dt = 5e-2
 Results_512_21.np    = NPS;
-Results_512_21.time  = [789   422   206   106    81    72    66    82]; %tmax 10
+Results_512_21.time  = [799   422   206   106    81    72    66    82]; %tmax 10
 % Results_512_21.time  = [0162, 0108, 0055, 0032,  0030, 0045, 0061, 0084]; %tmax 2
 
-% Handwritten results for 512x256, P,J=3,2, Tmax = 5, mu=0, dt 5e-2
+% Handwritten results for 512x256, P,J=3,2, Tmax = 10, mu=0, dt 5e-2
 Results_512_32.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
-Results_512_32.time  = [1221, 0608, 0307, 0163,  0130, 0127, 0194, 0260];
+Results_512_32.time  = [2421  1241  0598  0309   0221  0188  0159  0148];
 
 % Handwritten results for 1024x512, P,J=2,1, Tmax = 5 dt = 0.05, mu = 0
 Results_1024_21.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
 % Results_1024_21.time  = [1920, 0000, 0563, 0306,  0247, 0240, 0000, 0000];
-Results_1024_21.time  = [3808, 0000, 1108, 0586,  0465, 0443, 0483, 0496];
+Results_1024_21.time  = [1920, 1260, 0556, 0289,  0219, 0189, 0172, 0167];
 
 % Handwritten results for 1024x512, P,J=3,2, Tmax = 2 dt = 0.05, mu = 0
 Results_1024_32.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
-Results_1024_32.time  = [0000, 0000, 0000, 0000,  0000, 0000, 0000, 0000];
+Results_1024_32.time  = [2199, 1424, 0623, 0324,  0243, 0208, 0188, 0180];
 
 % 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];
diff --git a/wk/parameters_ZP.m b/wk/parameters_ZP.m
index 7991aa72404769f3f19cec2765246e24b6d214dc..41d2dd46a59ef0dd9d6f38bc31ed8f03f5cdf44b 100644
--- a/wk/parameters_ZP.m
+++ b/wk/parameters_ZP.m
@@ -12,19 +12,19 @@ ETAT    = 0.0;    % Temperature gradient
 MU      = 5e-4;   % Hyper diffusivity coefficient
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 256;     % Frequency gridpoints (Nkr = N/2)
-L       = 66;     % 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    = 400;  % Maximal time unit
+N       = 128;     % Frequency gridpoints (Nkr = N/2)
+L       = 33;     % Size of the squared frequency domain
+PMAXE   = 0;     % Highest electron Hermite polynomial degree
+JMAXE   = 0;     % Highest ''       Laguerre ''
+PMAXI   = 0;     % Highest ion      Hermite polynomial degree
+JMAXI   = 0;     % Highest ''       Laguerre ''
+%% TIME PARAMETERS
+TMAX    = 10;  % Maximal time unit
 DT      = 1e-2;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
 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
+RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 1;
 %% OPTIONS
 SIMID   = 'ZP';  % Name of the simulation
@@ -32,14 +32,11 @@ CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Do
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% unused
-KR0KH   = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
-NO_E    = 0;  % Remove electrons dynamic
 % DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
 KREQ0   = 0;      % put kr = 0
 KPAR    = 0.0;    % Parellel wave vector component
-LAMBDAD = 0.0; 
+LAMBDAD = 0.0;
 NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 LOAD_MARCONI = 0;
-CANCEL_ODD_P = 0;% Cancels the odd polynomials degree
 
 setup
diff --git a/wk/profiler.m b/wk/profiler.m
index 30aa7bfc8e92e29262583d8a19a448de62a1b524..4ca0045721d43681bbfcf726c2b66a7a68338484 100644
--- a/wk/profiler.m
+++ b/wk/profiler.m
@@ -1,7 +1,7 @@
 %% load profiling
 filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
 
-CPUTIME   = h5readatt(filename,'/data/input','cpu_time');
+CPUTIME   = double(h5readatt(filename,'/data/input','cpu_time'));
 DT_SIM    = h5readatt(filename,'/data/input','dt');
 
 
@@ -17,21 +17,41 @@ Ts0D         = h5read(filename,'/profiler/time');
 missing_Tc   = step_Tc - rhs_Tc - adv_field_Tc -...
                poisson_Tc - Sapj_Tc -diag_Tc -checkfield_Tc;
 total_Tc     = step_Tc;
+
+TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(poisson_Tc);...
+    diff(Sapj_Tc); diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)];
+TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/7,7]);
+
+TIME_PER_STEP = sum(TIME_PER_FCT,2);
+TIME_PER_CPU  = trapz(Ts0D(2:end),TIME_PER_STEP);
+
 %% Plots
-Y = [diff(rhs_Tc); diff(adv_field_Tc); diff(poisson_Tc);...
+TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(poisson_Tc);...
     diff(Sapj_Tc); diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)];
-Y = reshape(Y,[numel(Y)/7,7]);
+TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/7,7]);
 fig = figure;
 
-p1 = area(Ts0D(2:end),100*Y./diff(total_Tc),'LineStyle','none');
+p1 = area(Ts0D(2:end),TIME_PER_FCT,'LineStyle','none');
 legend('Compute RHS','Adv. fields','Poisson','Sapj','Check+Sym','Diag','Missing')
-xlabel('Sim. Time'); ylabel('Step Comp. Time [\%]')
-ylim([0,100]); xlim([Ts0D(2),Ts0D(end)]);
+xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('Step Comp. Time [s]')
+xlim([Ts0D(2),Ts0D(end)]);
+title(sprintf('Proc. #1, total sim. time  ~%.0f [h]',CPUTIME/3600))
 hold on
-yyaxis right
-p2 = plot(Ts0D(2:end),diff(total_Tc),'--r','LineWidth',1.0);
-ylabel('Step Comp. Time [s]')
-ylim([0,1.1*max(diff(total_Tc))])
-set(gca,'ycolor','r') 
 FIGNAME = 'profiler';
-save_figure
\ No newline at end of file
+save_figure
+
+%% Plots
+% fig = figure;
+% 
+% p1 = area(Ts0D(2:end),100*TIME_PER_FCT./diff(total_Tc),'LineStyle','none');
+% legend('Compute RHS','Adv. fields','Poisson','Sapj','Check+Sym','Diag','Missing')
+% xlabel('Sim. Time'); ylabel('Step Comp. Time [\%]')
+% ylim([0,100]); xlim([Ts0D(2),Ts0D(end)]);
+% hold on
+% yyaxis right
+% p2 = plot(Ts0D(2:end),diff(total_Tc),'--r','LineWidth',1.0);
+% ylabel('Step Comp. Time [s]')
+% ylim([0,1.1*max(diff(total_Tc))])
+% set(gca,'ycolor','r') 
+% FIGNAME = 'profiler';
+% save_figure
\ No newline at end of file
diff --git a/wk/setup.m b/wk/setup.m
index 1d84c183a2a4c43dc59d6a274158ae8722369cfa..0636b57da3df91f5bc76ecf1fc9a4546bff50730 100644
--- a/wk/setup.m
+++ b/wk/setup.m
@@ -2,7 +2,7 @@
 SIMDIR = ['../results/',SIMID,'/'];
 % Grid parameters
 GRID.pmaxe = PMAXE;  % Electron Hermite moments
-GRID.jmaxe = JMAXE;  % Electron Laguerre moments 
+GRID.jmaxe = JMAXE;  % Electron Laguerre moments
 GRID.pmaxi = PMAXI;  % Ion Hermite moments
 GRID.jmaxi = JMAXI;  % Ion Laguerre moments
 GRID.Nr    = N * (1-KREQ0) + KREQ0; % r grid resolution
@@ -10,11 +10,7 @@ GRID.Lr    = L * (1-KREQ0); % r length
 GRID.Nz    = N; % z ''
 GRID.Lz    = L; % z ''
 GRID.kpar  = KPAR;
-if CANCEL_ODD_P
-GRID.CANCEL_ODD_P = '.true.';
-else
-GRID.CANCEL_ODD_P = '.false.';
-end
+
 % Model parameters
 MODEL.CO      = CO;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 if 0;      MODEL.DK      = '.true.'; else; MODEL.DK      = '.false.';end;
@@ -36,9 +32,6 @@ MODEL.eta_n   = ETAN;        % source term kappa for HW
 MODEL.eta_T   = ETAT;        % Temperature
 MODEL.eta_B   = ETAB;        % Magnetic
 MODEL.lambdaD = LAMBDAD;
-% background phi drive for Kelvin-Helmholtz instability
-MODEL.kr0KH   = KR0KH;
-MODEL.A0KH    = A0KH;
 if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end;
 % Time integration and intialization parameters
 TIME_INTEGRATION.numerical_scheme  = '''RK4''';
@@ -76,8 +69,11 @@ BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
 BASIC.PARAMS = PARAMS;
 BASIC.SIMID  = SIMID;
 BASIC.nrun       = 1e8;
-BASIC.dt         = DT;   
+BASIC.dt         = DT;
 BASIC.tmax       = TMAX;    %time normalized to 1/omega_pe
+BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ...
+                   + str2num(CLUSTER.TIME(4:5))*60 ...
+                   + str2num(CLUSTER.TIME(7:8));
 % Outputs parameters
 if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end;
 OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT);
diff --git a/wk/test_parallel.m b/wk/test_parallel.m
index 70eadc474a6c2e72b7ddc8c2c115866cacdc9832..4a6afc8c86dcac22e63dba3a124f16a5b1d3e66a 100644
--- a/wk/test_parallel.m
+++ b/wk/test_parallel.m
@@ -2,6 +2,7 @@ clear all;
 addpath(genpath('../matlab')) % ... add
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
+CLUSTER.TIME  = '00:00:10'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
 NU      = 1e-1;   % Collision frequency
@@ -18,16 +19,17 @@ 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    = 40;  % Maximal time unit
+%% TIME PARAMETERS
+TMAX    = 20;  % 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
+SPSCP   = 1/10;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
-JOB2LOAD= 1;
+JOB2LOAD= 0;
 %% OPTIONS
-SIMID   = 'test_sapj_opt';  % Name of the simulation
+SIMID   = 'test_runtime';  % Name of the simulation
 CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -37,9 +39,8 @@ NO_E    = 0;  % Remove electrons dynamic
 % DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
 KREQ0   = 0;      % put kr = 0
 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
+LAMBDAD = 0.0;
+NON_LIN = 0 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 LOAD_MARCONI = 0;
 
 setup