diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 958c08597497a341cb9a49ab7f2e659b60713b92..c3be68638cef60832c2f3790c9d5c0e162fb791c 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -95,20 +95,18 @@ SUBROUTINE diagnose_full(kstep)
   USE grid,            ONLY: &
     local_nj,local_nky,local_nkx,local_nz,ngj,ngz,&
     parray_full,pmax,jarray_full,jmax,&
-    kyarray_full,kxarray_full,zarray_full, total_na
+    kyarray_full,kxarray_full,zarray_full
   USE diagnostics_par
   USE futils,          ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf, putarrnd
-  USE species,         ONLY: name
   USE array
   USE model,           ONLY: EM
   USE parallel,        ONLY: my_id, comm0
   USE collision,       ONLY: coll_outputinputs
   USE geometry,        ONLY: gxx,gxy,gyy,gxz,gyz,gzz,hatR,hatZ,hatB,dBdx,dBdy,dBdz,Jacobian,gradz_coeff,Ckxky
   IMPLICIT NONE
-  CHARACTER           :: letter_a
   INTEGER, INTENT(in) :: kstep
   INTEGER, parameter  :: BUFSIZE = 2
-  INTEGER :: rank = 0, ierr, ia
+  INTEGER :: rank = 0, ierr
   INTEGER :: dims(1) = (/0/)
   !____________________________________________________________________________
   !                   1.   Initial diagnostics
@@ -256,15 +254,12 @@ END SUBROUTINE diagnose_full
 !!-------------- Auxiliary routines -----------------!!
 SUBROUTINE diagnose_0d
   USE basic
-  USE futils, ONLY: append, attach, getatt
+  USE futils, ONLY: append, attach, getatt, creatd
   USE diagnostics_par
   USE prec_const
   USE processing
   USE model,   ONLY: Na
-  USE species, ONLY: name
   IMPLICIT NONE
-  CHARACTER :: letter_a
-  INTEGER   :: ia
   ! Time measurement data
   CALL append(fidres, "/profiler/Tc_rhs",       REAL(chrono_mrhs%ttot,dp),ionode=0)
   CALL append(fidres, "/profiler/Tc_adv_field", REAL(chrono_advf%ttot,dp),ionode=0)
@@ -278,6 +273,10 @@ SUBROUTINE diagnose_0d
   CALL append(fidres, "/profiler/Tc_grad",      REAL(chrono_grad%ttot,dp),ionode=0)
   CALL append(fidres, "/profiler/Tc_nadiab",    REAL(chrono_napj%ttot,dp),ionode=0)
   CALL append(fidres, "/profiler/Tc_step",      REAL(chrono_step%ttot,dp),ionode=0)
+#ifdef TEST_SVD
+  CALL creatd(fidres, 0, (/0/), "/profiler/Tc_DLRA", "cumulative total DLRA computation time")
+  CALL append(fidres, "/profiler/Tc_DLRA",      REAL(chrono_DLRA%ttot,dp),ionode=0) 
+#endif
   CALL append(fidres, "/profiler/time",                REAL(time,dp),ionode=0)
   ! Processing data
   CALL append(fidres,  "/data/var0d/time",      REAL(time,dp),ionode=0)
@@ -311,13 +310,9 @@ SUBROUTINE diagnose_3d
   USE prec_const
   USE processing, ONLY: compute_fluid_moments, compute_Napjz_spectrum
   USE model,      ONLY: EM
-  USE species,    ONLY: name
   USE parallel,   ONLY: manual_3D_bcast
   IMPLICIT NONE
-  CHARACTER :: letter_a
-  INTEGER   :: ia
   COMPLEX(xp), DIMENSION(local_na,local_nky,local_nkx,local_nz) :: Na00_
-  COMPLEX(xp), DIMENSION(local_na,local_np, local_nj, local_nz) :: Napjz_
   ! add current time, cstep and frame
   CALL append(fidres,  "/data/var3d/time",           REAL(time,dp),ionode=0)
   CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
@@ -337,8 +332,7 @@ SUBROUTINE diagnose_3d
     ENDIF
     CALL write_field3da_kykxz(Na00_, 'Na00')
     ! <<Napj>x>y spectrum
-    Napjz_ = Napjz
-    CALL write_field3da_pjz(Napjz_, 'Napjz')
+    CALL write_field3da_pjz(Napjz, 'Napjz')
   ENDIF
   !! Fuid moments
   IF (write_dens .OR. write_fvel .OR. write_temp) &
@@ -351,15 +345,15 @@ SUBROUTINE diagnose_3d
     CALL write_field3da_kykxz(uper, 'uper')
   ENDIF
   IF (write_temp) THEN
