diff --git a/Makefile b/Makefile
index d776c7daa5678172ff7fdf2d5f0461d7d9dc4133..21a929d8622fd9d4d8ba50dd9b76850caa8acce1 100644
--- a/Makefile
+++ b/Makefile
@@ -52,9 +52,9 @@ $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
 $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \
 $(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o \
 $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o \
-$(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/readinputs.o \
-$(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o $(OBJDIR)/stepon.o \
-$(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
+$(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o \
+$(OBJDIR)/readinputs.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o \
+$(OBJDIR)/stepon.o $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
 
  $(EXEC): $(FOBJ)
 	$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
@@ -86,10 +86,7 @@ $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
  $(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/control.F90 -o $@
 
- $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
-	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@
-
- $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/diagnose.o : src/diagnose.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnose.F90 -o $@
 
  $(OBJDIR)/diagnostics_par_mod.o : src/diagnostics_par_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
@@ -101,6 +98,9 @@ $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
  $(OBJDIR)/fields_mod.o : src/fields_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@
 
+ $(OBJDIR)/fourier_mod.o : src/fourier_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fourier_mod.F90 -o $@
+
  $(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/ppinit.o  $(OBJDIR)/time_integration_mod.o 
 		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@
 
@@ -137,6 +137,9 @@ $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
  $(OBJDIR)/prec_const_mod.o : src/prec_const_mod.F90
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/prec_const_mod.F90 -o $@
 
+ $(OBJDIR)/processing_mod.o : src/processing_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/basic_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/processing_mod.F90 -o $@
+
  $(OBJDIR)/readinputs.o : src/readinputs.F90  $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/readinputs.F90 -o $@
 
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 74bd6379d21b14918c9af2c19b277a48cd5394ef..44065bbc397a96a3c49595241f0ecebd930ce784 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -14,35 +14,51 @@ while(CONTINUE)
         % Load results of simulation #JOBNUM
         load_results
         % Check polynomials degrees
-        sz = size(Nepj); Pe_new= sz(1); Je_new= sz(2);
-        sz = size(Nipj); Pi_new= sz(1); Ji_new= sz(2);
+        Pe_new= numel(Pe); Je_new= numel(Je);
+        Pi_new= numel(Pi); Ji_new= numel(Ji);
         % If a degree is larger than previous job, put them in a larger array
         if (sum([Pe_new, Je_new, Pi_new, Ji_new]>[Pe_old, Je_old, Pi_old, Ji_old]) >= 1)
-            tmp = Nipj_; sz = size(tmp);
-            Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
-            Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
-            tmp = Nepj_; sz = size(tmp);
-            Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
-            Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
-            tmp = Sipj_; sz = size(tmp);
-            Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
-            Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
-            tmp = Sepj_; sz = size(tmp);
-            Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
-            Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
+            if W_NAPJ
+                tmp = Nipj_; sz = size(tmp);
+                Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
+                Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
+                tmp = Nepj_; sz = size(tmp);
+                Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
+                Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
+            end
+            if W_SAPJ
+                tmp = Sipj_; sz = size(tmp);
+                Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
+                Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
+                tmp = Sepj_; sz = size(tmp);
+                Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
+                Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
+            end
         end
         
-        
-        Nipj_ = cat(5,Nipj_,Nipj);
-        Nepj_ = cat(5,Nepj_,Nepj);
-        Ni00_ = cat(3,Ni00_,Ni00);
-        Ne00_ = cat(3,Ne00_,Ne00);
-        PHI_  = cat(3,PHI_,PHI);
-        Ts2D_   = cat(1,Ts2D_,Ts2D);
-        Ts5D_   = cat(1,Ts5D_,Ts5D);
-        
-        Sipj_ = cat(5,Sipj_,Sipj);
-        Sepj_ = cat(5,Sepj_,Sepj);
+
+        if W_PHI || W_NA00
+        	Ts2D_   = cat(1,Ts2D_,Ts2D);
+        end
+        if W_PHI
+            PHI_  = cat(3,PHI_,PHI);
+        end
+        if W_NA00
+            Ni00_ = cat(3,Ni00_,Ni00);
+            Ne00_ = cat(3,Ne00_,Ne00);
+        end
+
+        if W_NAPJ || W_SAPJ
+            Ts5D_   = cat(1,Ts5D_,Ts5D);
+        end
+        if W_NAPJ
+            Nipj_ = cat(5,Nipj_,Nipj);
+            Nepj_ = cat(5,Nepj_,Nepj);
+        end
+        if W_SAPJ
+     	  Sipj_ = cat(5,Sipj_,Sipj);
+          Sepj_ = cat(5,Sepj_,Sepj);
+        end
 
         JOBFOUND = JOBFOUND + 1;
         LASTJOB  = JOBNUM;
diff --git a/matlab/load_0D_data.m b/matlab/load_0D_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..7c2182f1281be5310dd2245144305eb005bbb50b
--- /dev/null
+++ b/matlab/load_0D_data.m
@@ -0,0 +1,9 @@
+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');
+    dt    = h5readatt(filename,'/data/input','dt');
+    cstart= h5readatt(filename,'/data/input','start_iframe0d'); 
+    data     = h5read(filename,['/data/var0d/',variablename]);
+
+end
+
diff --git a/matlab/load_2D_data.m b/matlab/load_2D_data.m
index f03bbaa3a0e23e209ebe32cd137b3b2ef9aa76a2..47ef91a645345a1f29fb59805dd52e90a33e9c8c 100644
--- a/matlab/load_2D_data.m
+++ b/matlab/load_2D_data.m
@@ -1,11 +1,8 @@
-function [ data, kr, kz, time, dt ] = load_2D_data( filename, variablename )
+function [ data, time, dt ] = load_2D_data( filename, variablename )
 %LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ
     time     = h5read(filename,'/data/var2d/time');
     kr       = h5read(filename,'/data/grid/coordkr');
     kz       = h5read(filename,'/data/grid/coordkz');
-%     kr       = h5read(filename,['/data/var2d/',variablename,'/coordkr']);
-%     kz       = h5read(filename,['/data/var2d/',variablename,'/coordkz']);
-
     dt    = h5readatt(filename,'/data/input','dt');
     cstart= h5readatt(filename,'/data/input','start_iframe2d'); 
     
diff --git a/matlab/load_5D_data.m b/matlab/load_5D_data.m
index 87edc9b9446f37900b4f016610902967822fd301..e36d767550e9bf69000364153c4563ff4410363a 100644
--- a/matlab/load_5D_data.m
+++ b/matlab/load_5D_data.m
@@ -1,14 +1,15 @@
-function [ data, p, j, kr, kz, time, dt ] = load_5D_data( filename, variablename )
+function [ data, time, dt ] = load_5D_data( filename, variablename )
 %LOAD_5D_DATA load a 5D variable stored in a hdf5 result file from HeLaZ
     time  = h5read(filename,'/data/var5d/time');
-    p     = h5read(filename,'/data/grid/coordp');
-    j     = h5read(filename,'/data/grid/coordj');
+    if strcmp(variablename,'moments_e') || strcmp(variablename,'Sepj')
+        p     = h5read(filename,'/data/grid/coordp_e');
+        j     = h5read(filename,'/data/grid/coordj_e');
+    else
+        p     = h5read(filename,'/data/grid/coordp_i');
+        j     = h5read(filename,'/data/grid/coordj_i');
+    end
     kr    = h5read(filename,'/data/grid/coordkr');
     kz    = h5read(filename,'/data/grid/coordkz');
-%     p     = h5read(filename,['/data/var5d/',variablename,'/coordp']);
-%     j     = h5read(filename,['/data/var5d/',variablename,'/coordj']);
-%     kr    = h5read(filename,['/data/var5d/',variablename,'/coordkr']);
-%     kz    = h5read(filename,['/data/var5d/',variablename,'/coordkz']);
 
     dt    = h5readatt(filename,'/data/input','dt');
     cstart= h5readatt(filename,'/data/input','start_iframe5d'); 
diff --git a/matlab/load_grid_data.m b/matlab/load_grid_data.m
new file mode 100644
index 0000000000000000000000000000000000000000..431815777c1b684f2cabe20512dfb3467eb40a3d
--- /dev/null
+++ b/matlab/load_grid_data.m
@@ -0,0 +1,9 @@
+function [ pe, je, pi, ji, kr, kz ] = load_grid_data( filename )
+%LOAD_GRID_DATA stored in a hdf5 result file from HeLaZ
+    pe    = h5read(filename,'/data/grid/coordp_e');
+    je    = h5read(filename,'/data/grid/coordj_e');
+    pi    = h5read(filename,'/data/grid/coordp_i');
+    ji    = h5read(filename,'/data/grid/coordj_i');
+    kr    = h5read(filename,'/data/grid/coordkr');
+    kz    = h5read(filename,'/data/grid/coordkz');
+end
\ No newline at end of file
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 0e47f493d87e33f4fa216db843facb4fd8c4e61b..c176b707a958f90b19e47f1168f3be2ef426a3d0 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -13,6 +13,12 @@ NZ      = h5readatt(filename,'/data/input','nz');
 L       = h5readatt(filename,'/data/input','Lr');
 CLOS    = h5readatt(filename,'/data/input','CLOS');
 MU = str2num(filename(end-18:end-14));
+W_GAMMA   = h5readatt(filename,'/data/input','write_gamma') == 'y';
+W_PHI     = h5readatt(filename,'/data/input','write_phi')   == 'y';
+W_NA00    = h5readatt(filename,'/data/input','write_Na00')  == 'y';
+W_NAPJ    = h5readatt(filename,'/data/input','write_Napj')  == 'y';
+W_SAPJ    = h5readatt(filename,'/data/input','write_Sapj')  == 'y';
+
 if NON_LIN == 'y'
     NON_LIN = 1;
 else
diff --git a/matlab/load_results.m b/matlab/load_results.m
index 1cb6be5de7e342181193f0dc03e3f53128a61dc2..fc41f9c4f13a013dc996d153051a3ca1d61c0dac 100644
--- a/matlab/load_results.m
+++ b/matlab/load_results.m
@@ -4,12 +4,30 @@ disp(['Loading ',filename])
 % Loading from output file
 CPUTIME   = h5readatt(filename,'/data/input','cpu_time');
 DT_SIM    = h5readatt(filename,'/data/input','dt');
-[Nipj, Pi, Ji, kr, kz, Ts5D, dt5D] = load_5D_data(filename, 'moments_i');
-[Nepj, Pe, Je                    ] = load_5D_data(filename, 'moments_e');
-[Ni00, kr, kz, Ts2D, dt2D]         = load_2D_data(filename, 'Ni00');
- Ne00                              = load_2D_data(filename, 'Ne00');
 
-PHI                            = load_2D_data(filename, 'phi');
+[Pe, Je, Pi, Ji, kr, kz] = load_grid_data(filename);
 
-% Sipj    = load_5D_data(filename, 'Sipj');
-% Sepj    = load_5D_data(filename, 'Sepj');
+if W_GAMMA
+    [ GGAMMA_RI, Ts0D, dt0D] = load_0D_data(filename, 'gflux_ri');
+      PGAMMA_RI              = load_0D_data(filename, 'pflux_ri');
+end
+
+if W_PHI
+    [ PHI, Ts2D, dt2D] = load_2D_data(filename, 'phi');
+end
+
+if W_NA00
+    [Ni00, Ts2D, dt2D] = load_2D_data(filename, 'Ni00');
+     Ne00              = load_2D_data(filename, 'Ne00');
+end
+
+
+if W_NAPJ
+    [Nipj, Ts5D, dt5D] = load_5D_data(filename, 'moments_i');
+    [Nepj            ] = load_5D_data(filename, 'moments_e');
+end
+
+if W_SAPJ
+    [Sipj, Ts5D, dt5D] = load_5D_data(filename, 'Sipj');
+     Sepj              = load_5D_data(filename, 'Sepj');
+end
diff --git a/matlab/setup.m b/matlab/setup.m
index 6f42c2b2fed67f53adf8ea9068bffc55151a54e8..a0da7b856f71d639f93a4e6076366aa5b7b052e4 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -97,10 +97,15 @@ OUTPUTS.nsave_1d = -1;
 OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT);
 OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT);
 OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT);
-OUTPUTS.write_doubleprecision = '.false.';
-OUTPUTS.resfile0      = '''outputs''';
-OUTPUTS.rstfile0      = '''checkpoint''';
-OUTPUTS.job2load      = JOB2LOAD;
+if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end;
+if W_GAMMA;  OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end;
+if W_PHI;    OUTPUTS.write_phi   = '.true.'; else; OUTPUTS.write_phi   = '.false.';end;
+if W_NA00;   OUTPUTS.write_Na00  = '.true.'; else; OUTPUTS.write_Na00  = '.false.';end;
+if W_NAPJ;   OUTPUTS.write_Napj  = '.true.'; else; OUTPUTS.write_Napj  = '.false.';end;
+if W_SAPJ;   OUTPUTS.write_Sapj  = '.true.'; else; OUTPUTS.write_Sapj  = '.false.';end;
+OUTPUTS.resfile0    = '''outputs''';
+OUTPUTS.rstfile0    = '''checkpoint''';
+OUTPUTS.job2load    = JOB2LOAD;
 %% Create directories
 if ~exist(SIMDIR, 'dir')
    mkdir(SIMDIR)
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 1e2f98da2a8fb75ab77763c20f26970f1a45884b..b5045f92035cfeebf40536d0158cf99a11b5dfaf 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -30,6 +30,11 @@ fprintf(fid,['  nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']);
 fprintf(fid,['  nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']);
 fprintf(fid,['  nsave_cp = ', num2str(OUTPUTS.nsave_cp),'\n']);
 fprintf(fid,['  write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
+fprintf(fid,['  write_gamma = ', OUTPUTS.write_gamma,'\n']);
+fprintf(fid,['  write_phi   = ', OUTPUTS.write_phi,'\n']);
+fprintf(fid,['  write_Na00 = ', OUTPUTS.write_Na00,'\n']);
+fprintf(fid,['  write_Napj = ', OUTPUTS.write_Napj,'\n']);
+fprintf(fid,['  write_Sapj = ', OUTPUTS.write_Sapj,'\n']);
 fprintf(fid,['  resfile0      = ', OUTPUTS.resfile0,'\n']);
 fprintf(fid,['  rstfile0      = ', OUTPUTS.rstfile0,'\n']);
 fprintf(fid,['  job2load      = ', num2str(OUTPUTS.job2load),'\n']);
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 15ca6353722e2d84501a5bfe7b305a0c226eeb74..ecfb4165852ec6a37b07f6561e28396160c4c89e 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -13,6 +13,7 @@ MODULE basic
 
   INTEGER :: comm0                 ! Default communicator with a topology
   INTEGER :: commp, commr          ! Communicators for 1-dim cartesian subgrids of comm0
+  INTEGER :: commr_p0              ! Communicators along kr for only rank 0 on p
 
   INTEGER :: jobnum  = 0           ! Job number
   INTEGER :: step    = 0           ! Calculation step of this run
@@ -24,11 +25,13 @@ MODULE basic
   INTEGER :: ierr                  ! flag for MPI error
   INTEGER :: my_id                 ! Rank in COMM_WORLD
   INTEGER :: num_procs             ! number of MPI processes
-  INTEGER :: num_procs_p, num_procs_kr              ! Number of processes in p and r
+  INTEGER :: num_procs_p           ! Number of processes in p
+  INTEGER :: num_procs_kr          ! Number of processes in r
   INTEGER :: rank_0, rank_p, rank_r! Ranks in comm0, commp, commr
-  INTEGER :: nbr_L, nbr_R ! Left and right neighbours (along p)
-  INTEGER :: nbr_T, nbr_B ! Top and bottom neighbours (along kr)
+  INTEGER :: nbr_L, nbr_R          ! Left and right neighbours (along p)
+  INTEGER :: nbr_T, nbr_B          ! Top and bottom neighbours (along kr)
 
+  INTEGER :: iframe0d              ! counting the number of times 0d datasets are outputed (for diagnose)
   INTEGER :: iframe1d              ! counting the number of times 1d datasets are outputed (for diagnose)
   INTEGER :: iframe2d              ! counting the number of times 2d datasets are outputed (for diagnose)
   INTEGER :: iframe3d              ! counting the number of times 3d datasets are outputed (for diagnose)
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index f4d9a2802236afaf410c6c918d06e3e2b78cea98..40094b6d813c60689ae2e314e151d02301ee1ec8 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -17,7 +17,8 @@ SUBROUTINE diagnose(kstep)
 
   INTEGER, INTENT(in) :: kstep
   INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank, dims(1) = (/0/)
+  INTEGER :: rank = 0
+  INTEGER :: dims(1) = (/0/)
   INTEGER :: cp_counter = 0
   CHARACTER(len=256) :: str, fname,test_
   ! putarr(...,pardim=1) does not work for 2D domain decomposition
@@ -44,18 +45,6 @@ SUBROUTINE diagnose(kstep)
 
      !  Data group
      CALL creatg(fidres, "/data", "data")
-     CALL creatg(fidres, "/data/var2d", "2d profiles")
-     CALL creatg(fidres, "/data/var5d", "5d profiles")
-
-     ! Initialize counter of number of saves for each category
-     IF (cstep==0) THEN
-         iframe2d=0
-     ENDIF
-     CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
-     IF (cstep==0) THEN
-         iframe5d=0
-     END IF
-     CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
 
      !  File group
      CALL creatg(fidres, "/files", "files")
@@ -96,31 +85,71 @@ SUBROUTINE diagnose(kstep)
       endif
      ENDIF
 
-     !  var2d group (electro. pot., Ni00 moment)
-     rank = 0
-     CALL creatd(fidres, rank, dims,  "/data/var2d/time",     "Time t*c_s/R")
-     CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
+     ! Grid info
      CALL creatg(fidres, "/data/grid", "Grid data")
      CALL putarr(fidres, "/data/grid/coordkr", krarray_full(1:nkr),"kr*rho_s0", ionode=0)
      CALL putarr(fidres, "/data/grid/coordkz", kzarray_full(1:nkz),"kz*rho_s0", ionode=0)
-     CALL putarr(fidres, "/data/grid/coordp" , parray_e_full(1:pmaxe+1), "p_e", ionode=0)
-     CALL putarr(fidres, "/data/grid/coordj" , jarray_e_full(1:jmaxe+1), "j_e", ionode=0)
+     CALL putarr(fidres, "/data/grid/coordp_e" , parray_e_full(1:pmaxe+1), "p_e", ionode=0)
+     CALL putarr(fidres, "/data/grid/coordj_e" , jarray_e_full(1:jmaxe+1), "j_e", ionode=0)
+     CALL putarr(fidres, "/data/grid/coordp_i" , parray_i_full(1:pmaxe+1), "p_i", ionode=0)
+     CALL putarr(fidres, "/data/grid/coordj_i" , jarray_i_full(1:jmaxe+1), "j_i", ionode=0)
+
+     !  var0d group (gyro transport)
+     IF (nsave_0d .GT. 0) THEN
+      CALL creatg(fidres, "/data/var0d", "0d profiles")
+      CALL creatd(fidres, rank, dims,  "/data/var0d/time",     "Time t*c_s/R")
+      CALL creatd(fidres, rank, dims, "/data/var0d/cstep", "iteration number")
+
+      IF (write_gamma) THEN
+        CALL creatd(fidres, rank, dims, "/data/var0d/gflux_ri", "Radial gyro ion transport")
+        CALL creatd(fidres, rank, dims, "/data/var0d/pflux_ri", "Radial part ion transport")
+      ENDIF
+      IF (cstep==0) THEN
+        iframe0d=0
+      ENDIF
+      CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
+     END IF
 
+
+     !  var2d group (electro. pot., Ni00 moment)
      IF (nsave_2d .GT. 0) THEN
+      CALL creatg(fidres, "/data/var2d", "2d profiles")
+      CALL creatd(fidres, rank, dims,  "/data/var2d/time",     "Time t*c_s/R")
+      CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
+
+      IF (write_phi) CALL creatg(fidres, "/data/var2d/phi", "phi")
+
+      IF (write_Na00) THEN
        CALL creatg(fidres, "/data/var2d/Ne00", "Ne00")
        CALL creatg(fidres, "/data/var2d/Ni00", "Ni00")
-       CALL creatg(fidres, "/data/var2d/phi", "phi")
+      ENDIF
+
+      IF (cstep==0) THEN
+        iframe2d=0
+      ENDIF
+      CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
      END IF
 
      !  var5d group (moments)
-     rank = 0
-     CALL creatd(fidres, rank, dims,  "/data/var5d/time",     "Time t*c_s/R")
-     CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
      IF (nsave_5d .GT. 0) THEN
-       CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
-       CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
-      !  CALL creatg(fidres, "/data/var5d/Sepj", "Sepj")
-      !  CALL creatg(fidres, "/data/var5d/Sipj", "Sipj")
+       CALL creatg(fidres, "/data/var5d", "5d profiles")
+       CALL creatd(fidres, rank, dims,  "/data/var5d/time",     "Time t*c_s/R")
+       CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
+
+       IF (write_Napj) THEN
+        CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
+        CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
+       ENDIF
+
+       IF (write_Sapj) THEN
+        CALL creatg(fidres, "/data/var5d/moments_e", "Sipj")
+        CALL creatg(fidres, "/data/var5d/moments_i", "Sepj")
+       ENDIF
+
+       IF (cstep==0) THEN
+        iframe5d=0
+       END IF
+       CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
      END IF
 
      !  Add input namelist variables as attributes of /data/input, defined in srcinfo.h
@@ -130,8 +159,6 @@ SUBROUTINE diagnose(kstep)
      IF (my_id .EQ. 0) WRITE(*,*)    'HOST=', HOST
 
      WRITE(str,'(a,i2.2)') "/data/input"
-
-     rank=0
      CALL creatd(fidres, rank,dims,TRIM(str),'Input parameters')
      CALL attach(fidres, TRIM(str),     "version",  VERSION) !defined in srcinfo.h
      CALL attach(fidres, TRIM(str),      "branch",   BRANCH) !defined in srcinfo.h
@@ -140,6 +167,7 @@ SUBROUTINE diagnose(kstep)
      CALL attach(fidres, TRIM(str),        "host",     HOST) !defined in srcinfo.h
      CALL attach(fidres, TRIM(str),  "start_time",     time)
      CALL attach(fidres, TRIM(str), "start_cstep",    cstep-1)
+     CALL attach(fidres, TRIM(str), "start_iframe0d", iframe0d)
      CALL attach(fidres, TRIM(str), "start_iframe2d", iframe2d)
      CALL attach(fidres, TRIM(str), "start_iframe5d", iframe5d)
      CALL attach(fidres, TRIM(str),          "dt",       dt)
@@ -149,6 +177,11 @@ SUBROUTINE diagnose(kstep)
      CALL attach(fidres, TRIM(str),       "Nproc",   num_procs)
      CALL attach(fidres, TRIM(str),       "Np_p" , num_procs_p)
      CALL attach(fidres, TRIM(str),       "Np_kr",num_procs_kr)
+     CALL attach(fidres, TRIM(str), "write_gamma",write_gamma)
+     CALL attach(fidres, TRIM(str),   "write_phi",write_phi)
+     CALL attach(fidres, TRIM(str),  "write_Na00",write_Na00)
+     CALL attach(fidres, TRIM(str),  "write_Napj",write_Napj)
+     CALL attach(fidres, TRIM(str),  "write_Sapj",write_Sapj)
 
      CALL grid_outputinputs(fidres, str)
 
@@ -242,12 +275,13 @@ END SUBROUTINE diagnose
 SUBROUTINE diagnose_0d
 
   USE basic
-  USE futils, ONLY: append
+  USE futils, ONLY: append, attach, getatt
   USE diagnostics_par
   USE prec_const
+  USE processing
 
   IMPLICIT NONE
-
+  ! Time measurement data
   CALL append(fidres, "/profiler/Tc_rhs",              tc_rhs,ionode=0)
   CALL append(fidres, "/profiler/Tc_adv_field",  tc_adv_field,ionode=0)
   CALL append(fidres, "/profiler/Tc_poisson",      tc_poisson,ionode=0)
@@ -257,6 +291,18 @@ SUBROUTINE diagnose_0d
   CALL append(fidres, "/profiler/Tc_comm",            tc_comm,ionode=0)
   CALL append(fidres, "/profiler/Tc_step",            tc_step,ionode=0)
   CALL append(fidres, "/profiler/time",                  time,ionode=0)
+  ! Processing data
+  CALL append(fidres,  "/data/var0d/time",           time,ionode=0)
+  CALL append(fidres, "/data/var0d/cstep", real(cstep,dp),ionode=0)
+  CALL getatt(fidres,      "/data/var0d/",       "frames",iframe2d)
+  iframe0d=iframe0d+1
+  CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
+  ! Ion transport data
+  IF (write_gamma) THEN
+    CALL compute_radial_ion_transport
+    CALL append(fidres, "/data/var0d/gflux_ri",gflux_ri,ionode=0)
+    CALL append(fidres, "/data/var0d/pflux_ri",pflux_ri,ionode=0)
+  ENDIF
 
 END SUBROUTINE diagnose_0d
 
@@ -282,78 +328,79 @@ SUBROUTINE diagnose_2d
   iframe2d=iframe2d+1
   CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
 
-  CALL write_field2d(phi (:,:), 'phi')
+  IF (write_phi) CALL write_field2d(phi (:,:), 'phi')
 
-  IF ( (ips_e .EQ. 1) .AND. (ips_i .EQ. 1) ) THEN
-    Ne00(ikrs:ikre,ikzs:ikze) = moments_e(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel)
-    Ni00(ikrs:ikre,ikzs:ikze) = moments_i(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel)
-  ENDIF
-
-  root = 0
-
-  !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
-  CALL MPI_COMM_RANK(commp,world_rank,ierr)
-  CALL MPI_COMM_SIZE(commp,world_size,ierr)
+  IF (write_Na00) THEN
+    IF ( (ips_e .EQ. 1) .AND. (ips_i .EQ. 1) ) THEN
+      Ne00(ikrs:ikre,ikzs:ikze) = moments_e(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel)
+      Ni00(ikrs:ikre,ikzs:ikze) = moments_i(ips_e,1,ikrs:ikre,ikzs:ikze,updatetlevel)
+    ENDIF
 
-  IF (world_size .GT. 1) THEN
-    !! Broadcast phi to the other processes on the same k range (communicator along p)
-    IF (world_rank .EQ. root) THEN
-      ! Fill the buffer
-      DO ikr = ikrs,ikre
-        DO ikz = ikzs,ikze
-          buffer(ikr,ikz) = Ne00(ikr,ikz)
+    root = 0
+    !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
+    CALL MPI_COMM_RANK(commp,world_rank,ierr)
+    CALL MPI_COMM_SIZE(commp,world_size,ierr)
+
+    IF (world_size .GT. 1) THEN
+      !! Broadcast phi to the other processes on the same k range (communicator along p)
+      IF (world_rank .EQ. root) THEN
+        ! Fill the buffer
+        DO ikr = ikrs,ikre
+          DO ikz = ikzs,ikze
+            buffer(ikr,ikz) = Ne00(ikr,ikz)
+          ENDDO
         ENDDO
-      ENDDO
-      ! Send it to all the other processes
-      DO i_ = 0,num_procs_p-1
-        IF (i_ .NE. world_rank) &
-        CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr)
-      ENDDO
-    ELSE
-      ! Recieve buffer from root
-      CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr)
-      ! Write it in phi
-      DO ikr = ikrs,ikre
-        DO ikz = ikzs,ikze
-          Ne00(ikr,ikz) = buffer(ikr,ikz)
+        ! Send it to all the other processes
+        DO i_ = 0,num_procs_p-1
+          IF (i_ .NE. world_rank) &
+          CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr)
+        ENDDO
+      ELSE
+        ! Recieve buffer from root
+        CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr)
+        ! Write it in phi
+        DO ikr = ikrs,ikre
+          DO ikz = ikzs,ikze
+            Ne00(ikr,ikz) = buffer(ikr,ikz)
+          ENDDO
         ENDDO
-      ENDDO
+      ENDIF
     ENDIF
-  ENDIF
 
-  CALL write_field2d(Ne00(ikrs:ikre,ikzs:ikze), 'Ne00')
+    CALL write_field2d(Ne00(ikrs:ikre,ikzs:ikze), 'Ne00')
 
-    !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
-  CALL MPI_COMM_RANK(commp,world_rank,ierr)
-  CALL MPI_COMM_SIZE(commp,world_size,ierr)
-
-  IF (world_size .GT. 1) THEN
-    !! Broadcast phi to the other processes on the same k range (communicator along p)
-    IF (world_rank .EQ. root) THEN
-      ! Fill the buffer
-      DO ikr = ikrs,ikre
-        DO ikz = ikzs,ikze
-          buffer(ikr,ikz) = Ni00(ikr,ikz)
+      !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
+    CALL MPI_COMM_RANK(commp,world_rank,ierr)
+    CALL MPI_COMM_SIZE(commp,world_size,ierr)
+
+    IF (world_size .GT. 1) THEN
+      !! Broadcast phi to the other processes on the same k range (communicator along p)
+      IF (world_rank .EQ. root) THEN
+        ! Fill the buffer
+        DO ikr = ikrs,ikre
+          DO ikz = ikzs,ikze
+            buffer(ikr,ikz) = Ni00(ikr,ikz)
+          ENDDO
         ENDDO
-      ENDDO
-      ! Send it to all the other processes
-      DO i_ = 0,num_procs_p-1
-        IF (i_ .NE. world_rank) &
-        CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr)
-      ENDDO
-    ELSE
-      ! Recieve buffer from root
-      CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr)
-      ! Write it in phi
-      DO ikr = ikrs,ikre
-        DO ikz = ikzs,ikze
-          Ni00(ikr,ikz) = buffer(ikr,ikz)
+        ! Send it to all the other processes
+        DO i_ = 0,num_procs_p-1
+          IF (i_ .NE. world_rank) &
+          CALL MPI_SEND(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, i_, 0, commp, ierr)
+        ENDDO
+      ELSE
+        ! Recieve buffer from root
+        CALL MPI_RECV(buffer, local_nkr * nkz , MPI_DOUBLE_COMPLEX, root, 0, commp, MPI_STATUS_IGNORE, ierr)
+        ! Write it in phi
+        DO ikr = ikrs,ikre
+          DO ikz = ikzs,ikze
+            Ni00(ikr,ikz) = buffer(ikr,ikz)
+          ENDDO
         ENDDO
