diff --git a/matlab/load_marconi.m b/matlab/load_marconi.m
new file mode 100644
index 0000000000000000000000000000000000000000..a8e687dbdca443c5d18faba5aac0137be5152fb6
--- /dev/null
+++ b/matlab/load_marconi.m
@@ -0,0 +1,10 @@
+function [ RESDIR ] = load_marconi( outfilename )
+%UNTITLED2 Summary of this function goes here
+%   Detailed explanation goes here
+    hostfolder  = outfilename(1:end-8);
+    hostfile    = [hostfolder,'/outputs*'];
+    RESDIR      = ['../',outfilename(46:end-8),'/'];
+    localfolder = [RESDIR,'.'];
+    system(['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',localfolder])
+end
+
diff --git a/matlab/load_params.m b/matlab/load_params.m
new file mode 100644
index 0000000000000000000000000000000000000000..d1eee972ddda5a92ec3e73fa1f2c45a59c2ed955
--- /dev/null
+++ b/matlab/load_params.m
@@ -0,0 +1,18 @@
+CONAME = h5readatt(filename,'/data/input','CO');
+if CONAME == -1
+    CONAME = 'FC';
+elseif CONAME == 0
+    CONAME = 'LB';
+elseif CONAME == 1
+    CONAME = 'DG';
+end
+ETAB = h5readatt(filename,'/data/input','eta_B');
+ETAN = h5readatt(filename,'/data/input','eta_n');
+ETAT = h5readatt(filename,'/data/input','eta_T');
+
+NON_LIN = h5readatt(filename,'/data/input','NON_LIN');
+if NON_LIN == 'y'
+    NON_LIN = 1;
+else
+    NON_LIN = 0;
+end
diff --git a/matlab/write_sbash.m b/matlab/write_sbash.m
index 7590654231b6c6421f1d12b600ee86371fd7cff2..ba785d725051314556609535cb7cbec36faa1728 100644
--- a/matlab/write_sbash.m
+++ b/matlab/write_sbash.m
@@ -25,6 +25,7 @@ fid = fopen(INPUT,'wt');
 
 fprintf(fid,[...
 '#!/bin/bash\n',...
+'#SBATCH --job-name=',CLUSTER.JNAME,'\n',...
 '#SBATCH --time=', CLUSTER.TIME,'\n',...
 '#SBATCH --nodes=', CLUSTER.NODES,'\n',...
 '#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',...
diff --git a/src/auxval.F90 b/src/auxval.F90
index 0ef6890bd1b1f0ec7cf68fa0a247f6098f56f582..07f7acde91e4074425cbd213350ed23a7b5f5391 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -4,11 +4,11 @@ subroutine auxval
   USE basic
   USE grid
   USE array
-  USE fourier, ONLY: init_grid_distr_and_plans
+  USE fourier, ONLY: init_grid_distr_and_plans, alloc_local_1, alloc_local_2
   use prec_const
   IMPLICIT NONE
 
-  INTEGER :: irows,irowe, irow, icol
+  INTEGER :: irows,irowe, irow, icol, i_
   IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ==='
 
   CALL init_grid_distr_and_plans(Nr,Nz)
@@ -21,4 +21,12 @@ subroutine auxval
 
   CALL memory ! Allocate memory for global arrays
 
+  DO i_ = 0,num_procs-1
+    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+    IF (my_id .EQ. i_) WRITE (*,'(I2,A9,I3,A8,I3,A8,I3,A8,I3,A15,I6)') &
+    i_,': ikrs = ', ikrs, ' ikre = ', ikre, 'ikzs = ', ikzs, ' ikze = ', ikze, &
+    'alloc_local = ', alloc_local_1+alloc_local_2
+    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+  ENDDO
+
 END SUBROUTINE auxval
diff --git a/src/compute_Sapj_old.F90 b/src/compute_Sapj_old.F90
deleted file mode 100644
index b4d02513667ab91c92554c6c62f5bc8d89a95336..0000000000000000000000000000000000000000
--- a/src/compute_Sapj_old.F90
+++ /dev/null
@@ -1,207 +0,0 @@
-SUBROUTINE compute_Sapj
-  ! This routine is meant to compute the non linear term for each specie and degree
-  !! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
-  !! Sapj = Sum_n (ikr Kn phi)#(ikz Sum_s d_njs Naps) - (ikz Kn phi)#(ikr Sum_s d_njs Naps)
-  !! where # denotes the convolution.
-  USE array, ONLY : dnjs, Sepj, Sipj
-  USE basic
-  USE fourier
-  USE fields!, ONLY : phi, moments_e, moments_i
-  USE grid
-  USE model
-  USE prec_const
-  USE time_integration!, ONLY : updatetlevel
-  IMPLICIT NONE
-
-  COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze) :: F_, G_, CONV
-  INTEGER :: in, is
-  REAL(dp):: kr, kz, kerneln, be_2, bi_2, factn
-  REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
-
-  ! Execution time start
-  CALL cpu_time(t0_Sapj)
-
-  !!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
-  sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the kerneln argument
-
-  ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
-    jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
-      Sepj(ip,ij,:,:)  = 0._dp
-      factn = 1
-
-      nloope: DO in = 1,jmaxe+1 ! Loop over laguerre for the sum
-
-        krloope1: DO ikr = ikrs,ikre ! Loop over kr
-          kzloope1: DO ikz = ikzs,ikze ! Loop over kz
-            kr     = krarray(ikr)
-            kz     = kzarray(ikz)
-            be_2    = sigmae2_taue_o2 * (kr**2 + kz**2)
-            kerneln = be_2**(in-1)/factn * EXP(-be_2)
-            ! First convolution term
-            IF ( NON_LIN ) THEN
-              F_(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
-            ELSE
-              F_(ikr,ikz) = 0._dp
-            ENDIF
-            ! Background sinusoidal electrostatic potential phi_0 for KH inst.
-            IF ( q_e .NE. 0._dp ) THEN ! If electron are not removed
-              IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0
-                IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0
-                  F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz)
-                ENDIF
-              ENDIF
-            ENDIF
-            ! Second convolution term
-            G_(ikr,ikz) = 0._dp ! initialization of the sum
-            DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
-              G_(ikr,ikz) = G_(ikr,ikz) + &
-               dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
-            ENDDO
-            G_(ikr,ikz) = imagu*kz*G_(ikr,ikz)
-
-          ENDDO kzloope1
-        ENDDO krloope1
-
-        CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve Fourier to Fourier
-
-        DO ikr = ikrs,ikre ! Loop over kr
-          DO ikz = ikzs,ikze ! Loop over kz
-            Sepj(ip,ij,ikr,ikz) = Sepj(ip,ij,ikr,ikz) + CONV(ikr,ikz) ! Add it to Sepj
-          ENDDO
-        ENDDO
-
-        IF ( NON_LIN ) THEN ! Fully non linear term ikz phi * ikr Napj
-          krloope2: DO ikr = ikrs,ikre ! Loop over kr
-            kzloope2: DO ikz = ikzs,ikze ! Loop over kz
-              kr     = krarray(ikr)
-              kz     = kzarray(ikz)
-              be_2    = sigmae2_taue_o2 * (kr**2 + kz**2)
-              kerneln = be_2**(in-1)/factn * EXP(-be_2)
-              ! First convolution term
-              F_(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
-              ! Second convolution term
-              G_(ikr,ikz) = 0._dp ! initialization of the sum
-              DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
-                G_(ikr,ikz) = G_(ikr,ikz) + &
-                 dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
-              ENDDO
-              G_(ikr,ikz) = imagu*kr*G_(ikr,ikz)
-
-            ENDDO kzloope2
-          ENDDO krloope2
-
-          CALL convolve_2D_F2F( F_, G_, CONV )
-
-          DO ikr = ikrs,ikre
-            DO ikz = ikzs,ikze
-              Sepj(ip,ij,ikr,ikz) = Sepj(ip,ij,ikr,ikz) - CONV(ikr,ikz) ! substract it to Sepj
-            ENDDO
-          ENDDO
-
-        ENDIF
-
-        IF ( in + 1 .LE. jmaxe+1 ) THEN
-          factn = real(in,dp) * factn ! compute (n+1)!
-        ENDIF
-      ENDDO nloope
-
-    ENDDO jloope
-  ENDDO ploope
-  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
-  sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp
-
-  ploopi: DO ip = ips_i,ipe_i ! Loop over Hermite moments
-    jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
-      Sipj(ip,ij,:,:)  = 0._dp
-      factn = 1
-
-      nloopi: DO in = 1,jmaxi+1 ! Loop over laguerre for the sum
-
-        krloopi1: DO ikr = ikrs,ikre ! Loop over kr
-          kzloopi1: DO ikz = ikzs,ikze ! Loop over kz
-            kr      = krarray(ikr)
-            kz      = kzarray(ikz)
-            bi_2    = sigmai2_taui_o2 * (kr**2 + kz**2)
-            kerneln = bi_2**(in-1)/factn * EXP(-bi_2)
-
-            ! First convolution term
-            IF ( NON_LIN ) THEN
-              F_(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
-            ELSE
-              F_(ikr,ikz) = 0._dp
-            ENDIF
-            ! Background sinusoidal electrostatic potential phi_0 for KH inst.
-            IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0
-              IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=kr0
-                F_(ikr,ikz) = -A0KH/2._dp + F_(ikr,ikz)
-              ENDIF
-            ENDIF
-            ! Second convolution term
-            G_(ikr,ikz) = 0._dp ! initialization of the sum
-            DO is = 1, MIN( in+ij-1, jmaxi+1 )
-              G_(ikr,ikz) = G_(ikr,ikz) + &
-               dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
-            ENDDO
-            G_(ikr,ikz) = imagu*kz*G_(ikr,ikz)
-
-          ENDDO kzloopi1
-        ENDDO krloopi1
-
-        CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve Fourier to Fourier
-        DO ikr = ikrs,ikre
-          DO ikz = ikzs,ikze
-            Sipj(ip,ij,ikr,ikz) = Sipj(ip,ij,ikr,ikz) + CONV(ikr,ikz) ! add it to Sipj
-          ENDDO
-        ENDDO
-
-        krloopi2: DO ikr = ikrs,ikre ! Loop over kr
-          kzloopi2: DO ikz = ikzs,ikze ! Loop over kz
-            kr      = krarray(ikr)
-            kz      = kzarray(ikz)
-            bi_2    = sigmai2_taui_o2 * (kr**2 + kz**2)
-            kerneln = bi_2**(in-1)/factn * EXP(-bi_2)
-            ! First convolution term
-            F_(ikr,ikz) = imagu*kz*phi(ikr,ikz) * kerneln
-            ! Second convolution term
-            G_(ikr,ikz) = 0._dp ! initialization of the sum
-            IF ( NON_LIN ) THEN
-              DO is = 1, MIN( in+ij-1, jmaxi+1 )
-                G_(ikr,ikz) = G_(ikr,ikz) + &
-                 dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
-              ENDDO
-              G_(ikr,ikz) = imagu*kr*G_(ikr,ikz)
-            ENDIF
-            ! Background sinusoidal ion denstiy n_i0 for KH inst.
-            IF ( kz .EQ. 0._dp ) THEN ! Kronecker kz=0
-              IF ( ABS(kr) .EQ. ABS(kr0KH) ) THEN ! kronecker kr=+-kr0
-                G_(ikr,ikz) = -A0KH*kr0KH**2/2._dp + G_(ikr,ikz)
-              ENDIF
-            ENDIF
-
-          ENDDO kzloopi2
-        ENDDO krloopi2
-
-        CALL convolve_2D_F2F( F_, G_, CONV ) ! Convolve and back to Fourier
-        DO ikr = ikrs,ikre
-          DO ikz = ikzs,ikze
-            Sipj(ip,ij,ikr,ikz) = Sipj(ip,ij,ikr,ikz) - CONV(ikr,ikz) ! substract it to Sepj
-          ENDDO
-        ENDDO
-
-        IF ( in + 1 .LE. jmaxi+1 ) THEN
-          factn = real(in,dp) * factn ! factorial(n+1)
-        ENDIF
-
-      ENDDO nloopi
-
-    ENDDO jloopi
-  ENDDO ploopi
-  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  ! Execution time end
-  CALL cpu_time(t1_Sapj)
-  tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
-
-END SUBROUTINE compute_Sapj
diff --git a/src/compute_Sapj_opt.F90 b/src/compute_Sapj_opt.F90
deleted file mode 100644
index 3fed22876f9a3fd3e060899d62049b95eea37be5..0000000000000000000000000000000000000000
--- a/src/compute_Sapj_opt.F90
+++ /dev/null
@@ -1,189 +0,0 @@
-SUBROUTINE compute_Sapj
-  ! This routine is meant to compute the non linear term for each specie and degree
-  !! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
-  !! Sapj = Sum_n (ikr Kn phi)#(ikz Sum_s d_njs Naps) - (ikz Kn phi)#(ikr Sum_s d_njs Naps)
-  !! where # denotes the convolution.
-  USE array, ONLY : dnjs, Sepj, Sipj
-  USE basic
-  USE fourier
-  USE fields!, ONLY : phi, moments_e, moments_i
-  USE grid
-  USE model
-  USE prec_const
-  USE time_integration!, ONLY : updatetlevel
-  IMPLICIT NONE
-  INCLUDE 'fftw3-mpi.f03'
-
-  COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze) :: Fr_cmpx, Gz_cmpx
-  COMPLEX(dp), DIMENSION(ikrs:ikre,ikzs:ikze) :: Fz_cmpx, Gr_cmpx, F_conv_G
-  REAL(dp),    DIMENSION(irs:ire,izs:ize)     :: fr_real, gz_real
-  REAL(dp),    DIMENSION(irs:ire,izs:ize)     :: fz_real, gr_real, f_times_g
-
-  INTEGER :: in, is
-  REAL(dp):: kr, kz, kerneln, be_2, bi_2, factn
-  REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
-
-  ! Execution time start
-  CALL cpu_time(t0_Sapj)
-
-  !!!!!!!!!!!!!!!!!!!! ELECTRON non linear term computation (Sepj)!!!!!!!!!!
-  sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the kerneln argument
-
-  ploope: DO ip = ips_e,ipe_e ! Loop over Hermite moments
-    jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
-      real_data_c = 0._dp ! initialize sum over real nonlinear term
-      factn = 1
-
-      nloope: DO in = 1,jmaxe+1 ! Loop over laguerre for the sum
-
-        krloope: DO ikr = ikrs,ikre ! Loop over kr
-          kzloope: DO ikz = ikzs,ikze ! Loop over kz
-            kr     = krarray(ikr)
-            kz     = kzarray(ikz)
-            be_2    = sigmae2_taue_o2 * (kr**2 + kz**2)
-            kerneln = be_2**(in-1)/factn * EXP(-be_2)
-            ! First convolution terms
-            Fr_cmpx(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
-            Fz_cmpx(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
-            ! Second convolution terms
-            Gz_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
-            Gr_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
-            DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
-              Gz_cmpx(ikr,ikz) = Gz_cmpx(ikr,ikz) + &
-               dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
-              Gr_cmpx(ikr,ikz) = Gr_cmpx(ikr,ikz) + &
-               dnjs(in,ij,is) * moments_e(ip,is,ikr,ikz,updatetlevel)
-            ENDDO
-            Gz_cmpx(ikr,ikz) = imagu*kz*Gz_cmpx(ikr,ikz)
-            Gr_cmpx(ikr,ikz) = imagu*kr*Gr_cmpx(ikr,ikz)
-          ENDDO kzloope
-        ENDDO krloope
-
-        ! First term drphi x dzf
-        DO ikr = ikrs, ikre
-          DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fr_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gz_cmpx(ikr,ikz)
-          ENDDO
-        ENDDO
-
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
-
-        real_data_c = real_data_c + real_data_f/Nz/Nr  * real_data_g/Nz/Nr
-
-        ! Second term -dzphi x drf
-        DO ikr = ikrs, ikre
-          DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fz_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gr_cmpx(ikr,ikz)
-          ENDDO
-        ENDDO
-
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
-
-        real_data_c = real_data_c - real_data_f/Nz/Nr  * real_data_g/Nz/Nr
-
-        IF ( in + 1 .LE. jmaxe+1 ) THEN
-          factn = real(in,dp) * factn ! compute (n+1)!
-        ENDIF
-      ENDDO nloope
-
-      ! Put the real nonlinear product into k-space
-      call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
-
-      ! Retrieve convolution in input format
-      DO ikr = ikrs, ikre
-        DO ikz = ikzs, ikze
-          Sepj(ip,ij,ikr,ikz) = cmpx_data_c(ikz,ikr-local_nkr_offset)*AA_r(ikr)*AA_z(ikz)
-        ENDDO
-      ENDDO
-
-    ENDDO jloope
-  ENDDO ploope
-  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  !!!!!!!!!!!!!!!!!!!! ION non linear term computation (Sipj)!!!!!!!!!!
-  sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! factor of the kerneln argument
-
-  ploopi: DO ip = ips_e,ipe_e ! Loop over Hermite moments
-    jloopi: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
-      real_data_c = 0._dp ! initialize sum over real nonlinear term
-      factn = 1
-
-      nloopi: DO in = 1,jmaxi+1 ! Loop over laguerre for the sum
-
-        krloopi: DO ikr = ikrs,ikre ! Loop over kr
-          kzloopi: DO ikz = ikzs,ikze ! Loop over kz
-            kr      = krarray(ikr)
-            kz      = kzarray(ikz)
-            bi_2    = sigmai2_taui_o2 * (kr**2 + kz**2)
-            kerneln = bi_2**(in-1)/factn * EXP(-bi_2)
-            ! First convolution terms
-            Fr_cmpx(ikr,ikz) = imagu*kr* phi(ikr,ikz) * kerneln
-            Fz_cmpx(ikr,ikz) = imagu*kz* phi(ikr,ikz) * kerneln
-            ! Second convolution terms
-            Gz_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
-            Gr_cmpx(ikr,ikz) = 0._dp ! initialization of the sum
-            DO is = 1, MIN( in+ij-1, jmaxi+1 ) ! sum truncation on number of moments
-              Gz_cmpx(ikr,ikz) = Gz_cmpx(ikr,ikz) + &
-               dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
-              Gr_cmpx(ikr,ikz) = Gr_cmpx(ikr,ikz) + &
-               dnjs(in,ij,is) * moments_i(ip,is,ikr,ikz,updatetlevel)
-            ENDDO
-            Gz_cmpx(ikr,ikz) = imagu*kz*Gz_cmpx(ikr,ikz)
-            Gr_cmpx(ikr,ikz) = imagu*kr*Gr_cmpx(ikr,ikz)
-          ENDDO kzloopi
-        ENDDO krloopi
-
-        ! First term drphi x dzf
-        DO ikr = ikrs, ikre
-          DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fr_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gz_cmpx(ikr,ikz)
-          ENDDO
-        ENDDO
-
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
-
-        real_data_c = real_data_c + real_data_f/Nz/Nr  * real_data_g/Nz/Nr
-
-        ! Second term -dzphi x drf
-        DO ikr = ikrs, ikre
-          DO ikz = ikzs, ikze
-            cmpx_data_f(ikz,ikr-local_nkr_offset) = Fz_cmpx(ikr,ikz)
-            cmpx_data_g(ikz,ikr-local_nkr_offset) = Gr_cmpx(ikr,ikz)
-          ENDDO
-        ENDDO
-
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
-        call fftw_mpi_execute_dft_c2r(planb, cmpx_data_g, real_data_g)
-
-        real_data_c = real_data_c - real_data_f/Nz/Nr  * real_data_g/Nz/Nr
-
-        IF ( in + 1 .LE. jmaxe+1 ) THEN
-          factn = real(in,dp) * factn ! compute (n+1)!
-        ENDIF
-      ENDDO nloopi
-
-      ! Put the real nonlinear product into k-space
-      call fftw_mpi_execute_dft_r2c(planf, real_data_c, cmpx_data_c)
-
-      ! Retrieve convolution in input format
-      DO ikr = ikrs, ikre
-        DO ikz = ikzs, ikze
-          Sipj(ip,ij,ikr,ikz) = cmpx_data_c(ikz,ikr-local_nkr_offset)*AA_r(ikr)*AA_z(ikz)
-        ENDDO
-      ENDDO
-
-    ENDDO jloopi
-  ENDDO ploopi
-  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  ! Execution time END
-  CALL cpu_time(t1_Sapj)
-  tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
-
-END SUBROUTINE compute_Sapj
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index f1d4eb1bd7b806416d57e317869f76c4cf2bdbe2..528e21ddede3301b8f2b3a8b4baa928f19c2b2c9 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -18,7 +18,8 @@ MODULE fourier
   type(C_PTR)                            :: cdatar_f, cdatar_g, cdatar_c
   type(C_PTR)                            :: cdatac_f, cdatac_g, cdatac_c
   type(C_PTR) , PUBLIC                   :: planf, planb
-  integer(C_INTPTR_T)                    :: i, ix, iy, alloc_local_1, alloc_local_2
+  integer(C_INTPTR_T)                    :: i, ix, iy
+  integer(C_INTPTR_T), PUBLIC            :: alloc_local_1, alloc_local_2
   integer(C_INTPTR_T)                    :: NR_, NZ_
 
   CONTAINS
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 77b6119bc79e41cffe5104f5e311abc2cc4045c2..8f5be3dfa155bf88a363ed3100a3c52a28e18984 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -94,12 +94,6 @@ CONTAINS
       ikrs = local_nkr_offset + 1
       ikre = ikrs + local_nkr - 1
 
-      DO i_ = 0,num_procs-1
-        CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-        IF (my_id .EQ. i_) print *, i_,': ikrs = ', ikrs, ' ikre = ', ikre
-        CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-      ENDDO
-
       IF (my_id .EQ. num_procs-1) ikre = Nkr
       !WRITE(*,*) 'ID = ',my_id,' ikrs = ', ikrs, ' ikre = ', ikre
       ! Grid spacings
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index dc5198904d5581f68e0b5ac458cb358db1a61010..4e01df312d990de8bb22a7dddf1d8e7eec44e60c 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,12 +1,14 @@
 %% Load results
-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])
+if 0
+    %%
+    outfile = '/marconi_scratch/userexternal/ahoffman/HeLaZ/results/Marconi/512x256_L_100_Pe_6_Je_3_Pi_6_Ji_3_nB_0.8_nN_1_nu_1e-01_FC_mu_5e-04/out.txt';
+    BASIC.RESDIR = load_marconi(outfile);
 end
+%%
 % JOBNUM = 0; load_results;
 % JOBNUM = 1; load_results;
 compile_results
+load_params
 
 %% Retrieving max polynomial degree and sampling info
 Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
@@ -16,10 +18,8 @@ Ns2D      = numel(Ts2D);
 % renaming and reshaping quantity of interest
 Ts5D      = Ts5D';
 Ts2D      = Ts2D';
-if strcmp(OUTPUTS.write_non_lin,'.true.')
-    Si00      = squeeze(Sipj(1,1,:,:,:));
-    Se00      = squeeze(Sepj(1,1,:,:,:));
-end
+Si00      = squeeze(Sipj(1,1,:,:,:));
+Se00      = squeeze(Sepj(1,1,:,:,:));
 %% Build grids
 Nkr = numel(kr); Nkz = numel(kz);
 [KZ,KR] = meshgrid(kz,kr);
@@ -53,11 +53,10 @@ for it = 1:numel(Ts2D)
     dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
 end
 
-if strcmp(OUTPUTS.write_non_lin,'.true.')
 for it = 1:numel(Ts5D)
     si00(:,:,it)      = real(fftshift(ifft2(squeeze(Si00(:,:,it)),Nr,Nz)));
 end
-end
+
 % Post processing
 disp('- post processing')
 E_pot    = zeros(1,Ns2D);      % Potential energy n^2
@@ -69,12 +68,11 @@ Flux_re  = zeros(1,Ns2D);      % Particle flux Gamma = <ne drphi>
 Flux_ze  = zeros(1,Ns2D);      % Particle flux Gamma = <ne dzphi>
 Ne_norm  = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Napj
 Ni_norm  = zeros(Npi,Nji,Ns5D);% .
-if strcmp(OUTPUTS.write_non_lin,'.true.')
 Se_norm  = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Sapj
 Si_norm  = zeros(Npi,Nji,Ns5D);% .
 Sne00_norm = zeros(1,Ns2D);    % Time evol. of the amp of e nonlin term
 Sni00_norm = zeros(1,Ns2D);    %
-end
+
 
 Ddr = 1i*KR; Ddz = 1i*KZ; lapl   = Ddr.^2 + Ddz.^2; 
 
@@ -96,14 +94,11 @@ dEdt     = diff(E_pot+E_kin)./dt2D;
 for it = 1:numel(Ts5D) % Loop over 5D arrays
     Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
     Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz;
-if strcmp(OUTPUTS.write_non_lin,'.true.')   
     Se_norm(:,:,it)= sum(sum(abs(Sepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
     Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4)/Nkr/Nkz;
     Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))))/Nkr/Nkz;
     Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))))/Nkr/Nkz;
 end
