From fafef78fc6581dfa131af0007babd21e8c744604 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Thu, 12 May 2022 17:33:54 +0200
Subject: [PATCH] First commit of a version that can create multiple outputs

---
 matlab/compile_results.m                |  14 +-
 matlab/compute/compute_fa.m             |  10 +-
 matlab/plot/show_napjz.m                |  10 +-
 matlab/process_field.m                  |  14 +-
 src/basic_mod.F90                       |  29 +-
 src/diag_headers/diagnose_fields.h      | 157 ++++++
 src/diag_headers/diagnose_full.h        | 478 ++++++++++++++++
 src/diag_headers/diagnose_gridgeom.h    |  50 ++
 src/diag_headers/diagnose_moments.h     | 136 +++++
 src/diag_headers/diagnose_momspectrum.h | 105 ++++
 src/diag_headers/diagnose_profiler.h    |  54 ++
 src/diag_headers/diagnose_timetraces.h  |  76 +++
 src/diagnose.F90                        | 696 +++---------------------
 src/diagnostics_par_mod.F90             |  43 +-
 src/nonlinear_mod.F90                   |  10 +-
 src/ppinit.F90                          |   4 +
 wk/analysis_3D.m                        |  28 +-
 wk/analysis_header.m                    |   3 +-
 wk/gene_analysis_3D.m                   |   4 +-
 wk/quick_run.m                          |  48 +-
 20 files changed, 1282 insertions(+), 687 deletions(-)
 create mode 100644 src/diag_headers/diagnose_fields.h
 create mode 100644 src/diag_headers/diagnose_full.h
 create mode 100644 src/diag_headers/diagnose_gridgeom.h
 create mode 100644 src/diag_headers/diagnose_moments.h
 create mode 100644 src/diag_headers/diagnose_momspectrum.h
 create mode 100644 src/diag_headers/diagnose_profiler.h
 create mode 100644 src/diag_headers/diagnose_timetraces.h

diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 50d383b3..d0892abf 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -71,18 +71,18 @@ while(CONTINUE)
             if W_NAPJ
                 tmp = Nipj_; sz = size(tmp);
                 Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
-                Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
+                Nipj_(1:Pi_old,1:Ji_old,:,:,:,:) = tmp;
                 tmp = Nepj_; sz = size(tmp);
-                Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
-                Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
+%                 Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
+%                 Nepj_(1:Pe_old,1:Je_old,:,:,:,:) = tmp;
             end
             if W_SAPJ
                 tmp = Sipj_; sz = size(tmp);
                 Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
-                Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
+                Sipj_(1:Pi_old,1:Ji_old,:,:,:,:) = tmp;
                 tmp = Sepj_; sz = size(tmp);
-                Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
-                Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
+%                 Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
+%                 Sepj_(1:Pe_old,1:Je_old,:,:,:,:) = tmp;
             end
         % If a degree is smaller than previous job, put zero to add. deg.
         elseif (sum([Pe_new, Je_new, Pi_new, Ji_new]<[Pe_old, Je_old, Pi_old, Ji_old]) >= 1 && Pe_old ~= 1e9)
@@ -258,7 +258,7 @@ DATA.ikx0 = ikx0; DATA.iky0 = iky0;
 DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky; 
 DATA.Pmaxe = numel(Pe); DATA.Pmaxi = numel(Pi); DATA.Jmaxe = numel(Je); DATA.Jmaxi = numel(Ji);
 DATA.dir      = DIRECTORY;
-DATA.localdir = ['..',DIRECTORY(20:end)];
+DATA.localdir = DIRECTORY;
 DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ...
     ', $\kappa_N=$',num2str(DATA.K_N),', $L=',num2str(DATA.L),'$, $N=',...
     num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',...
diff --git a/matlab/compute/compute_fa.m b/matlab/compute/compute_fa.m
index 3d8e4ec4..a580a2ee 100644
--- a/matlab/compute/compute_fa.m
+++ b/matlab/compute/compute_fa.m
@@ -46,7 +46,7 @@ Napj_ = squeeze(Napj_);
 
 Np = numel(parray); Nj = numel(jarray);
 
-FF = zeros(data.Nkx,data.Nky,numel(options.XPERP),numel(options.SPAR));
+FF = zeros(data.Nky,data.Nkx,numel(options.XPERP),numel(options.SPAR));
 % FF = zeros(numel(options.XPERP),numel(options.SPAR));
 
 FAM = FaM(SS,XX);
@@ -60,14 +60,14 @@ for ip_ = 1:Np
         HLF = HH.*LL.*FAM;
         for ikx = 1:data.Nkx
             for iky = 1:data.Nky