-      ENDDO
+      ENDIF
     ENDIF
-  ENDIF
 
-  CALL write_field2d(Ni00(ikrs:ikre,ikzs:ikze), 'Ni00')
+    CALL write_field2d(Ni00(ikrs:ikre,ikzs:ikze), 'Ni00')
+  ENDIF
 
 CONTAINS
 
@@ -405,11 +452,15 @@ SUBROUTINE diagnose_5d
    iframe5d=iframe5d+1
    CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
 
-   CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,:,:,updatetlevel), 'moments_e')
-   CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,:,:,updatetlevel), 'moments_i')
+   IF (write_Napj) THEN
+    CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,:,:,updatetlevel), 'moments_e')
+    CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,:,:,updatetlevel), 'moments_i')
+   ENDIF
 
-  !  CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,:,:), 'Sepj')
-  !  CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,:,:), 'Sipj')
+   IF (write_Sapj) THEN
+     CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,:,:), 'Sepj')
+     CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,:,:), 'Sipj')
+   ENDIF
 
  CONTAINS
 
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index a0fb7a60afa14084d1aa0dd480b9e8fe9c1e44d4..0f743896ec319531b828ae241127f4b5cd69187b 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -5,7 +5,10 @@ MODULE diagnostics_par
   IMPLICIT NONE
   PRIVATE
 