-end
-
 
 %% Compute growth rate
 if NON_LIN == 0
@@ -241,7 +236,7 @@ save_figure
 end
 
 %%
-t0    = 0;
+t0    = 1;
 skip_ = 1; 
 DELAY = 0.01*skip_;
 FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
@@ -310,56 +305,29 @@ fig = figure; FIGNAME = 'space_time_drphi';set(gcf, 'Position',  [100, 100, 1200
 save_figure
 end
 
-if 0
-%% Flow
-skip = 10;
-plt  = @(x) x(1:skip:end,1:skip:end);
-tf = 120; [~,it] = min(abs(Ts2D-tf));
-
-fig = figure; FIGNAME = 'space_time_drphi';
-    quiver(plt(RR),plt(ZZ),plt(dzphi(:,:,it)),plt(-drphi(:,:,it)),0);
-    xlabel('$r$'), ylabel('$z$')
-    title(sprintf('$v_E$, $t=%.0f c_s/R$',Ts2D(it)));
-end
-
-if 0
-%% Growth rate
-fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)];
-    [~,ikr0KH] = min(abs(kr-KR0KH));
-    plot(kz/kr(ikr0KH),gkr0kz_Ni00/(kr(ikr0KH)*A0KH),'DisplayName',['J = ',num2str(JMAXI)])
-    xlabel('$k_z/k_{r0}$'); ylabel('$\gamma_{Ni}/(A_0k_{r0})$');
-    xlim([0. 1.0]); ylim([0. 0.6]);
-save_figure
-
-end
 %%
 if 0
 %% Mode time evolution