-                FF(ikx,iky,:,:) = squeeze(FF(ikx,iky,:,:)) + Napj_(ip_,ij_,ikx,iky)*HLF;
+                FF(iky,ikx,:,:) = squeeze(FF(iky,ikx,:,:)) + Napj_(ip_,ij_,iky,ikx)*HLF;
             end
         end
     end
 end
-% FF = (FF.*conj(FF)); %|f_a|^2
-FF = abs(FF); %|f_a|
-FF = sqrt(squeeze(mean(mean(FF,1),2))); %sqrt(<|f_a|>kx,ky)
+FF = (FF.*conj(FF)); %|f_a|^2
+% FF = abs(FF); %|f_a|
+FF = sqrt(squeeze(mean(mean(FF,1),2))); %sqrt(<|f_a|^2>kx,ky)
 FF = FF./max(max(FF));
 FF = FF';
 % FF = FF.*conj(FF);
diff --git a/matlab/plot/show_napjz.m b/matlab/plot/show_napjz.m
index e521b944..002bc2eb 100644
--- a/matlab/plot/show_napjz.m
+++ b/matlab/plot/show_napjz.m
@@ -56,15 +56,7 @@ switch OPTIONS.PLOT_TYPE
     imagesc(DATA.Ts5D,p2ja,plt(Na_ST,1:numel(p2ja))); 
     set(gca,'YDir','normal')        
     xlabel('$t$'); ylabel('$p+2j$')
-    title('$\langle\sum_k |',name,'|\rangle_{p+2j=const}$')
-    if DATA.K_E
-    subplot(2,1,2)
-        imagesc(DATA.Ts5D,p2je,plt(Ne_ST,1:numel(p2ja))); 
-        set(gca,'YDir','normal')
-    xlabel('$t$'); ylabel('$p+2j$')
-    title('$\langle\sum_k |N_e^{pj}|\rangle_{p+2j=const}$')
-    suptitle(DATA.param_title);
-    end
+    title(['$\langle\sum_k |',name,'|^2\rangle_{p+2j=const}$'])
     
     case 'Tavg-1D'
     t0 = OPTIONS.TIME(1);
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 57096536..201d0a49 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -43,6 +43,10 @@ switch OPTIONS.PLAN
         XNAME = '$v_\parallel$'; YNAME = '$\mu$';
         [Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR);
         REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 0;
+        for i = 1:numel(OPTIONS.TIME)
+            [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts5D));
+        end
+        FRAMES = unique(FRAMES);
 end
 Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
 Lmin_ = min([Xmax_,Ymax_]);
@@ -440,13 +444,9 @@ switch OPTIONS.NAME
     case 'f_i'
         shift_x = @(x) x; shift_y = @(y) y;
         NAME = 'fi'; OPTIONS.SPECIE = 'i';