-  LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision=.FALSE.
+  LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision = .FALSE.
+  LOGICAL, PUBLIC, PROTECTED :: write_gamma
+  LOGICAL, PUBLIC, PROTECTED :: write_phi,  write_Na00 
+  LOGICAL, PUBLIC, PROTECTED :: write_Napj, write_Sapj
 
   INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d, nsave_cp
 
@@ -31,7 +34,8 @@ CONTAINS
     IMPLICIT NONE
 
     NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d, nsave_cp
-    NAMELIST /OUTPUT_PAR/ write_doubleprecision
+    NAMELIST /OUTPUT_PAR/ write_doubleprecision, write_gamma, write_phi
+    NAMELIST /OUTPUT_PAR/ write_Na00, write_Napj, write_Sapj
     NAMELIST /OUTPUT_PAR/ resfile0, rstfile0, job2load
 
     READ(lu_in,output_par)
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..a99c56719b222c9bbffc722c125a5577facd9d33
--- /dev/null
+++ b/src/processing_mod.F90
@@ -0,0 +1,65 @@
+MODULE processing
+    ! contains the Hermite-Laguerre collision operators. Solved using COSOlver.
+    USE basic
+    USE prec_const
+    USE grid
+    implicit none
+
+    REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri
+    
+    PUBLIC :: compute_radial_ion_transport
+    
+CONTAINS
+
+SUBROUTINE compute_radial_ion_transport
+    USE fields,           ONLY : moments_i, phi
+    USE array,            ONLY : kernel_i
+    USE time_integration, ONLY : updatetlevel
+    IMPLICIT NONE
+    COMPLEX(dp) :: pflux_local, gflux_local
+    REAL(dp)    :: kz_, buffer(1:2)
+    INTEGER     :: i_, world_rank, world_size, root
+
+    pflux_local = 0._dp
+    gflux_local = 0._dp
+    IF(ips_i .EQ. 1) THEN
+        ! Loop to compute gamma_kr = sum_kz sum_j -i k_z Kernel_j Ni00 * phi
+        DO ikz = ikzs,ikze
+            kz_ = kzarray(ikz)
+            DO ikr = ikrs,ikre
+                gflux_local = gflux_local - &
+                    imagu * kz_ * moments_i(1,1,ikr,ikz,updatetlevel) * CONJG(phi(ikr,ikz))
+                DO ij = ijs_i, ije_i
+                    pflux_local = pflux_local - &
+                        imagu * kz_ * kernel_i(ij,ikr,ikz) * moments_i(1,ij,ikr,ikz,updatetlevel) * CONJG(phi(ikr,ikz))
+                ENDDO
+            ENDDO
+        ENDDO
+
+        buffer(1) = REAL(gflux_local)
+        buffer(2) = REAL(pflux_local)
+
+        root = 0
+        !Gather manually among the rank_p=0 processes and perform the sum
+        gflux_ri = 0
+        pflux_ri = 0
+        IF (num_procs_kr .GT. 1) THEN
+            !! Everyone sends its local_sum to root = 0
+            IF (rank_r .NE. root) THEN
+                CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, commr, ierr)
+            ELSE
+                ! Recieve from all the other processes
+                DO i_ = 0,num_procs_kr-1
+                    IF (i_ .NE. rank_r) &
+                        CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, commr, MPI_STATUS_IGNORE, ierr)
+                        gflux_ri = gflux_ri + buffer(1)
+                        pflux_ri = pflux_ri + buffer(2)
+                ENDDO
+            ENDIF
+        ENDIF
+    ENDIF
+
+
+END SUBROUTINE compute_radial_ion_transport
+
+END MODULE processing
\ No newline at end of file
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index a9b122602be791eb2ec67710b7e80174fa8bf692..bd544c9959a2d23a0708e33dacb398ead31f7795 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -104,7 +104,7 @@ PFlux_ri  = zeros(1,Ns5D);      % Particle   flux
 GFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Ni00.*conj(PHI),1),2)))*(2*pi/Nr/Nz)^2;
 PFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Np_i.*conj(PHI_Ts5D),1),2)))*(2*pi/Nr/Nz)^2;
 
