diff --git a/wk/batch_script.sh b/wk/batch_script.sh
new file mode 100644
index 0000000000000000000000000000000000000000..50f9506474867f2cfd4d3a047f9d6aa925bdf028
--- /dev/null
+++ b/wk/batch_script.sh
@@ -0,0 +1,16 @@
+#!/bin/bash
+#SBATCH --time=24:00:00
+#SBATCH --nodes=1
+#SBATCH --cpus-per-task=1
+#SBATCH --ntasks-per-node=16
+#SBATCH --mem=32GB
+#SBATCH --error=../results/Marconi/512x256_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/err.txt
+#SBATCH --output=../results/Marconi/512x256_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/out.txt
+#SBATCH --account=FUA34_GBSedge
+#SBATCH --partition=skl_fua_prod
+
+module load intel
+module load intelmpi
+module load autoload hdf5/1.10.4--intelmpi--2018--binary
+module load fftw
+srun --cpu-bind=cores ./../bin/helaz
\ No newline at end of file
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..1a9988ae332849b59a7d52897d827c6b0f99f14d
--- /dev/null
+++ b/wk/marconi_run.m
@@ -0,0 +1,59 @@
+%clear all;
+addpath(genpath('../matlab')) % ... add
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% CLUSTER PARAMETERS
+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.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
+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 = 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 = 5e-3; % Time step
+SPS0D = 1/DT/4; % 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 = 0; % To restart from last checkpoint
+JOB2LOAD= 0;
+%% OPTIONS
+SIMID = 'Marconi'; % Name of the simulation
+CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% 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;
+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
+
+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
diff --git a/wk/marconi_scaling.m b/wk/marconi_scaling.m
new file mode 100644
index 0000000000000000000000000000000000000000..cd75f6415ea02f28692436208e822d63d5a845d7
--- /dev/null
+++ b/wk/marconi_scaling.m
@@ -0,0 +1,56 @@
+%clear all;
+addpath(genpath('../matlab')) % ... add
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% CLUSTER PARAMETERS
+CLUSTER.TIME = '12:00:00'; % allocation time hh:mm:ss
+CLUSTER.NODES = '1'; % MPI process
+CLUSTER.CPUPT = '1'; % CPU per task
+CLUSTER.NTPN = '20'; % N tasks per node
+CLUSTER.PART = 'prod'; % dbg or prod
+CLUSTER.MEM = '10GB'; % Memory
+%% PHYSICAL PARAMETERS
+NU = 1e-1; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+ETAB = 0.0; % Magnetic gradient
+ETAN = 1.0; % Density gradient
+ETAT = 0.0; % Temperature gradient
+MU = 0e-4; % Hyper diffusivity coefficient
+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 ''
+%% TIME PARAMETERS
+TMAX = 1; % Maximal time unit
+DT = 1e-2; % Time step
+SPS0D = 0; % 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
+JOB2LOAD= 0;
+%% OPTIONS
+SIMID = ['Scaling_np',num2str(CLUSTER.NTPN)]; % Name of the simulation
+CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% 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;
+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
diff --git a/wk/parallel_scaling.m b/wk/parallel_scaling.m
new file mode 100644
index 0000000000000000000000000000000000000000..206d117e17cea629e8acf9a85809a7c76fedb13b
--- /dev/null
+++ b/wk/parallel_scaling.m
@@ -0,0 +1,89 @@
+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
+Results_512_32.np = [ 1, 2, 4, 8, 16, 20];
+Results_512_32.time = [4450, 2267, 1136, 0595, 0363, 0323];
+
+% 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];
+
+% 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_22.np = [ 1, 2, 4, 6, 8, 10, 12, 16, 20];
+Results_1024_22.time = [2391, 1373, 0654, 0457, 0343, 0297, 0274, 0219, 0206];
+
+%
+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),'o-','DisplayName','$512\times256$, $P,J=2,1$');
+res = Results_512_32;
+plot(res.np,res.time(1)./(res.time),'o-','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),'o-','DisplayName','$1024\times512$, $P,J=2,2$');xlim([1,max(res.np)]);
+xlabel('$N_p$'); ylabel('speedup')
+legend('show')
+title('Strong scaling')
+grid on
+
+
+FIGNAME = [SIMDIR,'strong_scaling.png'];
+saveas(fig,FIGNAME);
+disp(['Figure saved @ : ',FIGNAME])
+
+%% 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 = [SIMDIR,'weak_scaling.png'];
+saveas(fig,FIGNAME);
+disp(['Figure saved @ : ',FIGNAME])
diff --git a/wk/run_helaz.sh b/wk/run_helaz.sh
new file mode 100644
index 0000000000000000000000000000000000000000..5703af9c44c22ff91c89fdabd4fbc6660b5d45f1
--- /dev/null
+++ b/wk/run_helaz.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+#SBATCH --time=01:00:00
+#SBATCH --nodes=1
+#SBATCH --cpus-per-task=1
+#SBATCH --ntasks-per-node=1
+#SBATCH --mem=10GB
+#SBATCH --error=../results/Scaling/1024x512_L_100_Pe_1_Je_1_Pi_1_Ji_1_nB_0_nN_1_nu_1e-01_DG_mu_0e+00/err.txt
+#SBATCH --output=../results/Scaling/1024x512_L_100_Pe_1_Je_1_Pi_1_Ji_1_nB_0_nN_1_nu_1e-01_DG_mu_0e+00/out.txt
+#SBATCH --account=FUA34_GBSedge
+#SBATCH --partition=skl_fua_dbg
+
+#SBATCH --job-name=1024x512_L_100_Pe_1_Je_1_Pi_1_Ji_1_nB_0_nN_1_nu_1e-01_DG_mu_0e+00
+
+module load intel
+module load intelmpi
+module load autoload hdf5/1.10.4--intelmpi--2018--binary
+module load fftw
+srun --cpu-bind=cores ./../bin/helaz
\ No newline at end of file
diff --git a/wk/setup_and_run.sh b/wk/setup_and_run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..52c28773a9a5cb2befa54a2f8604a72a31995606
--- /dev/null
+++ b/wk/setup_and_run.sh
@@ -0,0 +1,12 @@
+#!/bin/bash
+mkdir -p $CINECA_SCRATCH/HeLaZ/wk
+mkdir -p $CINECA_SCRATCH/HeLaZ/bin
+
+cd $CINECA_SCRATCH/HeLaZ/wk/
+cp $HOME/HeLaZ/wk/fort.90 .
+cp $HOME/HeLaZ/wk/batch_script.sh .
+cp -r $HOME/HeLaZ/iCa ..
+
+mkdir -p ../results/Marconi/512x256_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/
+sbatch batch_script.sh
+echo $CINECA_SCRATCH/HeLaZ/results/Marconi/512x256_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-01_FC_mu_5e-04/out.txt
\ No newline at end of file