-        [~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D));
-        [~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D));
-        dit_ = 1;%ceil((it1_-it0_+1)/10); 
-        FRAMES = it0_:dit_:it1_;
-        iz = 1;
-        for it = 1:numel(FRAMES)
-            OPTIONS.T = DATA.Ts5D(FRAMES(it));
+        for it = 1:numel(OPTIONS.TIME)
+            [~,it0_] =min(abs(OPTIONS.TIME(it)-DATA.Ts5D));
+            OPTIONS.T = DATA.Ts5D(it0_);
             [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
         end  
     case 'f_e'
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index ab7869f8..f1d67f07 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -16,7 +16,8 @@ MODULE basic
   ! Communicators for 1-dim cartesian subgrids of comm0
   INTEGER :: comm_p, comm_ky, comm_z
   INTEGER :: rank_p, rank_ky, rank_z! Ranks
-  INTEGER :: commr_p0              ! Communicators along kx for only rank 0 on p
+  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 :: jobnum  = 0           ! Job number
   INTEGER :: step    = 0           ! Calculation step of this run
@@ -80,6 +81,32 @@ CONTAINS
     tc_Sapj      = 0.; tc_diag     = 0.; tc_checkfield = 0.
 
   END SUBROUTINE basic_data
+
+
+  SUBROUTINE basic_outputinputs(fid, str)
+    !
+    !    Write the input parameters to the results_xx.h5 file
+    !
+    USE prec_const
+    USE futils, ONLY: attach
+    IMPLICIT NONE
+    INTEGER, INTENT(in) :: fid
+    CHARACTER(len=256), INTENT(in) :: str
+    CALL attach(fid, TRIM(str), "start_iframe0d", iframe0d)
+    CALL attach(fid, TRIM(str), "start_iframe2d", iframe2d)
+    CALL attach(fid, TRIM(str), "start_iframe3d", iframe3d)
+    CALL attach(fid, TRIM(str), "start_iframe5d", iframe5d)
+    CALL attach(fid, TRIM(str),  "start_time",     time)
+    CALL attach(fid, TRIM(str), "start_cstep",    cstep-1)
+    CALL attach(fid, TRIM(str),          "dt",       dt)
+    CALL attach(fid, TRIM(str),        "tmax",     tmax)
+    CALL attach(fid, TRIM(str),        "nrun",     nrun)
+    CALL attach(fid, TRIM(str),    "cpu_time",       -1)
+    CALL attach(fid, TRIM(str),       "Nproc",   num_procs)
+    CALL attach(fid, TRIM(str),       "Np_p" , num_procs_p)
+    CALL attach(fid, TRIM(str),       "Np_kx",num_procs_ky)
+    CALL attach(fid, TRIM(str),        "Np_z", num_procs_z)
+  END SUBROUTINE basic_outputinputs
   !================================================================================
   SUBROUTINE find_input_file
     IMPLICIT NONE
diff --git a/src/diag_headers/diagnose_fields.h b/src/diag_headers/diagnose_fields.h
new file mode 100644
index 00000000..a68cb8fb
--- /dev/null
+++ b/src/diag_headers/diagnose_fields.h
@@ -0,0 +1,157 @@
+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
new file mode 100644
index 00000000..5f9cc79d
--- /dev/null
+++ b/src/diag_headers/diagnose_full.h
@@ -0,0 +1,478 @@
+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
new file mode 100644
index 00000000..9ef1855c
--- /dev/null
+++ b/src/diag_headers/diagnose_gridgeom.h
@@ -0,0 +1,50 @@
+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
new file mode 100644
index 00000000..58582260
--- /dev/null
+++ b/src/diag_headers/diagnose_moments.h
@@ -0,0 +1,136 @@
+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
new file mode 100644
index 00000000..36befdda
--- /dev/null
+++ b/src/diag_headers/diagnose_momspectrum.h
@@ -0,0 +1,105 @@
+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
new file mode 100644
index 00000000..b99f67dd
--- /dev/null
+++ b/src/diag_headers/diagnose_profiler.h
@@ -0,0 +1,54 @@
+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
new file mode 100644
index 00000000..3c688aee
--- /dev/null
+++ b/src/diag_headers/diagnose_timetraces.h
@@ -0,0 +1,76 @@
+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 ced695af..894263d0 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -2,635 +2,115 @@ SUBROUTINE diagnose(kstep)
   !   Diagnostics, writing simulation state to disk
 
   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
 
-  INCLUDE 'srcinfo.h'
-
   INTEGER, INTENT(in) :: kstep
-  INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank = 0
-  INTEGER :: dims(1) = (/0/)
-  INTEGER :: cp_counter = 0
-  CHARACTER(len=256) :: str, fname,test_
 
   CALL cpu_time(t0_diag) ! Measuring time
-  !_____________________________________________________________________________
-  !                   1.   Initial diagnostics
 
+  !! Basic diagnose loop for reading input file, displaying advancement and ending
   IF ((kstep .EQ. 0)) THEN
-    ! Writing output filename
-     WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
-     !                      1.1   Initial run
-     ! Main output file creation
-     IF (write_doubleprecision) THEN
-        CALL creatf(resfile, fidres, real_prec='d', mpicomm=comm0)
-     ELSE
-        CALL creatf(resfile, fidres, mpicomm=comm0)
-     END IF
-     IF (my_id .EQ. 0) WRITE(*,'(3x,a,a)') TRIM(resfile), ' created'
-
-     !  Data group
-     CALL creatg(fidres, "/data", "data")
-
-     !  File group
-     CALL creatg(fidres, "/files", "files")
-     CALL attach(fidres, "/files",  "jobnum", jobnum)
-
-     ! 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
-
-     !  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
-
-     WRITE(str,'(a,i2.2)') "/data/input"
-     CALL creatd(fidres, rank,dims,TRIM(str),'Input parameters')
-     CALL attach(fidres, TRIM(str),     "version",  VERSION) !defined in srcinfo.h
-     CALL attach(fidres, TRIM(str),      "branch",   BRANCH) !defined in srcinfo.h
-     CALL attach(fidres, TRIM(str),      "author",   AUTHOR) !defined in srcinfo.h
-     CALL attach(fidres, TRIM(str),    "execdate", EXECDATE) !defined in srcinfo.h
-     CALL attach(fidres, TRIM(str),        "host",     HOST) !defined in srcinfo.h
-     CALL attach(fidres, TRIM(str),  "start_time",     time)
-     CALL attach(fidres, TRIM(str), "start_cstep",    cstep-1)
-     CALL attach(fidres, TRIM(str), "start_iframe0d", iframe0d)
-     CALL attach(fidres, TRIM(str), "start_iframe2d", iframe2d)
-     CALL attach(fidres, TRIM(str), "start_iframe3d", iframe3d)
-     CALL attach(fidres, TRIM(str), "start_iframe5d", iframe5d)
-     CALL attach(fidres, TRIM(str),          "dt",       dt)
-     CALL attach(fidres, TRIM(str),        "tmax",     tmax)
-     CALL attach(fidres, TRIM(str),        "nrun",     nrun)
-     CALL attach(fidres, TRIM(str),    "cpu_time",       -1)
-     CALL attach(fidres, TRIM(str),       "Nproc",   num_procs)
-     CALL attach(fidres, TRIM(str),       "Np_p" , num_procs_p)
-     CALL attach(fidres, TRIM(str),       "Np_kx",num_procs_ky)
-     CALL attach(fidres, TRIM(str), "write_gamma", write_gamma)
-     CALL attach(fidres, TRIM(str),    "write_hf",    write_hf)
-     CALL attach(fidres, TRIM(str),   "write_phi",   write_phi)
-     CALL attach(fidres, TRIM(str),  "write_Na00",  write_Na00)
-     CALL attach(fidres, TRIM(str),  "write_Napj",  write_Napj)
-     CALL attach(fidres, TRIM(str),  "write_Sapj",  write_Sapj)
-     CALL attach(fidres, TRIM(str),  "write_dens",  write_dens)
-     CALL attach(fidres, TRIM(str),  "write_fvel",  write_fvel)
-     CALL attach(fidres, TRIM(str),  "write_temp",  write_temp)
-
-     CALL grid_outputinputs(fidres, str)
-     CALL geometry_outputinputs(fidres, str)
-     CALL diag_par_outputinputs(fidres, str)
-     CALL model_outputinputs(fidres, str)
-     CALL coll_outputinputs(fidres, str)
-     CALL initial_outputinputs(fidres, str)
-     CALL time_integration_outputinputs(fidres, str)
-
-
-     !  Save STDIN (input file) of this run
-     IF(jobnum .LE. 99) THEN
-        WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
-     ELSE
-        WRITE(str,'(a,i3.2)') "/files/STDIN.",jobnum
-     END IF
-
-     INQUIRE(unit=lu_in, name=fname)
-     CLOSE(lu_in)
-
-     CALL putfile(fidres, TRIM(str), TRIM(fname),ionode=0)
-
-   ELSEIF((kstep .EQ. 0)) THEN
-
-      IF(jobnum .LE. 99) THEN
-         WRITE(resfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
-      ELSE
-         WRITE(resfile,'(a,a1,i3.2,a3)') TRIM(resfile0),'_',jobnum,'.h5'
-      END IF
-
-      CALL openf(resfile,fidres, 'D');
-
-   ENDIF
-
-  !_____________________________________________________________________________
-  !                   2.   Periodic diagnostics
-  !
+    INQUIRE(unit=lu_in, name=input_fname)
+    CLOSE(lu_in)
+  ENDIF
   IF (kstep .GE. 0) THEN
-
     ! Terminal info
     IF (MOD(cstep, INT(1.0/dt)) == 0 .AND. (my_id .EQ. 0)) THEN
      WRITE(*,"(F6.0,A,F6.0)") time,"/",tmax
     ENDIF
-
-     !                       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 cpu_time(finish)
-     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
+    CALL cpu_time(finish)
      ! Display computational time cost
      IF (my_id .EQ. 0) CALL display_h_min_s(finish-start)
-
-     !   Close all diagnostic files
-     CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-     CALL closef(fidres)
-
   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)
 
-END SUBROUTINE diagnose
-
-
-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_2d
-
-  USE basic
-  USE futils, ONLY: append, getatt, attach, putarrnd
-  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
-
-  COMPLEX(dp) :: buffer(ikys:ikye,ikxs:ikxe)
-  INTEGER     :: i_, root, world_rank, world_size
-
-  CALL append(fidres,  "/data/var2d/time",           time,ionode=0)
-  CALL append(fidres, "/data/var2d/cstep", real(cstep,dp),ionode=0)
-  CALL getatt(fidres,      "/data/var2d/",       "frames",iframe2d)
-  iframe2d=iframe2d+1
-  CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
-
-CONTAINS
-
-  SUBROUTINE write_field2d(field, text)
-    USE futils, ONLY: attach, putarr
-    USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx
-    USE prec_const
-    USE basic, ONLY : comm_ky, num_procs_p, rank_p
-    IMPLICIT NONE
-
-    COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe), INTENT(IN) :: field
-    CHARACTER(*), INTENT(IN) :: text
-    COMPLEX(dp) :: buffer_dist(ikys:ikye,ikxs:ikxe)
-    COMPLEX(dp) :: buffer_full(1:Nky,1:Nkx)
-    INTEGER     :: scount, rcount
-    CHARACTER(LEN=50) :: dset_name
-
-    scount = (ikye-ikys+1) * (ikxe-ikxs+1)
-    rcount = scount
-
-    WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var2d", TRIM(text), iframe2d
-    IF (num_procs .EQ. 1) THEN ! no data distribution
-      CALL putarr(fidres, dset_name, field(ikys:ikye,ikxs:ikxe), ionode=0)
-    ELSE
-      CALL putarrnd(fidres, dset_name, field(ikys:ikye,ikxs:ikxe),  (/1, 1/))
-    ENDIF
-    CALL attach(fidres, dset_name, "time", time)
-
-  END SUBROUTINE write_field2d
-
-END SUBROUTINE diagnose_2d
-
-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
+END SUBROUTINE diagnose
 