-phi_avg  = zeros(1,Ns2D);        % Time evol. of the norm of phi
+phi_max  = zeros(1,Ns2D);        % Time evol. of the norm of phi
 Ne_norm  = zeros(Npe,Nje,Ns5D);  % Time evol. of the norm of Napj
 Ni_norm  = zeros(Npi,Nji,Ns5D);  % .
 
@@ -112,7 +112,7 @@ Ddr = 1i*KR; Ddz = 1i*KZ; lapl   = Ddr.^2 + Ddz.^2;
 
 for it = 1:numel(Ts2D) % Loop over 2D arrays
     NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
-    phi_avg(it)   = real(squeeze(PH_(ikr0,ikz0)))/2;
+    phi_max(it)   = max(max(squeeze(phi(:,:,it))));
     ExB(it)       = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz))));
     GFlux_ri(it)  = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
     GFlux_zi(it)  = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
@@ -151,7 +151,7 @@ if 1
 %% Time evolutions and growth rate
 fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS];
 set(gcf, 'Position',  [100, 100, 900, 800])
-    subplot(221); 
+    subplot(421); 
     for ip = 1:Npe
         for ij = 1:Nje
             plt      = @(x) squeeze(x(ip,ij,:));
@@ -163,7 +163,7 @@ set(gcf, 'Position',  [100, 100, 900, 800])
         end
     end
     grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