+[~,ikr ] = min(abs(kr-dkr));
 [~,ik00] = min(abs(kz));
 [~,idk]  = min(abs(kz-dkz));
-[~,ik50] = min(abs(kz-0.5*KR0KH));
-[~,ik75] = min(abs(kz-0.7*KR0KH));
-[~,ik10] = min(abs(kz-1.00*KR0KH));
+[~,ik50] = min(abs(kz-0.1*max(kz)));
+[~,ik75] = min(abs(kz-0.2*max(kz)));
+[~,ik10] = min(abs(kz-0.3*max(kz)));
 plt = @(x) abs(squeeze(x));
 fig = figure; FIGNAME = ['mode_time_evolution',sprintf('_%.2d',JOBNUM)];
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,:)),'DisplayName', ...
-            ['$k_z/k_{r0} = $',num2str(kz(ik00)/KR0KH)]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr0KH,idk,:)),'DisplayName', ...
-            ['$k_z/k_{r0} = $',num2str(kz(idk)/KR0KH)]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,:)),'DisplayName', ...
-            ['$k_z/k_{r0} = $',num2str(kz(ik50)/KR0KH)]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,:)),'DisplayName', ...
-            ['$k_z/k_{r0} = $',num2str(kz(ik75)/KR0KH)]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,:)),'DisplayName', ...
-            ['$k_z/k_{r0} = $',num2str(kz(ik10)/KR0KH)]); hold on
+        semilogy(Ts2D,plt(Ni00(ikr,ik00,:)),'DisplayName', ...
+            ['$k_z = $',num2str(kz(ik00))]); hold on
+        semilogy(Ts2D,plt(Ni00(ikr,idk,:)),'DisplayName', ...
+            ['$k_z = $',num2str(kz(idk))]); hold on
+        semilogy(Ts2D,plt(Ni00(ikr,ik50,:)),'DisplayName', ...
+            ['$k_z = $',num2str(kz(ik50))]); hold on
+        semilogy(Ts2D,plt(Ni00(ikr,ik75,:)),'DisplayName', ...
+            ['$k_z = $',num2str(kz(ik75))]); hold on
+        semilogy(Ts2D,plt(Ni00(ikr,ik10,:)),'DisplayName', ...
+            ['$k_z = $',num2str(kz(ik10))]); hold on
         xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show');
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik00,end))*exp(gkr0kz_Ni00(ik00)*(Ts2D-Ts2D(end))),'k--')
-        semilogy(Ts2D,plt(Ni00(ikr0KH,idk,end))*exp(gkr0kz_Ni00(idk)*(Ts2D-Ts2D(end))),'k--')
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik50,end))*exp(gkr0kz_Ni00(ik50)*(Ts2D-Ts2D(end))),'k--')
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik75,end))*exp(gkr0kz_Ni00(ik75)*(Ts2D-Ts2D(end))),'k--')
-        semilogy(Ts2D,plt(Ni00(ikr0KH,ik10,end))*exp(gkr0kz_Ni00(ik10)*(Ts2D-Ts2D(end))),'k--')
-        plot(tstart*[1 1],ylim,'k-','LineWidth',0.5);
-        plot(tend*[1 1],ylim,'k-','LineWidth',0.5);
+title(sprintf('$k_r=$ %1.1f',kr(ikr)))
 save_figure
 end
 