+SUBROUTINE init_outfile(comm,file0,file,fid)
+  USE diagnostics_par, ONLY : write_doubleprecision, diag_par_outputinputs, input_fname
+  USE basic,           ONLY : my_id, jobnum, lu_in, basic_outputinputs
+  USE grid,            ONLY : grid_outputinputs
+  USE geometry,        ONLY : geometry_outputinputs
+  USE model,           ONLY : model_outputinputs
+  USE collision,       ONLY : coll_outputinputs
+  USE initial_par,     ONLY : initial_outputinputs
+  USE time_integration,ONLY : time_integration_outputinputs
+  USE futils,          ONLY : creatf, creatg, creatd, attach, putfile
   IMPLICIT NONE
+  !input
+  INTEGER,            INTENT(IN)    :: comm
+  CHARACTER(len=256), INTENT(IN)    :: file0
+  CHARACTER(len=256), INTENT(OUT)   :: file
+  INTEGER,            INTENT(OUT)   :: fid
+  CHARACTER(len=256)                :: str,fname
+  INCLUDE 'srcinfo.h'
 
-  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
+  ! Writing output filename
+  WRITE(file,'(a,a1,i2.2,a3)') TRIM(file0)   ,'_',jobnum,'.h5'
+  !                      1.1   Initial run
+  ! Main output file creation
+  IF (write_doubleprecision) THEN
+    CALL creatf(file, fid, real_prec='d', mpicomm=comm)
+  ELSE
+    CALL creatf(file, fid, mpicomm=comm)
+  END IF
+  IF (my_id .EQ. 0) WRITE(*,'(3x,a,a)')  TRIM(file), ' created'
+  !  basic data group
+  CALL creatg(fid, "/data", "data")
+  !  File group
+  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')
+  CALL attach(fid, TRIM(str),     "version",  VERSION) !defined in srcinfo.h
+  CALL attach(fid, TRIM(str),      "branch",   BRANCH) !defined in srcinfo.h
+  CALL attach(fid, TRIM(str),      "author",   AUTHOR) !defined in srcinfo.h
+  CALL attach(fid, TRIM(str),    "execdate", EXECDATE) !defined in srcinfo.h
+  CALL attach(fid, TRIM(str),        "host",     HOST) !defined in srcinfo.h
+
+  CALL basic_outputinputs(fid,str)
+  CALL grid_outputinputs(fid, str)
+  CALL geometry_outputinputs(fid, str)
+  CALL diag_par_outputinputs(fid, str)
+  CALL model_outputinputs(fid, str)
+  CALL coll_outputinputs(fid, str)
+  CALL initial_outputinputs(fid, str)
+  CALL time_integration_outputinputs(fid, str)
+
+  !  Save STDIN (input file) of this run
+  IF(jobnum .LE. 99) THEN
+     WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
+  ELSE
+     WRITE(str,'(a,i3.2)') "/files/STDIN.",jobnum
+  END IF
+  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'
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index 59a7578f..9cbe257d 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -14,12 +14,19 @@ MODULE diagnostics_par
   INTEGER, PUBLIC, PROTECTED :: nsave_0d, nsave_1d, nsave_2d, nsave_3d, nsave_5d
 
   !  HDF5 file