-    subplot(222)
+    subplot(423)
     for ip = 1:Npi
         for ij = 1:Nji
             plt      = @(x) squeeze(x(ip,ij,:));
@@ -174,7 +174,14 @@ set(gcf, 'Position',  [100, 100, 900, 800])
                 'Color',clr,'LineStyle',lstyle{1}); hold on;
         end
     end
-    grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$');
+    grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$')
+    subplot(222)
+        semilogy(Ts0D,GGAMMA_RI*(2*pi/Nr/Nz)^2); hold on;
+%         plot(Ts2D,GFLUX_RI)
+        plot(Ts0D,PGAMMA_RI*(2*pi/Nr/Nz)^2);
+%         plot(Ts5D,PFLUX_RI,'--');
+        legend(['Gyro. flux';'Part. flux']);
+        grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$')
     subplot(223)
         plot(kz,g_(1,:),'-','DisplayName','$\gamma$'); hold on;
         grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); %legend('show');
@@ -185,11 +192,11 @@ set(gcf, 'Position',  [100, 100, 900, 800])
             plotname = '$\langle\phi\rangle_{r,z}(t)$';
             clr      = line_colors(min(ip,numel(line_colors(:,1))),:);
             lstyle   = line_styles(min(ij,numel(line_styles)));
