From a65d383b011f73d0c32d95c1d80fe4f398fe007a Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Thu, 19 May 2022 13:21:59 +0200
Subject: [PATCH] tentative of a single process output version of the code

---
 Makefile                         | 112 +++++++++++++++-------
 src/array_mod.F90                |   9 +-
 src/auxval.F90                   |   3 +
 src/basic_mod.F90                |   7 ++
 src/calculus_mod.F90             |   2 +-
 src/diag_headers/diagnose_full.h |  19 +++-
 src/diagnose.F90                 |   6 --
 src/grid_mod.F90                 |  14 ++-
 src/memory.F90                   |   1 +
 src/parallel_mod.F90             | 155 +++++++++++++++++++++++++++++++
 src/poisson.F90                  |   2 +-
 src/ppinit.F90                   |  21 ++++-
 src/processing_mod.F90           |  15 ---
 src/utility_mod.F90              | 118 +----------------------
 wk/analysis_3D.m                 |   8 +-
 wk/analysis_header.m             |   2 +-
 16 files changed, 304 insertions(+), 190 deletions(-)
 create mode 100644 src/parallel_mod.F90

diff --git a/Makefile b/Makefile
index 382dd631..80ae8671 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/src/array_mod.F90 b/src/array_mod.F90
index 96485a88..0b0db612 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -63,7 +63,14 @@ MODULE array
   ! Poisson operator (ikx,iky,iz)
   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_poisson_op
 
-  !! Diagnostics
+  !! Diagnostics (full arrays to gather on process 0 for output)
+  ! moments
+  COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_e_full
+  COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_i_full
+  
+  ! ES potential
+  COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: phi_full
+
   ! 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 db695e67..7a9d9852 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -9,6 +9,7 @@ subroutine auxval
   use prec_const
   USE numerics
   USE geometry
+  USE parallel, ONLY: init_parallel_var
   IMPLICIT NONE
 
   INTEGER :: irows,irowe, irow, icol, i_
@@ -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
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index f1d67f07..786901f2 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)
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index 76e22a86..8cb9c02b 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, &
diff --git a/src/diag_headers/diagnose_full.h b/src/diag_headers/diagnose_full.h
index 5f9cc79d..93b45ca0 100644
--- a/src/diag_headers/diagnose_full.h
+++ b/src/diag_headers/diagnose_full.h
@@ -8,11 +8,12 @@ USE model
 USE initial_par
 USE fields
 USE time_integration
-USE utility
+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
@@ -106,6 +107,9 @@ IF ((kstep .EQ. 0)) THEN
    CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
 
    IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
+   IF (write_phi) THEN
+    CALL creatg(fidres, "/data/var3d/phi_gatherv", "phi_gatherv")
+   ENDIF
 
    IF (write_Na00) THEN
     IF(KIN_E)&
@@ -280,12 +284,13 @@ SUBROUTINE diagnose_3d
   USE diagnostics_par
   USE prec_const
   USE processing
+  USE parallel, ONLY : gather_xyz
   USE model, ONLY: KIN_E
 
   IMPLICIT NONE
 
-  INTEGER     :: i_, root, world_rank, world_size
-
+  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)
@@ -293,6 +298,12 @@ SUBROUTINE diagnose_3d
   CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
 
   IF (write_phi) CALL write_field3d_kykxz(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
+  IF (write_phi) THEN
+    CALL gather_xyz(phi(ikys:ikye,1:Nkx,izs:ize),phi_full(1:Nky,1:Nkx,1:Nz))
+    WRITE(dset_name, "(A, '/', i6.6)") "/data/var3d/phi_gatherv", iframe3d
+    CALL putarr(fidres, dset_name, phi_full(1:Nky,1:Nkx,1:Nz), ionode=0)
+    CALL attach(fidres, dset_name, "time", time)
+  ENDIF
 
   IF (write_Na00) THEN
     IF(KIN_E)THEN
@@ -349,7 +360,7 @@ SUBROUTINE diagnose_3d
     IMPLICIT NONE
     COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: field
     CHARACTER(*), INTENT(IN) :: text
-    CHARACTER(LEN=50) :: dset_name
+    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)
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index c1b81f25..4296d0eb 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -108,9 +108,3 @@ 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'
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 36585c36..c4d29cb3 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -48,10 +48,12 @@ 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
@@ -275,7 +277,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 +299,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
diff --git a/src/memory.F90 b/src/memory.F90
index a3837b7c..6806e644 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -14,6 +14,7 @@ SUBROUTINE memory
 
   ! Electrostatic potential
   CALL allocate_array(           phi, ikys,ikye, ikxs,ikxe, izgs,izge)