-  CHARACTER(len=256), PUBLIC :: resfile0 = "outputs"   ! Head of main result file name
-  CHARACTER(len=256), PUBLIC :: resfile                ! Main result file
-  CHARACTER(len=256), PUBLIC :: rstfile                ! restart result file
-  INTEGER, PUBLIC            :: job2load               ! jobnum of the checkpoint to load
-  INTEGER, PUBLIC            :: fidres                 ! FID for resfile
-  INTEGER, PUBLIC            :: fidrst                 ! FID for restart file
+  CHARACTER(len=256), PUBLIC :: resfile,resfile0 = "outputs"            ! Head of main result file name
+  CHARACTER(len=256), PUBLIC :: momfile,momfile0 = "moments"   ! Head of the moment spectrum file (N_a(p,j,z))
+  CHARACTER(len=256), PUBLIC :: mspfile,mspfile0 = "moments_spectrum"   ! Head of the moment spectrum file (N_a(p,j,z))
+  CHARACTER(len=256), PUBLIC :: fldfile,fldfile0 = "fields"             ! Head of field (phi,A)
+  CHARACTER(len=256), PUBLIC :: ttrfile,ttrfile0 = "time_traces"        ! Head of time traces (gamma_x,Q_x)
+  CHARACTER(len=256), PUBLIC :: ggmfile,ggmfile0 = "grid_geometry"        ! Head of time traces (gamma_x,Q_x)
+  CHARACTER(len=256), PUBLIC :: prffile,prffile0 = "profiler"        ! Head of time traces (gamma_x,Q_x)
+  CHARACTER(len=256), PUBLIC :: input_fname
+  CHARACTER(len=256), PUBLIC :: rstfile                     ! restart result file
+  INTEGER, PUBLIC            :: job2load                    ! jobnum of the checkpoint to load
+  INTEGER, PUBLIC            :: fidres,fidmsp,fidfld,fidttr ! FID for output
+  INTEGER, PUBLIC            :: fidmom,fidggm, fidprf
+  INTEGER, PUBLIC            :: fidrst                      ! FID for restart file
 
   PUBLIC :: diag_par_readinputs, diag_par_outputinputs
 