-            plot(Ts2D,phi_avg,'DisplayName',plotname,...
+            plot(Ts2D,phi_max,'DisplayName',plotname,...
                 'Color',clr,'LineStyle',lstyle{1}); hold on;
         end
     end
-    grid on; xlabel('$t c_s/R$'); ylabel('$\tilde\phi_{00}/2$'); %legend('show');
+    grid on; xlabel('$t c_s/R$'); ylabel('$\max_{r,z}(\phi)$'); %legend('show');
 % suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]);
 save_figure
 end
diff --git a/wk/local_run.m b/wk/local_run.m
index dfe9ccea17e3d27ad348abbf003d771ff9339be5..9d897ea7e55bb261a8cf79dcd5ce8044759e3c9e 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -18,13 +18,13 @@ JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS
 TMAX    = 100;  % Maximal time unit
 DT      = 2e-2;   % Time step
-SPS0D   = 1/DT;    % Sampling per time unit for profiler
+SPS0D   = 1;    % Sampling per time unit for profiler
 SPS2D   = 1/2;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/4;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints/10
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
-%% OPTIONS
+%% OPTIONS AND NAMING
 % SIMID   = ['local_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
 % SIMID   = sprintf(SIMID,NU);
 % SIMID   = 'test_init_phi';  % Name of the simulation
@@ -33,6 +33,13 @@ CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Do
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
 KERN    = 0;   % Kernel model (0 : GK)
 INIT_PHI= 1;   % Start simulation with a noisy phi and moments
+%% OUTPUTS
+W_DOUBLE = 0;
+W_GAMMA  = 1;
+W_PHI    = 1;
+W_NA00   = 1;
+W_NAPJ   = 1;
+W_SAPJ   = 0;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% unused
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 92e35314a4783a068a9be1494fd9193a0f663cab..48bab70fcd2d9c2a6a35147a3e86f3f2d46f843a 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -16,29 +16,36 @@ NU      = 0.1;   % Collision frequency
 ETAB    = 0.6;   % Magnetic gradient
 NU_HYP  = 0.1;   % Hyperdiffusivity coefficient
 %% GRID PARAMETERS
-N       = 050;   % Frequency gridpoints (Nkr = N/2)
-L       = 100;   % Size of the squared frequency domain
-P       = 10;    % Electron and Ion highest Hermite polynomial degree
-J       = 01;    % Electron and Ion highest Laguerre polynomial degree
+N       = 200;   % Frequency gridpoints (Nkr = N/2)
+L       = 120;   % Size of the squared frequency domain
+P       = 20;    % Electron and Ion highest Hermite polynomial degree
+J       = 08;    % Electron and Ion highest Laguerre polynomial degree
 MU_P    = 0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARAMETERS
-TMAX    = 50;  % Maximal time unit
-DT      = 2e-2;  % Time step
+TMAX    = 200;  % Maximal time unit
+DT      = 1e-2;  % Time step
 SPS0D   = 1;      % Sampling per time unit for profiler
-SPS2D   = 1/2;   % Sampling per time unit for 2D arrays
-SPS5D   = 1/10;  % Sampling per time unit for 5D arrays
+SPS2D   = 1/10;   % Sampling per time unit for 2D arrays
+SPS5D   = 0/50;  % 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   = ['HeLaZ_v2.4_2_12_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
+SIMID   = ['HeLaZ_v2.4_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
 % SIMID   = 'Marconi_test';  % Name of the simulation
 SIMID   = sprintf(SIMID,NU);
 CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
 KERN    = 0;   % Kernel model (0 : GK)
 INIT_PHI= 1;   % Start simulation with a noisy phi and moments
+%% OUTPUTS
+W_DOUBLE = 0;
+W_GAMMA  = 0;
+W_PHI    = 1;
+W_NA00   = 1;
+W_NAPJ   = 0;
+W_SAPJ   = 0;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% fixed parameters (for current study)