diff --git a/.gitignore b/.gitignore
index c03a314529fbe9cbacdbedcf7bb6102da2ceac79..cbb062382e2da7a4de5caef39e883b00da2938fe 100644
--- a/.gitignore
+++ b/.gitignore
@@ -38,3 +38,4 @@ src/srcinfo/srcinfo.h
 local/
 */*.sh
 Gallery/
+.vscode/settings.json
diff --git a/Makefile b/Makefile
index 382dd6310315ee132848398e4273aa513478dab4..80ae86717d3d9a6419f5f2b8f0560935a45693f9 100644
--- a/Makefile
+++ b/Makefile
@@ -68,7 +68,7 @@ $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR
 $(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/geometry_mod.o \
 $(OBJDIR)/main.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs_mod.o \
 $(OBJDIR)/numerics_mod.o $(OBJDIR)/poisson.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)/parallel_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
 
@@ -81,40 +81,57 @@ $(OBJDIR)/utility_mod.o
  $(EDBG): $(FOBJ)
 	$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
 
- $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o
+ $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o \
+   $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o \
+	 $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@
 
  $(OBJDIR)/array_mod.o : src/array_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/array_mod.F90 -o $@
 
- $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/geometry_mod.o  $(OBJDIR)/grid_mod.o $(OBJDIR)/numerics_mod.o
+ $(OBJDIR)/auxval.o : src/auxval.F90 $(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o \
+   $(OBJDIR)/model_mod.o $(OBJDIR)/geometry_mod.o  $(OBJDIR)/grid_mod.o $(OBJDIR)/numerics_mod.o\
+   $(OBJDIR)/parallel_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
 
  $(OBJDIR)/basic_mod.o : src/basic_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/basic_mod.F90 -o $@
 
- $(OBJDIR)/calculus_mod.o : src/calculus_mod.F90  $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
+ $(OBJDIR)/calculus_mod.o : src/calculus_mod.F90  $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
+   $(OBJDIR)/grid_mod.o $(OBJDIR)/parallel_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/calculus_mod.F90 -o $@
 
- $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o
-		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@
+ $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@
 
- $(OBJDIR)/closure_mod.o : src/closure_mod.F90 $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o
-	 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/closure_mod.F90 -o $@
+ $(OBJDIR)/closure_mod.o : src/closure_mod.F90 $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/closure_mod.F90 -o $@
 
- $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
-	 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@
+ $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o \
+	 $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@
 
- $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
-		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -o $@
+ $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o\
+	 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -o $@
 
- $(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/geometry_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
+ $(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/geometry_mod.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)/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
+ $(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\
+	 $(OBJDIR)/parallel_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
+ $(OBJDIR)/diagnostics_par_mod.o : src/diagnostics_par_mod.F90 $(OBJDIR)/prec_const_mod.o \
+   $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/diagnostics_par_mod.F90 -o $@
 
  $(OBJDIR)/endrun.o : src/endrun.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
@@ -123,68 +140,93 @@ $(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
+ $(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)/geometry_mod.o : src/geometry_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o
+ $(OBJDIR)/geometry_mod.o : src/geometry_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o\
+   $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/utility_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_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 $@
+ $(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 $@
 
- $(OBJDIR)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o
-		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@
+ $(OBJDIR)/grid_mod.o : src/grid_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o\
+   $(OBJDIR)/prec_const_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/grid_mod.F90 -o $@
 
- $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
+ $(OBJDIR)/inital.o : src/inital.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/fields_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/restarts_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/inital.F90 -o $@
 
- $(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/initial_par_mod.o : src/initial_par_mod.F90 $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/initial_par_mod.F90 -o $@
 
  $(OBJDIR)/main.o : src/main.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/main.F90 -o $@
 
- $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/collision_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/grid_mod.o
+ $(OBJDIR)/memory.o : src/memory.F90 $ $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/collision_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
+	 $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
 
  $(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
 
- $(OBJDIR)/moments_eq_rhs_mod.o : src/moments_eq_rhs_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/moments_eq_rhs_mod.o : src/moments_eq_rhs_mod.F90 $(OBJDIR)/array_mod.o \
+   $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs_mod.F90 -o $@
 
- $(OBJDIR)/numerics_mod.o : src/numerics_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/utility_mod.o
+ $(OBJDIR)/numerics_mod.o : src/numerics_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o \
+  $(OBJDIR)/coeff_mod.o $(OBJDIR)/utility_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/numerics_mod.F90 -o $@
 
- $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o \
+  $(OBJDIR)/grid_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o \
+	$(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/parallel_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/poisson.F90 -o $@
 
- $(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o
+ $(OBJDIR)/parallel_mod.o : src/parallel_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
+   $(OBJDIR)/grid_mod.o
+	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/parallel_mod.F90 -o $@
+
+ $(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/coeff_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppexit.F90 -o $@
 
- $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/ppinit.o : src/ppinit.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o \
+   $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o \
+	 $(OBJDIR)/basic_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ppinit.F90 -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
+ $(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
+ $(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 $@
 
- $(OBJDIR)/restarts_mod.o : src/restarts_mod.F90  $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
+ $(OBJDIR)/restarts_mod.o : src/restarts_mod.F90  $(OBJDIR)/diagnostics_par_mod.o \
+   $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/restarts_mod.F90 -o $@
 
- $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
+ $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o \
+   $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
 
  $(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/tesend.F90 -o $@
 
- $(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
+ $(OBJDIR)/time_integration_mod.o : src/time_integration_mod.F90 $(OBJDIR)/basic_mod.o \
+   $(OBJDIR)/prec_const_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/time_integration_mod.F90 -o $@
 
- $(OBJDIR)/utility_mod.o : src/utility_mod.F90  $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o
+ $(OBJDIR)/utility_mod.o : src/utility_mod.F90  $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o \
+   $(OBJDIR)/grid_mod.o
 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/utility_mod.F90 -o $@
diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m
index 4c7ee6716b8194d336649283d2fc54a8a5d3c25b..58b1baef595e674b081bce2c28b24c75f693ab3f 100644
--- a/matlab/plot/create_film.m
+++ b/matlab/plot/create_film.m
@@ -1,6 +1,6 @@
 function create_film(DATA,OPTIONS,format)
 %% Plot options
-FPS = 30; DELAY = 1/FPS;
+FPS = 16; DELAY = 1/FPS;
 BWR = 1; NORMALIZED = 1;
 if ~strcmp(OPTIONS.PLAN,'sx')
     T = DATA.Ts3D;
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 0e81e20b28b18b9b2a30cba59819db08ec3ffe31..93e0c5cb2ea39338bdead42162b7b7d422339701 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -83,9 +83,9 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
     toplot = process_field(DATA,OPTIONS);
     f2plot = toplot.FIELD;
     
-    clim = max(max(max(abs(plt(f2plot(:,:,its3D:ite3D))))));
+    clim = max(max(max(abs(plt(f2plot(:,:,:))))));
     subplot(313)
-        [TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
+        [TY,TX] = meshgrid(DATA.x,DATA.Ts3D(toplot.FRAMES));
         pclr = pcolor(TX,TY,squeeze(plt(f2plot))'); 
         set(pclr, 'edgecolor','none'); 
         legend(['$\langle ',OPTIONS.ST_FIELD,' \rangle_{y,z}$']) %colorbar;
@@ -103,14 +103,14 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
 %% Zonal vs NZonal energies    
     subplot(312)
     it0 = 1; itend = Ns3D;
-    trange = it0:itend;
+    trange = toplot.FRAMES;
     plt1 = @(x) x;%-x(1);
-    plt2 = @(x) x./max((x(its3D:ite3D)));
+    plt2 = @(x) x./max((x(:)));
     toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
 %     plty = @(x) x(500:end)./max(squeeze(x(500:end)));
         yyaxis left
 %         plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
-        plot(plt1(DATA.Ts3D(trange)),plt2(toplot(trange)),'DisplayName','Sum $A^2$');
+        plot(plt1(DATA.Ts3D(trange)),plt2(toplot(:)),'DisplayName','Sum $A^2$');
         ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
         yyaxis right
         plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 755272047c1a3f45eceb293f9d9b6ce8bcf67c32..3dca4e309433da22a41cbeae65d4477232297d87 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -63,7 +63,6 @@ MODULE array
   ! Poisson operator (ikx,iky,iz)
   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
 
-  !! Diagnostics
   ! Gyrocenter density for electron and ions (ikx,iky,iz)
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ne00
   COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ni00
diff --git a/src/auxval.F90 b/src/auxval.F90
index 34164c591b32c44fc2ccb4902bad0c4a1948b364..8ade415977c1a31c36045127a2c0dd79127039f8 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -9,17 +9,18 @@ subroutine auxval
   use prec_const
   USE numerics
   USE geometry
+  USE parallel, ONLY: init_parallel_var
   IMPLICIT NONE
 
   INTEGER :: irows,irowe, irow, icol, i_
   IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ==='
 
   IF (LINEARITY .NE. 'linear') THEN
-    write(*,*) 'FFTW3 y-grid distribution'
+    IF (my_id .EQ. 0) write(*,*) 'FFTW3 y-grid distribution'
     CALL init_grid_distr_and_plans(Nx,Ny)
   ELSE
     CALL init_1Dgrid_distr
-    write(*,*) 'Manual y-grid distribution'
+    IF (my_id .EQ. 0) write(*,*) 'Manual y-grid distribution'
   ENDIF
   ! Init the grids
   CALL set_pgrid ! parallel kin (MPI distributed)
@@ -35,6 +36,8 @@ subroutine auxval
 
   CALL memory ! Allocate memory for global arrays
 
+  CALL init_parallel_var
+
   CALL eval_magnetic_geometry ! precompute coeff for lin equation
 
   CALL compute_lin_coeff ! precompute coeff for lin equation and geometry
@@ -78,16 +81,10 @@ subroutine auxval
   ENDDO
   CALL mpi_barrier(MPI_COMM_WORLD, ierr)
 
-  ! IF((my_id.EQ.0) .AND. (CLOS .EQ. 1)) THEN
-  ! IF(KIN_E) &
-  ! write(*,*) 'Closure = 1 -> Maximal Nepj degree is min(Pmaxe,2*Jmaxe+1): De = ', dmaxi
-  ! write(*,*) 'Closure = 1 -> Maximal Nipj degree is min(Pmaxi,2*Jmaxi+1): Di = ', dmaxi
-  ! ENDIF
-  ! DO ip = ips_i,ipe_i
-  !   DO ij = ijs_i,ije_i
-  !     IF((parray_i(ip)+2*jarray_i(ij) .LE. dmaxi) .AND. (rank_ky + rank_z .EQ. 0))&
-  !     print*, '(',parray_i(ip),',',jarray_i(ij),')'
-  !   ENDDO
-  ! ENDDO
+  IF((my_id.EQ.0) .AND. (CLOS .EQ. 1)) THEN
+  IF(KIN_E) &
+  write(*,*) 'Closure = 1 -> Maximal Nepj degree is min(Pmaxe,2*Jmaxe+1): De = ', dmaxi
+  write(*,*) 'Closure = 1 -> Maximal Nipj degree is min(Pmaxi,2*Jmaxi+1): Di = ', dmaxi
+  ENDIF
 
 END SUBROUTINE auxval
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index f1d67f07028219981f609f23d3cd5ccd8b529d07..bc79d1b238f85b376028ec3eeb8c349f126a032b 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -12,12 +12,17 @@ MODULE basic
   real(dp) :: time   = 0           ! Current simulation time (Init from restart file)
 
   INTEGER :: comm0                 ! Default communicator with a topology
+  INTEGER :: group0                ! Default group with a topology
   INTEGER :: rank_0                ! Ranks in comm0
   ! Communicators for 1-dim cartesian subgrids of comm0
   INTEGER :: comm_p, comm_ky, comm_z
   INTEGER :: rank_p, rank_ky, rank_z! Ranks
   INTEGER :: comm_pz,  rank_pz      ! 2D comm for N_a(p,j,z) output (mspfile)
   INTEGER :: comm_kyz, rank_kyz     ! 2D comm for N_a(p,j,z) output (mspfile)
+  INTEGER :: comm_ky0, rank_ky0     ! comm along ky with p=0
+  INTEGER :: comm_z0,  rank_z0      ! comm along z  with p=0
+
+  INTEGER :: group_ky0, group_z0
 
   INTEGER :: jobnum  = 0           ! Job number
   INTEGER :: step    = 0           ! Calculation step of this run
@@ -31,6 +36,8 @@ MODULE basic
   INTEGER :: num_procs_p           ! Number of processes in p
   INTEGER :: num_procs_ky          ! Number of processes in r
   INTEGER :: num_procs_z           ! Number of processes in z
+  INTEGER :: num_procs_pz          ! Number of processes in pz comm
+  INTEGER :: num_procs_kyz         ! Number of processes in kyz comm
   INTEGER :: nbr_L, nbr_R          ! Left and right neighbours (along p)
   INTEGER :: nbr_T, nbr_B          ! Top and bottom neighbours (along kx)
   INTEGER :: nbr_U, nbr_D          ! Upstream and downstream neighbours (along z)
@@ -55,6 +62,8 @@ MODULE basic
               tc_step, tc_clos, tc_ghost, tc_coll, tc_process
   real(dp) :: maxruntime = 1e9 ! Maximum simulation CPU time
 
+  LOGICAL :: GATHERV_OUTPUT = .true.
+
   INTERFACE allocate_array
     MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4, allocate_array_dp5, allocate_array_dp6
     MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3,allocate_array_dc4, allocate_array_dc5, allocate_array_dc6
@@ -159,22 +168,22 @@ CONTAINS
       hours = (time/3600./24. - days) * 24
       mins  = (time/3600. - days*24. - hours) * 60
       secs  = (time/60. - days*24.*60 - hours*60 - mins) * 60
-      WRITE(*,*) 'CPU Time = ', days, '[day]', hours, '[h]', mins, '[min]', secs, '[s]'
-      WRITE(*,*) '(',time,'[s])'
+      IF (my_id .EQ. 0) WRITE(*,*) 'CPU Time = ', days, '[day]', hours, '[h]', mins, '[min]', secs, '[s]'
+      IF (my_id .EQ. 0) WRITE(*,*) '(',time,'[s])'
 
     ELSEIF ( hours .GT. 0 ) THEN !display h min s
       mins  = (time/3600. - hours) * 60
       secs  = (time/60. - hours*60 - mins) * 60
-      WRITE(*,*) 'CPU Time = ', hours, '[h]', mins, '[min]', secs, '[s]'
-      WRITE(*,*) '(',time,'[s])'
+      IF (my_id .EQ. 0) WRITE(*,*) 'CPU Time = ', hours, '[h]', mins, '[min]', secs, '[s]'
+      IF (my_id .EQ. 0) WRITE(*,*) '(',time,'[s])'
 
     ELSEIF ( mins .GT. 0 ) THEN !display min s
       secs  = (time/60. - mins) * 60
-      WRITE(*,*) 'CPU Time = ', mins, '[min]', secs, '[s]'
-      WRITE(*,*) '(',time,'[s])'
+      IF (my_id .EQ. 0) WRITE(*,*) 'CPU Time = ', mins, '[min]', secs, '[s]'
+      IF (my_id .EQ. 0) WRITE(*,*) '(',time,'[s])'
 
     ELSE ! display s
-      WRITE(*,*) 'CPU Time = ', FLOOR(time), '[s]'
+      IF (my_id .EQ. 0) WRITE(*,*) 'CPU Time = ', FLOOR(time), '[s]'
 
     ENDIF
   END SUBROUTINE display_h_min_s
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index 466df5318d10a40205c64d8842d26e2980ee1d9b..a724845d354b87fdcacf5a11f3628781427b24ee 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -3,7 +3,7 @@ MODULE calculus
   USE basic
   USE prec_const
   USE grid
-  USE utility
+  USE parallel, ONLY: manual_0D_bcast
   IMPLICIT NONE
   REAL(dp), dimension(-2:2) :: dz_usu = &
    (/  onetwelfth, -twothird, &
@@ -16,7 +16,7 @@ MODULE calculus
    REAL(dp), dimension(-2:2) :: dz2_usu = &
    (/-1.0_dp/12.0_dp, 4.0_dp/3.0_dp, -5.0_dp/2.0_dp, 4.0_dp/3.0_dp, -1.0_dp/12.0_dp /)! 2th derivative, 4th order (for parallel hypdiff)
    REAL(dp), dimension(-2:2) :: dz4_usu = &
-   (/ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp /) ! 4th derivative, 2nd order (for parallel hypdiff)
+   (/  1._dp, -4._dp, 6._dp, -4._dp, 1._dp /) ! 4th derivative, 2nd order (for parallel hypdiff)
   PUBLIC :: simpson_rule_z, interp_z, grad_z, grad_z4
 
 CONTAINS
@@ -99,7 +99,7 @@ SUBROUTINE grad_z4(f,ddz4f)
   IF(Nz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
       DO iz = izs,ize
        ddz4f(iz) = dz4_usu(-2)*f(iz-2) + dz4_usu(-1)*f(iz-1) &
-                  +dz4_usu( 0)*f(iz  )&
+                  +dz4_usu( 0)*f(iz)&
                   +dz4_usu( 1)*f(iz+1) + dz4_usu( 2)*f(iz+2)
       ENDDO
   ELSE
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index ea43b1d3fc804b528b8e74fe5b54cd657ca7bfc3..f2ca2956cfe02df87cf24e46ee46bcbda6f953ac 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -499,7 +499,7 @@ CONTAINS
                   ! Sum up all the sub collision terms on root 0
                   CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
                   ! distribute the sum over the process among p
-                  CALL MPI_SCATTERV(buffer_e, counts_np_e, displs_np_e, MPI_DOUBLE_COMPLEX,&
+                  CALL MPI_SCATTERV(buffer_e, rcv_p_e, dsp_p_e, MPI_DOUBLE_COMPLEX,&
                                     TColl_distr_e, local_np_e, MPI_DOUBLE_COMPLEX,&
                                     0, comm_p, ierr)
                 ELSE
@@ -522,7 +522,7 @@ CONTAINS
               CALL MPI_REDUCE(local_sum_i, buffer_i, total_np_i, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
               ! buffer contains the entire collision term along p, we scatter it between
               ! the other processes (use of scatterv since Pmax/Np is not an integer)
-              CALL MPI_SCATTERV(buffer_i, counts_np_i, displs_np_i, MPI_DOUBLE_COMPLEX,&
+              CALL MPI_SCATTERV(buffer_i, rcv_p_i, dsp_p_i, MPI_DOUBLE_COMPLEX,&
                                 TColl_distr_i, local_np_i, MPI_DOUBLE_COMPLEX, &
                                 0, comm_p, ierr)
             ELSE
diff --git a/src/diag_headers/diagnose_fields.h b/src/diag_headers/diagnose_fields.h
deleted file mode 100644
index a68cb8fbb596044e23eef10df7e648be9760fe26..0000000000000000000000000000000000000000
--- a/src/diag_headers/diagnose_fields.h
+++ /dev/null
@@ -1,157 +0,0 @@
-SUBROUTINE diagnose_fields(kstep)
-  USE basic
-  USE grid
-  USE diagnostics_par
-  USE futils
-  USE array
-  USE model
-  USE initial_par
-  USE fields
-  USE time_integration
-  USE utility
-  USE prec_const
-  USE collision, ONLY: coll_outputinputs
-  USE geometry
-  USE processing
-  IMPLICIT NONE
-  INTEGER, INTENT(in) :: kstep
-  INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank = 0
-  INTEGER :: dims(1) = (/0/)
-  !____________________________________________________________________________!
-  IF ((kstep .EQ. 0)) THEN
-   !  fields result file (electro. pot., Ni00 moment, fluid moments etc.)
-   IF (nsave_3d .GT. 0) THEN
-    CALL init_outfile(comm_kyz,fldfile0,fldfile,fidfld)
-    CALL creatd(fidfld, rank, dims,  "/data/time",     "Time t*c_s/R")
-    CALL creatd(fidfld, rank, dims, "/data/cstep", "iteration number")
-
-    IF (write_phi) CALL creatg(fidfld, "/data/phi", "phi")
-
-    IF (write_Na00) THEN
-     IF(KIN_E)&
-     CALL creatg(fidfld, "/data/Ne00", "Ne00")
-     CALL creatg(fidfld, "/data/Ni00", "Ni00")
-    ENDIF
-
-    IF (write_dens) THEN
-      IF(KIN_E)&
-     CALL creatg(fidfld, "/data/dens_e", "dens_e")
-     CALL creatg(fidfld, "/data/dens_i", "dens_i")
-    ENDIF
-
-    IF (write_fvel) THEN
-      IF(KIN_E) THEN
-      CALL creatg(fidfld, "/data/upar_e", "upar_e")
-      CALL creatg(fidfld, "/data/uper_e", "uper_e")
-      ENDIF
-      CALL creatg(fidfld, "/data/upar_i", "upar_i")
-      CALL creatg(fidfld, "/data/uper_i", "uper_i")
-    ENDIF
-
-    IF (write_temp) THEN
-      IF(KIN_E) THEN
-      CALL creatg(fidfld, "/data/Tper_e", "Tper_e")
-      CALL creatg(fidfld, "/data/Tpar_e", "Tpar_e")
-      CALL creatg(fidfld, "/data/temp_e", "temp_e")
-      ENDIF
-      CALL creatg(fidfld, "/data/Tper_i", "Tper_i")
-      CALL creatg(fidfld, "/data/Tpar_i", "Tpar_i")
-      CALL creatg(fidfld, "/data/temp_i", "temp_i")
-    ENDIF
-
-    IF (cstep==0) THEN
-      iframe3d=0
-    ENDIF
-    CALL attach(fidfld,"/data/" , "frames", iframe3d)
-   END IF
-  ENDIF
-
-  !_____________________________________________________________________________
-  !                   2.   Periodic diagnostics
-  !
-  IF (kstep .GE. 0) THEN
-
-  !                       2.2   3d profiles
-  IF (nsave_3d .GT. 0) THEN
-     IF (MOD(cstep, nsave_3d) == 0) THEN
-
-       CALL append(fidfld,  "/data/time",           time,ionode=0)
-       CALL append(fidfld, "/data/cstep", real(cstep,dp),ionode=0)
-       CALL getatt(fidfld,      "/data/",       "frames",iframe3d)
-       iframe3d=iframe3d+1
-       CALL attach(fidfld,"/data/" , "frames", iframe3d)
-
-       IF (write_phi) CALL write_field3d_kykxz(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
-
-       IF (write_Na00) THEN
-         IF(KIN_E)THEN
-         IF (ips_e .EQ. 1) THEN
-           Ne00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_e(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
-         ENDIF
-         CALL write_field3d_kykxz(Ne00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ne00')
-         ENDIF
-         IF (ips_i .EQ. 1) THEN
-           Ni00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_i(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
-         ENDIF
-         CALL write_field3d_kykxz(Ni00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ni00')
-       ENDIF
-
-       !! Fuid moments
-       IF (write_dens .OR. write_fvel .OR. write_temp) &
-       CALL compute_fluid_moments
-
-       IF (write_dens) THEN
-         IF(KIN_E)&
-         CALL write_field3d_kykxz(dens_e(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_e')
-         CALL write_field3d_kykxz(dens_i(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_i')
-       ENDIF
-
-       IF (write_fvel) THEN
-         IF(KIN_E)&
-         CALL write_field3d_kykxz(upar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_e')
-         CALL write_field3d_kykxz(upar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_i')
-         IF(KIN_E)&
-         CALL write_field3d_kykxz(uper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_e')
-         CALL write_field3d_kykxz(uper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_i')
-       ENDIF
-
-       IF (write_temp) THEN
-         IF(KIN_E)&
-         CALL write_field3d_kykxz(Tpar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_e')
-         CALL write_field3d_kykxz(Tpar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_i')
-         IF(KIN_E)&
-         CALL write_field3d_kykxz(Tper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_e')
-         CALL write_field3d_kykxz(Tper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_i')
-         IF(KIN_E)&
-         CALL write_field3d_kykxz(temp_e(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_e')
-         CALL write_field3d_kykxz(temp_i(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_i')
-       ENDIF
-
-     END IF
-  END IF
-  !_____________________________________________________________________________
-  !                   3.   Final diagnostics
-  ELSEIF (kstep .EQ. -1) THEN
-    !   Close diagnostic files
-    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-    CALL closef(fidfld)
-
-  END IF
-
-  CONTAINS
-    SUBROUTINE write_field3d_kykxz(field, text)
-      IMPLICIT NONE
-      COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: field
-      CHARACTER(*), INTENT(IN) :: text
-      CHARACTER(LEN=50) :: dset_name
-      WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data", TRIM(text), iframe3d
-      IF (num_procs .EQ. 1) THEN ! no data distribution
-        CALL putarr(fidfld, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), ionode=0)
-      ELSE
-        CALL putarrnd(fidfld, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize),  (/1, 3/))
-      ENDIF
-      CALL attach(fidfld, dset_name, "time", time)
-    END SUBROUTINE write_field3d_kykxz
-
-END SUBROUTINE diagnose_fields
diff --git a/src/diag_headers/diagnose_full.h b/src/diag_headers/diagnose_full.h
deleted file mode 100644
index 5f9cc79de522ecb1a82861832dec3cf877c6c19d..0000000000000000000000000000000000000000
--- a/src/diag_headers/diagnose_full.h
+++ /dev/null
@@ -1,478 +0,0 @@
-SUBROUTINE diagnose_full(kstep)
-USE basic
-USE grid
-USE diagnostics_par
-USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf, putarrnd
-USE array
-USE model
-USE initial_par
-USE fields
-USE time_integration
-USE utility
-USE prec_const
-USE collision, ONLY: coll_outputinputs
-USE geometry
-IMPLICIT NONE
-INTEGER, INTENT(in) :: kstep
-INTEGER, parameter  :: BUFSIZE = 2
-INTEGER :: rank = 0
-INTEGER :: dims(1) = (/0/)
-INTEGER :: cp_counter = 0
-CHARACTER(len=256) :: str,test_
-!____________________________________________________________________________
-!                   1.   Initial diagnostics
-
-IF ((kstep .EQ. 0)) THEN
-  CALL init_outfile(comm0,   resfile0,resfile,fidres)
-
-  ! Profiler time measurement
-  CALL creatg(fidres, "/profiler", "performance analysis")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_rhs",        "cumulative rhs computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_adv_field",  "cumulative adv. fields computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_clos",       "cumulative closure computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_ghost",       "cumulative communication time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_coll",       "cumulative collision computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_poisson",    "cumulative poisson computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_Sapj",       "cumulative Sapj computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_checkfield", "cumulative checkfield computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_diag",       "cumulative sym computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_process",    "cumulative process computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/Tc_step",       "cumulative total step computation time")
-  CALL creatd(fidres, 0, dims, "/profiler/time",          "current simulation time")
-
-  ! Grid info
-  CALL creatg(fidres, "/data/grid", "Grid data")
-  CALL putarr(fidres, "/data/grid/coordkx",   kxarray_full,  "kx*rho_s0", ionode=0)
-  CALL putarr(fidres, "/data/grid/coordky",   kyarray_full,  "ky*rho_s0", ionode=0)
-  CALL putarr(fidres, "/data/grid/coordz",    zarray_full,   "z/R", ionode=0)
-  CALL putarr(fidres, "/data/grid/coordp_e" , parray_e_full, "p_e", ionode=0)
-  CALL putarr(fidres, "/data/grid/coordj_e" , jarray_e_full, "j_e", ionode=0)
-  CALL putarr(fidres, "/data/grid/coordp_i" , parray_i_full, "p_i", ionode=0)
-  CALL putarr(fidres, "/data/grid/coordj_i" , jarray_i_full, "j_i", ionode=0)
-
-  ! Metric info
-  CALL   creatg(fidres, "/data/metric", "Metric data")
-  CALL putarrnd(fidres, "/data/metric/gxx",            gxx(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gxy",            gxy(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gyy",            gyy(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gyz",            gyz(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gzz",            gzz(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/hatR",          hatR(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/hatZ",          hatZ(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/hatB",          hatB(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gradxB",      gradxB(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gradyB",      gradyB(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gradzB",      gradzB(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/Jacobian",    Jacobian(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1, 1/))
-  CALL putarrnd(fidres, "/data/metric/Ckxky",       Ckxky(ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 3/))
-  CALL putarrnd(fidres, "/data/metric/kernel_i",    kernel_i(ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/ 1, 2, 4/))
-
-  !  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 (write_hf) THEN
-     CALL creatd(fidres, rank, dims, "/data/var0d/hflux_x", "Radial part ion heat flux")
-   ENDIF
-   IF (cstep==0) THEN
-     iframe0d=0
-   ENDIF
-   CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
-  END IF
-
-
-  !  var2d group (??)
-  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 (cstep==0) THEN
-     iframe2d=0
-   ENDIF
-   CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
-  END IF
-
-  !  var3d group (electro. pot., Ni00 moment)
-  IF (nsave_3d .GT. 0) THEN
-   CALL creatg(fidres, "/data/var3d", "3d profiles")
-   CALL creatd(fidres, rank, dims,  "/data/var3d/time",     "Time t*c_s/R")
-   CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
-
-   IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
-
-   IF (write_Na00) THEN
-    IF(KIN_E)&
-    CALL creatg(fidres, "/data/var3d/Ne00", "Ne00")
-    CALL creatg(fidres, "/data/var3d/Ni00", "Ni00")
-    IF(KIN_E)&
-    CALL creatg(fidres, "/data/var3d/Nepjz", "Nepjz")
-    CALL creatg(fidres, "/data/var3d/Nipjz", "Nipjz")
-   ENDIF
-
-   IF (write_dens) THEN
-     IF(KIN_E)&
-    CALL creatg(fidres, "/data/var3d/dens_e", "dens_e")
-    CALL creatg(fidres, "/data/var3d/dens_i", "dens_i")
-   ENDIF
-
-   IF (write_fvel) THEN
-     IF(KIN_E) THEN
-     CALL creatg(fidres, "/data/var3d/upar_e", "upar_e")
-     CALL creatg(fidres, "/data/var3d/uper_e", "uper_e")
-     ENDIF
-     CALL creatg(fidres, "/data/var3d/upar_i", "upar_i")
-     CALL creatg(fidres, "/data/var3d/uper_i", "uper_i")
-   ENDIF
-
-   IF (write_temp) THEN
-     IF(KIN_E) THEN
-     CALL creatg(fidres, "/data/var3d/Tper_e", "Tper_e")
-     CALL creatg(fidres, "/data/var3d/Tpar_e", "Tpar_e")
-     CALL creatg(fidres, "/data/var3d/temp_e", "temp_e")
-     ENDIF
-     CALL creatg(fidres, "/data/var3d/Tper_i", "Tper_i")
-     CALL creatg(fidres, "/data/var3d/Tpar_i", "Tpar_i")
-     CALL creatg(fidres, "/data/var3d/temp_i", "temp_i")
-   ENDIF
-
-   IF (cstep==0) THEN
-     iframe3d=0
-   ENDIF
-   CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
-  END IF
-
-  !  var5d group (moments)
-  IF (nsave_5d .GT. 0) THEN
-    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
-      IF(KIN_E)&
-     CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
-     CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
-    ENDIF
-
-    IF (write_Sapj) THEN
-      IF(KIN_E)&
-     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
-ENDIF
-
-!_____________________________________________________________________________
-!                   2.   Periodic diagnostics
-!
-IF (kstep .GE. 0) THEN
-
-   !                       2.1   0d history arrays
-   IF (nsave_0d .GT. 0) THEN
-      IF ( MOD(cstep, nsave_0d) == 0 ) THEN
-         CALL diagnose_0d
-      END IF
-   END IF
-
-   !                       2.2   1d profiles
-   ! empty in our case
-
-   !                       2.3   2d profiles
-   ! empty in our case
-
-
-   !                       2.3   2d profiles
-   IF (nsave_3d .GT. 0) THEN
-      IF (MOD(cstep, nsave_3d) == 0) THEN
-        CALL diagnose_3d
-      END IF
-   END IF
-
-   !                       2.4   3d profiles
-   IF (nsave_5d .GT. 0 .AND. cstep .GT. 0) THEN
-      IF (MOD(cstep, nsave_5d) == 0) THEN
-         CALL diagnose_5d
-      END IF
-   END IF
-
-!_____________________________________________________________________________
-!                   3.   Final diagnostics
-
-ELSEIF (kstep .EQ. -1) THEN
-   CALL attach(fidres, "/data/input","cpu_time",finish-start)
-
-   ! make a checkpoint at last timestep if not crashed
-   IF(.NOT. crashed) THEN
-     IF(my_id .EQ. 0) write(*,*) 'Saving last state'
-     IF (nsave_5d .GT. 0) &
-     CALL diagnose_5d
-   ENDIF
-
-   !   Close all diagnostic files
-   CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-   CALL closef(fidres)
-
-END IF
-
-END SUBROUTINE diagnose_full
-
-!!-------------- Auxiliary routines -----------------!!
-SUBROUTINE diagnose_0d
-
-  USE basic
-  USE futils, ONLY: append, attach, getatt
-  USE diagnostics_par
-  USE prec_const
-  USE processing
-  USE model, ONLY: KIN_E
-
-  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_clos",            tc_clos,ionode=0)
-  CALL append(fidres, "/profiler/Tc_ghost",          tc_ghost,ionode=0)
-  CALL append(fidres, "/profiler/Tc_coll",            tc_coll,ionode=0)
-  CALL append(fidres, "/profiler/Tc_poisson",      tc_poisson,ionode=0)
-  CALL append(fidres, "/profiler/Tc_Sapj",            tc_Sapj,ionode=0)
-  CALL append(fidres, "/profiler/Tc_checkfield",tc_checkfield,ionode=0)
-  CALL append(fidres, "/profiler/Tc_diag",            tc_diag,ionode=0)
-  CALL append(fidres, "/profiler/Tc_process",      tc_process,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
-  IF (write_hf) THEN
-    CALL compute_radial_heatflux
-    CALL append(fidres, "/data/var0d/hflux_x",hflux_x,ionode=0)
-  ENDIF
-END SUBROUTINE diagnose_0d
-
-SUBROUTINE diagnose_3d
-
-  USE basic
-  USE futils, ONLY: append, getatt, attach, putarrnd, putarr
-  USE fields
-  USE array
-  USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx, ikx, iky, ips_e, ips_i
-  USE time_integration
-  USE diagnostics_par
-  USE prec_const
-  USE processing
-  USE model, ONLY: KIN_E
-
-  IMPLICIT NONE
-
-  INTEGER     :: i_, root, world_rank, world_size
-
-  CALL append(fidres,  "/data/var3d/time",           time,ionode=0)
-  CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
-  CALL getatt(fidres,      "/data/var3d/",       "frames",iframe3d)
-  iframe3d=iframe3d+1
-  CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
-
-  IF (write_phi) CALL write_field3d_kykxz(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
-
-  IF (write_Na00) THEN
-    IF(KIN_E)THEN
-    IF (ips_e .EQ. 1) THEN
-      Ne00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_e(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
-    ENDIF
-    CALL write_field3d_kykxz(Ne00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ne00')
-    ENDIF
-    IF (ips_i .EQ. 1) THEN
-      Ni00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_i(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
-    ENDIF
-    CALL write_field3d_kykxz(Ni00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ni00')
-
-    CALL compute_Napjz_spectrum
-    IF(KIN_E) &
-    CALL write_field3d_pjz_e(Nepjz(ips_e:ipe_e,ijs_e:ije_e,izs:ize), 'Nepjz')
-    CALL write_field3d_pjz_i(Nipjz(ips_i:ipe_i,ijs_i:ije_i,izs:ize), 'Nipjz')
-  ENDIF
-
-  !! Fuid moments
-  IF (write_dens .OR. write_fvel .OR. write_temp) &
-  CALL compute_fluid_moments
-
-  IF (write_dens) THEN
-    IF(KIN_E)&
-    CALL write_field3d_kykxz(dens_e(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_e')
-    CALL write_field3d_kykxz(dens_i(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_i')
-  ENDIF
-
-  IF (write_fvel) THEN
-    IF(KIN_E)&
-    CALL write_field3d_kykxz(upar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_e')
-    CALL write_field3d_kykxz(upar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_i')
-    IF(KIN_E)&
-    CALL write_field3d_kykxz(uper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_e')
-    CALL write_field3d_kykxz(uper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_i')
-  ENDIF
-
-  IF (write_temp) THEN
-    IF(KIN_E)&
-    CALL write_field3d_kykxz(Tpar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_e')
-    CALL write_field3d_kykxz(Tpar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_i')
-    IF(KIN_E)&
-    CALL write_field3d_kykxz(Tper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_e')
-    CALL write_field3d_kykxz(Tper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_i')
-    IF(KIN_E)&
-    CALL write_field3d_kykxz(temp_e(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_e')
-    CALL write_field3d_kykxz(temp_i(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_i')
-  ENDIF
-
-  CONTAINS
-
-  SUBROUTINE write_field3d_kykxz(field, text)
-    IMPLICIT NONE
-    COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: field
-    CHARACTER(*), INTENT(IN) :: text
-    CHARACTER(LEN=50) :: dset_name
-    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
-    IF (num_procs .EQ. 1) THEN ! no data distribution
-      CALL putarr(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), ionode=0)
-    ELSE
-      CALL putarrnd(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize),  (/1, 1, 3/))
-    ENDIF
-    CALL attach(fidres, dset_name, "time", time)
-  END SUBROUTINE write_field3d_kykxz
-
-  SUBROUTINE write_field3d_pjz_i(field, text)
-    IMPLICIT NONE
-    REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
-    CHARACTER(*), INTENT(IN) :: text
-    CHARACTER(LEN=50) :: dset_name
-    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
-    IF (num_procs .EQ. 1) THEN ! no data distribution
-      CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize), ionode=0)
-    ELSE
-      CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize),  (/1, 0, 3/))
-    ENDIF
-    CALL attach(fidres, dset_name, "time", time)
-  END SUBROUTINE write_field3d_pjz_i
-
-  SUBROUTINE write_field3d_pjz_e(field, text)
-    IMPLICIT NONE
-    REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
-    CHARACTER(*), INTENT(IN) :: text
-    CHARACTER(LEN=50) :: dset_name
-    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
-    IF (num_procs .EQ. 1) THEN ! no data distribution
-      CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize), ionode=0)
-    ELSE
-      CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize),  (/1, 0, 3/))
-    ENDIF
-    CALL attach(fidres, dset_name, "time", time)
-  END SUBROUTINE write_field3d_pjz_e
-
-END SUBROUTINE diagnose_3d
-
-SUBROUTINE diagnose_5d
-
-  USE basic
-  USE futils, ONLY: append, getatt, attach, putarrnd
-  USE fields
-  USE array!, ONLY: Sepj, Sipj
-  USE grid, ONLY: ips_e,ipe_e, ips_i, ipe_i, &
-                 ijs_e,ije_e, ijs_i, ije_i, &
-                 ikxs,ikxe,ikys,ikye,izs,ize
-  USE time_integration
-  USE diagnostics_par
-  USE prec_const
-  USE model, ONLY: KIN_E
-  IMPLICIT NONE
-
-  CALL append(fidres,  "/data/var5d/time",           time,ionode=0)
-  CALL append(fidres, "/data/var5d/cstep", real(cstep,dp),ionode=0)
-  CALL getatt(fidres,      "/data/var5d/",       "frames",iframe5d)
-  iframe5d=iframe5d+1
-  CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
-
-  IF (write_Napj) THEN
-   IF(KIN_E)&
-  CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_e')
-  CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_i')
-  ENDIF
-
-  IF (write_Sapj) THEN
-   IF(KIN_E)&
-   CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), 'Sepj')
-   CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), 'Sipj')
-  ENDIF
-
-  CONTAINS
-
-  SUBROUTINE write_field5d_e(field, text)
-    USE futils, ONLY: attach, putarr, putarrnd
-    USE grid,   ONLY: ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize
-    USE prec_const
-    IMPLICIT NONE
-
-    COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
-    CHARACTER(*), INTENT(IN) :: text
-
-    CHARACTER(LEN=50) :: dset_name
-
-    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
-    IF (num_procs .EQ. 1) THEN
-     CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
-    ELSE
-     CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize),  (/1,3,5/))
-    ENDIF
-    CALL attach(fidres, dset_name, 'cstep', cstep)
-    CALL attach(fidres, dset_name, 'time', time)
-    CALL attach(fidres, dset_name, 'jobnum', jobnum)
-    CALL attach(fidres, dset_name, 'dt', dt)
-    CALL attach(fidres, dset_name, 'iframe2d', iframe2d)
-    CALL attach(fidres, dset_name, 'iframe5d', iframe5d)
-
-  END SUBROUTINE write_field5d_e
-
-  SUBROUTINE write_field5d_i(field, text)
-    USE futils, ONLY: attach, putarr, putarrnd
-    USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize
-    USE prec_const
-    IMPLICIT NONE
-
-    COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
-    CHARACTER(*), INTENT(IN) :: text
-
-    CHARACTER(LEN=50) :: dset_name
-
-    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
-    IF (num_procs .EQ. 1) THEN
-      CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
-    ELSE
-      CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize),  (/1,3,5/))
-    ENDIF
-    CALL attach(fidres, dset_name, 'cstep', cstep)
-    CALL attach(fidres, dset_name, 'time', time)
-    CALL attach(fidres, dset_name, 'jobnum', jobnum)
-    CALL attach(fidres, dset_name, 'dt', dt)
-    CALL attach(fidres, dset_name, 'iframe2d', iframe2d)
-    CALL attach(fidres, dset_name, 'iframe5d', iframe5d)
-
-  END SUBROUTINE write_field5d_i
-END SUBROUTINE diagnose_5d
diff --git a/src/diag_headers/diagnose_gridgeom.h b/src/diag_headers/diagnose_gridgeom.h
deleted file mode 100644
index 9ef1855c928957eeb795102eb70a28f7269b0c07..0000000000000000000000000000000000000000
--- a/src/diag_headers/diagnose_gridgeom.h
+++ /dev/null
@@ -1,50 +0,0 @@
-SUBROUTINE diagnose_gridgeom(kstep)
-  USE basic
-  USE grid
-  USE diagnostics_par
-  USE futils, ONLY: putarr, creatg, putarrnd, closef
-  USE time_integration
-  USE prec_const
-  USE geometry
-  USE array
-  IMPLICIT NONE
-  INTEGER, INTENT(in) :: kstep
-  INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank = 0
-  INTEGER :: dims(1) = (/0/)
-  IF (kstep .EQ. 0) THEN
-    ! Grid info
-    CALL init_outfile(comm0,ggmfile0,ggmfile,fidggm)
-    CALL creatg(fidggm, "/data/grid", "Grid data")
-    CALL putarr(fidggm, "/data/grid/coordkx",   kxarray_full,  "kx*rho_s0", ionode=0)
-    CALL putarr(fidggm, "/data/grid/coordky",   kyarray_full,  "ky*rho_s0", ionode=0)
-    CALL putarr(fidggm, "/data/grid/coordz",    zarray_full,   "z/R", ionode=0)
-    CALL putarr(fidggm, "/data/grid/coordp_e" , parray_e_full, "p_e", ionode=0)
-    CALL putarr(fidggm, "/data/grid/coordj_e" , jarray_e_full, "j_e", ionode=0)
-    CALL putarr(fidggm, "/data/grid/coordp_i" , parray_i_full, "p_i", ionode=0)
-    CALL putarr(fidggm, "/data/grid/coordj_i" , jarray_i_full, "j_i", ionode=0)
-
-    ! Metric info
-    CALL   creatg(fidggm, "/data/metric", "Metric data")
-    CALL putarrnd(fidggm, "/data/metric/gxx",            gxx(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gxy",            gxy(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gyy",            gyy(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gyz",            gyz(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gzz",            gzz(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/hatR",          hatR(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/hatZ",          hatZ(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/hatB",          hatB(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gradxB",      gradxB(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gradyB",      gradyB(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gradzB",      gradzB(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/Jacobian",    Jacobian(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1, 1/))
-    CALL putarrnd(fidggm, "/data/metric/Ckxky",       Ckxky(ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 3/))
-    CALL putarrnd(fidggm, "/data/metric/kernel_i",    kernel_i(ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/ 1, 2, 4/))
-  ENDIF
-  IF (kstep .EQ. -1) THEN
-    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-    CALL closef(fidggm)
-  ENDIF
-
-END SUBROUTINE diagnose_gridgeom
diff --git a/src/diag_headers/diagnose_moments.h b/src/diag_headers/diagnose_moments.h
deleted file mode 100644
index 585822602281cbdeff42e29dbc0d3e91c958c795..0000000000000000000000000000000000000000
--- a/src/diag_headers/diagnose_moments.h
+++ /dev/null
@@ -1,136 +0,0 @@
-SUBROUTINE diagnose_moments(kstep)
-  USE basic
-  USE grid
-  USE diagnostics_par
-  USE futils,         ONLY: creatg,creatd,append,putarr,putarrnd,attach,closef,getatt
-  USE model,          ONLY: KIN_E
-  USE array,          ONLY: Sepj,Sipj
-  USE fields,         ONLY: moments_e, moments_i
-  USE time_integration
-  USE utility
-  USE prec_const
-  IMPLICIT NONE
-  INTEGER, INTENT(in) :: kstep
-  INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank = 0
-  INTEGER :: dims(1) = (/0/)
-  !____________________________________________________________________________!
-  IF ((kstep .EQ. 0)) THEN
-
-  !  var5d group (moments)
-  IF (nsave_5d .GT. 0) THEN
-    CALL init_outfile(comm0,   momfile0,momfile,fidmom)
-
-    CALL creatd(fidmom, rank, dims,  "/data/time",     "Time t*c_s/R")
-    CALL creatd(fidmom, rank, dims, "/data/cstep", "iteration number")
-
-    IF (write_Napj) THEN
-      IF(KIN_E)&
-     CALL creatg(fidmom, "/data/moments_e", "moments_e")
-     CALL creatg(fidmom, "/data/moments_i", "moments_i")
-    ENDIF
-
-    IF (write_Sapj) THEN
-      IF(KIN_E)&
-     CALL creatg(fidmom, "/data/moments_e", "Sipj")
-     CALL creatg(fidmom, "/data/moments_i", "Sepj")
-    ENDIF
-
-    IF (cstep==0) THEN
-     iframe5d=0
-    END IF
-    CALL attach(fidmom,"/data/" , "frames", iframe5d)
-  END IF
-
-  ENDIF
-
-  !_____________________________________________________________________________
-  !                   2.   Periodic diagnostics
-  !
-  IF ((kstep .GE. 0) .OR. ((kstep .EQ. -1) .AND. (.NOT. CRASHED))) THEN
-
-  IF((kstep .EQ. -1) .AND. (.NOT. CRASHED .AND. (my_id .EQ. 0))) &
-    write(*,*) 'Saving last state'
-
-  !                       2.3   5d profiles
-    IF (MOD(cstep, nsave_5d) == 0) THEN
-     CALL append(fidmom,  "/data/time",           time,ionode=0)
-     CALL append(fidmom, "/data/cstep", real(cstep,dp),ionode=0)
-     CALL getatt(fidmom,      "/data/",       "frames",iframe5d)
-     iframe5d=iframe5d+1
-     CALL attach(fidmom,"/data/" , "frames", iframe5d)
-
-     IF (write_Napj) THEN
-      IF(KIN_E)&
-     CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_e')
-     CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_i')
-     ENDIF
-
-     IF (write_Sapj) THEN
-      IF(KIN_E)&
-      CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), 'Sepj')
-      CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), 'Sipj')
-     ENDIF
-    END IF
-  END IF
-  !_____________________________________________________________________________
-  !                   3.   Final diagnostics
-  IF (kstep .EQ. -1) THEN
-    !   Close diagnostic files
-    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-    CALL closef(fidmom)
-  END IF
-
-  CONTAINS
-    SUBROUTINE write_field5d_e(field, text)
-      USE futils, ONLY: attach, putarr, putarrnd
-      USE grid,   ONLY: ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize
-      USE prec_const
-      IMPLICIT NONE
-
-      COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
-      CHARACTER(*), INTENT(IN) :: text
-
-      CHARACTER(LEN=50) :: dset_name
-
-      WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data", TRIM(text), iframe5d
-      IF (num_procs .EQ. 1) THEN
-       CALL putarr(fidmom, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
-      ELSE
-       CALL putarrnd(fidmom, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize),  (/1,3,5/))
-      ENDIF
-      CALL attach(fidmom, dset_name, 'cstep', cstep)
-      CALL attach(fidmom, dset_name, 'time', time)
-      CALL attach(fidmom, dset_name, 'jobnum', jobnum)
-      CALL attach(fidmom, dset_name, 'dt', dt)
-      CALL attach(fidmom, dset_name, 'iframe2d', iframe2d)
-      CALL attach(fidmom, dset_name, 'iframe5d', iframe5d)
-
-    END SUBROUTINE write_field5d_e
-
-    SUBROUTINE write_field5d_i(field, text)
-      USE futils, ONLY: attach, putarr, putarrnd
-      USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize
-      USE prec_const
-      IMPLICIT NONE
-
-      COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
-      CHARACTER(*), INTENT(IN) :: text
-
-      CHARACTER(LEN=50) :: dset_name
-
-      WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data", TRIM(text), iframe5d
-      IF (num_procs .EQ. 1) THEN
-        CALL putarr(fidmom, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
-      ELSE
-        CALL putarrnd(fidmom, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize),  (/1,3,5/))
-      ENDIF
-      CALL attach(fidmom, dset_name, 'cstep', cstep)
-      CALL attach(fidmom, dset_name, 'time', time)
-      CALL attach(fidmom, dset_name, 'jobnum', jobnum)
-      CALL attach(fidmom, dset_name, 'dt', dt)
-      CALL attach(fidmom, dset_name, 'iframe2d', iframe2d)
-      CALL attach(fidmom, dset_name, 'iframe5d', iframe5d)
-
-    END SUBROUTINE write_field5d_i
-END SUBROUTINE diagnose_moments
diff --git a/src/diag_headers/diagnose_momspectrum.h b/src/diag_headers/diagnose_momspectrum.h
deleted file mode 100644
index 36befdda16321d1172fb86b89760da2e3a2e36ee..0000000000000000000000000000000000000000
--- a/src/diag_headers/diagnose_momspectrum.h
+++ /dev/null
@@ -1,105 +0,0 @@
-SUBROUTINE diagnose_momspectrum(kstep)
-  USE basic
-  USE grid
-  USE diagnostics_par
-  USE futils
-  USE array
-  USE model
-  USE initial_par
-  USE fields
-  USE time_integration
-  USE utility
-  USE prec_const
-  USE collision, ONLY: coll_outputinputs
-  USE geometry
-  USE processing
-  IMPLICIT NONE
-  INTEGER, INTENT(in) :: kstep
-  INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank = 0, frame
-  INTEGER :: dims(1) = (/0/)
-  !____________________________________________________________________________!
-  IF ((kstep .EQ. 0)) THEN
-    !  var3d group (pjz, spectrum of moments)
-    IF (nsave_3d .GT. 0) THEN
-     CALL init_outfile(comm_pz, mspfile0,mspfile,fidmsp)
-     CALL creatd(fidmsp, rank, dims,  "/data/time",     "Time t*c_s/R")
-     CALL creatd(fidmsp, rank, dims, "/data/cstep", "iteration number")
-
-     IF (write_Na00) THEN
-      IF(KIN_E)&
-      CALL creatg(fidmsp, "/data/Nepjz", "Nepjz")
-      CALL creatg(fidmsp, "/data/Nipjz", "Nipjz")
-     ENDIF
-
-     IF (cstep==0) THEN
-       iframe3d=0
-     ENDIF
-     CALL attach(fidmsp,"/data/" , "frames", iframe3d)
-    END IF
-
-  ENDIF
-
-  !_____________________________________________________________________________
-  !                   2.   Periodic diagnostics
-  !
-  IF (kstep .GE. 0) THEN
-  !                       2.2   3d profiles
-    IF (nsave_3d .GT. 0) THEN
-       IF (MOD(cstep, nsave_3d) == 0) THEN
-         CALL append(fidmsp,  "/data/time",           time,ionode=0)
-         CALL append(fidmsp, "/data/cstep", real(cstep,dp),ionode=0)
-         CALL getatt(fidmsp,      "/data/",       "frames",frame)
-         frame=frame+1
-         CALL attach(fidmsp,"/data/" , "frames", frame)
-
-         IF (write_Na00) THEN
-           CALL compute_Napjz_spectrum
-           IF(KIN_E) &
-           CALL write_field3d_pjz_e(Nepjz(ips_e:ipe_e,ijs_e:ije_e,izs:ize), 'Nepjz', frame)
-           CALL write_field3d_pjz_i(Nipjz(ips_i:ipe_i,ijs_i:ije_i,izs:ize), 'Nipjz', frame)
-         ENDIF
-       END IF
-    END IF
-
-  !_____________________________________________________________________________
-  !                   3.   Final diagnostics
-  ELSEIF (kstep .EQ. -1) THEN
-      !   Close diagnostic files
-      CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-      CALL closef(fidmsp)
-
-  END IF
-
-    CONTAINS
-      SUBROUTINE write_field3d_pjz_i(field, text, frame)
-        IMPLICIT NONE
-        REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
-        CHARACTER(*), INTENT(IN) :: text
-        INTEGER, INTENT(IN) :: frame
-        CHARACTER(LEN=50) :: dset_name
-        WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data", TRIM(text), frame
-        IF (num_procs .EQ. 1) THEN ! no data distribution
-          CALL putarr(fidmsp, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize), ionode=0)
-        ELSE
-          CALL putarrnd(fidmsp, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize),  (/1, 3/))
-        ENDIF
-        CALL attach(fidmsp, dset_name, "time", time)
-      END SUBROUTINE write_field3d_pjz_i
-
-      SUBROUTINE write_field3d_pjz_e(field, text, frame)
-        IMPLICIT NONE
-        REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
-        CHARACTER(*), INTENT(IN) :: text
-        INTEGER, INTENT(IN) :: frame
-        CHARACTER(LEN=50) :: dset_name
-        WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data", TRIM(text), frame
-        IF (num_procs .EQ. 1) THEN ! no data distribution
-          CALL putarr  (fidmsp, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize), ionode=0)
-        ELSE
-          CALL putarrnd(fidmsp, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize),  (/1, 3/))
-        ENDIF
-        CALL attach(fidmsp, dset_name, "time", time)
-      END SUBROUTINE write_field3d_pjz_e
-
-END SUBROUTINE diagnose_momspectrum
diff --git a/src/diag_headers/diagnose_profiler.h b/src/diag_headers/diagnose_profiler.h
deleted file mode 100644
index b99f67dda7826fc97b4e3a49ecea405a8b6baca4..0000000000000000000000000000000000000000
--- a/src/diag_headers/diagnose_profiler.h
+++ /dev/null
@@ -1,54 +0,0 @@
-SUBROUTINE diagnose_profiler(kstep)
-  USE basic
-  USE diagnostics_par
-  USE futils, ONLY: creatd,creatg,append,closef
-  USE time_integration
-  USE utility
-  USE prec_const
-  IMPLICIT NONE
-  INTEGER, INTENT(in) :: kstep
-  INTEGER :: dims(1) = (/0/)
-  !____________________________________________________________________________!
-  IF ((kstep .EQ. 0)) THEN
-    CALL init_outfile(comm0,   prffile0,prffile,fidprf)
-
-    ! data time measurement
-    CALL creatd(fidprf, 0, dims, "/data/Tc_rhs",        "cumulative rhs computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_adv_field",  "cumulative adv. fields computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_clos",       "cumulative closure computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_ghost",       "cumulative communication time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_coll",       "cumulative collision computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_poisson",    "cumulative poisson computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_Sapj",       "cumulative Sapj computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_checkfield", "cumulative checkfield computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_diag",       "cumulative sym computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_process",    "cumulative process computation time")
-    CALL creatd(fidprf, 0, dims, "/data/Tc_step",       "cumulative total step computation time")
-    CALL creatd(fidprf, 0, dims, "/data/time",          "current simulation time")
-  ENDIF
-  !_____________________________________________________________________________
-  !                   2.   Periodic diagnostics
-  !
-  IF (kstep .GE. 0) THEN
-  ! Time measurement data
-  CALL append(fidprf, "/data/Tc_rhs",              tc_rhs,ionode=0)
-  CALL append(fidprf, "/data/Tc_adv_field",  tc_adv_field,ionode=0)
-  CALL append(fidprf, "/data/Tc_clos",            tc_clos,ionode=0)
-  CALL append(fidprf, "/data/Tc_ghost",          tc_ghost,ionode=0)
-  CALL append(fidprf, "/data/Tc_coll",            tc_coll,ionode=0)
-  CALL append(fidprf, "/data/Tc_poisson",      tc_poisson,ionode=0)
-  CALL append(fidprf, "/data/Tc_Sapj",            tc_Sapj,ionode=0)
-  CALL append(fidprf, "/data/Tc_checkfield",tc_checkfield,ionode=0)
-  CALL append(fidprf, "/data/Tc_diag",            tc_diag,ionode=0)
-  CALL append(fidprf, "/data/Tc_process",      tc_process,ionode=0)
-  CALL append(fidprf, "/data/Tc_step",            tc_step,ionode=0)
-  CALL append(fidprf, "/data/time",                  time,ionode=0)
-  !_____________________________________________________________________________
-  !                   3.   Final diagnostics
-  ELSEIF (kstep .EQ. -1) THEN
-    !   Close diagnostic files
-    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-    CALL closef(fidprf)
-
-  END IF
-END SUBROUTINE diagnose_profiler
diff --git a/src/diag_headers/diagnose_timetraces.h b/src/diag_headers/diagnose_timetraces.h
deleted file mode 100644
index 3c688aee8c1230f1fc3292f4a4e1111dee0de34a..0000000000000000000000000000000000000000
--- a/src/diag_headers/diagnose_timetraces.h
+++ /dev/null
@@ -1,76 +0,0 @@
-SUBROUTINE diagnose_timetraces(kstep)
-  USE basic
-  USE grid
-  USE diagnostics_par
-  USE futils
-  USE array
-  USE model
-  USE initial_par
-  USE fields
-  USE time_integration
-  USE utility
-  USE prec_const
-  USE collision, ONLY: coll_outputinputs
-  USE geometry
-  USE processing
-  IMPLICIT NONE
-  INTEGER, INTENT(in) :: kstep
-  INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank = 0
-  INTEGER :: dims(1) = (/0/)
-  !____________________________________________________________________________!
-  IF ((kstep .EQ. 0)) THEN
-    CALL init_outfile(comm0,   ttrfile0,ttrfile,fidttr)
-
-    !  var0d group (gyro transport)
-    CALL creatd(fidttr, rank, dims,  "/data/time",     "Time t*c_s/R")
-    CALL creatd(fidttr, rank, dims, "/data/cstep", "iteration number")
-
-    IF (write_gamma) THEN
-     CALL creatd(fidttr, rank, dims, "/data/gflux_ri", "Radial gyro ion transport")
-     CALL creatd(fidttr, rank, dims, "/data/pflux_ri", "Radial part ion transport")
-    ENDIF
-    IF (write_hf) THEN
-     CALL creatd(fidttr, rank, dims, "/data/hflux_x", "Radial part ion heat flux")
-    ENDIF
-    IF (cstep==0) THEN
-     iframe0d=0
-    ENDIF
-    CALL attach(fidttr,"/data/" , "frames", iframe0d)
-  ENDIF
-
-  !_____________________________________________________________________________
-  !                   2.   Periodic diagnostics
-  !
-  IF (kstep .GE. 0) THEN
-
-  !                       0d time traces arrays
-  IF ( MOD(cstep, nsave_0d) == 0 ) THEN
-   ! Processing data
-   CALL append(fidttr,  "/data/time",           time,ionode=0)
-   CALL append(fidttr, "/data/cstep", real(cstep,dp),ionode=0)
-   CALL getatt(fidttr,      "/data/",       "frames",iframe0d)
-   iframe0d=iframe0d+1
-   CALL attach(fidttr,"/data/" , "frames", iframe0d)
-   ! Ion transport data
-   IF (write_gamma) THEN
-     CALL compute_radial_ion_transport
-     CALL append(fidttr, "/data/gflux_ri",gflux_ri,ionode=0)
-     CALL append(fidttr, "/data/pflux_ri",pflux_ri,ionode=0)
-   ENDIF
-   IF (write_hf) THEN
-     CALL compute_radial_heatflux
-     CALL append(fidttr, "/data/hflux_x",hflux_x,ionode=0)
-   ENDIF
-  END IF
-  !_____________________________________________________________________________
-  !                   3.   Final diagnostics
-  ELSEIF (kstep .EQ. -1) THEN
-    !   Close diagnostic files
-    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-    CALL closef(fidttr)
-
-  END IF
-END SUBROUTINE diagnose_timetraces
-!____________________________________________________________________________!
-!                       AUXILIARY ROUTINES
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index c1b81f2558712578a2e368633aab61ab1726c6be..1451b5e95c6fc454e9b604deb7d74203d4ef0080 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -6,7 +6,6 @@ SUBROUTINE diagnose(kstep)
   IMPLICIT NONE
 
   INTEGER, INTENT(in) :: kstep
-
   CALL cpu_time(t0_diag) ! Measuring time
 
   !! Basic diagnose loop for reading input file, displaying advancement and ending
@@ -26,12 +25,6 @@ SUBROUTINE diagnose(kstep)
   END IF
   !! Specific diagnostic calls
   CALL diagnose_full(kstep)
-  ! IF(nsave_5d .GT. 0) CALL diagnose_moments(kstep)
-  ! IF(nsave_3d .GT. 0) CALL diagnose_momspectrum(kstep)
-  ! IF(nsave_3d .GT. 0) CALL diagnose_fields(kstep)
-  ! IF(nsave_0d .GT. 0) CALL diagnose_profiler(kstep)
-  ! IF(nsave_0d .GT. 0) CALL diagnose_gridgeom(kstep)
-  ! IF(nsave_0d .GT. 0) CALL diagnose_timetraces(kstep)
 
   CALL cpu_time(t1_diag); tc_diag = tc_diag + (t1_diag - t0_diag)
 
@@ -73,12 +66,6 @@ SUBROUTINE init_outfile(comm,file0,file,fid)
   CALL creatg(fid, "/files", "files")
   CALL attach(fid, "/files",  "jobnum", jobnum)
 
-  !  Add input namelist variables as attributes of /data/input, defined in srcinfo.h
-  ! IF (my_id .EQ. 0) WRITE(*,*) 'VERSION=', VERSION
-  ! IF (my_id .EQ. 0) WRITE(*,*)  'BRANCH=', BRANCH
-  ! IF (my_id .EQ. 0) WRITE(*,*)  'AUTHOR=', AUTHOR
-  ! IF (my_id .EQ. 0) WRITE(*,*)    'HOST=', HOST
-
   ! Add the code info and parameters to the file
   WRITE(str,'(a,i2.2)') "/data/input"
   CALL creatd(fid, 0,(/0/),TRIM(str),'Input parameters')
@@ -106,11 +93,501 @@ SUBROUTINE init_outfile(comm,file0,file,fid)
   CALL putfile(fid, TRIM(str), TRIM(input_fname),ionode=0)
 END SUBROUTINE init_outfile
 
-!! Auxiliary routines hidden in headers
-INCLUDE 'diag_headers/diagnose_full.h'
-! INCLUDE 'diag_headers/diagnose_moments.h'
-! INCLUDE 'diag_headers/diagnose_momspectrum.h'
-! INCLUDE 'diag_headers/diagnose_fields.h'
-! INCLUDE 'diag_headers/diagnose_profiler.h'
-! INCLUDE 'diag_headers/diagnose_gridgeom.h'
-! INCLUDE 'diag_headers/diagnose_timetraces.h'
+SUBROUTINE diagnose_full(kstep)
+  USE basic
+  USE grid
+  USE diagnostics_par
+  USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf, putarrnd
+  USE array
+  USE model
+  USE initial_par
+  USE fields
+  USE time_integration
+  USE parallel
+  USE prec_const
+  USE collision, ONLY: coll_outputinputs
+  USE geometry
+  IMPLICIT NONE
+
+  INTEGER, INTENT(in) :: kstep
+  INTEGER, parameter  :: BUFSIZE = 2
+  INTEGER :: rank = 0
+  INTEGER :: dims(1) = (/0/)
+  INTEGER :: cp_counter = 0
+  CHARACTER(len=256) :: str,test_
+  !____________________________________________________________________________
+  !                   1.   Initial diagnostics
+
+  IF ((kstep .EQ. 0)) THEN
+    CALL init_outfile(comm0,   resfile0,resfile,fidres)
+
+    ! Profiler time measurement
+    CALL creatg(fidres, "/profiler", "performance analysis")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_rhs",        "cumulative rhs computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_adv_field",  "cumulative adv. fields computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_clos",       "cumulative closure computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_ghost",       "cumulative communication time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_coll",       "cumulative collision computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_poisson",    "cumulative poisson computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_Sapj",       "cumulative Sapj computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_checkfield", "cumulative checkfield computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_diag",       "cumulative sym computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_process",    "cumulative process computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/Tc_step",       "cumulative total step computation time")
+    CALL creatd(fidres, 0, dims, "/profiler/time",          "current simulation time")
+
+    ! Grid info
+    CALL creatg(fidres, "/data/grid", "Grid data")
+    CALL putarr(fidres, "/data/grid/coordkx",   kxarray_full,  "kx*rho_s0", ionode=0)
+    CALL putarr(fidres, "/data/grid/coordky",   kyarray_full,  "ky*rho_s0", ionode=0)
+    CALL putarr(fidres, "/data/grid/coordz",    zarray_full,   "z/R", ionode=0)
+    CALL putarr(fidres, "/data/grid/coordp_e" , parray_e_full, "p_e", ionode=0)
+    CALL putarr(fidres, "/data/grid/coordj_e" , jarray_e_full, "j_e", ionode=0)
+    CALL putarr(fidres, "/data/grid/coordp_i" , parray_i_full, "p_i", ionode=0)
+    CALL putarr(fidres, "/data/grid/coordj_i" , jarray_i_full, "j_i", ionode=0)
+
+    ! Metric info
+    CALL   creatg(fidres, "/data/metric", "Metric data")
+    CALL putarrnd(fidres, "/data/metric/gxx",            gxx(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gxy",            gxy(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gyy",            gyy(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gyz",            gyz(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gzz",            gzz(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/hatR",          hatR(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/hatZ",          hatZ(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/hatB",          hatB(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gradxB",      gradxB(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gradyB",      gradyB(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gradzB",      gradzB(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/Jacobian",    Jacobian(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1, 1/))
+    CALL putarrnd(fidres, "/data/metric/Ckxky",       Ckxky(ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 3/))
+    CALL putarrnd(fidres, "/data/metric/kernel_i",    kernel_i(ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/ 1, 2, 4/))
+
+    !  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 (write_hf) THEN
+       CALL creatd(fidres, rank, dims, "/data/var0d/hflux_x", "Radial part ion heat flux")
+     ENDIF
+     IF (cstep==0) THEN
+       iframe0d=0
+     ENDIF
+     CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
+    END IF
+
+
+    !  var2d group (??)
+    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 (cstep==0) THEN
+       iframe2d=0
+     ENDIF
+     CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
+    END IF
+
+    !  var3d group (electro. pot., Ni00 moment)
+    IF (nsave_3d .GT. 0) THEN
+     CALL creatg(fidres, "/data/var3d", "3d profiles")
+     CALL creatd(fidres, rank, dims,  "/data/var3d/time",     "Time t*c_s/R")
+     CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
+
+     IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
+
+     IF (write_Na00) THEN
+      IF(KIN_E)&
+      CALL creatg(fidres, "/data/var3d/Ne00", "Ne00")
+      CALL creatg(fidres, "/data/var3d/Ni00", "Ni00")
+      IF(KIN_E)&
+      CALL creatg(fidres, "/data/var3d/Nepjz", "Nepjz")
+      CALL creatg(fidres, "/data/var3d/Nipjz", "Nipjz")
+     ENDIF
+
+     IF (write_dens) THEN
+       IF(KIN_E)&
+      CALL creatg(fidres, "/data/var3d/dens_e", "dens_e")
+      CALL creatg(fidres, "/data/var3d/dens_i", "dens_i")
+     ENDIF
+
+     IF (write_fvel) THEN
+       IF(KIN_E) THEN
+       CALL creatg(fidres, "/data/var3d/upar_e", "upar_e")
+       CALL creatg(fidres, "/data/var3d/uper_e", "uper_e")
+       ENDIF
+       CALL creatg(fidres, "/data/var3d/upar_i", "upar_i")
+       CALL creatg(fidres, "/data/var3d/uper_i", "uper_i")
+     ENDIF
+
+     IF (write_temp) THEN
+       IF(KIN_E) THEN
+       CALL creatg(fidres, "/data/var3d/Tper_e", "Tper_e")
+       CALL creatg(fidres, "/data/var3d/Tpar_e", "Tpar_e")
+       CALL creatg(fidres, "/data/var3d/temp_e", "temp_e")
+       ENDIF
+       CALL creatg(fidres, "/data/var3d/Tper_i", "Tper_i")
+       CALL creatg(fidres, "/data/var3d/Tpar_i", "Tpar_i")
+       CALL creatg(fidres, "/data/var3d/temp_i", "temp_i")
+     ENDIF
+
+     IF (cstep==0) THEN
+       iframe3d=0
+     ENDIF
+     CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
+    END IF
+
+    !  var5d group (moments)
+    IF (nsave_5d .GT. 0) THEN
+      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
+        IF(KIN_E)&
+       CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
+       CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
+      ENDIF
+
+      IF (write_Sapj) THEN
+        IF(KIN_E)&
+        CALL creatg(fidres, "/data/var5d/Sepj", "Sepj")
+        CALL creatg(fidres, "/data/var5d/Sipj", "Sipj")
+      ENDIF
+
+      IF (cstep==0) THEN
+       iframe5d=0
+      END IF
+      CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
+    END IF
+  ENDIF
+
+  !_____________________________________________________________________________
+  !                   2.   Periodic diagnostics
+  !
+  IF (kstep .GE. 0) THEN
+
+     !                       2.1   0d history arrays
+     IF (nsave_0d .GT. 0) THEN
+        IF ( MOD(cstep, nsave_0d) == 0 ) THEN
+           CALL diagnose_0d
+        END IF
+     END IF
+
+     !                       2.2   1d profiles
+     ! empty in our case
+
+     !                       2.3   2d profiles
+     ! empty in our case
+
+
+     !                       2.3   2d profiles
+     IF (nsave_3d .GT. 0) THEN
+        IF (MOD(cstep, nsave_3d) == 0) THEN
+          CALL diagnose_3d
+        END IF
+     END IF
+
+     !                       2.4   3d profiles
+     IF (nsave_5d .GT. 0 .AND. cstep .GT. 0) THEN
+        IF (MOD(cstep, nsave_5d) == 0) THEN
+           CALL diagnose_5d
+        END IF
+     END IF
+
+  !_____________________________________________________________________________
+  !                   3.   Final diagnostics
+
+  ELSEIF (kstep .EQ. -1) THEN
+     CALL attach(fidres, "/data/input","cpu_time",finish-start)
+
+     ! make a checkpoint at last timestep if not crashed
+     IF(.NOT. crashed) THEN
+       IF(my_id .EQ. 0) write(*,*) 'Saving last state'
+       IF (nsave_5d .GT. 0) &
+       CALL diagnose_5d
+     ENDIF
+
+     !   Close all diagnostic files
+     CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+     CALL closef(fidres)
+
+  END IF
+
+END SUBROUTINE diagnose_full
+
+!!-------------- Auxiliary routines -----------------!!
+SUBROUTINE diagnose_0d
+
+  USE basic
+  USE futils, ONLY: append, attach, getatt
+  USE diagnostics_par
+  USE prec_const
+  USE processing
+  USE model, ONLY: KIN_E
+
+  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_clos",            tc_clos,ionode=0)
+  CALL append(fidres, "/profiler/Tc_ghost",          tc_ghost,ionode=0)
+  CALL append(fidres, "/profiler/Tc_coll",            tc_coll,ionode=0)
+  CALL append(fidres, "/profiler/Tc_poisson",      tc_poisson,ionode=0)
+  CALL append(fidres, "/profiler/Tc_Sapj",            tc_Sapj,ionode=0)
+  CALL append(fidres, "/profiler/Tc_checkfield",tc_checkfield,ionode=0)
+  CALL append(fidres, "/profiler/Tc_diag",            tc_diag,ionode=0)
+  CALL append(fidres, "/profiler/Tc_process",      tc_process,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
+  IF (write_hf) THEN
+    CALL compute_radial_heatflux
+    CALL append(fidres, "/data/var0d/hflux_x",hflux_x,ionode=0)
+  ENDIF
+END SUBROUTINE diagnose_0d
+
+SUBROUTINE diagnose_3d
+
+  USE basic
+  USE futils, ONLY: append, getatt, attach, putarrnd, putarr
+  USE fields
+  USE array
+  USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx, ikx, iky, ips_e, ips_i
+  USE time_integration
+  USE diagnostics_par
+  USE prec_const
+  USE processing
+  USE model, ONLY: KIN_E
+
+  IMPLICIT NONE
+
+  INTEGER        :: i_, root, world_rank, world_size
+  CHARACTER(256) :: dset_name
+  CALL append(fidres,  "/data/var3d/time",           time,ionode=0)
+  CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
+  CALL getatt(fidres,      "/data/var3d/",       "frames",iframe3d)
+  iframe3d=iframe3d+1
+  CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
+
+  IF (write_phi) CALL write_field3d_kykxz(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
+
+  IF (write_Na00) THEN
+    IF(KIN_E)THEN
+    IF (ips_e .EQ. 1) &
+      Ne00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_e(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
+    CALL write_field3d_kykxz(Ne00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ne00')
+    ENDIF
+    IF (ips_i .EQ. 1) &
+      Ni00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_i(ips_i,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
+    CALL write_field3d_kykxz(Ni00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ni00')
+
+    ! CALL compute_Napjz_spectrum
+    ! IF(KIN_E) &
+    ! CALL write_field3d_pjz_e(Nepjz(ips_e:ipe_e,ijs_e:ije_e,izs:ize), 'Nepjz')
+    ! CALL write_field3d_pjz_i(Nipjz(ips_i:ipe_i,ijs_i:ije_i,izs:ize), 'Nipjz')
+  ENDIF
+
+  !! Fuid moments
+  IF (write_dens .OR. write_fvel .OR. write_temp) &
+  CALL compute_fluid_moments
+
+  IF (write_dens) THEN
+    IF(KIN_E)&
+    CALL write_field3d_kykxz(dens_e(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_e')
+    CALL write_field3d_kykxz(dens_i(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_i')
+  ENDIF
+
+  IF (write_fvel) THEN
+    IF(KIN_E)&
+    CALL write_field3d_kykxz(upar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_e')
+    CALL write_field3d_kykxz(upar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_i')
+    IF(KIN_E)&
+    CALL write_field3d_kykxz(uper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_e')
+    CALL write_field3d_kykxz(uper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_i')
+  ENDIF
+
+  IF (write_temp) THEN
+    IF(KIN_E)&
+    CALL write_field3d_kykxz(Tpar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_e')
+    CALL write_field3d_kykxz(Tpar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_i')
+    IF(KIN_E)&
+    CALL write_field3d_kykxz(Tper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_e')
+    CALL write_field3d_kykxz(Tper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_i')
+    IF(KIN_E)&
+    CALL write_field3d_kykxz(temp_e(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_e')
+    CALL write_field3d_kykxz(temp_i(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_i')
+  ENDIF
+
+  CONTAINS
+
+  SUBROUTINE write_field3d_kykxz(field, text)
+    USE parallel, ONLY : gather_xyz
+    IMPLICIT NONE
+    COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: field
+    CHARACTER(*), INTENT(IN) :: text
+    COMPLEX(dp), DIMENSION(1:Nky,1:Nkx,1:Nz) :: field_full
+    CHARACTER(256) :: dset_name
+    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
+
+    IF (num_procs .EQ. 1) THEN ! no data distribution
+      CALL putarr(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), ionode=0)
+
+    ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
+      CALL gather_xyz(field(ikys:ikye,1:Nkx,izs:ize),field_full(1:Nky,1:Nkx,1:Nz))
+      CALL putarr(fidres, dset_name, field_full(1:Nky,1:Nkx,1:Nz), ionode=0)
+    ELSE ! output using putarrnd (very slow on marconi)
+      CALL putarrnd(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize),  (/1, 1, 3/))
+    ENDIF
+    CALL attach(fidres, dset_name, "time", time)
+  END SUBROUTINE write_field3d_kykxz
+
+  SUBROUTINE write_field3d_pjz_i(field, text)
+    IMPLICIT NONE
+    REAL(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,izs:ize), INTENT(IN) :: field
+    CHARACTER(*), INTENT(IN) :: text
+    CHARACTER(LEN=50) :: dset_name
+    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
+    IF (num_procs .EQ. 1) THEN ! no data distribution
+      CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize), ionode=0)
+    ELSE
+      CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,izs:ize),  (/1, 0, 3/))
+    ENDIF
+    CALL attach(fidres, dset_name, "time", time)
+  END SUBROUTINE write_field3d_pjz_i
+
+  SUBROUTINE write_field3d_pjz_e(field, text)
+    IMPLICIT NONE
+    REAL(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,izs:ize), INTENT(IN) :: field
+    CHARACTER(*), INTENT(IN) :: text
+    CHARACTER(LEN=50) :: dset_name
+    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
+    IF (num_procs .EQ. 1) THEN ! no data distribution
+      CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize), ionode=0)
+    ELSE
+      CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,izs:ize),  (/1, 0, 3/))
+    ENDIF
+    CALL attach(fidres, dset_name, "time", time)
+  END SUBROUTINE write_field3d_pjz_e
+
+END SUBROUTINE diagnose_3d
+
+SUBROUTINE diagnose_5d
+
+  USE basic
+  USE futils, ONLY: append, getatt, attach, putarrnd, putarr
+  USE fields
+  USE array!, ONLY: Sepj, Sipj
+  USE grid, ONLY: ips_e,ipe_e, ips_i, ipe_i, &
+                 ijs_e,ije_e, ijs_i, ije_i, &
+                 Np_i, Nj_i, Np_e, Nj_e, Nky, Nkx, Nz, &
+                 ikxs,ikxe,ikys,ikye,izs,ize
+  USE time_integration
+  USE diagnostics_par
+  USE prec_const
+  USE model, ONLY: KIN_E
+  IMPLICIT NONE
+  CHARACTER(LEN=50) :: dset_name
+
+  CALL append(fidres,  "/data/var5d/time",           time,ionode=0)
+  CALL append(fidres, "/data/var5d/cstep", real(cstep,dp),ionode=0)
+  CALL getatt(fidres,      "/data/var5d/",       "frames",iframe5d)
+  iframe5d=iframe5d+1
+  CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
+
+  IF (write_Napj) THEN
+   IF(KIN_E)&
+  CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_e')
+  CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_i')
+  ENDIF
+
+  IF (write_Sapj) THEN
+   IF(KIN_E)&
+   CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), 'Sepj')
+   CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), 'Sipj')
+  ENDIF
+
+  CONTAINS
+
+  SUBROUTINE write_field5d_e(field, text)
+    USE futils, ONLY: attach, putarr, putarrnd
+    USE parallel, ONLY: gather_pjxyz_e
+    USE grid,   ONLY: ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize
+    USE prec_const
+    IMPLICIT NONE
+
+    COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
+    CHARACTER(*), INTENT(IN) :: text
+
+    COMPLEX(dp), DIMENSION(1:Np_e,1:Nj_e,1:Nky,1:Nkx,1:Nz) :: field_full
+    CHARACTER(LEN=50) :: dset_name
+
+    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
+    IF (num_procs .EQ. 1) THEN
+     CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
+   ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
+     CALL gather_pjxyz_e(field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize),&
+                         field_full(1:Np_e,1:Nj_e,1:Nky,1:Nkx,1:Nz))
+     CALL putarr(fidres, dset_name, field_full(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz), ionode=0)
+   ELSE
+     CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize),  (/1,3,5/))
+    ENDIF
+    CALL attach(fidres, dset_name, 'cstep', cstep)
+    CALL attach(fidres, dset_name, 'time', time)
+    CALL attach(fidres, dset_name, 'jobnum', jobnum)
+    CALL attach(fidres, dset_name, 'dt', dt)
+    CALL attach(fidres, dset_name, 'iframe2d', iframe2d)
+    CALL attach(fidres, dset_name, 'iframe5d', iframe5d)
+
+  END SUBROUTINE write_field5d_e
+
+  SUBROUTINE write_field5d_i(field, text)
+    USE futils, ONLY: attach, putarr, putarrnd
+    USE parallel, ONLY: gather_pjxyz_i
+    USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize
+    USE prec_const
+    IMPLICIT NONE
+
+    COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
+    CHARACTER(*), INTENT(IN) :: text
+
+    COMPLEX(dp), DIMENSION(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz) :: field_full
+    CHARACTER(LEN=50) :: dset_name
+
+    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
+    IF (num_procs .EQ. 1) THEN
+      CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
+    ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
+      CALL gather_pjxyz_i(field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize),&
+                        field_full(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz))
+      CALL putarr(fidres, dset_name, field_full(1:Np_i,1:Nj_i,1:Nky,1:Nkx,1:Nz), ionode=0)
+    ELSE
+      CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize),  (/1,3,5/))
+    ENDIF
+    CALL attach(fidres, dset_name, 'cstep', cstep)
+    CALL attach(fidres, dset_name, 'time', time)
+    CALL attach(fidres, dset_name, 'jobnum', jobnum)
+    CALL attach(fidres, dset_name, 'dt', dt)
+    CALL attach(fidres, dset_name, 'iframe2d', iframe2d)
+    CALL attach(fidres, dset_name, 'iframe5d', iframe5d)
+
+  END SUBROUTINE write_field5d_i
+END SUBROUTINE diagnose_5d
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index cdc2fefa87e990ffae94a1d50baaf3460b622e7c..9eec0d8c7c3c5c685d768ebc05a91b3b2f4ce58b 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -68,7 +68,7 @@ MODULE fourier
     planb = fftw_mpi_plan_dft_c2r_2D(NX_, NY_, cmpx_data_f, real_data_f, communicator,  ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN))
 
    if ((.not. c_associated(planf)) .OR. (.not. c_associated(planb))) then
-      write(*,*) "plan creation error!!"
+      IF (my_id .EQ. 0) write(*,*) "plan creation error!!"
       stop
    end if
 
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 01d7e1ab9ebacad2b095a2da9a1c265ac1c8ca3e..77ea54e5313f6d0f3be9f0a5ac11efd1be01a220 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -48,16 +48,18 @@ MODULE grid
   INTEGER,  PUBLIC, PROTECTED  ::  izgs, izge ! ghosts
   LOGICAL,  PUBLIC, PROTECTED  ::  SG = .true.! shifted grid flag
   INTEGER,  PUBLIC :: ir,iz ! counters
-  ! Data about parallel distribution for kx
+  ! Data about parallel distribution for ky.kx
   integer(C_INTPTR_T), PUBLIC :: local_nkx, local_nky
   integer(C_INTPTR_T), PUBLIC :: local_nkx_offset, local_nky_offset
   INTEGER,             PUBLIC :: local_nkp
+  INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_nkx, counts_nky
+  INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_nkx, displs_nky
   ! "" for p
   INTEGER,             PUBLIC :: local_np_e, local_np_i
-  INTEGER,             PUBLIC :: total_np_e, total_np_i
+  INTEGER,             PUBLIC :: total_np_e, total_np_i, Np_e, Np_i
   integer(C_INTPTR_T), PUBLIC :: local_np_e_offset, local_np_i_offset
-  INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_np_e, counts_np_i
-  INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_np_e, displs_np_i
+  INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: rcv_p_e, rcv_p_i
+  INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: dsp_p_e, dsp_p_i
   ! "" for z
   INTEGER,             PUBLIC :: local_nz
   INTEGER,             PUBLIC :: total_nz
@@ -65,7 +67,7 @@ MODULE grid
   INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_nz
   INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_nz
   ! "" for j (not parallelized)
-  INTEGER,             PUBLIC :: local_nj_e, local_nj_i
+  INTEGER,             PUBLIC :: local_nj_e, local_nj_i, Nj_e, Nj_i
   ! Grids containing position in fourier space
   REAL(dp), DIMENSION(:),     ALLOCATABLE, PUBLIC :: kxarray, kxarray_full
   REAL(dp), DIMENSION(:),     ALLOCATABLE, PUBLIC :: kyarray, kyarray_full
@@ -161,6 +163,8 @@ CONTAINS
     ! Total number of Hermite polynomials we will evolve
     total_np_e = (Pmaxe/deltape) + 1
     total_np_i = (Pmaxi/deltapi) + 1
+    Np_e = total_np_e ! Reduced names (redundant)
+    Np_i = total_np_i
     ! Build the full grids on process 0 to diagnose it without comm
     ALLOCATE(parray_e_full(1:total_np_e))
     ALLOCATE(parray_i_full(1:total_np_i))
@@ -177,17 +181,17 @@ CONTAINS
     ipgs_e = ips_e - 2/deltape; ipge_e = ipe_e + 2/deltape;
     ipgs_i = ips_i - 2/deltapi; ipge_i = ipe_i + 2/deltapi;
     ! List of shift and local numbers between the different processes (used in scatterv and gatherv)
-    ALLOCATE(counts_np_e (1:num_procs_p))
-    ALLOCATE(counts_np_i (1:num_procs_p))
-    ALLOCATE(displs_np_e (1:num_procs_p))
-    ALLOCATE(displs_np_i (1:num_procs_p))
+    ALLOCATE(rcv_p_e (1:num_procs_p))
+    ALLOCATE(rcv_p_i (1:num_procs_p))
+    ALLOCATE(dsp_p_e (1:num_procs_p))
+    ALLOCATE(dsp_p_i (1:num_procs_p))
     DO in = 0,num_procs_p-1
       CALL decomp1D(total_np_e, num_procs_p, in, istart, iend)
-      counts_np_e(in+1) = iend-istart+1
-      displs_np_e(in+1) = istart-1
+      rcv_p_e(in+1) = iend-istart+1
+      dsp_p_e(in+1) = istart-1
       CALL decomp1D(total_np_i, num_procs_p, in, istart, iend)
-      counts_np_i(in+1) = iend-istart+1
-      displs_np_i(in+1) = istart-1
+      rcv_p_i(in+1) = iend-istart+1
+      dsp_p_i(in+1) = istart-1
     ENDDO
 
     ! local grid computation
@@ -244,7 +248,9 @@ CONTAINS
     USE prec_const
     IMPLICIT NONE
     INTEGER :: ij
-
+    ! Total number of J degrees
+    Nj_e = jmaxe+1
+    Nj_i = jmaxi+1
     ! Build the full grids on process 0 to diagnose it without comm
     ALLOCATE(jarray_e_full(1:jmaxe+1))
     ALLOCATE(jarray_i_full(1:jmaxi+1))
@@ -275,7 +281,7 @@ CONTAINS
     USE prec_const
     USE model, ONLY: LINEARITY
     IMPLICIT NONE
-    INTEGER :: i_
+    INTEGER :: i_, in, istart, iend
     Nky = Ny/2+1 ! Defined only on positive kx since fields are real
     ! Grid spacings
     IF (Ny .EQ. 1) THEN
@@ -297,6 +303,14 @@ CONTAINS
     ikye = ikys + local_nky - 1
     ALLOCATE(kyarray(ikys:ikye))
     local_kymax = 0._dp
+    ! List of shift and local numbers between the different processes (used in scatterv and gatherv)
+    ALLOCATE(counts_nky (1:num_procs_ky))
+    ALLOCATE(displs_nky (1:num_procs_ky))
+    DO in = 0,num_procs_ky-1
+      CALL decomp1D(Nky, num_procs_ky, in, istart, iend)
+      counts_nky(in+1) = iend-istart+1
+      displs_nky(in+1) = istart-1
+    ENDDO
     ! Creating a grid ordered as dk*(0 1 2 3)
     DO iky = ikys,ikye
       kyarray(iky) = REAL(iky-1,dp) * deltaky
@@ -395,6 +409,7 @@ CONTAINS
 
   SUBROUTINE set_zgrid
     USE prec_const
+    USE model, ONLY: mu_z
     IMPLICIT NONE
     INTEGER :: i_, fid
     REAL    :: grid_shift, Lz
@@ -405,9 +420,12 @@ CONTAINS
     ! Z stepping (#interval = #points since periodic)
     deltaz        = Lz/REAL(Nz,dp)
     inv_deltaz    = 1._dp/deltaz
-    diff_dz_coeff = -(deltaz/2._dp)**4 ! adaptive fourth derivative
-    ! non adaptive
-    ! diff_dz_coeff = -1._dp
+    ! Parallel hyperdiffusion coefficient
+    IF(mu_z .GT. 0) THEN
+      diff_dz_coeff = -(deltaz/2._dp)**4 ! adaptive fourth derivative (~GENE)
+    ELSE
+      diff_dz_coeff = 1._dp    ! non adaptive (positive sign to compensate mu_z neg)
+    ENDIF
     IF (SG) THEN
       grid_shift = deltaz/2._dp
     ELSE
diff --git a/src/memory.F90 b/src/memory.F90
index 3b12ca1a5ae9a99eb69eb00e370c44ae744b6c03..c9509278094ed98800c478c2e1469bf47a0728db 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -124,4 +124,5 @@ SUBROUTINE memory
       ENDIF
       CALL allocate_array(  Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), 1,1, 1,1, 1,1)
  ENDIF
+
 END SUBROUTINE memory
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index db861edde3c92bba30ea5d711579507fdf8df9d2..2931da9aa3dc7cd72313502b8d70a2e9b35dbec7 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -59,7 +59,7 @@ SUBROUTINE moments_eq_rhs_e
             eo    = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
             kperp2= kparray(iky,ikx,iz,eo)**2
 
-          IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN
+          IF((CLOS .EQ. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN
             !! Compute moments mixing terms
             Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
             ! Perpendicular dynamic
@@ -179,7 +179,7 @@ SUBROUTINE moments_eq_rhs_i
             p_int = parray_i(ip)    ! Hermite degree
             eo    = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
             kperp2= kparray(iky,ikx,iz,eo)**2
-            IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN
+            IF((CLOS .EQ. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN
               !! Compute moments mixing terms
               Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
               ! Perpendicular dynamic
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..96d8c27749b423139c06c832d946f945d7885917
--- /dev/null
+++ b/src/parallel_mod.F90
@@ -0,0 +1,410 @@
+MODULE parallel
+  USE basic
+  USE grid, ONLY: Nkx,Nky,Nz, ikys,ikye, izs,ize, local_nky, local_nz, &
+                  local_np_e, local_np_i, Np_e, Np_i, Nj_e, Nj_i, &
+                  pmaxi, pmaxe, ips_e, ipe_e, ips_i, ipe_i, &
+                  jmaxi, jmaxe, ijs_e, ije_e, ijs_i, ije_i, &
+                  rcv_p_e, rcv_p_i, dsp_p_e, dsp_p_i
+  use prec_const
+  IMPLICIT NONE
+  ! recieve and displacement counts for gatherv
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_y, dsp_y, rcv_zy, dsp_zy
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_zp_e,  dsp_zp_e
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_yp_e,  dsp_yp_e
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_zyp_e, dsp_zyp_e
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_zp_i,  dsp_zp_i
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_yp_i,  dsp_yp_i
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_zyp_i, dsp_zyp_i
+
+  PUBLIC :: manual_0D_bcast, manual_3D_bcast, init_parallel_var, &
+            gather_xyz, gather_pjz_i, gather_pjxyz_e, gather_pjxyz_i
+
+CONTAINS
+
+  SUBROUTINE init_parallel_var
+    INTEGER :: i_
+
+    !!!!!! XYZ gather variables
+    !! Y reduction at constant x and z
+    ! number of points recieved and displacement for the y reduction
+    ALLOCATE(rcv_y(0:num_procs_ky-1),dsp_y(0:num_procs_ky-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nky,1,MPI_INTEGER,rcv_y,1,MPI_INTEGER,comm_ky,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_y(0)=0
+    DO i_=1,num_procs_ky-1
+       dsp_y(i_) =dsp_y(i_-1) + rcv_y(i_-1)
+    END DO
+    !! Z reduction for full slices of y data but constant x
+    ! number of points recieved and displacement for the z reduction
+    ALLOCATE(rcv_zy(0:num_procs_z-1),dsp_zy(0:num_procs_z-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nz*Nky,1,MPI_INTEGER,rcv_zy,1,MPI_INTEGER,comm_z,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_zy(0)=0
+    DO i_=1,num_procs_z-1
+       dsp_zy(i_) =dsp_zy(i_-1) + rcv_zy(i_-1)
+    END DO
+
+    !!!!! PJZ gather variables
+    ! IONS
+    !! Z reduction for full slices of p data but constant j
+    ! number of points recieved and displacement for the z reduction
+    ALLOCATE(rcv_zp_i(0:num_procs_z-1),dsp_zp_i(0:num_procs_z-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nz*Np_i,1,MPI_INTEGER,rcv_zp_i,1,MPI_INTEGER,comm_z,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_zp_i(0)=0
+    DO i_=1,num_procs_z-1
+       dsp_zp_i(i_) =dsp_zp_i(i_-1) + rcv_zp_i(i_-1)
+    END DO
+
+    !!!!! PJXYZ gather variables
+    !! Y reduction for full slices of p data but constant j
+    ! number of points recieved and displacement for the y reduction
+    ALLOCATE(rcv_yp_i(0:num_procs_ky-1),dsp_yp_i(0:num_procs_ky-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nky*Np_i,1,MPI_INTEGER,rcv_yp_i,1,MPI_INTEGER,comm_ky,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_yp_i(0)=0
+    DO i_=1,num_procs_ky-1
+       dsp_yp_i(i_) =dsp_yp_i(i_-1) + rcv_yp_i(i_-1)
+    END DO
+    !! Z reduction for full slices of py data but constant j
+    ! number of points recieved and displacement for the z reduction
+    ALLOCATE(rcv_zyp_i(0:num_procs_z-1),dsp_zyp_i(0:num_procs_z-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nz*Np_i*Nky,1,MPI_INTEGER,rcv_zyp_i,1,MPI_INTEGER,comm_z,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_zyp_i(0)=0
+    DO i_=1,num_procs_z-1
+       dsp_zyp_i(i_) =dsp_zyp_i(i_-1) + rcv_zyp_i(i_-1)
+    END DO
+
+    ! ELECTONS
+    !! Z reduction for full slices of p data but constant j
+    ! number of points recieved and displacement for the z reduction
+    ALLOCATE(rcv_zp_e(0:num_procs_z-1),dsp_zp_e(0:num_procs_z-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nz*Np_e,1,MPI_INTEGER,rcv_zp_e,1,MPI_INTEGER,comm_z,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_zp_e(0)=0
+    DO i_=1,num_procs_z-1
+       dsp_zp_e(i_) =dsp_zp_e(i_-1) + rcv_zp_e(i_-1)
+    END DO
+
+    !!!!! PJXYZ gather variables
+    !! Y reduction for full slices of p data but constant j
+    ! number of points recieved and displacement for the y reduction
+    ALLOCATE(rcv_yp_e(0:num_procs_ky-1),dsp_yp_e(0:num_procs_ky-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nky*Np_e,1,MPI_INTEGER,rcv_yp_e,1,MPI_INTEGER,comm_ky,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_yp_e(0)=0
+    DO i_=1,num_procs_ky-1
+       dsp_yp_e(i_) =dsp_yp_e(i_-1) + rcv_yp_e(i_-1)
+    END DO
+    !! Z reduction for full slices of py data but constant j
+    ! number of points recieved and displacement for the z reduction
+    ALLOCATE(rcv_zyp_e(0:num_procs_z-1),dsp_zyp_e(0:num_procs_z-1)) !Displacement sizes for balance diagnostic
+    ! all processes share their local number of points
+    CALL MPI_ALLGATHER(local_nz*Np_e*Nky,1,MPI_INTEGER,rcv_zyp_e,1,MPI_INTEGER,comm_z,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_zyp_e(0)=0
+    DO i_=1,num_procs_z-1
+       dsp_zyp_e(i_) =dsp_zyp_e(i_-1) + rcv_zyp_e(i_-1)
+    END DO
+
+  END SUBROUTINE init_parallel_var
+
+  !!!!! Gather a field in spatial coordinates on rank 0 !!!!!
+  SUBROUTINE gather_xyz(field_sub,field_full)
+    COMPLEX(dp), DIMENSION(ikys:ikye, 1:Nkx, izs:ize), INTENT(IN)    :: field_sub
+    COMPLEX(dp), DIMENSION(   1:Nky,  1:Nkx,   1:Nz),  INTENT(INOUT) :: field_full
+    COMPLEX(dp), DIMENSION(ikys:ikye)           :: buffer_ly_cz !local y, constant z
+    COMPLEX(dp), DIMENSION(   1:Nky )           :: buffer_fy_cz !full  y, constant z
+    COMPLEX(dp), DIMENSION(   1:Nky,  izs:ize ) :: buffer_fy_lz !full  y, local z
+    COMPLEX(dp), DIMENSION(   1:Nky,    1:Nz  ) :: buffer_fy_fz !full  y, full  z
+    INTEGER :: snd_y, snd_z, root_p, root_z, root_ky, ix, iz
+
+    snd_y  = local_nky    ! Number of points to send along y (per z)
+    snd_z  = Nky*local_nz ! Number of points to send along z (full y)
+
+    root_p = 0; root_z = 0; root_ky = 0
+    IF(rank_p .EQ. root_p) THEN
+      DO ix = 1,Nkx
+        DO iz = izs,ize
+          ! fill a buffer to contain a slice of data at constant kx and z
+          buffer_ly_cz(ikys:ikye) = field_sub(ikys:ikye,ix,iz)
+          CALL MPI_GATHERV(buffer_ly_cz, snd_y,        MPI_DOUBLE_COMPLEX, &
+                           buffer_fy_cz, rcv_y, dsp_y, MPI_DOUBLE_COMPLEX, &
+                           root_ky, comm_ky, ierr)
+          buffer_fy_lz(1:Nky,iz) = buffer_fy_cz(1:Nky)
+        ENDDO
+
+        ! send the full line on y contained by root_kyas
+        IF(rank_ky .EQ. 0) THEN
+          CALL MPI_GATHERV(buffer_fy_lz, snd_z,        MPI_DOUBLE_COMPLEX, &
+                           buffer_fy_fz, rcv_zy, dsp_zy, MPI_DOUBLE_COMPLEX, &
+                           root_z, comm_z, ierr)
+        ENDIF
+        ! ID 0 (the one who output) rebuild the whole array
+        IF(my_id .EQ. 0) &
+          field_full(1:Nky,ix,1:Nz) = buffer_fy_fz(1:Nky,1:Nz)
+      ENDDO
+    ENDIF
+  END SUBROUTINE gather_xyz
+
+  !!!!! Gather a field in kinetic + z coordinates on rank 0 !!!!!
+  SUBROUTINE gather_pjz_i(field_sub,field_full)
+    COMPLEX(dp), DIMENSION(ips_i:ipe_i, ijs_i:ije_i, izs:ize), INTENT(IN)    :: field_sub
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1,  1:jmaxi+1,   1:Nz),  INTENT(INOUT) :: field_full
+    COMPLEX(dp), DIMENSION(ips_i:ipe_i)         :: buffer_lp_cz !local p, constant z
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1 )       :: buffer_fp_cz !full  p, constant z
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1,  izs:ize ) :: buffer_fp_lz !full  p, local z
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1,    1:Nz  ) :: buffer_fp_fz !full  p, full  z
+    INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz
+
+    snd_p  = local_np_i   ! Number of points to send along y (per z)
+    snd_z  = Np_i*local_nz ! Number of points to send along z (full y)
+
+    root_p = 0; root_z = 0; root_ky = 0
+    IF(rank_ky .EQ. root_ky) THEN
+      DO ij = 1,jmaxi+1
+        DO iz = izs,ize
+          ! fill a buffer to contain a slice of data at constant kx and z
+          buffer_lp_cz(ips_i:ipe_i) = field_sub(ips_i:ipe_i,ij,iz)
+          CALL MPI_GATHERV(buffer_lp_cz, snd_p,            MPI_DOUBLE_COMPLEX, &
+                           buffer_fp_cz, rcv_p_i, dsp_p_i, MPI_DOUBLE_COMPLEX, &
+                           root_p, comm_p, ierr)
+          buffer_fp_lz(1:Np_i,iz) = buffer_fp_cz(1:Np_i)
+        ENDDO
+
+        ! send the full line on y contained by root_kyas
+        IF(rank_p .EQ. 0) THEN
+          CALL MPI_GATHERV(buffer_fp_lz, snd_z,              MPI_DOUBLE_COMPLEX, &
+                           buffer_fp_fz, rcv_zp_i, dsp_zp_i, MPI_DOUBLE_COMPLEX, &
+                           root_z, comm_z, ierr)
+        ENDIF
+        ! ID 0 (the one who output) rebuild the whole array
+        IF(my_id .EQ. 0) &
+          field_full(1:Np_i,ij,1:Nz) = buffer_fp_fz(1:Np_i,1:Nz)
+      ENDDO
+    ENDIF
+  END SUBROUTINE gather_pjz_i
+
+  SUBROUTINE gather_pjz_e(field_sub,field_full)
+    COMPLEX(dp), DIMENSION(ips_e:ipe_e, ijs_e:ije_e, izs:ize), INTENT(IN)    :: field_sub
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1,  1:jmaxi+1,   1:Nz),  INTENT(INOUT) :: field_full
+    COMPLEX(dp), DIMENSION(ips_e:ipe_e)         :: buffer_lp_cz !local p, constant z
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1 )       :: buffer_fp_cz !full  p, constant z
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1,  izs:ize ) :: buffer_fp_lz !full  p, local z
+    COMPLEX(dp), DIMENSION(   1:pmaxi+1,    1:Nz  ) :: buffer_fp_fz !full  p, full  z
+    INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz
+
+    snd_p  = local_np_e   ! Number of points to send along y (per z)
+    snd_z  = Np_e*local_nz ! Number of points to send along z (full y)
+
+    root_p = 0; root_z = 0; root_ky = 0
+    IF(rank_ky .EQ. root_ky) THEN
+      DO ij = 1,jmaxi+1
+        DO iz = izs,ize
+          ! fill a buffer to contain a slice of data at constant kx and z
+          buffer_lp_cz(ips_e:ipe_e) = field_sub(ips_e:ipe_e,ij,iz)
+          CALL MPI_GATHERV(buffer_lp_cz, snd_p,            MPI_DOUBLE_COMPLEX, &
+                           buffer_fp_cz, rcv_p_e, dsp_p_e, MPI_DOUBLE_COMPLEX, &
+                           root_p, comm_p, ierr)
+          buffer_fp_lz(1:Np_e,iz) = buffer_fp_cz(1:Np_e)
+        ENDDO
+
+        ! send the full line on y contained by root_kyas
+        IF(rank_p .EQ. 0) THEN
+          CALL MPI_GATHERV(buffer_fp_lz, snd_z,              MPI_DOUBLE_COMPLEX, &
+                           buffer_fp_fz, rcv_zp_e, dsp_zp_e, MPI_DOUBLE_COMPLEX, &
+                           root_z, comm_z, ierr)
+        ENDIF
+        ! ID 0 (the one who output) rebuild the whole array
+        IF(my_id .EQ. 0) &
+          field_full(1:Np_e,ij,1:Nz) = buffer_fp_fz(1:Np_e,1:Nz)
+      ENDDO
+    ENDIF
+  END SUBROUTINE gather_pjz_e
+
+  !!!!! Gather a field in kinetic + spatial coordinates on rank 0 !!!!!
+  !!!!! Gather a field in spatial coordinates on rank 0 !!!!!
+  SUBROUTINE gather_pjxyz_i(field_sub,field_full)
+    COMPLEX(dp), DIMENSION( ips_i:ipe_i, 1:Nj_i, ikys:ikye, 1:Nkx, izs:ize), INTENT(IN)    :: field_sub
+    COMPLEX(dp), DIMENSION(     1:Np_i,  1:Nj_i,    1:Nky,  1:Nkx,   1:Nz),  INTENT(INOUT) :: field_full
+    COMPLEX(dp), DIMENSION(ips_i:ipe_i)       :: buffer_lp_cy_cz     !local p, constant y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_i)            :: buffer_fp_cy_cz     ! full p, constant y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_i, ikys:ikye)    :: buffer_fp_ly_cz     ! full p,    local y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_i,    1:Nky )    :: buffer_fp_fy_cz     ! full p,     full y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_i,    1:Nky, izs:ize ) :: buffer_fp_fy_lz ! full p,     full y, local    z
+    COMPLEX(dp), DIMENSION(1:Np_i,    1:Nky,   1:Nz  ) :: buffer_fp_fy_fz ! full p,     full y, full     z
+    INTEGER :: snd_p, snd_y, snd_z, root_p, root_z, root_ky, iy, ix, iz, ij
+
+    snd_p  =     local_np_i    ! Number of points to send along p (per z,y)
+    snd_y  = Np_i*local_nky    ! Number of points to send along y (per z, full p)
+    snd_z  = Np_i*Nky*local_nz ! Number of points to send along z (full y, full p)
+
+    root_p = 0; root_z = 0; root_ky = 0
+
+    j: DO ij = 1,Nj_i
+      x: DO ix = 1,Nkx
+        z: DO iz = izs,ize
+          y: DO iy = ikys,ikye
+            ! fill a buffer to contain a slice of p data at constant j, ky, kx and z
+            buffer_lp_cy_cz(ips_i:ipe_i) = field_sub(ips_i:ipe_i,ij,iy,ix,iz)
+            CALL MPI_GATHERV(buffer_lp_cy_cz, snd_p,            MPI_DOUBLE_COMPLEX, &
+                             buffer_fp_cy_cz, rcv_p_i, dsp_p_i, MPI_DOUBLE_COMPLEX, &
+                             root_p, comm_p, ierr)
+            buffer_fp_ly_cz(1:Np_i,iy) = buffer_fp_cy_cz(1:Np_i)
+          ENDDO y
+          ! send the full line on p contained by root_p
+          IF(rank_p .EQ. 0) THEN
+            CALL MPI_GATHERV(buffer_fp_ly_cz, snd_y,        MPI_DOUBLE_COMPLEX, &
+                             buffer_fp_fy_cz, rcv_yp_i, dsp_yp_i, MPI_DOUBLE_COMPLEX, &
+                             root_ky, comm_ky, ierr)
+            buffer_fp_fy_lz(1:Np_i,1:Nky,iz) = buffer_fp_fy_cz(1:Np_i,1:Nky)
+          ENDIF
+        ENDDO z
+        ! send the full line on y contained by root_kyas
+        IF(rank_ky .EQ. 0) THEN
+          CALL MPI_GATHERV(buffer_fp_fy_lz, snd_z,                MPI_DOUBLE_COMPLEX, &
+                           buffer_fp_fy_fz, rcv_zyp_i, dsp_zyp_i, MPI_DOUBLE_COMPLEX, &
+                           root_z, comm_z, ierr)
+        ENDIF
+        ! ID 0 (the one who output) rebuild the whole array
+        IF(my_id .EQ. 0) &
+          field_full(1:Np_i,ij,1:Nky,ix,1:Nz) = buffer_fp_fy_fz(1:Np_i,1:Nky,1:Nz)
+      ENDDO x
+    ENDDO j
+
+  END SUBROUTINE gather_pjxyz_i
+
+  SUBROUTINE gather_pjxyz_e(field_sub,field_full)
+    COMPLEX(dp), DIMENSION( ips_e:ipe_e, 1:Nj_e, ikys:ikye, 1:Nkx, izs:ize), INTENT(IN)    :: field_sub
+    COMPLEX(dp), DIMENSION(     1:Np_e,  1:Nj_e,    1:Nky,  1:Nkx,   1:Nz),  INTENT(INOUT) :: field_full
+    COMPLEX(dp), DIMENSION(ips_e:ipe_e)       :: buffer_lp_cy_cz     !local p, constant y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_e)            :: buffer_fp_cy_cz     ! full p, constant y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_e, ikys:ikye)    :: buffer_fp_ly_cz     ! full p,    local y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_e,    1:Nky )    :: buffer_fp_fy_cz     ! full p,     full y, constant z
+    COMPLEX(dp), DIMENSION(1:Np_e,    1:Nky, izs:ize ) :: buffer_fp_fy_lz ! full p,     full y, local    z
+    COMPLEX(dp), DIMENSION(1:Np_e,    1:Nky,   1:Nz  ) :: buffer_fp_fy_fz ! full p,     full y, full     z
+    INTEGER :: snd_p, snd_y, snd_z, root_p, root_z, root_ky, iy, ix, iz, ij
+
+    snd_p  =     local_np_e    ! Number of points to send along p (per z,y)
+    snd_y  = Np_e*local_nky    ! Number of points to send along y (per z, full p)
+    snd_z  = Np_e*Nky*local_nz ! Number of points to send along z (full y, full p)
+
+    root_p = 0; root_z = 0; root_ky = 0
+
+    j: DO ij = 1,Nj_e
+      x: DO ix = 1,Nkx
+        z: DO iz = izs,ize
+          y: DO iy = ikys,ikye
+            ! fill a buffer to contain a slice of p data at constant j, ky, kx and z
+            buffer_lp_cy_cz(ips_e:ipe_e) = field_sub(ips_e:ipe_e,ij,iy,ix,iz)
+            CALL MPI_GATHERV(buffer_lp_cy_cz, snd_p,            MPI_DOUBLE_COMPLEX, &
+                             buffer_fp_cy_cz, rcv_p_e, dsp_p_e, MPI_DOUBLE_COMPLEX, &
+                             root_p, comm_p, ierr)
+            buffer_fp_ly_cz(1:Np_e,iy) = buffer_fp_cy_cz(1:Np_e)
+          ENDDO y
+          ! send the full line on p contained by root_p
+          IF(rank_p .EQ. 0) THEN
+            CALL MPI_GATHERV(buffer_fp_ly_cz, snd_y,        MPI_DOUBLE_COMPLEX, &
+                             buffer_fp_fy_cz, rcv_yp_e, dsp_yp_e, MPI_DOUBLE_COMPLEX, &
+                             root_ky, comm_ky, ierr)
+            buffer_fp_fy_lz(1:Np_e,1:Nky,iz) = buffer_fp_fy_cz(1:Np_e,1:Nky)
+          ENDIF
+        ENDDO z
+        ! send the full line on y contained by root_kyas
+        IF(rank_ky .EQ. 0) THEN
+          CALL MPI_GATHERV(buffer_fp_fy_lz, snd_z,                MPI_DOUBLE_COMPLEX, &
+                           buffer_fp_fy_fz, rcv_zyp_e, dsp_zyp_e, MPI_DOUBLE_COMPLEX, &
+                           root_z, comm_z, ierr)
+        ENDIF
+        ! ID 0 (the one who output) rebuild the whole array
+        IF(my_id .EQ. 0) &
+          field_full(1:Np_e,ij,1:Nky,ix,1:Nz) = buffer_fp_fy_fz(1:Np_e,1:Nky,1:Nz)
+      ENDDO x
+    ENDDO j
+
+  END SUBROUTINE gather_pjxyz_e
+
+  !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
+  SUBROUTINE manual_3D_bcast(field_)
+    USE grid
+    IMPLICIT NONE
+    COMPLEX(dp), INTENT(INOUT) :: field_(ikys:ikye,ikxs:ikxe,izs:ize)
+    COMPLEX(dp) :: buffer(ikys:ikye,ikxs:ikxe,izs:ize)
+    INTEGER     :: i_, root, world_rank, world_size, count
+    root = 0;
+    count = (ikye-ikys+1) * (ikxe-ikxs+1) * (ize-izs+1);
+
+    CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
+    CALL MPI_COMM_SIZE(comm_p,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 iz = izs,ize
+          DO ikx = ikxs,ikxe
+            DO iky = ikys,ikye
+                buffer(iky,ikx,iz) = field_(iky,ikx,iz)
+              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, count, MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
+        ENDDO
+      ELSE
+        ! Recieve buffer from root
+        CALL MPI_RECV(buffer, count, MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
+        ! Write it in phi
+        DO iz = izs,ize
+          DO ikx = ikxs,ikxe
+            DO iky = ikys,ikye
+              field_(iky,ikx,iz) = buffer(iky,ikx,iz)
+            ENDDO
+          ENDDO
+        ENDDO
+      ENDIF
+    ENDIF
+  END SUBROUTINE manual_3D_bcast
+
+  !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
+  SUBROUTINE manual_0D_bcast(v)
+    USE grid
+    IMPLICIT NONE
+    COMPLEX(dp), INTENT(INOUT) :: v
+    COMPLEX(dp) :: buffer
+    INTEGER     :: i_, root, world_rank, world_size, count
+    root = 0;
+    count = 1;
+
+    CALL MPI_COMM_RANK(comm_z,world_rank,ierr)
+    CALL MPI_COMM_SIZE(comm_z,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
+        buffer = v
+        ! Send it to all the other processes
+        DO i_ = 0,num_procs_z-1
+          IF (i_ .NE. world_rank) &
+          CALL MPI_SEND(buffer, count, MPI_DOUBLE_COMPLEX, i_, 0, comm_z, ierr)
+        ENDDO
+      ELSE
+        ! Recieve buffer from root
+        CALL MPI_RECV(buffer, count, MPI_DOUBLE_COMPLEX, root, 0, comm_z, MPI_STATUS_IGNORE, ierr)
+        ! Write it in phi
+        v = buffer
+      ENDIF
+    ENDIF
+  END SUBROUTINE manual_0D_bcast
+
+
+END MODULE parallel
diff --git a/src/poisson.F90 b/src/poisson.F90
index 3936324eee0339353c29520b0486c7fd79ef04ee..75e2b887a670919ba094a6c80df19f307761d65f 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -7,7 +7,7 @@ SUBROUTINE poisson
   USE fields
   USE grid
   USE calculus, ONLY : simpson_rule_z
-  USE utility,  ONLY : manual_3D_bcast
+  USE parallel,  ONLY : manual_3D_bcast
   use model, ONLY : qe2_taue, qi2_taui, q_e, q_i, lambdaD, KIN_E
   USE processing, ONLY : compute_density
   USE prec_const
diff --git a/src/ppinit.F90 b/src/ppinit.F90
index 55db146013a68d8f93ed35fe625249bdbb27d4c2..1863dcf9852c2e85ab2a9fe87e1627bb8de368c4 100644
--- a/src/ppinit.F90
+++ b/src/ppinit.F90
@@ -11,7 +11,8 @@ SUBROUTINE ppinit
   INTEGER, DIMENSION(ndims) :: dims=0, coords=0, coords_L=0, coords_R=0
   LOGICAL :: periods(ndims) = .FALSE., reorder=.FALSE.
   CHARACTER(len=32) :: str
-  INTEGER :: nargs, i, l
+  INTEGER :: nargs, i, l, r_
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rank2exclude
 
   CALL MPI_INIT(ierr)
 
@@ -49,9 +50,8 @@ SUBROUTINE ppinit
   periods(3)=.TRUE.
 
   CALL MPI_CART_CREATE(MPI_COMM_WORLD, ndims, dims, periods, reorder, comm0, ierr)
-
+  CALL MPI_COMM_GROUP(comm0,group0, ierr)
   CALL MPI_COMM_RANK(comm0, rank_0,  ierr)
-
   CALL MPI_CART_COORDS(comm0,rank_0,ndims,coords,ierr)
 
   !
@@ -68,10 +68,25 @@ SUBROUTINE ppinit
   ! 2D communicator
   CALL MPI_CART_SUB (comm0, (/.TRUE.,.FALSE.,.TRUE./),  comm_pz,  ierr)
   CALL MPI_CART_SUB (comm0, (/.FALSE.,.TRUE.,.TRUE./),  comm_kyz, ierr)
+  ! Count the number of processes in 2D comms
+  CALL MPI_COMM_SIZE(comm_pz, num_procs_pz, ierr)
+  CALL MPI_COMM_SIZE(comm_kyz,num_procs_kyz,ierr)
+  ! Find id inside the 1d-sub communicators
+  CALL MPI_COMM_RANK(comm_pz,  rank_pz,   ierr)
+  CALL MPI_COMM_RANK(comm_kyz, rank_kyz,  ierr)
 
   ! Find neighbours
   CALL MPI_CART_SHIFT(comm0, 0, 1, nbr_L, nbr_R, ierr) !left   right neighbours
   CALL MPI_CART_SHIFT(comm0, 1, 1, nbr_B, nbr_T, ierr) !bottom top   neighbours
   CALL MPI_CART_SHIFT(comm0, 2, 1, nbr_D, nbr_U, ierr) !down   up    neighbours
 
+  ! Create the communicator for groups used in gatherv
+  ! CALL MPI_COMM_GROUP(comm0,group_ky0)
+  ! ALLOCATE(rank2include(0:num_procs_ky))
+  ! DO r_ = 0,rank_0
+  !   IF(rank_y .EQ. 0) &
+  !     rank2exclude
+  ! ENDDO
+  ! CALL MPI_COMM_CREATE_GROUPE(comm0, group_p0, comm_ky0)
+
 END SUBROUTINE ppinit
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index ffa5ac89ebc4615c31eb8ccf78aebe17c0a01370..048266cb9aa2f658284ddd979b7b0e8490b2d760 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -375,9 +375,6 @@ SUBROUTINE compute_density
         ENDDO
       ENDDO
   ENDIF
-  ! IF(KIN_E)&
-  ! CALL manual_3D_bcast(dens_e(ikxs:ikxe,ikys:ikye,izs:ize))
-  ! CALL manual_3D_bcast(dens_i(ikxs:ikxe,ikys:ikye,izs:ize))
 END SUBROUTINE compute_density
 
 ! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre)
@@ -412,9 +409,6 @@ SUBROUTINE compute_uperp
         ENDDO
       ENDDO
   ENDIF
-  ! IF(KIN_E)&
-  ! CALL manual_3D_bcast(uper_e(ikxs:ikxe,ikys:ikye,izs:ize))
-  ! CALL manual_3D_bcast(uper_i(ikxs:ikxe,ikys:ikye,izs:ize))
 END SUBROUTINE compute_uperp
 
 ! Compute the 2D particle fluid par velocity for electron and ions (sum over Laguerre)
@@ -450,9 +444,6 @@ SUBROUTINE compute_upar
     upar_e = 0
     upar_i = 0
   ENDIF
-  ! IF(KIN_E)&
-  ! CALL manual_3D_bcast(upar_e(ikxs:ikxe,ikys:ikye,izs:ize))
-  ! CALL manual_3D_bcast(upar_i(ikxs:ikxe,ikys:ikye,izs:ize))
 END SUBROUTINE compute_upar
 
 ! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
@@ -495,9 +486,6 @@ SUBROUTINE compute_tperp
         ENDDO
       ENDDO
   ENDIF
-  ! IF(KIN_E)&
-  ! CALL manual_3D_bcast(Tper_e(ikxs:ikxe,ikys:ikye,izs:ize))
-  ! CALL manual_3D_bcast(Tper_i(ikxs:ikxe,ikys:ikye,izs:ize))
 END SUBROUTINE compute_Tperp
 
 ! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
@@ -538,9 +526,6 @@ SUBROUTINE compute_Tpar
         ENDDO
       ENDDO
   ENDIF
-  ! IF(KIN_E)&
-  ! CALL manual_3D_bcast(Tpar_e(ikxs:ikxe,ikys:ikye,izs:ize))
-  ! CALL manual_3D_bcast(Tpar_i(ikxs:ikxe,ikys:ikye,izs:ize))
 END SUBROUTINE compute_Tpar
 
 ! Compute the 2D particle fluid moments for electron and ions (sum over Laguerre)
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index ea1697c311180786225e5ce74880114b9108f669..581b5d271e938fd3a20c7e2d8c475ef9337c8ff0 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -39,13 +39,11 @@ CONTAINS
         ENDIF
         CALL getatt(fidrst,"/data/input/" , "pmaxi", pmaxi_cp)
         CALL getatt(fidrst,"/data/input/" , "jmaxi", jmaxi_cp)
-        IF (my_id .EQ. 0) WRITE(*,*) "Pi_cp = ", pmaxi_cp
-        IF (my_id .EQ. 0) WRITE(*,*) "Ji_cp = ", jmaxi_cp
         CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
 
         IF ((KIN_E .AND. ((pmaxe_cp .NE. pmaxe) .OR. (jmaxe_cp .NE. jmaxe))) .OR.&
          (pmaxi_cp .NE. pmaxi) .OR. (jmaxi_cp .NE. jmaxi)) THEN
-         IF(my_id.EQ.0)WRITE(*,*) '! Extending the polynomials basis !'
+         IF(my_id.EQ.0) WRITE(*,*) '! Extending the polynomials basis !'
          CALL load_output_adapt_pj
         ELSE
 
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index b6e989a746f8dedcdf7f634ee7d769835060af71..b962a4c4c5c4f6113f7e8b92aa0b1735c7b2c96f 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -1,11 +1,9 @@
 MODULE utility
 
   USE basic
-
+  USE grid, ONLY: Nkx,Nky,Nz, ikys,ikye, izs,ize, local_nky, local_nz
   use prec_const
   IMPLICIT NONE
-  PUBLIC :: manual_0D_bcast, manual_2D_bcast, manual_3D_bcast
-
 CONTAINS
 
   FUNCTION is_nan(x,str) RESULT(isn)
@@ -69,118 +67,4 @@ CONTAINS
     mlend= is_nan( REAL(sumfield),str).OR.is_inf( REAL(sumfield),str) &
       .OR. is_nan(AIMAG(sumfield),str).OR.is_inf(AIMAG(sumfield),str)
   END FUNCTION checkfield
-
-
-    !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
-  SUBROUTINE manual_2D_bcast(field_)
-    USE grid
-    IMPLICIT NONE
-    COMPLEX(dp), INTENT(INOUT) :: field_(ikys:ikye,ikxs:ikxe)
-    COMPLEX(dp) :: buffer(ikys:ikye,ikxs:ikxe)
-    INTEGER     :: i_, root, world_rank, world_size
-    root = 0;
-    CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
-    CALL MPI_COMM_SIZE(comm_p,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 ikx = ikxs,ikxe
-          DO iky = ikys,ikye
-            buffer(iky,ikx) = field_(iky,ikx)
-          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_nky * Nkx , MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
-        ENDDO
-      ELSE
-        ! Recieve buffer from root
-        CALL MPI_RECV(buffer, local_nky * Nkx , MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
-        ! Write it in phi
-        DO ikx = ikxs,ikxe
-          DO iky = ikys,ikye
-            field_(iky,ikx) = buffer(iky,ikx)
-          ENDDO
-        ENDDO
-      ENDIF
-    ENDIF
-  END SUBROUTINE manual_2D_bcast
-
-  !!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
-SUBROUTINE manual_3D_bcast(field_)
-  USE grid
-  IMPLICIT NONE
-  COMPLEX(dp), INTENT(INOUT) :: field_(ikys:ikye,ikxs:ikxe,izs:ize)
-  COMPLEX(dp) :: buffer(ikys:ikye,ikxs:ikxe,izs:ize)
-  INTEGER     :: i_, root, world_rank, world_size, count
-  root = 0;
-  count = (ikye-ikys+1) * (ikxe-ikxs+1) * (ize-izs+1);
-
-  CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
-  CALL MPI_COMM_SIZE(comm_p,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 iz = izs,ize
-        DO ikx = ikxs,ikxe
-          DO iky = ikys,ikye
-              buffer(iky,ikx,iz) = field_(iky,ikx,iz)
-            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, count, MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
-      ENDDO
-    ELSE
-      ! Recieve buffer from root
-      CALL MPI_RECV(buffer, count, MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
-      ! Write it in phi
-      DO iz = izs,ize
-        DO ikx = ikxs,ikxe
-          DO iky = ikys,ikye
-            field_(iky,ikx,iz) = buffer(iky,ikx,iz)
-          ENDDO
-        ENDDO
-      ENDDO
-    ENDIF
-  ENDIF
-END SUBROUTINE manual_3D_bcast
-
-!!!!! This is a manual way to do MPI_BCAST !!!!!!!!!!!
-SUBROUTINE manual_0D_bcast(v)
-USE grid
-IMPLICIT NONE
-COMPLEX(dp), INTENT(INOUT) :: v
-COMPLEX(dp) :: buffer
-INTEGER     :: i_, root, world_rank, world_size, count
-root = 0;
-count = 1;
-
-CALL MPI_COMM_RANK(comm_z,world_rank,ierr)
-CALL MPI_COMM_SIZE(comm_z,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
-    buffer = v
-    ! Send it to all the other processes
-    DO i_ = 0,num_procs_z-1
-      IF (i_ .NE. world_rank) &
-      CALL MPI_SEND(buffer, count, MPI_DOUBLE_COMPLEX, i_, 0, comm_z, ierr)
-    ENDDO
-  ELSE
-    ! Recieve buffer from root
-    CALL MPI_RECV(buffer, count, MPI_DOUBLE_COMPLEX, root, 0, comm_z, MPI_STATUS_IGNORE, ierr)
-    ! Write it in phi
-    v = buffer
-  ENDIF
-ENDIF
-END SUBROUTINE manual_0D_bcast
-
-
 END MODULE utility
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index 1d0382518b1aac8bd7a56a7129af36367271a619..f2fb0d33c4e9aa67b9f752c4b613420d8564f3e7 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -10,9 +10,9 @@ system(['mkdir -p ',MISCDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20;
+JOBNUMMIN = 10; JOBNUMMAX = 20;
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
-data.FIGDIR = LOCALDIR;
+data.localdir = LOCALDIR;
 
 %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 default_plots_options
@@ -21,7 +21,7 @@ FMT = '.fig';
 
 if 0
 %% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0   = 0.8*data.Ts3D(end); 
+options.TAVG_0   = 0.98*data.Ts3D(end);
 options.TAVG_1   = data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
 % options.ST_FIELD = '\Gamma_x';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -48,12 +48,12 @@ options.NAME      = '\phi';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'yz';
+options.PLAN      = 'xz';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 % options.TIME      = dat.Ts5D;
-options.TIME      = 700:1:960;
+options.TIME      = 1250:1:1500;
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -71,11 +71,11 @@ options.NAME      = '\phi';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'yz';
+options.PLAN      = 'kxky';
 % options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
-options.COMP      = 1;
-options.TIME      =[800 900 950];
+options.COMP      = 'avg';
+options.TIME      = [1200 1300 1400 1500];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 save_figure(data,fig)
@@ -97,8 +97,8 @@ if 0
 % options.XPERP     = linspace( 0,6,64);
 options.SPAR      = gene_data.vp';
 options.XPERP     = gene_data.mu';
-options.iz        = 'avg';
-options.T         = 1000;
+options.iz        = 9;
+options.T         = 30;
 options.PLT_FCT   = 'contour';
 options.ONED      = 0;
 options.non_adiab = 1;
@@ -111,12 +111,12 @@ end
 if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
-options.P2J        = 0;
-options.ST         = 1;
+options.P2J        = 1;
+options.ST         = 0;
 options.PLOT_TYPE  = 'space-time';
 options.NORMALIZED = 1;
 options.JOBNUM     = 0;
-options.TIME       = [900 950];
+options.TIME       = [1300 1500];
 options.specie     = 'i';
 options.compz      = 'avg';
 fig = show_moments_spectrum(data,options);
@@ -162,7 +162,7 @@ if 0
 %% Mode evolution
 options.NORMALIZED = 1;
 options.K2PLOT = 1;
-options.TIME   = [0.5 1]*data.Ts3D(end);
+options.TIME   = 5:30;
 options.NMA    = 1;
 options.NMODES = 15;
 options.iz     = 9;
@@ -211,4 +211,4 @@ options.PLT_MTOPO = 1;
 data.rho_o_R      = 2e-3; % Sound larmor radius over Machine size ratio
 fig = show_geometry(data,options);
 save_figure(data,fig)
-end
\ No newline at end of file
+end
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 0fe05469c77ecd4444c4c8922bc0ddaa687fc4d8..281885dec2b0d23f1c5263e81273452e2035fa17 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -1,8 +1,7 @@
 % folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
-folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.2/';
+folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/linear_s_alpha_CBC_100/';
-% folder = '/misc/gene_results/shearless_cyclone/allmodes_CBC_100/';
 % folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
 % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
 % folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
@@ -12,7 +11,7 @@ gene_data = load_gene_data(folder);
 gene_data = invert_kxky_to_kykx_gene_results(gene_data);
 if 1
 %% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0   = 0.8*gene_data.Ts3D(end); 
+options.TAVG_0   = 0.8*gene_data.Ts3D(end);
 options.TAVG_1   = gene_data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
 options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
@@ -46,14 +45,14 @@ end
 if 0
 %% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Options
-options.INTERP    = 0;
+options.INTERP    = 1;
 options.POLARPLOT = 0;
 options.NAME      = '\phi';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'xy';
+options.PLAN      = 'xz';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
@@ -73,7 +72,7 @@ subplot(311)
     plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
     end
     xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry');
-    
+
 subplot(312)
     for i = 7:10
     plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
diff --git a/wk/debug_script.m b/wk/debug_script.m
new file mode 100644
index 0000000000000000000000000000000000000000..1d53f6cb1691fe60287931ae0b3850490a570851
--- /dev/null
+++ b/wk/debug_script.m
@@ -0,0 +1,30 @@
+% system('cd ../results/dev/test_diag; ./helaz3_dbg; cd $HOME/HeLaZ/wk');
+% system('cd ../results/dev/test_diag; mpirun -np 2 ./helaz3_dbg 1 2 1; cd $HOME/HeLaZ/wk');
+system('cd ../results/dev/test_diag; mpirun -np 8 ./helaz3_dbg 2 2 2; cd $HOME/HeLaZ/wk');
+
+
+filename = '../results/dev/test_diag/outputs_00.h5';
+
+% test phi
+f_        = h5read(filename,'/data/var3d/phi/000004/');
+f_gatherv = h5read(filename,'/data/var3d/phi_gatherv/000004/');
+f_        = f_.real        + 1i*f_.imaginary;
+f_gatherv = f_gatherv.real + 1i*f_gatherv.imaginary;
+err = sum(sum(sum(abs(f_-f_gatherv))))/sum(sum(sum(abs(f_))));
+disp(['error_phi = ',sprintf('%2.2e',err)]);    
+
+% test Ni00
+f_        = h5read(filename,'/data/var3d/Ni00/000004/');
+f_gatherv = h5read(filename,'/data/var3d/Ni00_gatherv/000004/');
+f_        = f_.real        + 1i*f_.imaginary;
+f_gatherv = f_gatherv.real + 1i*f_gatherv.imaginary;
+err = sum(sum(sum(abs(f_-f_gatherv))))/sum(sum(sum(abs(f_))));
+disp(['error_Ni00 = ',sprintf('%2.2e',err)]);    
+
+% test Nipj
+f_        = h5read(filename,'/data/var5d/moments_i/000004/');
+f_gatherv = h5read(filename,'/data/var5d/moments_i_gatherv/000004/');
+f_        = f_.real        + 1i*f_.imaginary;
+f_gatherv = f_gatherv.real + 1i*f_gatherv.imaginary;
+err = sum(sum(sum(sum(sum(abs(f_-f_gatherv))))))/sum(sum(sum(sum(sum(abs(f_))))));
+disp(['error_Nipj = ',sprintf('%2.2e',err)]);    
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 113d91144a967916ca693fcd1b9683d5cbaee6b2..05a07248e807018a60b0b34822d654440f0369b5 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -5,9 +5,9 @@ helazdir = '/home/ahoffman/HeLaZ/';
 outfile ='';
 outfile ='';
 outfile ='';
-% outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
+outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
 % outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
-outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0';
+% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0';
 % outfile ='quick_run/CLOS_1_64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
 % outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1';
 % outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK';
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 8d437c867bd14636a22241ce8e5604aa2b1e381d..8165d782312639b7731bca73b44d2b68bb8ad44e 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -113,7 +113,7 @@ end
 
 if 1
 %% Ballooning plot
-options.time_2_plot = data.Ts3D(end);
+options.time_2_plot = [0.9 1]*data.Ts3D(end);
 options.kymodes     = [0.5];
 options.normalized  = 1;
 options.sheared     = 0;