@@ -44,21 +51,31 @@ CONTAINS
   END SUBROUTINE diag_par_readinputs
 
 
-  SUBROUTINE diag_par_outputinputs(fidres, str)
+  SUBROUTINE diag_par_outputinputs(fid, str)
     !
     !    Write the input parameters to the results_xx.h5 file
     !
     USE prec_const
     USE futils, ONLY: attach
     IMPLICIT NONE
-    INTEGER, INTENT(in) :: fidres
+    INTEGER, INTENT(in) :: fid
     CHARACTER(len=256), INTENT(in) :: str
 
-    CALL attach(fidres, TRIM(str), "write_doubleprecision", write_doubleprecision)
-    CALL attach(fidres, TRIM(str), "nsave_0d", nsave_0d)
-    CALL attach(fidres, TRIM(str), "nsave_1d", nsave_1d)
-    CALL attach(fidres, TRIM(str), "nsave_2d", nsave_2d)
-    CALL attach(fidres, TRIM(str), "nsave_5d", nsave_5d)
+    CALL attach(fid, TRIM(str), "write_doubleprecision", write_doubleprecision)
+    CALL attach(fid, TRIM(str), "nsave_0d", nsave_0d)
+    CALL attach(fid, TRIM(str), "nsave_1d", nsave_1d)
+    CALL attach(fid, TRIM(str), "nsave_2d", nsave_2d)
+    CALL attach(fid, TRIM(str), "nsave_3d", nsave_3d)
+    CALL attach(fid, TRIM(str), "nsave_5d", nsave_5d)
+    CALL attach(fid, TRIM(str), "write_gamma", write_gamma)
+    CALL attach(fid, TRIM(str),    "write_hf",    write_hf)
+    CALL attach(fid, TRIM(str),   "write_phi",   write_phi)
+    CALL attach(fid, TRIM(str),  "write_Na00",  write_Na00)
+    CALL attach(fid, TRIM(str),  "write_Napj",  write_Napj)
+    CALL attach(fid, TRIM(str),  "write_Sapj",  write_Sapj)
+    CALL attach(fid, TRIM(str),  "write_dens",  write_dens)
+    CALL attach(fid, TRIM(str),  "write_fvel",  write_fvel)
+    CALL attach(fid, TRIM(str),  "write_temp",  write_temp)
 
   END SUBROUTINE diag_par_outputinputs
 
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 4d83dfc1..0fd98081 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -74,14 +74,14 @@ SUBROUTINE compute_nonlinear
     p_int = parray_e(ip)
     jloope: DO ij = ijs_e, ije_e ! Loop over Laguerre moments
       j_int=jarray_e(ij)
-      IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxe)) THEN !compute
+      IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxe)) THEN !compute
         ! Set non linear sum truncation
         IF (NL_CLOS .EQ. -2) THEN
           nmax = Jmaxe
         ELSEIF (NL_CLOS .EQ. -1) THEN
           nmax = Jmaxe-j_int
         ELSE
-          nmax = NL_CLOS
+          nmax = min(NL_CLOS,Jmaxe-j_int)
         ENDIF
         bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
         nloope: DO in = 1,nmax+1 ! Loop over laguerre for the sum