-    CALL write_field3d_kykxz(Tpar, 'Tpar')
-    CALL write_field3d_kykxz(Tper, 'Tper')
-    CALL write_field3d_kykxz(temp, 'temp')
+    CALL write_field3da_kykxz(Tpar, 'Tpar')
+    CALL write_field3da_kykxz(Tper, 'Tper')
+    CALL write_field3da_kykxz(temp, 'temp')
   ENDIF
   CONTAINS
     SUBROUTINE write_field3d_kykxz(field, text)
       USE parallel, ONLY : gather_xyz, num_procs
       IMPLICIT NONE
-      COMPLEX(xp), DIMENSION(local_nky,total_nkx,local_nz), INTENT(IN) :: field
+      COMPLEX(xp), DIMENSION(:,:,:), INTENT(IN) :: field
       CHARACTER(*), INTENT(IN) :: text
       COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: field_full
       CHARACTER(50) :: dset_name
@@ -376,16 +370,23 @@ SUBROUTINE diagnose_3d
       SUBROUTINE write_field3da_kykxz(field, text)
         USE parallel, ONLY : gather_xyz, num_procs
         IMPLICIT NONE
-        COMPLEX(xp), DIMENSION(local_na,local_nky,total_nkx,local_nz), INTENT(IN) :: field
+        COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(IN) :: field
         CHARACTER(*), INTENT(IN) :: text
         COMPLEX(xp), DIMENSION(total_na,total_nky,total_nkx,total_nz) :: field_full
+        COMPLEX(xp), DIMENSION(local_nky,total_nkx,local_nz) :: buff_local
+        COMPLEX(xp), DIMENSION(total_nky,total_nkx,total_nz) :: buff_full
         CHARACTER(50) :: dset_name
+        INTEGER :: ia
         field_full = 0;
         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, ionode=0)
         ELSE