+  CALL allocate_array(      phi_full,     1,Nky,     1,Nkx,      1,Nz)
   CALL allocate_array(        phi_ZF, ikxs,ikxe, izs,ize)
   CALL allocate_array(        phi_EM, ikys,ikye, izs,ize)
   CALL allocate_array(inv_poisson_op, ikys,ikye, ikxs,ikxe, izs,ize)
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
new file mode 100644
index 00000000..1598bf33
--- /dev/null
+++ b/src/parallel_mod.F90
@@ -0,0 +1,155 @@
+MODULE parallel
+  USE basic
+  USE grid, ONLY: Nkx,Nky,Nz, ikys,ikye, izs,ize, local_nky, local_nz
+  use prec_const
+  IMPLICIT NONE
+  ! recieve and displacement counts for gatherv
+  INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_z, rcv_y, dsp_z, dsp_y
+
+  PUBLIC :: manual_0D_bcast, manual_3D_bcast, init_parallel_var, &
+            gather_xyz!, gather_pjz, gather_pjxyz
+
+CONTAINS
+
+  SUBROUTINE init_parallel_var
+    INTEGER :: i_
+    !! 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
+    print*, rcv_y, dsp_y
+    !! Z reduction for full slices of y data but constant x
+    ! number of points recieved and displacement for the z reduction
+    ALLOCATE(rcv_z(0:num_procs_z-1),dsp_z(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_z,1,MPI_INTEGER,comm_z,ierr)
+    ! the displacement array can be build from each local_np as
+    dsp_z(0)=0
+    DO i_=1,num_procs_z-1
+       dsp_z(i_) =dsp_z(i_-1) + rcv_z(i_-1)
+    END DO
+    print*, rcv_z, dsp_z
+  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 :: send_count, root_p, root_z, root_ky, ix, iz
+
+    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
+          print*,iz
+          buffer_ly_cz(ikys:ikye) = field_sub(ikys:ikye,ix,iz)
+          CALL MPI_GATHERV(buffer_ly_cz,local_nky,MPI_DOUBLE_COMPLEX, &
+                           buffer_fy_cz, rcv_y, dsp_y, MPI_DOUBLE_COMPLEX, &
+                           root_ky, comm_ky)
+          buffer_fy_lz(1:Nky,iz) = buffer_fy_cz(1:Nky)
+        ENDDO
+
+        ! send the full line on y contained by root_ky
+        IF(rank_ky .EQ. 0) THEN
+          CALL MPI_GATHERV(buffer_fy_lz,Nky*local_nz,MPI_DOUBLE_COMPLEX, &
+                           buffer_fy_fz, rcv_z, dsp_z, MPI_DOUBLE_COMPLEX, &
+                           root_z, comm_z)
+        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 !!!!!
+
+  !!!!! Gather a field in kinetic + spatial coordinates on rank 0 !!!!!
+
+  !!!!! 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 3936324e..75e2b887 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 55db1460..f103cdd9 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)
   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 5a47ed6f..d30b795c 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/utility_mod.F90 b/src/utility_mod.F90
index b6e989a7..b962a4c4 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_3D.m b/wk/analysis_3D.m
index e9ab6b17..84d0bd3a 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -99,9 +99,9 @@ if 0
 % options.XPERP     = linspace( 0,6,64);
 options.SPAR      = gene_data.vp';
 options.XPERP     = gene_data.mu';
-options.iz        = 13;
-options.T         = 20;
-options.PLT_FCT   = 'pcolor';
+options.iz        = 9;
+options.T         = 30;
+options.PLT_FCT   = 'contour';
 options.ONED      = 0;
 options.non_adiab = 1;
 options.SPECIE    = 'i';
@@ -164,7 +164,7 @@ if 0
 %% Mode evolution
 options.NORMALIZED = 1;
 options.K2PLOT = 1;
-options.TIME   = 2:20;
+options.TIME   = 5:30;
 options.NMA    = 1;
 options.NMODES = 15;
 options.iz     = 9;
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index 3ff76bdc..28d53d9f 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -7,7 +7,7 @@ outfile ='';
 outfile ='';
 % outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
 % outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
-outfile ='shearless_cyclone/linear_CBC_120';
+outfile ='shearless_cyclone/linear_CBC_100_D_16';
 % 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';
-- 
GitLab