@@ -89,7 +89,7 @@ SUBROUTINE compute_nonlinear
           F_cmpx(ikxs:ikxe,ikxs:ikxe) = phi(ikxs:ikxe,ikxs:ikxe,iz) * kernel_e(in, ikxs:ikxe,ikxs:ikxe, iz, eo)
           ! Second convolution terms
           G_cmpx(ikxs:ikxe,ikxs:ikxe) = 0._dp ! initialization of the sum
-          smax = MIN( (in-1)+(ij-1), jmaxe );
+          smax = MIN( (in-1)+(ij-1), Jmaxe );
           DO is = 1, smax+1 ! sum truncation on number of moments
             G_cmpx(ikxs:ikxe,ikxs:ikxe)  = G_cmpx(ikxs:ikxe,ikxs:ikxe) + &
               dnjs(in,ij,is) * moments_e(ip,is,ikxs:ikxe,ikxs:ikxe,iz,updatetlevel)
@@ -120,14 +120,14 @@ zloopi: DO iz = izs,ize
     eo = MODULO(parray_i(ip),2)
     jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
       j_int=jarray_i(ij)
-      IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxi)) THEN !compute
+      IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute
         ! Set non linear sum truncation
         IF (NL_CLOS .EQ. -2) THEN
           nmax = Jmaxi
         ELSEIF (NL_CLOS .EQ. -1) THEN
           nmax = Jmaxi-j_int
         ELSE
-          nmax = NL_CLOS
+          nmax = min(NL_CLOS,Jmaxi-j_int)
         ENDIF
         bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
         nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
diff --git a/src/ppinit.F90 b/src/ppinit.F90
index 1a662bab..55db1460 100644
--- a/src/ppinit.F90
+++ b/src/ppinit.F90
@@ -65,6 +65,10 @@ SUBROUTINE ppinit
   CALL MPI_COMM_RANK(comm_ky, rank_ky, ierr)
   CALL MPI_COMM_RANK(comm_z,  rank_z,  ierr)
 