-          CALL gather_xyz(field,field_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
+          DO ia = 1,total_Na
+            buff_local = field(ia,:,:,:)
+            CALL gather_xyz(buff_local,buff_full,local_nky,total_nky,total_nkx,local_nz,total_nz)
+            field_full(ia,:,:,:) = buff_full
+          ENDDO
           CALL putarr(fidres, dset_name, field_full, ionode=0)
         ENDIF
         END SUBROUTINE write_field3da_kykxz
@@ -393,7 +394,7 @@ SUBROUTINE diagnose_3d
     SUBROUTINE write_field3da_pjz(field, text)
       USE parallel, ONLY : gather_pjz, num_procs
       IMPLICIT NONE
-      COMPLEX(xp), DIMENSION(local_na,local_np,local_nj,local_nz), INTENT(IN) :: field
+      COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(IN) :: field
       CHARACTER(*), INTENT(IN) :: text
       COMPLEX(xp), DIMENSION(total_na,total_np,total_nj,total_nz) :: field_full
       CHARACTER(LEN=50) :: dset_name
@@ -402,7 +403,7 @@ SUBROUTINE diagnose_3d
       IF (num_procs .EQ. 1) THEN ! no data distribution
         CALL putarr(fidres, dset_name, field, ionode=0)
       ELSE
-        CALL gather_pjz(field,field_full,local_np,total_np,total_nj,local_nz,total_nz)
+        CALL gather_pjz(field,field_full,total_na,local_np,total_np,total_nj,local_nz,total_nz)
         CALL putarr(fidres, dset_name, field_full, ionode=0)
       ENDIF
       CALL attach(fidres, dset_name, "time", time)
@@ -418,7 +419,6 @@ SUBROUTINE diagnose_5d
   USE grid,   ONLY:total_np, total_nj, total_nky, total_nkx, total_nz, &
                    local_np, local_nj, local_nky, local_nkx, local_nz, &
                    ngp, ngj, ngz, total_na
-  USE time_integration, ONLY: updatetlevel, ntimelevel
   USE diagnostics_par
   USE prec_const, ONLY: xp,dp
   IMPLICIT NONE
@@ -436,18 +436,18 @@ SUBROUTINE diagnose_5d
   CONTAINS
 
   SUBROUTINE write_field5d(field, text)
-    USE basic,    ONLY: GATHERV_OUTPUT, jobnum, dt
-    USE futils,   ONLY: attach, putarr, putarrnd
-    USE parallel, ONLY: gather_pjxyz, num_procs
-    USE prec_const, ONLY: xp
+    USE basic,            ONLY: GATHERV_OUTPUT, jobnum, dt
+    USE futils,           ONLY: attach, putarr, putarrnd
+    USE parallel,         ONLY: gather_pjxyz, num_procs
+    USE prec_const,       ONLY: xp
     IMPLICIT NONE
-    COMPLEX(xp), DIMENSION(total_na,local_np+ngp,local_nj+ngj,local_nky,local_nkx,local_nz+ngz,ntimelevel), INTENT(IN) :: field
+    COMPLEX(xp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: field
     CHARACTER(*), INTENT(IN) :: text
     COMPLEX(xp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
     COMPLEX(xp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
     CHARACTER(LEN=50) :: dset_name
-    field_sub  = field(1:total_na,(1+ngp/2):(local_np+ngp/2),((1+ngj/2)):((local_nj+ngj/2)),&
-                          1:local_nky,1:local_nkx,  ((1+ngz/2)):((local_nz+ngz/2)),updatetlevel)
+    field_sub  = field(1:total_na,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
+                          1:local_nky,1:local_nkx,(1+ngz/2):(local_nz+ngz/2),1)
     field_full = 0;
     WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
     IF (num_procs .EQ. 1) THEN
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
index e132fd04c4787c82b3fa9b6c64c88ebf10f0a08a..b4e398b6e7e13b67ae36524cd8b29279bd02166d 100644
--- a/src/parallel_mod.F90
+++ b/src/parallel_mod.F90
@@ -239,8 +239,8 @@ CONTAINS
   SUBROUTINE gather_xyz(field_loc,field_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot
-    COMPLEX(xp), DIMENSION(nky_loc,nkx_tot,nz_loc), INTENT(IN)  :: field_loc
-    COMPLEX(xp), DIMENSION(nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
+    COMPLEX(xp), DIMENSION(:,:,:), INTENT(IN)  :: field_loc
+    COMPLEX(xp), DIMENSION(:,:,:), INTENT(OUT) :: field_tot
     COMPLEX(xp), DIMENSION(nky_tot,nz_loc) :: buffer_yt_zl !full  y, local z
     COMPLEX(xp), DIMENSION(nky_tot,nz_tot) :: buffer_yt_zt !full  y, full  z
     COMPLEX(xp), DIMENSION(nky_loc):: buffer_yl_zc !local y, constant z
@@ -272,42 +272,45 @@ CONTAINS
       ENDDO
     ENDIF
   END SUBROUTINE gather_xyz
+  
 
   !!!!! Gather a field in kinetic + z coordinates on rank 0 !!!!!
-  SUBROUTINE gather_pjz(field_loc,field_tot,np_loc,np_tot,nj_tot,nz_loc,nz_tot)
+  SUBROUTINE gather_pjz(field_loc,field_tot,na_tot,np_loc,np_tot,nj_tot,nz_loc,nz_tot)
     IMPLICIT NONE
-    INTEGER, INTENT(IN) :: np_loc,np_tot, nj_tot, nz_loc,nz_tot
-    COMPLEX(xp), DIMENSION(np_loc,nj_tot,nz_loc), INTENT(IN)  :: field_loc
-    COMPLEX(xp), DIMENSION(np_tot,nj_tot,nz_tot), INTENT(OUT) :: field_tot
+    INTEGER, INTENT(IN) :: na_tot,np_loc,np_tot, nj_tot, nz_loc,nz_tot
+    COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(IN)  :: field_loc
+    COMPLEX(xp), DIMENSION(:,:,:,:), INTENT(OUT) :: field_tot
     COMPLEX(xp), DIMENSION(np_tot,nz_loc) :: buffer_pt_zl !full  p, local z
     COMPLEX(xp), DIMENSION(np_tot,nz_tot) :: buffer_pt_zt !full  p, full  z
     COMPLEX(xp), DIMENSION(np_loc) :: buffer_pl_zc !local p, constant z
     COMPLEX(xp), DIMENSION(np_tot) :: buffer_pt_zc !full  p, constant z
-    INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz
+    INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz, ia
 
     snd_p  = np_loc    ! Number of points to send along p (per z)
     snd_z  = np_tot*nz_loc ! Number of points to send along z (full p)
 
     root_p = 0; root_z = 0; root_ky = 0
     IF(rank_ky .EQ. root_ky) THEN
-      DO ij = 1,nj_tot
-        DO iz = 1,nz_loc
-          ! fill a buffer to contain a slice of data at constant j and z
-          buffer_pl_zc(1:np_loc) = field_loc(1:np_loc,ij,iz)
-          CALL MPI_GATHERV(buffer_pl_zc, snd_p,        mpi_xp_c, &
-                           buffer_pt_zc, rcv_p, dsp_p, mpi_xp_c, &
-                           root_p, comm_p, ierr)
-          buffer_pt_zl(1:np_tot,iz) = buffer_pt_zc(1:np_tot)
+      DO ia= 1,na_tot
+        DO ij = 1,nj_tot
+          DO iz = 1,nz_loc
+            ! fill a buffer to contain a slice of data at constant j and z
+            buffer_pl_zc(1:np_loc) = field_loc(ia,1:np_loc,ij,iz)
+            CALL MPI_GATHERV(buffer_pl_zc, snd_p,        mpi_xp_c, &
+                            buffer_pt_zc, rcv_p, dsp_p, mpi_xp_c, &
+                            root_p, comm_p, ierr)
+            buffer_pt_zl(1:np_tot,iz) = buffer_pt_zc(1:np_tot)
+          ENDDO
+          ! send the full line on y contained by root_p
+          IF(rank_p .EQ. root_p) THEN
+            CALL MPI_GATHERV(buffer_pt_zl, snd_z,         mpi_xp_c, &
+                            buffer_pt_zt, rcv_zp, dsp_zp, mpi_xp_c, &
+                            root_z, comm_z, ierr)
+          ENDIF
+          ! ID 0 (the one who output) rebuild the whole array
+          IF(my_id .EQ. 0) &
+            field_tot(ia,1:np_tot,ij,1:nz_tot) = buffer_pt_zt(1:np_tot,1:nz_tot)
         ENDDO
-        ! send the full line on y contained by root_p
-        IF(rank_p .EQ. root_p) THEN
-          CALL MPI_GATHERV(buffer_pt_zl, snd_z,          mpi_xp_c, &
-                           buffer_pt_zt, rcv_zp, dsp_zp, mpi_xp_c, &
-                           root_z, comm_z, ierr)
-        ENDIF
-        ! ID 0 (the one who output) rebuild the whole array
-        IF(my_id .EQ. 0) &
-          field_tot(1:np_tot,ij,1:nz_tot) = buffer_pt_zt(1:np_tot,1:nz_tot)
       ENDDO
     ENDIF
   END SUBROUTINE gather_pjz
@@ -317,8 +320,8 @@ CONTAINS
   SUBROUTINE gather_pjxyz(field_loc,field_tot,na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: na_tot,np_loc,np_tot,nj_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot
-    COMPLEX(xp), DIMENSION(np_loc,nj_tot,nky_loc,nkx_tot,nz_loc), INTENT(IN)  :: field_loc
-    COMPLEX(xp), DIMENSION(na_tot,np_tot,nj_tot,nky_tot,nkx_tot,nz_tot), INTENT(OUT) :: field_tot
+    COMPLEX(xp), DIMENSION(:,:,:,:,:,:), INTENT(IN)  :: field_loc
+    COMPLEX(xp), DIMENSION(:,:,:,:,:,:), INTENT(OUT) :: field_tot
     COMPLEX(xp), DIMENSION(np_tot,nky_tot,nz_loc) :: buffer_pt_yt_zl ! full p,     full y, local    z
     COMPLEX(xp), DIMENSION(np_tot,nky_tot,nz_tot) :: buffer_pt_yt_zt ! full p,     full y, full     z
     COMPLEX(xp), DIMENSION(np_tot,nky_loc) :: buffer_pt_yl_zc     ! full p,    local y, constant z
@@ -336,7 +339,7 @@ CONTAINS
           z: DO iz = 1,nz_loc
             y: DO iy = 1,nky_loc
               ! fill a buffer to contain a slice of p data at constant j, ky, kx and z
-              buffer_pl_cy_zc(1:np_loc) = field_loc(1:np_loc,ij,iy,ix,iz)
+              buffer_pl_cy_zc(1:np_loc) = field_loc(ia,1:np_loc,ij,iy,ix,iz)
               CALL MPI_GATHERV(buffer_pl_cy_zc, snd_p,       mpi_xp_c, &
                               buffer_pt_cy_zc, rcv_p, dsp_p, mpi_xp_c, &
                               root_p, comm_p, ierr)
@@ -368,7 +371,7 @@ CONTAINS
   SUBROUTINE manual_3D_bcast(field_,n1,n2,n3)
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: n1,n2,n3
-    COMPLEX(xp), DIMENSION(n1,n2,n3), INTENT(INOUT) :: field_
+    COMPLEX(xp), DIMENSION(:,:,:), INTENT(INOUT) :: field_
     COMPLEX(xp) :: buffer(n1,n2,n3)
     INTEGER     :: i_, root, world_rank, world_size, count, i1,i2,i3
     root = 0;
@@ -441,7 +444,7 @@ CONTAINS
   SUBROUTINE exchange_ghosts_1D(f,nbr_L,nbr_R,np,ng)
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: np,ng, nbr_L, nbr_R
-    COMPLEX(xp), DIMENSION(np+ng), INTENT(INOUT) :: f
+    COMPLEX(xp), DIMENSION(:), INTENT(INOUT) :: f
     INTEGER :: ierr, first, last, ig
     first = 1 + ng/2
     last  = np + ng/2