@@ -384,3 +352,37 @@ if strcmp(OUTPUTS.write_non_lin,'.true.')
 end
 save_figure
 end
+
+%% Phase space distribution function
+% M_ = 25;
+% spar = linspace(0,4,M_); xperp = spar;
+% 
+% PSDF_e = zeros(numel(kr),numel(kz),M_,M_,numel(Ts5D));
+% for ikr = 1:numel(kr)
+%     for ikz = 1:numel(kz)
+%         for it = 1:numel(Ts5D)
+%         PSDF_e(ikr,ikz,:,:,it) = compute_fa(Nepj(:,:,ikr,ikz,it), spar, xperp);
+%         end
+%     end
+% end
+% 
+% ktarget = 0.3*max(kz);
+% ttarget = 310;
+% 
+% [~,ik10] = min(abs(kz-ktarget));
+% [~,it]   = min(abs(Ts5D-ttarget));
+% if 0
+%     %%
+% plt = @(x) real(x(ik10,ik10,:,:,it));
+% pclr = pcolor(spar,xperp,plt(PSDF_e)); set(pclr, 'edgecolor','none'); colorbar;
+% xlabel('$s_\parallel$'); ylabel('$x_\perp$'); title(sprintf('$t c_s/R=%.0f$',Ts5D(it))); 
+% legend(['$Re(f_e),k_\perp \approx$',sprintf('%01.0f',norm([kr(ik10),kz(ik10)]))]);
+% end
+% if 0
+% %% Phase space distribution function time evolution
+% GIFNAME = ['f_e',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
+% FIELD = real(ni00+ne00); X = RR; Y = ZZ; T = Ts5D;
+% FIELDNAME = '$n_i-n_e$'; XNAME = '$r\rho_s$'; YNAME = '$z\rho_s$';
+% create_gif
+% end
+
diff --git a/wk/load_results.m b/wk/load_results.m
index 43bdbdeb1d508f8aef6373ca1528731c7b1e45f2..db7c69dc3571dc64f0bd6e82836f209510d61573 100644
--- a/wk/load_results.m
+++ b/wk/load_results.m
@@ -4,10 +4,8 @@ disp(['Loading ',filename])
 % Loading from output file
 CPUTIME   = h5readatt(filename,'/data/input','cpu_time');
 DT_SIM    = h5readatt(filename,'/data/input','dt');
-if strcmp(OUTPUTS.write_moments,'.true.') 
 [Nipj, Pi, Ji, kr, kz, Ts5D, dt5D] = load_5D_data(filename, 'moments_i');
 [Nepj, Pe, Je                    ] = load_5D_data(filename, 'moments_e');
-end
 [Ni00, kr, kz, Ts2D, dt2D] = load_2D_data(filename, 'Ni00');
  Ne00                      = load_2D_data(filename, 'Ne00');
 %Pi = [0]; Ji = Pi; Pe = Pi; Je = Pi;
@@ -16,7 +14,5 @@ end
 % Nipj(1,1,:,:,:) = Ni00; Nepj(1,1,:,:,:) = Ne00;
 PHI                            = load_2D_data(filename, 'phi');
 
-if strcmp(OUTPUTS.write_non_lin,'.true.') 
-    Sipj    = load_5D_data(filename, 'Sipj');
-    Sepj    = load_5D_data(filename, 'Sepj');
-end
\ No newline at end of file
+Sipj    = load_5D_data(filename, 'Sipj');
+Sepj    = load_5D_data(filename, 'Sepj');
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 4a2ae9df68cd32a7b06d353f2738715857497c4b..3fc57d14fbc0c9bf0a4804f72a20653b2bcf78bb 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -8,26 +8,27 @@ CLUSTER.TIME  = '24: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   = '32GB';     % Memory
+CLUSTER.PART  = 'prod';     % dbg or prod
+CLUSTER.MEM   = '16GB';     % Memory
+CLUSTER.JNAME = 'gamma_inf';     % Job name
 %% PHYSICAL PARAMETERS
 NU      = 1e-1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.4;    % Magnetic gradient
+ETAB    = 0.66;    % 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       = 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 ''
+L       = 100;     % Size of the squared frequency domain
+PMAXE   = 8;    % Highest electron Hermite polynomial degree
+JMAXE   = 4;     % Highest ''       Laguerre ''
+PMAXI   = 8;     % Highest ion      Hermite polynomial degree
+JMAXI   = 4;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS
-TMAX    = 200;  % Maximal time unit
-DT      = 5e-3;   % Time step
+TMAX    = 400;  % Maximal time unit
+DT      = 1e-2;   % Time step
 SPS0D   = 10;    % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/10;    % Sampling per time unit for 5D arrays
diff --git a/wk/marconi_scaling.m b/wk/marconi_scaling.m
index a65e63789f7fd01c817dfa2808d1ea78488e6246..fb087636e48d9402eba222268caac3aead745bc4 100644
--- a/wk/marconi_scaling.m
+++ b/wk/marconi_scaling.m
@@ -5,11 +5,12 @@ addpath(genpath('../matlab')) % ... add
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% CLUSTER PARAMETERS
 CLUSTER.TIME  = '01:00:00'; % allocation time hh:mm:ss
-CLUSTER.NODES = '01';       % MPI process
+CLUSTER.NODES = '02';       % MPI process
 CLUSTER.CPUPT = '01';       % CPU per task
-CLUSTER.NTPN  = '24';       % N tasks per node (openMP)
-CLUSTER.PART  = 'prod';      % dbg or prod
-CLUSTER.MEM   = '16GB';     % Memory
+CLUSTER.NTPN  = '14';       % N tasks per node (openMP)
+CLUSTER.PART  = 'dbg';      % dbg or prod
+CLUSTER.MEM   = '8GB';     % Memory
+CLUSTER.JNAME = 'scaling';     % Job name
 %% PHYSICAL PARAMETERS
 NU      = 1e-1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
@@ -19,14 +20,14 @@ ETAT    = 0.0;    % Temperature gradient
 MU      = 0e-4;   % Hyper diffusivity coefficient
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 128;     % Frequency gridpoints (Nkr = N/2)
+N       = 1024;     % Frequency gridpoints (Nkr = N/2)
 L       = 10;     % 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   = 2;     % Highest electron Hermite polynomial degree
+JMAXE   = 1;     % Highest ''       Laguerre ''
+PMAXI   = 2;     % Highest ion      Hermite polynomial degree
+JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS
-TMAX    = 2;  % Maximal time unit
+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
@@ -35,7 +36,7 @@ SPSCP   = -1;      % 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   = ['nn',num2str(CLUSTER.NODES),'_np',num2str(CLUSTER.NTPN)];  % Name of the simulation
 CO      = -2;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -53,4 +54,8 @@ setup
 
 write_sbash
 
+system(['mkdir -p ../results/Scaling']);
+system(['cp -r ../results/',SIMID,' ../results/Scaling/']);
+system(['rm -r ../results/',SIMID]);
+
 MARCONI = 1;
diff --git a/wk/parallel_scaling_new.m b/wk/parallel_scaling_new.m
index c9608265e7ed5bdc343f27a279464439d2bce4cb..1ad96fa61aedff5d11c5219dd668749cad167246 100644
--- a/wk/parallel_scaling_new.m
+++ b/wk/parallel_scaling_new.m
@@ -1,68 +1,103 @@
 default_plots_options
 
-NPS = [01 02 04 08 12 16 20 24];
+% NPS = [01 02 04 08 12 16 20 24];
+NPS = [02 04 10];
 TIMES = NPS;
 if 0
 %% Load times
 i_ = 1;
+nn = 2;
 for np = NPS
-    SIM_NAME = sprintf('Scaling_np%02d',np);
+    SIM_NAME = sprintf('nn%02d_np%02d',nn,np)
     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'];
+    localtarget= ['../results/',sprintf('Scaling/nn%02d_np%02d',nn,np),'/'];
+    system(['scp -r ahoffman@login.marconi.cineca.it:',hostfile,' ',localtarget]);
+    filename = ['../results/',sprintf('Scaling/nn%02d_np%02d',nn,np),'/',BASIC.PARAMS,'/outputs_00.h5'];
     TIMES(i_)   = h5readatt(filename,'/data/input','cpu_time');
     i_ = i_ + 1;
 end
-TIMES
+disp(PARAMS)
+disp(TIMES')
 end
 %% 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  = [799   422   206   106    81    72    66    82]; %tmax 10
-% Results_512_21.time  = [0162, 0108, 0055, 0032,  0030, 0045, 0061, 0084]; %tmax 2
+Results_512_21.tproc = [857   447   217   107    74    59    46    41]; %Nthreads
+Results_512_21.tnode = [854   442   199    96    75    64    57    55]; %Nnodes
 
 % 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  = [2421  1241  0598  0309   0221  0188  0159  0148];
+Results_512_32.tproc = [2546  1290  0630  0316   0220  0177  0143  0125];
+
+% Handwritten results for 512x256, P,J=6,4, Tmax = 5, mu=0, dt 5e-2
+Results_512_64.np    = [   1,    2,    4,    8,    12,   16,   20,   24];
+Results_512_64.tproc = [5801  2974  1467  0746   0507  0408  0320  0279];
 
 % 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  = [1920, 1260, 0556, 0289,  0219, 0189, 0172, 0167];
+Results_1024_21.tproc = [2051, 1296, 0583, 0296,  0205, 0163, 0132, 0115];
+Results_1024_21.tnode = [2052  1327  0511  0231   0157  0109  0099  0083]; %Nnodes
 
 % 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  = [2199, 1424, 0623, 0324,  0243, 0208, 0188, 0180];
+Results_1024_32.tproc = [2248, 1526, 0644, 0326,  0233, 0182, 0148, 0130];
 
 % 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];
+Results_1024_64.tproc = [10330,6548, 2901, 1497,  1001, 0769, 0610, 0520];
 
 %
 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_21;
-plot(res.np,res.time(1)./(res.time),'o-','DisplayName','$1024\times512$, $P,J=2,1$');
-res = Results_1024_32;
-plot(res.np,res.time(1)./(res.time),'s-','DisplayName','$1024\times512$, $P,J=3,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  
+    subplot(121)
+    plot(1:24,1:24,'-k','DisplayName','Ideal')
+    hold on
+    res = Results_512_21;
+    plot(res.np,res.tproc(1)./(res.tproc),'v-','Color',line_colors(1,:),...
+        'DisplayName','$512\times256$, $P,J=2,1$');
+
+    res = Results_512_32;
+    plot(res.np,res.tproc(1)./(res.tproc),'v-','Color',line_colors(2,:),...
+        'DisplayName','$512\times256$, $P,J=3,2$');
+
+    res = Results_512_64;
+    plot(res.np,res.tproc(1)./(res.tproc),'v-','Color',line_colors(3,:),...
+        'DisplayName','$512\times256$, $P,J=6,4$');
+
+    res = Results_1024_21;
+    plot(res.np,res.tproc(1)./(res.tproc),'^--','Color',line_colors(1,:),...
+        'DisplayName','$1024\times512$, $P,J=2,1$');
+
+    res = Results_1024_32;
+    plot(res.np,res.tproc(1)./(res.tproc),'^--','Color',line_colors(2,:),...
+        'DisplayName','$1024\times512$, $P,J=3,2$');
+
+    res = Results_1024_64;
+    plot(res.np,res.tproc(1)./(res.tproc),'^--','Color',line_colors(3,:),...
+        'DisplayName','$1024\times512$, $P,J=6,4$');
+
+    xlabel('$N_p$'); ylabel('speedup')
+    xlim([1,24]);
+    legend('show')
+    title('$N_n=01$')
     
+subplot(122)
+    plot(1:24,1:24,'-k','DisplayName','Ideal')
+    hold on
+    res = Results_512_21;
+    plot(res.np,res.tnode(1)./(res.tnode),'o-','Color',line_colors(1,:),...
+        'DisplayName','$512\times256$, $P,J=2,1$');
+
+    res = Results_1024_21;
+    plot(res.np,res.tnode(1)./(res.tnode),'o--','Color',line_colors(1,:),...
+        'DisplayName','$1024\times512$, $P,J=2,1$');
+
+    xlabel('$N_n$'); ylabel('speedup')
+    xlim([1,24]);
+    legend('show')
+    title('$N_p=01$')
+        
 
 FIGNAME = '/home/ahoffman/HeLaZ/results/strong_scaling_new.pdf';
 saveas(fig,FIGNAME);
diff --git a/wk/profiler.m b/wk/profiler.m
index 4ca0045721d43681bbfcf726c2b66a7a68338484..e578665159ac924839c9c4da56c09a0eb703807c 100644
--- a/wk/profiler.m
+++ b/wk/profiler.m
@@ -1,5 +1,5 @@
 %% load profiling
-filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
+% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
 
 CPUTIME   = double(h5readatt(filename,'/data/input','cpu_time'));
 DT_SIM    = h5readatt(filename,'/data/input','dt');
diff --git a/wk/test_parallel.m b/wk/test_parallel.m
index 4a6afc8c86dcac22e63dba3a124f16a5b1d3e66a..8ae627881de323155b03b22d8b7ee45d1d20b593 100644
--- a/wk/test_parallel.m
+++ b/wk/test_parallel.m
@@ -2,7 +2,7 @@ clear all;
 addpath(genpath('../matlab')) % ... add
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
-CLUSTER.TIME  = '00:00:10'; % allocation time hh:mm:ss
+CLUSTER.TIME  = '00:10:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
 NU      = 1e-1;   % Collision frequency
@@ -29,18 +29,17 @@ SPSCP   = 1/10;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
-SIMID   = 'test_runtime';  % Name of the simulation
-CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+SIMID   = 'debug';  % Name of the simulation
+CO      = 0;  % 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 = 0 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
+NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 LOAD_MARCONI = 0;
 
 setup