+  ! 2D communicator
+  CALL MPI_CART_SUB (comm0, (/.TRUE.,.FALSE.,.TRUE./),  comm_pz,  ierr)
+  CALL MPI_CART_SUB (comm0, (/.FALSE.,.TRUE.,.TRUE./),  comm_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
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index f8ea45ee..07f4a3d0 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -42,18 +42,18 @@ if 0
 % Options
 options.INTERP    = 0;
 options.POLARPLOT = 0;
-% options.NAME      = '\phi';
-options.NAME      = 'N_i^{00}';
+options.NAME      = '\phi';
+% options.NAME      = 'N_i^{00}';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'kxky';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 9;
 % options.TIME      = dat.Ts5D;
-options.TIME      = 200:1:600;
+options.TIME      = 300:1:4500;
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -62,7 +62,7 @@ end
 if 0
 %% 2D snapshots
 % Options
-options.INTERP    = 1;
+options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 options.NAME      = '\phi';
@@ -71,11 +71,11 @@ options.NAME      = '\phi';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
-options.PLAN      = 'xz';
-% options.NAME      = 'f_e';
+options.PLAN      = 'kxky';
+% options.NAME      = 'f_i';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [260:265];
+options.TIME      = [300 350 400];
 data.a = data.EPS * 2e3;
 fig = photomaton(data,options);
 save_figure(data,fig)
@@ -100,7 +100,7 @@ if 0
 options.SPAR      = gene_data.vp';
 options.XPERP     = gene_data.mu';
 options.Z         = 'avg';
-options.T         = 340;
+options.T         = 900;
 options.PLT_FCT   = 'contour';
 options.ONED      = 0;
 options.non_adiab = 1;
@@ -113,16 +113,16 @@ end
 if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
-options.P2J        = 1;
+options.P2J        = 0;
 options.ST         = 0;
 options.PLOT_TYPE  = 'space-time';
 options.NORMALIZED = 0;
-options.JOBNUM     = 03;
-options.TIME       = [200 500];
+options.JOBNUM     = 0;
+options.TIME       = [300 500];
 options.specie     = 'i';
 options.compz      = 'avg';
-% fig = show_moments_spectrum(data,options);
-fig = show_napjz(data,options);
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
 save_figure(data,fig)
 end
 
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index 5188876d..012bd79e 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -5,7 +5,8 @@ helazdir = '/home/ahoffman/HeLaZ/';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
+% outfile ='shearless_cyclone/128x128x16xdmax_6_L_120_CBC_1.0';
+outfile ='shearless_cyclone/128x128x16xdmax_8_L_120_CBC_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/gene_analysis_3D.m b/wk/gene_analysis_3D.m
index 5716569d..5f83d870 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -26,7 +26,7 @@ options.INTERP    = 0;
 options.POLARPLOT = 0;
 options.AXISEQUAL = 1;
 % options.NAME      = 'Q_x';
-options.NAME      = 'n_i';
+options.NAME      = '\phi';
 % options.NAME      = 'T_i';
 % options.NAME      = '\Gamma_x';
 % options.NAME      = 'k^2n_e';
@@ -34,7 +34,7 @@ options.PLAN      = 'xy';
 % options.NAME      ='f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
-options.TIME      = [200];
+options.TIME      = [500];
 gene_data.a = data.EPS * 2000;
 fig = photomaton(gene_data,options);
 save_figure(gene_data,fig)
diff --git a/wk/quick_run.m b/wk/quick_run.m
index b3f5829b..f5a8c35d 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -7,7 +7,7 @@ RUN = 1; % To run or just to load
 addpath(genpath('../matlab')) % ... add
 default_plots_options
 HELAZDIR = '/home/ahoffman/HeLaZ/';
-EXECNAME = 'helaz3.12';
+EXECNAME = 'helaz3';
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
@@ -21,25 +21,25 @@ K_E     = 0.0;   % Electrostat '''
 SIGMA_E = 0.0233380;   % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
 KIN_E   = 0;     % 1: kinetic electrons, 2: adiabatic electrons
 %% GRID PARAMETERS
-PMAXE   = 8;     % Hermite basis size of electrons
-JMAXE   = 4;     % Laguerre "
-PMAXI   = 8;     % " ions
-JMAXI   = 4;     % "
+PMAXE   = 4;     % Hermite basis size of electrons
+JMAXE   = 2;     % Laguerre "
+PMAXI   = 4;     % " ions
+JMAXI   = 2;     % "
 NX      = 1;    % real space x-gridpoints
-NY      = 20;     %     ''     y-gridpoints
-LX      = 120;   % Size of the squared frequency domain
-LY      = 120;     % Size of the squared frequency domain
-NZ      = 16;     % number of perpendicular planes (parallel grid)
+NY      = 24;     %     ''     y-gridpoints
+LX      = 100;   % Size of the squared frequency domain
+LY      = 100;     % Size of the squared frequency domain
+NZ      = 32;     % number of perpendicular planes (parallel grid)
 NPOL    = 1;
 SG      = 0;     % Staggered z grids option
 %% GEOMETRY
-GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-% GEOMETRY= 's-alpha';
+% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+GEOMETRY= 's-alpha';
 Q0      = 1.4;    % safety factor
 SHEAR   = 0.0;    % magnetic shear (Not implemented yet)
 EPS     = 0.18;    % inverse aspect ratio
 %% TIME PARMETERS
-TMAX    = 50;  % Maximal time unit
+TMAX    = 10;  % Maximal time unit
 DT      = 2e-2;   % Time step
 SPS0D   = 1;      % Sampling per time unit for 2D arrays
 SPS2D   = 0;      % Sampling per time unit for 2D arrays
@@ -90,7 +90,7 @@ system(['rm fort*.90']);
 % Run linear simulation
 if RUN
     system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',HELAZDIR,'bin/',EXECNAME,' 2 1 1 0; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
 %     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 3 2 0; cd ../../../wk'])
 end
@@ -104,7 +104,7 @@ JOBNUMMIN = 00; JOBNUMMAX = 00;
 data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 
 %% Short analysis
-if 1
+if 0
 %% linear growth rate (adapted for 2D zpinch and fluxtube)
 trange = [0.5 1]*data.Ts3D(end);
 nplots = 1;
@@ -114,12 +114,30 @@ end
 if 0
 %% Ballooning plot
 options.time_2_plot = data.Ts3D(end);
-options.kymodes     = [0.2 0.5 0.7];
+options.kymodes     = [0.4 0.5 0.6];
 options.normalized  = 1;
 options.field       = 'phi';
 fig = plot_ballooning(data,options);
 end
 
+if 1
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.P2J        = 1;
+options.ST         = 1;
+options.PLOT_TYPE  = 'space-time';
+% options.PLOT_TYPE  =   'Tavg-1D';
+% options.PLOT_TYPE  = 'Tavg-2D';
+options.NORMALIZED = 0;
+options.JOBNUM     = 0;
+options.TIME       = [0 50];
+options.specie     = 'i';
+options.compz      = 'avg';
+% fig = show_moments_spectrum(data,options);
+fig = show_napjz(data,options);
+save_figure(data,fig)
+end
+
 if 0
 %% linear growth rate for 3D Zpinch (kz fourier transform)
 trange = [0.5 1]*data.Ts3D(end);
-- 
GitLab