diff --git a/src/advance_field.F90 b/src/advance_field.F90
index c615d87..0c27b59 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -1,56 +1,56 @@
 MODULE advance_field_routine
 
 USE prec_const
 implicit none
 
 CONTAINS
 
   SUBROUTINE advance_time_level
 
     USE basic
     USE time_integration
     use prec_const
     IMPLICIT NONE
 
-    call set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
+    CALL set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
 
   END SUBROUTINE advance_time_level
 
 
 
   SUBROUTINE advance_field( f, f_rhs )
 
     USE basic
     USE time_integration
     USE array
     USE grid
     use prec_const
     IMPLICIT NONE
 
     COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f
     COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f_rhs
     INTEGER :: istage
 
     SELECT CASE (updatetlevel)
     CASE(1)
       DO ikz=ikzs,ikze
         DO ikr=ikrs,ikre
           DO istage=1,ntimelevel
             f(ikr,ikz,1) = f(ikr,ikz,1) + dt*b_E(istage)*f_rhs(ikr,ikz,istage)
           END DO
         END DO
       END DO
     CASE DEFAULT
       DO ikz=ikzs,ikze
         DO ikr=ikrs,ikre
           f(ikr,ikz,updatetlevel) = f(ikr,ikz,1);
           DO istage=1,updatetlevel-1
             f(ikr,ikz,updatetlevel) = f(ikr,ikz,updatetlevel) + &
                                   dt*A_E(updatetlevel,istage)*f_rhs(ikr,ikz,istage)
           END DO
         END DO
       END DO
     END SELECT
   END SUBROUTINE advance_field
 
 END MODULE advance_field_routine
diff --git a/src/advance_field_adapt.F90 b/src/advance_field_adapt.F90
index d2d18f2..2380904 100644
--- a/src/advance_field_adapt.F90
+++ b/src/advance_field_adapt.F90
@@ -1,80 +1,80 @@
 MODULE advance_field_routine
 
 USE prec_const
 implicit none
 
 CONTAINS
 
   SUBROUTINE advance_time_level
 
     USE basic
     USE time_integration
     use prec_const
     IMPLICIT NONE
 
-    call set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
+    CALL set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
 
   END SUBROUTINE advance_time_level
 
 
 
   SUBROUTINE advance_field( f, f_rhs )
 
     USE basic
     USE time_integration
     USE array
     USE grid
     use prec_const
     IMPLICIT NONE
 
     COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f
     COMPLEX(dp), DIMENSION ( ikrs:ikre, ikzs:ikze, ntimelevel ) :: f_rhs
     REAL(dp)    :: error
     INTEGER     :: istage
 
     SELECT CASE (updatetlevel)
 
     CASE(1)
       SELECT CASE (numerical_scheme)
 
       CASE ('DOPRI5_ADAPT')
         error = 0._dp
         DO ikz=ikzs,ikze
           DO ikr=ikrs,ikre
             fs = f(ikr,ikz,1)
             DO istage=1,ntimelevel
               f(ikr,ikz,1) = f(ikr,ikz,1) + dt*b_E(istage)*f_rhs(ikr,ikz,istage)
               fs = fs + dt*b_Es(istage)*f_rhs(ikr,ikz,istage)
             END DO
             IF ( ABS(f(ikr,ikz,1) - fs) .GT. error ) THEN
               error = ABS(f(ikr,ikz,1) - fs)
             ENDIF
           END DO
         END DO
         IF (error > TOL)
       CASE DEFAULT
         DO ikz=ikzs,ikze
           DO ikr=ikrs,ikre
             fs = f(ikr,ikz,1)
             DO istage=1,ntimelevel
             f(ikr,ikz,1) = f(ikr,ikz,1) + dt*b_E(istage)*f_rhs(ikr,ikz,istage)
             END DO
           END DO
         END DO
       END SELECT
 
     CASE DEFAULT
       DO ikz=ikzs,ikze
         DO ikr=ikrs,ikre
           f(ikr,ikz,updatetlevel) = f(ikr,ikz,1);
           DO istage=1,updatetlevel-1
             f(ikr,ikz,updatetlevel) = f(ikr,ikz,updatetlevel) + &
                                   dt*A_E(updatetlevel,istage)*f_rhs(ikr,ikz,istage)
           END DO
         END DO
       END DO
     END SELECT
 
   END SUBROUTINE advance_field
 
 END MODULE advance_field_routine
diff --git a/src/auxval.F90 b/src/auxval.F90
index f0c451c..c0cb1b6 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -1,21 +1,21 @@
 subroutine auxval
   !   Set auxiliary values, at beginning of simulation
 
   USE basic
   USE grid
   USE array
   ! use mumps_bsplines, only: putrow_mumps => putrow, updtmat_mumps => updtmat, factor_mumps => factor, bsolve_mumps => bsolve ! from BSPLINES
   use prec_const
   IMPLICIT NONE
 
   INTEGER :: irows,irowe, irow, icol
   WRITE(*,*) '=== Set auxiliary values ==='
 
-  call set_rgrid
-  call set_zgrid
-  call set_krgrid
-  call set_kzgrid
-  call set_pj
+  CALL set_rgrid
+  CALL set_zgrid
+  CALL set_krgrid
+  CALL set_kzgrid
+  CALL set_pj
 
   CALL memory ! Allocate memory for global arrays
 END SUBROUTINE auxval
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 735a3b4..eb0541e 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -1,238 +1,266 @@
 MODULE basic
   !   Basic module for time dependent problems
 
   use prec_const
   IMPLICIT none
 
   INTEGER  :: nrun   = 1           ! Number of time steps to run
   real(dp) :: tmax   = 100000.0    ! Maximum simulation time
   real(dp) :: dt     = 1.0         ! Time step
   real(dp) :: time   = 0           ! Current simulation time (Init from restart file)
 
   INTEGER :: jobnum  = 0           ! Job number
   INTEGER :: step    = 0           ! Calculation step of this run
   INTEGER :: cstep   = 0           ! Current step number (Init from restart file)
   LOGICAL :: RESTART = .FALSE.     ! Signal end of run
   LOGICAL :: nlend   = .FALSE.     ! Signal end of run
 
   INTEGER :: ierr                  ! flag for MPI error
   INTEGER :: iframe1d              ! counting the number of times 1d datasets are outputed (for diagnose)
   INTEGER :: iframe2d              ! counting the number of times 2d datasets are outputed (for diagnose)
   INTEGER :: iframe3d              ! counting the number of times 3d datasets are outputed (for diagnose)
   INTEGER :: iframe5d              ! counting the number of times 5d datasets are outputed (for diagnose)
 
   !  List of logical file units
   INTEGER :: lu_in   = 90              ! File duplicated from STDIN
   INTEGER :: lu_job  = 91              ! myjob file
 
   real :: start, finish ! To measure computation time
+  real :: start_est, finish_est, time_est
 
   INTERFACE allocate_array
     MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4
     MODULE PROCEDURE allocate_array_dc1,allocate_array_dc2,allocate_array_dc3,allocate_array_dc4, allocate_array_dc5
     MODULE PROCEDURE allocate_array_i1,allocate_array_i2,allocate_array_i3,allocate_array_i4
     MODULE PROCEDURE allocate_array_l1,allocate_array_l2,allocate_array_l3,allocate_array_l4
   END INTERFACE allocate_array
 
 CONTAINS
 
   !================================================================================
   SUBROUTINE basic_data
     !   Read basic data for input file
 
     use prec_const
     IMPLICIT NONE
 
     NAMELIST /BASIC/  nrun, dt, tmax, RESTART
 
     READ(lu_in,basic)
 
   END SUBROUTINE basic_data
   !================================================================================
   SUBROUTINE daytim(str)
     !   Print date and time
 
     use prec_const
     IMPLICIT NONE
 
     CHARACTER(len=*) , INTENT(in) :: str
     CHARACTER(len=16) :: d, t, dat, time
     !________________________________________________________________________________
     !
     CALL DATE_AND_TIME(d,t)
     dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4)
     time=t(1:2) // ':' // t(3:4) // ':' // t(5:10)
     WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), time(1:12)
     !
   END SUBROUTINE daytim
   !================================================================================
+  SUBROUTINE display_h_min_s(time)
+    real :: time
+    integer :: days, hours, mins, secs
+    days = FLOOR(time/24./3600.);
+    hours= FLOOR(time/3600.);
+    mins = FLOOR(time/60.);
+    secs = FLOOR(time);
+
+    IF ( days .GT. 0 ) THEN !display day h min s
+      hours = (time/3600./24. - days) * 24
+      mins  = (time/3600. - days*24. - hours) * 60
+      secs  = (time/60. - days*24.*60 - hours*60 - mins) * 60
+      WRITE(*,*) 'CPU Time = ', days, '[day]', hours, '[h]', mins, '[min]', secs, '[s]'
+
+    ELSEIF ( hours .GT. 0 ) THEN !display h min s
+      mins  = (time/3600. - hours) * 60
+      secs  = (time/60. - hours*60 - mins) * 60
+      WRITE(*,*) 'CPU Time = ', hours, '[h]', mins, '[min]', secs, '[s]'
+
+    ELSEIF ( mins .GT. 0 ) THEN !display min s
+      secs  = (time/60. - mins) * 60
+      WRITE(*,*) 'CPU Time = ', mins, '[min]', secs, '[s]'
+
+    ELSE ! display s
+      WRITE(*,*) 'CPU Time = ', FLOOR(time), '[s]'
+    ENDIF
+  END SUBROUTINE display_h_min_s
+  !================================================================================
 
 ! To allocate arrays of doubles, integers, etc. at run time
-
   SUBROUTINE allocate_array_dp1(a,is1,ie1)
     IMPLICIT NONE
     real(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1
     ALLOCATE(a(is1:ie1))
     a=0.0_dp
   END SUBROUTINE allocate_array_dp1
 
   SUBROUTINE allocate_array_dp2(a,is1,ie1,is2,ie2)
     IMPLICIT NONE
     real(dp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
     ALLOCATE(a(is1:ie1,is2:ie2))
     a=0.0_dp
   END SUBROUTINE allocate_array_dp2
 
   SUBROUTINE allocate_array_dp3(a,is1,ie1,is2,ie2,is3,ie3)
     IMPLICIT NONE
     real(dp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
     a=0.0_dp
   END SUBROUTINE allocate_array_dp3
 
   SUBROUTINE allocate_array_dp4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
     IMPLICIT NONE
     real(dp), DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
     a=0.0_dp
   END SUBROUTINE allocate_array_dp4
 
   SUBROUTINE allocate_array_dp5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
     IMPLICIT NONE
     real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
     a=0.0_dp
   END SUBROUTINE allocate_array_dp5
 
   !========================================
 
   SUBROUTINE allocate_array_dc1(a,is1,ie1)
     IMPLICIT NONE
     DOUBLE COMPLEX, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1
     ALLOCATE(a(is1:ie1))
     a=CMPLX(0.0_dp,0.0_dp)
   END SUBROUTINE allocate_array_dc1
 
   SUBROUTINE allocate_array_dc2(a,is1,ie1,is2,ie2)
     IMPLICIT NONE
     DOUBLE COMPLEX, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
     ALLOCATE(a(is1:ie1,is2:ie2))
     a=CMPLX(0.0_dp,0.0_dp)
   END SUBROUTINE allocate_array_dc2
 
   SUBROUTINE allocate_array_dc3(a,is1,ie1,is2,ie2,is3,ie3)
     IMPLICIT NONE
     DOUBLE COMPLEX, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
     a=CMPLX(0.0_dp,0.0_dp)
   END SUBROUTINE allocate_array_dc3
 
   SUBROUTINE allocate_array_dc4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
     IMPLICIT NONE
     DOUBLE COMPLEX, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
     a=CMPLX(0.0_dp,0.0_dp)
   END SUBROUTINE allocate_array_dc4
 
   SUBROUTINE allocate_array_dc5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
     IMPLICIT NONE
     DOUBLE COMPLEX, DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
     a=CMPLX(0.0_dp,0.0_dp)
   END SUBROUTINE allocate_array_dc5
 
   !========================================
 
   SUBROUTINE allocate_array_i1(a,is1,ie1)
     IMPLICIT NONE
     INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1
     ALLOCATE(a(is1:ie1))
     a=0
   END SUBROUTINE allocate_array_i1
 
   SUBROUTINE allocate_array_i2(a,is1,ie1,is2,ie2)
     IMPLICIT NONE
     INTEGER, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
     ALLOCATE(a(is1:ie1,is2:ie2))
     a=0
   END SUBROUTINE allocate_array_i2
 
   SUBROUTINE allocate_array_i3(a,is1,ie1,is2,ie2,is3,ie3)
     IMPLICIT NONE
     INTEGER, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
     a=0
   END SUBROUTINE allocate_array_i3
 
   SUBROUTINE allocate_array_i4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
     IMPLICIT NONE
     INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
     a=0
   END SUBROUTINE allocate_array_i4
 
   SUBROUTINE allocate_array_i5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
     IMPLICIT NONE
     real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
     a=0
   END SUBROUTINE allocate_array_i5
 
   !========================================
 
   SUBROUTINE allocate_array_l1(a,is1,ie1)
     IMPLICIT NONE
     LOGICAL, DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1
     ALLOCATE(a(is1:ie1))
     a=.false.
   END SUBROUTINE allocate_array_l1
 
   SUBROUTINE allocate_array_l2(a,is1,ie1,is2,ie2)
     IMPLICIT NONE
     LOGICAL, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2
     ALLOCATE(a(is1:ie1,is2:ie2))
     a=.false.
   END SUBROUTINE allocate_array_l2
 
   SUBROUTINE allocate_array_l3(a,is1,ie1,is2,ie2,is3,ie3)
     IMPLICIT NONE
     LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3))
     a=.false.
   END SUBROUTINE allocate_array_l3
 
   SUBROUTINE allocate_array_l4(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4)
     IMPLICIT NONE
     LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4))
     a=.false.
   END SUBROUTINE allocate_array_l4
 
   SUBROUTINE allocate_array_l5(a,is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5)
     IMPLICIT NONE
     real(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE, INTENT(INOUT) :: a
     INTEGER, INTENT(IN) :: is1,ie1,is2,ie2,is3,ie3,is4,ie4,is5,ie5
     ALLOCATE(a(is1:ie1,is2:ie2,is3:ie3,is4:ie4,is5:ie5))
     a=.false.
   END SUBROUTINE allocate_array_l5
 
 END MODULE basic
diff --git a/src/control.F90 b/src/control.F90
index 59abdcf..92638b0 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -1,74 +1,86 @@
 SUBROUTINE control
   !   Control the run
 
   use basic
   use prec_const
   IMPLICIT NONE
 
-  call cpu_time(start)
+  CALL cpu_time(start)
   !________________________________________________________________________________
   !              1.   Prologue
   !                   1.1     Initialize the parallel environment
   WRITE(*,*) 'Initialize MPI...'
   CALL ppinit
   WRITE(*,'(a/)') '...MPI initialized'
 
 
-  call daytim('Start at ')
+  CALL daytim('Start at ')
 
   !                   1.2     Define data specific to run
   WRITE(*,*) 'Load basic data...'
   CALL basic_data
   WRITE(*,'(a/)') '...basic data loaded.'
 
 
   !                   1.3   Read input parameters from input file
   WRITE(*,*) 'Read input parameters...'
   CALL readinputs
   WRITE(*,'(a/)') '...input parameters read'
 
   !                   1.4     Set auxiliary values (allocate arrays, set grid, ...)
   WRITE(*,*) 'Calculate auxval...'
   CALL auxval
   WRITE(*,'(a/)') '...auxval calculated'
   !
   !                   1.5     Initial conditions
   WRITE(*,*) 'Create initial state...'
   CALL inital
   WRITE(*,'(a/)') '...initial state created'
 
+  ! !                   1.5.1   Computation time estimation
+  ! WRITE(*,*) 'Estimation of computation time...'
+  ! CALL cpu_time(start_est)
+  ! CALL stepon
+  ! CALL cpu_time(finish_est)
+  ! time_est = 4.*tmax/dt*(finish_est-start_est)
+  ! CALL display_h_min_s(time_est)
+  ! WRITE(*,*) '... reinitializing'
+  ! CALL inital
+  ! WRITE(*,'(a/)') '... done'
+
+
   !                   1.6     Initial diagnostics
   WRITE(*,*) 'Initial diagnostics...'
   CALL diagnose(0)
   WRITE(*,'(a/)') '...initial diagnostics done'
   !
   CALL FLUSH(stdout)
   !________________________________________________________________________________
 
   WRITE(*,*) 'Time integration loop..'
   !________________________________________________________________________________
   !              2.   Main loop
   DO
      step  = step  + 1
      cstep = cstep + 1
 
      CALL stepon
      time  = time  + dt
 
      CALL tesend
      IF( nlend ) EXIT ! exit do loop
 
      CALL diagnose(step)
 
   END DO
   WRITE(*,'(a/)') '...time integration done'
   !________________________________________________________________________________
   !              9.   Epilogue
 
   CALL diagnose(-1)
   CALL endrun
 
   CALL daytim('Done at ')
   CALL ppexit
 
 END SUBROUTINE control
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index e975693..8c91e56 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -1,405 +1,399 @@
 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
   USE model
   USE initial_par
   USE fields
   USE time_integration
 
   USE prec_const
   IMPLICIT NONE
 
   INCLUDE 'srcinfo.h'
 
   INTEGER, INTENT(in) :: kstep
   INTEGER, parameter  :: BUFSIZE = 2
   INTEGER :: rank, dims(1) = (/0/)
   CHARACTER(len=256) :: str, fname
 
   !_____________________________________________________________________________
   !                   1.   Initial diagnostics
 
   IF ((kstep .EQ. 0)) THEN
      ! jobnum is either 0 from initialization or some integer values read from chkrst(0)
      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
 
      !                       1.1   Initial run
      IF (write_doubleprecision) THEN
         CALL creatf(resfile, fidres, real_prec='d')
      ELSE
         CALL creatf(resfile, fidres)
      END IF
 
 
      WRITE(*,'(3x,a,a)') TRIM(resfile), ' created'
-     call flush(6)
+     CALL flush(6)
 
 
      !  Data group
      CALL creatg(fidres, "/data", "data")
      CALL creatg(fidres, "/data/var2d", "2d profiles")
      CALL creatg(fidres, "/data/var5d", "5d profiles")
 
 
      ! Initialize counter of number of saves for each category
      IF (cstep==0) THEN
          iframe2d=0
      ENDIF
      CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
      IF (cstep==0) THEN
          iframe5d=0
      END IF
      CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
 
      !  File group
      CALL creatg(fidres, "/files", "files")
      CALL attach(fidres, "/files",  "jobnum", jobnum)
 
      !  var2d group (electro. pot., Ni00 moment)
      rank = 0
      CALL creatd(fidres, rank, dims,  "/data/var2d/time",     "Time t*c_s/R")
      CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
 
      IF (write_Na00) THEN
        CALL creatg(fidres, "/data/var2d/Ne00", "Ne00")
        CALL putarr(fidres, "/data/var2d/Ne00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
        CALL putarr(fidres, "/data/var2d/Ne00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
 
        CALL creatg(fidres, "/data/var2d/Ni00", "Ni00")
        CALL putarr(fidres, "/data/var2d/Ni00/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
        CALL putarr(fidres, "/data/var2d/Ni00/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      IF (write_phi) THEN
        CALL creatg(fidres, "/data/var2d/phi", "phi")
        CALL putarr(fidres, "/data/var2d/phi/coordkr", krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
        CALL putarr(fidres, "/data/var2d/phi/coordkz", kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      !  var5d group (moments)
      rank = 0
      CALL creatd(fidres, rank, dims,  "/data/var5d/time",     "Time t*c_s/R")
      CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
      IF (write_moments) THEN
        CALL creatg(fidres, "/data/var5d/moments_e", "moments_e")
        CALL putarr(fidres,  "/data/var5d/moments_e/coordp", parray_e(ips_e:ipe_e),       "p_e",ionode=0)
        CALL putarr(fidres,  "/data/var5d/moments_e/coordj", jarray_e(ijs_e:ije_e),       "j_e",ionode=0)
        CALL putarr(fidres, "/data/var5d/moments_e/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
        CALL putarr(fidres, "/data/var5d/moments_e/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
 
        CALL creatg(fidres, "/data/var5d/moments_i", "moments_i")
        CALL putarr(fidres,  "/data/var5d/moments_i/coordp", parray_i(ips_i:ipe_i),       "p_i",ionode=0)
        CALL putarr(fidres,  "/data/var5d/moments_i/coordj", jarray_i(ijs_i:ije_i),       "j_i",ionode=0)
        CALL putarr(fidres, "/data/var5d/moments_i/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
        CALL putarr(fidres, "/data/var5d/moments_i/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      IF (write_non_lin) THEN
        CALL creatg(fidres, "/data/var5d/Sepj", "Sepj")
        CALL putarr(fidres,  "/data/var5d/Sepj/coordp", parray_e(ips_e:ipe_e),       "p_e",ionode=0)
        CALL putarr(fidres,  "/data/var5d/Sepj/coordj", jarray_e(ijs_e:ije_e),       "j_e",ionode=0)
        CALL putarr(fidres, "/data/var5d/Sepj/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
        CALL putarr(fidres, "/data/var5d/Sepj/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
 
        CALL creatg(fidres, "/data/var5d/Sipj", "Sipj")
        CALL putarr(fidres,  "/data/var5d/Sipj/coordp", parray_i(ips_i:ipe_i),       "p_i",ionode=0)
        CALL putarr(fidres,  "/data/var5d/Sipj/coordj", jarray_i(ijs_i:ije_i),       "j_i",ionode=0)
        CALL putarr(fidres, "/data/var5d/Sipj/coordkr",    krarray(ikrs:ikre), "kr*rho_s0",ionode=0)
        CALL putarr(fidres, "/data/var5d/Sipj/coordkz",    kzarray(ikzs:ikze), "kz*rho_s0",ionode=0)
      END IF
 
      !  Add input namelist variables as attributes of /data/input, defined in srcinfo.h
      WRITE(*,*) 'VERSION=', VERSION
      WRITE(*,*)  'BRANCH=', BRANCH
      WRITE(*,*)  'AUTHOR=', AUTHOR
      WRITE(*,*)    'HOST=', HOST
 
      WRITE(str,'(a,i2.2)') "/data/input"
 
      rank=0
      CALL creatd(fidres, rank,dims,TRIM(str),'Input parameters')
      CALL attach(fidres, TRIM(str),     "version",  VERSION) !defined in srcinfo.h
      CALL attach(fidres, TRIM(str),      "branch",   BRANCH) !defined in srcinfo.h
      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_iframe2d", iframe2d)
      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 grid_outputinputs(fidres, str)
 
      CALL output_par_outputinputs(fidres, str)
 
      CALL model_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
   !
   IF (kstep .GE. 0) THEN
 
     ! Terminal info
     IF (MOD(cstep, INT(1.0/dt)) == 0) THEN
      WRITE(*,"(F5.0,A,F5.0)") time,"/",tmax
     ENDIF
 
      !                       2.1   0d history arrays
      IF (nsave_0d .NE. 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
      IF (nsave_2d .NE. 0) THEN
         IF (MOD(cstep, nsave_2d) == 0) THEN
            CALL diagnose_2d
         END IF
      END IF
 
      !                       2.4   3d profiles
      IF (nsave_5d .NE. 0) THEN
         IF (MOD(cstep, nsave_5d) == 0) THEN
            CALL diagnose_5d
         END IF
      END IF
 
      !                       2.5   Backups
      nsave_cp = INT(5/dt)
      IF (MOD(cstep, nsave_cp) == 0) THEN
        CALL checkpoint_save
      ENDIF
 
   !_____________________________________________________________________________
   !                   3.   Final diagnostics
 
   ELSEIF (kstep .EQ. -1) THEN
      CALL cpu_time(finish)
      CALL attach(fidres, "/data/input","cpu_time",finish-start)
 
      ! Display computational time cost
-     IF     ( FLOOR((finish-start)/3600.) .GT. 0 ) THEN !display h min s
-       WRITE(*,*) 'CPU Time = ', FLOOR((finish-start)/3600.), '[h]', &
-       FLOOR((finish-start)/3600. - FLOOR((finish-start)/3600.))*60., '[min]',&
-       FLOOR((finish-start)/60.   - FLOOR((finish-start)/3600. - FLOOR((finish-start)/3600.))*60.)*60, '[s]'
-     ELSEIF ( FLOOR((finish-start)/60.)   .GT. 0 ) THEN !display min s
-       WRITE(*,*) 'CPU Time = ', FLOOR((finish-start)/60.), '[min]', &
-       FLOOR((finish-start)/60. - FLOOR(finish-start)/60.)*60., '[s]'
-     ELSE ! display s
-       WRITE(*,*) 'CPU Time = ', FLOOR((finish-start)), '[s]'
-     ENDIF
+     CALL display_h_min_s(finish-start)
+     ! WRITE(*,*) 'Time estimated :'
+     ! CALL display_h_min_s(time_est)
+
      !   Close all diagnostic files
      CALL closef(fidres)
   END IF
 
 END SUBROUTINE diagnose
 
 
 SUBROUTINE diagnose_0d
 
   USE basic
   USE futils, ONLY: append
   USE diagnostics_par
   USE prec_const
 
   IMPLICIT NONE
   WRITE(*,'(a,1x,i7.7,a1,i7.7,20x,a,1pe10.3,10x,a,1pe10.3)') &
           '*** Timestep (this run/total) =', step, '/', cstep, 'Time =', time, 'dt =', dt
   WRITE(*,*)
 
-  ! flush stdout of all ranks. Usually ONLY rank 0 should write, but error messages might be written from other ranks as well
+  ! flush stdout of all ranks. Usually ONLY rank 0 should WRITE, but error messages might be written from other ranks as well
   CALL FLUSH(stdout)
 
 END SUBROUTINE diagnose_0d
 
 
 SUBROUTINE diagnose_2d
 
   USE basic
   USE futils, ONLY: append, getatt, attach, putarrnd
   USE fields
   USE time_integration
   USE diagnostics_par
   USE prec_const
   IMPLICIT NONE
 
   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)
 
   IF (write_phi) THEN
      CALL write_field2d(phi(:,:), 'phi')
   END IF
 
   IF (write_Na00) THEN
     CALL write_field2d(moments_e(1,1,:,:,updatetlevel), 'Ne00')
     CALL write_field2d(moments_i(1,1,:,:,updatetlevel), 'Ni00')
   END IF
 CONTAINS
 
   SUBROUTINE write_field2d(field, text)
     USE futils, ONLY: attach, putarr
     USE grid, ONLY: ikrs,ikre, ikzs,ikze
     USE prec_const
     IMPLICIT NONE
 
     COMPLEX(dp), DIMENSION(ikrs:ikre, ikzs:ikze), INTENT(IN) :: field
     CHARACTER(*), INTENT(IN) :: text
 
     CHARACTER(LEN=50) :: dset_name
 
     WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var2d", TRIM(text), iframe2d
     CALL putarr(fidres, dset_name, field(ikrs:ikre, ikzs:ikze),ionode=0)
 
     CALL attach(fidres, dset_name, "time", time)
 
   END SUBROUTINE write_field2d
 
 END SUBROUTINE diagnose_2d
 
 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
    USE time_integration
    USE diagnostics_par
    USE prec_const
    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_moments) THEN
       CALL write_field5d_e(moments_e(:,:,:,:,updatetlevel), 'moments_e')
       CALL write_field5d_i(moments_i(:,:,:,:,updatetlevel), 'moments_i')
    END IF
 
    IF (write_non_lin) THEN
       CALL write_field5d_e(Sepj, 'Sepj')
       CALL write_field5d_i(Sipj, 'Sipj')
    END IF
 
  CONTAINS
 
    SUBROUTINE write_field5d_e(field, text)
      USE futils, ONLY: attach, putarr
      USE grid,   ONLY: ips_e,ipe_e, ijs_e,ije_e, ikrs,ikre, ikzs,ikze
      USE prec_const
      IMPLICIT NONE
 
      COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
      CHARACTER(*), INTENT(IN) :: text
 
      CHARACTER(LEN=50) :: dset_name
 
      WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
      CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze),ionode=0)
 
      CALL attach(fidres, dset_name, "time", time)
 
    END SUBROUTINE write_field5d_e
 
    SUBROUTINE write_field5d_i(field, text)
       USE futils, ONLY: attach, putarr
       USE grid, ONLY: ips_i,ipe_i, ijs_i,ije_i, ikrs,ikre, ikzs,ikze
       USE prec_const
       IMPLICIT NONE
 
       COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze), INTENT(IN) :: field
       CHARACTER(*), INTENT(IN) :: text
 
       CHARACTER(LEN=50) :: dset_name
 
       WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
       CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze),ionode=0)
 
       CALL attach(fidres, dset_name, "time", time)
 
     END SUBROUTINE write_field5d_i
 
 END SUBROUTINE diagnose_5d
 
 SUBROUTINE checkpoint_save
   USE basic
   USE grid
   USE diagnostics_par
   USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf
   USE model
   USE initial_par
   USE fields
   USE time_integration
   IMPLICIT NONE
 
   WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',jobnum,'.h5'
 
   CALL creatf(rstfile, fidrst, real_prec='d')
   CALL creatg(fidrst, '/Basic', 'Basic data')
   CALL attach(fidrst, '/Basic', 'cstep', cstep)
   CALL attach(fidrst, '/Basic', 'time', time)
   CALL attach(fidrst, '/Basic', 'jobnum', jobnum)
   CALL attach(fidrst, '/Basic', 'dt', dt)
   CALL attach(fidrst, '/Basic', 'iframe2d', iframe2d)
   CALL attach(fidrst, '/Basic', 'iframe5d', iframe5d)
 
   ! Write state of system to restart file
   CALL putarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,&
                                                     ikrs:ikre,ikzs:ikze,1),ionode=0)
   CALL putarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,&
                                                     ikrs:ikre,ikzs:ikze,1),ionode=0)
   CALL closef(fidrst)
   WRITE(*,'(3x,a)') "Checkpoint file "//TRIM(rstfile)//" saved"
 
 END SUBROUTINE checkpoint_save
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index 6f065ec..30a6424 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -1,71 +1,71 @@
 MODULE diagnostics_par
   !   Module for diagnostic parameters
 
   USE prec_const
   IMPLICIT NONE
   PRIVATE
 
   LOGICAL, PUBLIC, PROTECTED :: write_Na00=.TRUE.
   LOGICAL, PUBLIC, PROTECTED :: write_moments=.TRUE.
   LOGICAL, PUBLIC, PROTECTED :: write_phi=.TRUE.
   LOGICAL, PUBLIC, PROTECTED :: write_non_lin=.TRUE.
   LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision=.FALSE.
 
   INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d
   INTEGER, PUBLIC :: nsave_cp = 1e4
 
   !  HDF5 file
-  CHARACTER(len=128), PUBLIC :: resfile0 = "results"   ! Head of main result file name
-  CHARACTER(len=128), PUBLIC :: resfile                ! Main result file
+  CHARACTER(len=256), PUBLIC :: resfile0 = "results"   ! Head of main result file name
+  CHARACTER(len=256), PUBLIC :: resfile                ! Main result file
   INTEGER, PUBLIC            :: job2load               ! jobnum of the checkpoint to load
   INTEGER, PUBLIC            :: fidres                 ! FID for resfile
-  CHARACTER(len=128), PUBLIC :: rstfile0 = "restart"   ! Head of restart file name
-  CHARACTER(len=128), PUBLIC :: rstfile                ! Full restart file
+  CHARACTER(len=256), PUBLIC :: rstfile0 = "restart"   ! Head of restart file name
+  CHARACTER(len=256), PUBLIC :: rstfile                ! Full restart file
   INTEGER, PUBLIC            :: fidrst                 ! FID for restart file
 
   PUBLIC :: output_par_readinputs, output_par_outputinputs
 
 CONTAINS
 
 
   SUBROUTINE output_par_readinputs
     !    Read the input parameters
 
     USE basic, ONLY : lu_in
     USE prec_const
     IMPLICIT NONE
 
     NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d
     NAMELIST /OUTPUT_PAR/ write_Na00, write_moments, write_phi, write_non_lin, write_doubleprecision
     NAMELIST /OUTPUT_PAR/ resfile0, rstfile0, job2load
 
     READ(lu_in,output_par)
 
   END SUBROUTINE output_par_readinputs
 
 
   SUBROUTINE output_par_outputinputs(fidres, str)
     !
     !    Write the input parameters to the results_xx.h5 file
     !
     USE prec_const
     USE futils, ONLY: attach
     IMPLICIT NONE
     INTEGER, INTENT(in) :: fidres
     CHARACTER(len=256), INTENT(in) :: str
 
     CALL attach(fidres, TRIM(str), "write_Na00", write_Na00)
     CALL attach(fidres, TRIM(str), "write_moments", write_moments)
     CALL attach(fidres, TRIM(str), "write_phi", write_phi)
     CALL attach(fidres, TRIM(str), "write_non_lin", write_non_lin)
     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(fidres, TRIM(str), "resfile0", resfile0)
 
   END SUBROUTINE output_par_outputinputs
 
 
 END MODULE diagnostics_par
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 48eaac9..7617963 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -1,128 +1,128 @@
 MODULE fourier
 
   USE prec_const
   USE grid
   use, intrinsic :: iso_c_binding
   implicit none
 
   INCLUDE 'fftw3.f03'
 
   PRIVATE
   PUBLIC :: fft_r2cc
   PUBLIC :: ifft_cc2r
   PUBLIC :: fft2_r2cc
   PUBLIC :: ifft2_cc2r
   PUBLIC :: convolve_2D_F2F
 
   CONTAINS
 
   SUBROUTINE fft_r2cc(fx_in, Fk_out)
 
     IMPLICIT NONE
     REAL(dp),    DIMENSION(Nr),  INTENT(IN) :: fx_in
     COMPLEX(dp), DIMENSION(Nkr), INTENT(OUT):: Fk_out
     integer*8 plan
 
-    call dfftw_plan_dft_r2c_1d(plan,Nr,fx_in,Fk_out,FFTW_FORWARD,FFTW_ESTIMATE)
-    call dfftw_execute_dft_r2c(plan, fx_in, Fk_out)
-    call dfftw_destroy_plan(plan)
+    CALL dfftw_plan_dft_r2c_1d(plan,Nr,fx_in,Fk_out,FFTW_FORWARD,FFTW_ESTIMATE)
+    CALL dfftw_execute_dft_r2c(plan, fx_in, Fk_out)
+    CALL dfftw_destroy_plan(plan)
 
   END SUBROUTINE fft_r2cc
 
 
 
   SUBROUTINE ifft_cc2r(Fk_in, fx_out)
 
     IMPLICIT NONE
     COMPLEX(dp), DIMENSION(Nkr),  INTENT(IN) :: Fk_in
     REAL(dp),    DIMENSION(Nr),   INTENT(OUT):: fx_out
     integer*8 plan
 
-    call dfftw_plan_dft_c2r_1d(plan,Nr,Fk_in,fx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
-    call dfftw_execute_dft_c2r(plan, Fk_in, fx_out)
-    call dfftw_destroy_plan(plan)
+    CALL dfftw_plan_dft_c2r_1d(plan,Nr,Fk_in,fx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
+    CALL dfftw_execute_dft_c2r(plan, Fk_in, fx_out)
+    CALL dfftw_destroy_plan(plan)
     fx_out = fx_out/Nr
 
   END SUBROUTINE ifft_cc2r
 
 
 
   SUBROUTINE fft2_r2cc( ffx_in, FFk_out )
 
       IMPLICIT NONE
       REAL(dp),    DIMENSION(Nr,Nz), INTENT(IN)   :: ffx_in
       COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT):: FFk_out
       integer*8 plan
 
       !!! 2D Forward FFT ________________________!
-      call dfftw_plan_dft_r2c_2d(plan,Nr,Nz,ffx_in,FFk_out,FFTW_FORWARD,FFTW_ESTIMATE)
-      call dfftw_execute_dft_r2c(plan,ffx_in,FFk_out)
-      call dfftw_destroy_plan(plan)
+      CALL dfftw_plan_dft_r2c_2d(plan,Nr,Nz,ffx_in,FFk_out,FFTW_FORWARD,FFTW_ESTIMATE)
+      CALL dfftw_execute_dft_r2c(plan,ffx_in,FFk_out)
+      CALL dfftw_destroy_plan(plan)
 
   END SUBROUTINE fft2_r2cc
 
 
 
   SUBROUTINE ifft2_cc2r( FFk_in, ffx_out )
 
       IMPLICIT NONE
       COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN) :: FFk_in
       REAL(dp),    DIMENSION(Nr,Nz), INTENT(OUT)  :: ffx_out
       COMPLEX(dp), DIMENSION(Nkr,Nkz) :: tmp_c
       integer*8 plan
 
       tmp_c = FFk_in
       !!! 2D Backward FFT ________________________!
-      call dfftw_plan_dft_c2r_2d(plan,Nr,Nz,tmp_c,ffx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
-      call dfftw_execute_dft_c2r(plan,tmp_c,ffx_out)
-      call dfftw_destroy_plan(plan)
+      CALL dfftw_plan_dft_c2r_2d(plan,Nr,Nz,tmp_c,ffx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
+      CALL dfftw_execute_dft_c2r(plan,tmp_c,ffx_out)
+      CALL dfftw_destroy_plan(plan)
       ffx_out = ffx_out/Nr/Nz
 
   END SUBROUTINE ifft2_cc2r
 
 
 
   !!! Convolution 2D Fourier to Fourier
   !   - Compute the convolution using the convolution theorem and MKL
   SUBROUTINE convolve_2D_F2F( F_2D, G_2D, C_2D )
 
     USE grid, ONLY : AA_r, AA_z, Lr, Lz
     IMPLICIT NONE
     COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(IN)  :: F_2D, G_2D  ! input fields
     COMPLEX(dp), DIMENSION(Nkr,Nkz), INTENT(OUT) :: C_2D  ! output convolutioned field
     REAL(dp), DIMENSION(Nr,Nz) :: ff, gg ! iFFT of input fields
     REAL(dp), DIMENSION(Nr,Nz) :: ffgg  ! will receive the product of f*g in Real
     INTEGER :: ikr, ikz
     REAL    :: a_r
 
     ! 2D inverse Fourier transform
     CALL ifft2_cc2r(F_2D,ff);
     CALL ifft2_cc2r(G_2D,gg);
 
     ! Product in physical space
     ffgg = ff * gg;
 
     ! 2D Fourier tranform
     CALL fft2_r2cc(ffgg,C_2D);
 
     ! Anti aliasing (2/3 rule)
     DO ikr = 1,Nkr
       a_r = AA_r(ikr)
       DO ikz = 1,Nkz
          C_2D(ikr,ikz) = C_2D(ikr,ikz) * a_r * AA_z(ikz)
       ENDDO
     ENDDO
 
 END SUBROUTINE convolve_2D_F2F
 
 
 ! Empty set/free routines to switch easily with MKL DFTI
 SUBROUTINE set_descriptors
 
 END SUBROUTINE set_descriptors
 
 SUBROUTINE free_descriptors
 
 END SUBROUTINE free_descriptors
 
 END MODULE fourier
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 125fa4a..9c98a1a 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -1,253 +1,253 @@
 MODULE grid
   ! Grid module for spatial discretization
   USE prec_const
   IMPLICIT NONE
   PRIVATE
 
   !   GRID Namelist
   INTEGER,  PUBLIC, PROTECTED :: pmaxe = 1      ! The maximal electron Hermite-moment computed
   INTEGER,  PUBLIC, PROTECTED :: jmaxe = 1      ! The maximal electron Laguerre-moment computed
   INTEGER,  PUBLIC, PROTECTED :: pmaxi = 1      ! The maximal ion Hermite-moment computed
   INTEGER,  PUBLIC, PROTECTED :: jmaxi = 1      ! The maximal ion Laguerre-moment computed
   INTEGER,  PUBLIC, PROTECTED :: maxj  = 1      ! The maximal ion Laguerre-moment computed
   INTEGER,  PUBLIC, PROTECTED :: Nr    = 16     ! Number of total internal grid points in r
   REAL(dp), PUBLIC, PROTECTED :: Lr    = 1._dp  ! horizontal length of the spatial box
   INTEGER,  PUBLIC, PROTECTED :: Nz    = 16     ! Number of total internal grid points in z
   REAL(dp), PUBLIC, PROTECTED :: Lz    = 1._dp  ! vertical length of the spatial box
   INTEGER,  PUBLIC, PROTECTED :: Nkr   = 8      ! Number of total internal grid points in kr
   REAL(dp), PUBLIC, PROTECTED :: Lkr   = 1._dp  ! horizontal length of the fourier box
   INTEGER,  PUBLIC, PROTECTED :: Nkz   = 16     ! Number of total internal grid points in kz
   REAL(dp), PUBLIC, PROTECTED :: Lkz   = 1._dp  ! vertical length of the fourier box
   REAL(dp), PUBLIC, PROTECTED :: kpar  = 0_dp   ! parallel wave vector component
 
   ! For Orszag filter
   REAL(dp), PUBLIC, PROTECTED :: two_third_krmax
   REAL(dp), PUBLIC, PROTECTED :: two_third_kzmax
 
   ! 1D Antialiasing arrays (2/3 rule)
   REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_r
   REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: AA_z
 
   ! Grids containing position in physical space
   REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: rarray
   REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray
   REAL(dp), PUBLIC, PROTECTED ::  deltar,  deltaz
   INTEGER,  PUBLIC, PROTECTED  ::  irs,  ire,  izs,  ize
   INTEGER,  PUBLIC :: ir,iz ! counters
 
   ! Grids containing position in fourier space
   REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: krarray
   REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kzarray
   REAL(dp), PUBLIC, PROTECTED ::  deltakr, deltakz
   INTEGER,  PUBLIC, PROTECTED ::  ikrs, ikre, ikzs, ikze, a
   INTEGER,  PUBLIC, PROTECTED :: ikr_0, ikz_0 ! Indices of k-grid origin
   INTEGER,  PUBLIC :: ikr, ikz, ip, ij ! counters
 
   ! Grid containing the polynomials degrees
   INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e
   INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_i
   INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_e
   INTEGER,  DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i
   INTEGER, PUBLIC, PROTECTED ::  ips_e,ipe_e, ijs_e,ije_e
   INTEGER, PUBLIC, PROTECTED ::  ips_i,ipe_i, ijs_i,ije_i
 
   ! Public Functions
   PUBLIC :: set_pj
   PUBLIC :: set_rgrid,  set_zgrid
   PUBLIC :: set_krgrid, set_kzgrid
   PUBLIC :: grid_readinputs, grid_outputinputs
   PUBLIC :: bare, bari
 CONTAINS
 
   SUBROUTINE set_pj
     USE prec_const
     IMPLICIT NONE
     INTEGER :: ip, ij
     ips_e = 1; ipe_e = pmaxe + 1
     ips_i = 1; ipe_i = pmaxi + 1
     ALLOCATE(parray_e(ips_e:ipe_e))
     ALLOCATE(parray_i(ips_i:ipe_i))
     DO ip = ips_e,ipe_e; parray_e(ip) = ip-1; END DO
     DO ip = ips_i,ipe_i; parray_i(ip) = ip-1; END DO
 
     ijs_e = 1; ije_e = jmaxe + 1
     ijs_i = 1; ije_i = jmaxi + 1
     ALLOCATE(jarray_e(ijs_e:ije_e))
     ALLOCATE(jarray_i(ijs_i:ije_i))
     DO ij = ijs_e,ije_e; jarray_e(ij) = ij-1; END DO
     DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO
     maxj  = MAX(jmaxi, jmaxe)
   END SUBROUTINE set_pj
 
   SUBROUTINE set_rgrid
     USE prec_const
     IMPLICIT NONE
     INTEGER :: ir
     ! Start and END indices of grid
     irs = 1
     ire = Nr
     ! Grid spacing
     IF ( Nr .GT. 1 ) THEN ! To avoid case with 0 intervals
       deltar = Lr / REAL(Nr-1,dp)
     ELSE
       deltar = 1;
     ENDIF
     ! Discretized r positions ordered as dx*(0 1 2 -3 -2 -1)
     ALLOCATE(rarray(irs:ire))
     IF (NR .GT. 1) THEN
       DO ir = irs,ire
         rarray(ir) = deltar*(MODULO(ir-1,Nr/2)-Nr/2*FLOOR(2.*real(ir-1)/real(Nr)))
       END DO
       rarray(Nr/2+1) = -rarray(Nr/2+1)
     ELSE
       rarray(1) = 0._dp
     ENDIF
 
   END SUBROUTINE set_rgrid
 
   SUBROUTINE set_zgrid
     USE prec_const
     IMPLICIT NONE
     INTEGER :: iz
     ! Start and END indices of grid
     izs = 1
     ize = Nr
     ! Grid spacing
     IF ( Nz .GT. 1 ) THEN ! To avoid case with 0 intervals
       deltaz = Lz / REAL(Nz-1,dp)
     ELSE
       deltaz = 1;
     ENDIF
     ! Discretized r positions ordered as dx*(0 1 2 -3 -2 -1)
     ALLOCATE(zarray(irs:ire))
     DO iz = izs,ize
       zarray(iz) = deltaz*(MODULO(iz-1,Nz/2)-Nz/2*FLOOR(2.*real(iz-1)/real(Nz)))
     END DO
     zarray(Nz/2+1) = -zarray(Nz/2+1)
 
   END SUBROUTINE set_zgrid
 
   SUBROUTINE set_krgrid
     USE prec_const
     IMPLICIT NONE
     INTEGER :: ikr
 
     Nkr = Nr/2+1 ! Defined only on positive kr since fields are real
     ! Start and END indices of grid
     ikrs = 1
     ikre = Nkr
     ! Grid spacings
     deltakr = 2._dp*PI/Nr/deltar
 
     ! Discretized kr positions ordered as dk*(0 1 2)
     ALLOCATE(krarray(ikrs:ikre))
     DO ikr = ikrs,ikre
       krarray(ikr) = REAL(ikr-1,dp) * deltakr
       if (krarray(ikr) .EQ. 0) THEN
         ikr_0 = ikr
       ENDIF
     END DO
 
     ! Orszag 2/3 filter
     two_third_krmax = 2._dp/3._dp*deltakr*Nkr
     ALLOCATE(AA_r(ikrs:ikre))
     DO ikr = ikrs,ikre
       IF ( (krarray(ikr) .GT. -two_third_krmax) .AND. (krarray(ikr) .LT. two_third_krmax) ) THEN
         AA_r(ikr) = 1._dp;
       ELSE
         AA_r(ikr) = 0._dp;
       ENDIF
     END DO
   END SUBROUTINE set_krgrid
 
   SUBROUTINE set_kzgrid
     USE prec_const
     USE model, ONLY : kr0KH, ikz0KH
     IMPLICIT NONE
 
     Nkz = Nz;
     ! Start and END indices of grid
     ikzs = 1
     ikze = Nkz
     ! Grid spacings
     deltakz = 2._dp*PI/Nkz/deltaz
 
     ! Discretized kz positions ordered as dk*(0 1 2 -3 -2 -1)
     ALLOCATE(kzarray(ikzs:ikze))
     DO ikz = ikzs,ikze
       kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz)))
       if (kzarray(ikz) .EQ. 0) THEN
         ikz_0 = ikz
       ENDIF
     END DO
     kzarray(Nz/2+1) = -kzarray(Nz/2+1)
 
     ! Orszag 2/3 filter
     two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2);
     ALLOCATE(AA_z(ikzs:ikze))
     DO ikz = ikzs,ikze
       IF ( (kzarray(ikz) .GT. -two_third_kzmax) .AND. (kzarray(ikz) .LT. two_third_kzmax) ) THEN
         AA_z(ikz) = 1._dp;
       ELSE
         AA_z(ikz) = 0._dp;
       ENDIF
     END DO
 
     ! Put kz0RT to the nearest grid point on kz
     ikz0KH = NINT(kr0KH/deltakr)+1
     kr0KH  = kzarray(ikz0KH)
-    write(*,*) 'ikz0KH = ', ikz0KH
-    write(*,*) 'kr0KH = ', kr0KH
+    WRITE(*,*) 'ikz0KH = ', ikz0KH
+    WRITE(*,*) 'kr0KH = ', kr0KH
 
   END SUBROUTINE set_kzgrid
 
   SUBROUTINE grid_readinputs
     ! Read the input parameters
     USE prec_const
     IMPLICIT NONE
     INTEGER :: lu_in   = 90              ! File duplicated from STDIN
 
     NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
                     Nr,  Lr,  Nz,  Lz, kpar
     READ(lu_in,grid)
 
   END SUBROUTINE grid_readinputs
 
 
   SUBROUTINE grid_outputinputs(fidres, str)
     ! Write the input parameters to the results_xx.h5 file
 
     USE futils, ONLY: attach
 
     USE prec_const
     IMPLICIT NONE
 
     INTEGER, INTENT(in) :: fidres
     CHARACTER(len=256), INTENT(in) :: str
     CALL attach(fidres, TRIM(str), "pmaxe", pmaxe)
     CALL attach(fidres, TRIM(str), "jmaxe", jmaxe)
     CALL attach(fidres, TRIM(str), "pmaxi", pmaxi)
     CALL attach(fidres, TRIM(str), "jmaxi", jmaxi)
     CALL attach(fidres, TRIM(str),   "nr",   nr)
     CALL attach(fidres, TRIM(str),   "Lr",   Lr)
     CALL attach(fidres, TRIM(str),   "nz",   nz)
     CALL attach(fidres, TRIM(str),   "Lz",   Lz)
     CALL attach(fidres, TRIM(str),   "nkr",   nkr)
     CALL attach(fidres, TRIM(str),   "Lkr",   Lkr)
     CALL attach(fidres, TRIM(str),   "nkz",   nkz)
     CALL attach(fidres, TRIM(str),   "Lkz",   Lkz)
   END SUBROUTINE grid_outputinputs
 
   FUNCTION bare(p_,j_)
     IMPLICIT NONE
     INTEGER :: bare, p_, j_
     bare = (jmaxe+1)*p_ + j_ + 1
   END FUNCTION
 
   FUNCTION bari(p_,j_)
     IMPLICIT NONE
     INTEGER :: bari, p_, j_
     bari = (jmaxi+1)*p_ + j_ + 1
   END FUNCTION
 
 END MODULE grid
diff --git a/src/inital.F90 b/src/inital.F90
index b3d9fc9..f1fa2ef 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -1,253 +1,237 @@
 !******************************************************************************!
 !!!!!! initialize the moments and load/build coeff tables
 !******************************************************************************!
 SUBROUTINE inital
 
   USE basic
   USE model, ONLY : CO, NON_LIN
   USE prec_const
   USE coeff
   USE time_integration
   USE array, ONLY : Sepj,Sipj
 
   implicit none
 
+  real :: start_init, end_init, time_estimation
+
   CALL set_updatetlevel(1)
   !!!!!! Set the moments arrays Nepj, Nipj !!!!!!
-  write(*,*) 'Init moments'
+  ! WRITE(*,*) 'Init moments'
   IF ( RESTART ) THEN
     CALL load_cp
   ELSE
     CALL init_moments
   ENDIF
   !!!!!! Set phi !!!!!!
-  write(*,*) 'Init phi'
+  ! WRITE(*,*) 'Init phi'
   CALL poisson
 
   !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
   IF ( NON_LIN .OR. (A0KH .NE. 0)) THEN;
-    write(*,*) 'Init Sapj'
+    ! WRITE(*,*) 'Init Sapj'
     CALL compute_Sapj
-    WRITE(*,*) 'Building Dnjs table'
+
+    ! WRITE(*,*) 'Building Dnjs table'
     CALL build_dnjs_table
   ENDIF
   !!!!!! Load the full coulomb collision operator coefficients !!!!!!
   IF (CO .EQ. -1) THEN
-    WRITE(*,*) '=== Load Full Coulomb matrix ==='
+    ! WRITE(*,*) '=== Load Full Coulomb matrix ==='
     CALL load_FC_mat
   ENDIF
+
 END SUBROUTINE inital
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Initialize the moments randomly
 !******************************************************************************!
 SUBROUTINE init_moments
   USE basic
   USE grid
   USE fields
   USE prec_const
   USE utility, ONLY: checkfield
   USE initial_par
   IMPLICIT NONE
 
   REAL(dp) :: noise
   REAL(dp) :: kr, kz, sigma, gain, kz_shift
   INTEGER, DIMENSION(12) :: iseedarr
 
   ! Seed random number generator
   iseedarr(:)=iseed
   CALL RANDOM_SEED(PUT=iseedarr)
 
   IF ( only_Na00 ) THEN ! Spike initialization on density only
 
     DO ikr=ikrs,ikre
       DO ikz=ikzs,ikze
         moments_e( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
         moments_i( 1,1, ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
       END DO
     END DO
 
   ELSE
-    sigma    = 5._dp     ! Gaussian sigma
-    gain     = 0.5_dp    ! Gaussian mean
-    kz_shift = 0.5_dp    ! Gaussian z shift
     !**** Gaussian initialization (Hakim 2017) *********************************
+    ! sigma    = 5._dp     ! Gaussian sigma
+    ! gain     = 0.5_dp    ! Gaussian mean
+    ! kz_shift = 0.5_dp    ! Gaussian z shift
     ! moments_i = 0; moments_e = 0;
     !   DO ikr=ikrs,ikre
     !     kr = krarray(ikr)
     !     DO ikz=ikzs,ikze
     !       kz = kzarray(ikz)
     !       moments_i( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+(kz-kz_shift)**2)*sigma**2/4._dp)
     !       moments_e( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+(kz-kz_shift)**2)*sigma**2/4._dp)
     !     END DO
     !   END DO
 
     !**** Broad noise initialization *******************************************
     DO ip=ips_e,ipe_e
       DO ij=ijs_e,ije_e
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_e( ip,ij,     ikr,    ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) &
-                                                    * AA_r(ikr) * AA_z(ikz)
+            moments_e( ip,ij,     ikr,    ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) !&
+                                                    ! * AA_r(ikr) * AA_z(ikz)
           END DO
         END DO
         DO ikz=2,Nkz/2 !symmetry at kr = 0
           CALL RANDOM_NUMBER(noise)
           moments_e( ip,ij,1,ikz, :) = moments_e( ip,ij,1,Nkz+2-ikz, :)
         END DO
       END DO
     END DO
     DO ip=ips_i,ipe_i
       DO ij=ijs_i,ije_i
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) &
-                                                    * AA_r(ikr) * AA_z(ikz)
+            moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) !&
+                                                    ! * AA_r(ikr) * AA_z(ikz)
           END DO
         END DO
         DO ikz=2,Nkz/2 !symmetry at kr = 0
           CALL RANDOM_NUMBER(noise)
           moments_i( ip,ij,1,ikz, :) = moments_i( ip,ij,1,Nkz+2-ikz, :)
         END DO
       END DO
     END DO
 
-    !**** Sinusoidal phi initialization for Kelvin-Helmholtz *******************
-  !   phi(1,FLOOR(0.5_dp/deltakr)) = 1._dp ! Trigger only mode kr = 1
-  !   ! moments_e( :,:, 1,FLOOR(0.5_dp/deltakr), :) = 1._dp
-  !   ! moments_i( :,:, 1,FLOOR(0.5_dp/deltakr), :) = 1._dp
-  !
-  !   DO ikr=ikrs,ikre
-  !     DO ikz=ikzs,ikze
-  !       CALL RANDOM_NUMBER(noise)
-  !       moments_e( :,:, ikr,ikz, :) = moments_e( :,:, ikr,ikz, :) &
-  !       + initnoise_moments*(noise-0.5_dp) ! adding noise
-  !
-  !       CALL RANDOM_NUMBER(noise)
-  !       moments_i( :,:, ikr,ikz, :) = moments_i( :,:, ikr,ikz, :) &
-  !       + initnoise_moments*(noise-0.5_dp) ! adding noise
-  !
-  !       CALL RANDOM_NUMBER(noise)
-  !       phi(ikr,ikz) = phi(ikr,ikz) + initnoise_moments*(noise-0.5_dp) ! adding noise
-  !     ENDDO
-  !   ENDDO
-  !   CALL moments_eq_rhs
   ENDIF
 END SUBROUTINE init_moments
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Load moments from a previous save
 !******************************************************************************!
 SUBROUTINE load_cp
   USE basic
   USE futils,          ONLY: openf, closef, getarr, getatt
   USE grid
   USE fields
   USE diagnostics_par
   USE time_integration
   IMPLICIT NONE
 
   WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',job2load,'.h5'
 
   WRITE(*,'(3x,a)') "Resume from previous run"
 
   CALL openf(rstfile, fidrst)
   CALL getatt(fidrst, '/Basic', 'cstep', cstep)
   CALL getatt(fidrst, '/Basic', 'time', time)
   CALL getatt(fidrst, '/Basic', 'jobnum', jobnum)
   jobnum = jobnum+1
   CALL getatt(fidrst, '/Basic', 'iframe2d',iframe2d)
   CALL getatt(fidrst, '/Basic', 'iframe5d',iframe5d)
   iframe2d = iframe2d-1; iframe5d = iframe5d-1
 
   ! Read state of system from restart file
   CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0)
   CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0)
   CALL closef(fidrst)
 
   WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!"
 
 END SUBROUTINE load_cp
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Build the Laguerre-Laguerre coupling coefficient table
 !******************************************************************************!
 SUBROUTINE build_dnjs_table
   USE basic
   USE array, Only : dnjs
   USE grid, Only : jmaxe, jmaxi
   USE coeff
   IMPLICIT NONE
 
   INTEGER :: in, ij, is, J
   INTEGER :: n_, j_, s_
 
   J = max(jmaxe,jmaxi)
 
   DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry
     n_ = in - 1
     DO ij = in,J+1
       j_ = ij - 1
       DO is = ij,J+1
         s_ = is - 1
 
         dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0))
         ! By symmetry
         dnjs(in,is,ij) = dnjs(in,ij,is)
         dnjs(ij,in,is) = dnjs(in,ij,is)
         dnjs(ij,is,in) = dnjs(in,ij,is)
         dnjs(is,ij,in) = dnjs(in,ij,is)
         dnjs(is,in,ij) = dnjs(in,ij,is)
       ENDDO
     ENDDO
   ENDDO
 END SUBROUTINE
 !******************************************************************************!
 
 !******************************************************************************!
 !!!!!!! Load the Full coulomb coefficient table from COSOlver results
 !******************************************************************************!
 SUBROUTINE load_FC_mat ! Works only for a partiular file for now with P,J <= 15,6
   USE basic
   USE grid
   USE initial_par
   USE array
   use model
   USE futils, ONLY : openf, getarr, closef
   IMPLICIT NONE
 
   INTEGER :: ip2,ij2
   INTEGER :: fid1, fid2, fid3, fid4
 
   !!!!!!!!!!!! Electron matrices !!!!!!!!!!!!
   ! get the self electron colision matrix
   CALL openf(selfmat_file,fid1, 'r', 'D');
   CALL getarr(fid1, '/Caapj/Ceepj', Ceepj) ! get array (moli format)
   CALL closef(fid1)
   ! get the Test and Back field electron ion collision matrix
   CALL openf(eimat_file,fid2, 'r', 'D');
   CALL getarr(fid2, '/Ceipj/CeipjT', CeipjT)
   CALL getarr(fid2, '/Ceipj/CeipjF', CeipjF)
   CALL closef(fid2)
 
   !!!!!!!!!!!!!!! Ion matrices !!!!!!!!!!!!!!
   ! get the self electron colision matrix
   CALL openf(selfmat_file, fid3, 'r', 'D');
   IF ( (pmaxe .EQ. pmaxi) .AND. (jmaxe .EQ. jmaxi) ) THEN ! if same degrees ion and electron matrices are alike so load Ceepj
     CALL getarr(fid3, '/Caapj/Ceepj', Ciipj) ! get array (moli format)
   ELSE
     CALL getarr(fid3, '/Caapj/Ciipj', Ciipj) ! get array (moli format)
   ENDIF
   CALL closef(fid3)
   ! get the Test and Back field ion electron collision matrix
   CALL openf(iemat_file, fid4, 'r', 'D');
   CALL getarr(fid4, '/Ciepj/CiepjT', CiepjT)
   CALL getarr(fid4, '/Ciepj/CiepjF', CiepjF)
   CALL closef(fid4)
 END SUBROUTINE load_FC_mat
 !******************************************************************************!
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index bdecea0..35cc155 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -1,444 +1,444 @@
 SUBROUTINE moments_eq_rhs
 
   USE basic
   USE time_integration
   USE array
   USE fields
   USE grid
   USE model
   USE prec_const
   USE utility, ONLY : is_nan
   IMPLICIT NONE
 
   INTEGER     :: ip2, ij2 ! loops indices
   REAL(dp)    :: ip_dp, ij_dp
   REAL(dp)    :: kr, kz, kperp2
   REAL(dp)    :: taue_qe_etaB, taui_qi_etaB
   REAL(dp)    :: sqrtTaue_qe, sqrtTaui_qi, qe_sigmae_sqrtTaue, qi_sigmai_sqrtTaui
   REAL(dp)    :: kernelj, kerneljp1, kerneljm1, be_2, bi_2 ! Kernel functions and variable
   REAL(dp)    :: factj, sigmae2_taue_o2, sigmai2_taui_o2 ! Auxiliary variables
   REAL(dp)    :: xNapj, xNapp1j, xNapm1j, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! Mom. factors depending on the pj loop
   REAL(dp)    :: xphij, xphijp1, xphijm1, xphijpar ! ESpot. factors depending on the pj loop
   REAL(dp)    :: xCapj,   xCa20,   xCa01, xCa10 ! Coll. factors depending on the pj loop
   COMPLEX(dp) :: TNapj, TNapp1j, TNapm1j, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
   COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
   COMPLEX(dp) :: test_nan
   REAL(dp)    :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency
 
   !Precompute species dependant factors
   IF( q_e .NE. 0._dp ) THEN
     taue_qe_etaB    = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
     sqrtTaue_qe     = sqrt(tau_e)/q_e   ! factor of parallel moment term
   ELSE
     taue_qe_etaB  = 0._dp
     sqrtTaue_qe   = 0._dp
   ENDIF
 
   taui_qi_etaB    = tau_i/q_i * eta_B ! factor of the magnetic moment coupling
   sqrtTaui_qi     = sqrt(tau_i)/q_i   ! factor of parallel moment term
   qe_sigmae_sqrtTaue = q_e/sigma_e/SQRT(tau_e) ! factor of parallel phi term
   qi_sigmai_sqrtTaui = q_i/sigma_i/SQRT(tau_i)
   sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
   sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp
   nu_e  = nu ! electron-ion collision frequency (where already multiplied by 0.532)
   nu_i  = nu * sigma_e * (tau_i)**(-3._dp/2._dp)/SQRT2 ! ion-ion collision frequ.
   nu_ee  = nu_e/SQRT2 ! e-e coll. frequ.
   nu_ie  = nu*sigma_e**2 ! i-e coll. frequ.
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!! Electrons moments RHS !!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1
     ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxe)
 
     ! N_e^{p+1,j} coeff
     xNapp1j = sqrtTaue_qe * SQRT(ip_dp + 1)
     ! N_e^{p-1,j} coeff
     xNapm1j = sqrtTaue_qe * SQRT(ip_dp)
 
     ! N_e^{p+2,j} coeff
     xNapp2j = taue_qe_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
     ! N_e^{p-2,j} coeff
     xNapm2j = taue_qe_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
 
     factj = 1.0 ! Start of the recursive factorial
 
     jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1
       ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxe)
 
       ! N_e^{p,j+1} coeff
       xNapjp1 = -taue_qe_etaB * (ij_dp + 1._dp)
       ! N_e^{p,j-1} coeff
       xNapjm1 = -taue_qe_etaB * ij_dp
       ! N_e^{pj} coeff
       xNapj   = taue_qe_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
 
       !! Collision operator pj terms
       xCapj = -nu_e*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
       ! Dougherty part
       IF ( CO .EQ. -2) THEN
         IF     ((ip .EQ. 3) .AND. (ij .EQ. 1)) THEN ! kronecker pj20
           xCa20 = nu_e * 2._dp/3._dp
           xCa01 = -SQRT2 * xCa20
           xCa10 = 0._dp
         ELSEIF ((ip .EQ. 1) .AND. (ij .EQ. 2)) THEN ! kronecker pj01
           xCa20 = -nu_e * SQRT2 * 2._dp/3._dp
           xCa01 = -SQRT2 * xCa20
           xCa10 = 0._dp
         ELSEIF ((ip .EQ. 2) .AND. (ij .EQ. 1)) THEN ! kronecker pj10
           xCa20 = 0._dp
           xCa01 = 0._dp
           xCa10 = nu_e
         ELSE
           xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp
         ENDIF
       ENDIF
 
       !! Electrostatic potential pj terms
       IF (ip .EQ. 1) THEN ! kronecker p0
         xphij   =  (eta_n + 2.*ij_dp*eta_T - 2._dp*eta_B*(ij_dp+1._dp) )
         xphijp1 = -(eta_T - eta_B)*(ij_dp+1._dp)
         xphijm1 = -(eta_T - eta_B)* ij_dp
         xphijpar= 0._dp
       ELSE IF (ip .EQ. 2) THEN ! kronecker p1
         xphijpar =  qe_sigmae_sqrtTaue
         xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp
       ELSE IF (ip .EQ. 3) THEN ! kronecker p2
         xphij   =  (eta_T/SQRT2 - SQRT2*eta_B)
         xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar= 0._dp
       ELSE
         xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar= 0._dp
       ENDIF
 
       ! Recursive factorial
       IF (ij_dp .GT. 0) THEN
         factj = factj * ij_dp
       ELSE
         factj = 1._dp
       ENDIF
 
       ! Loop on kspace
       krloope : DO ikr = ikrs,ikre
         kzloope : DO ikz = ikzs,ikze
           kr     = krarray(ikr)   ! Poloidal wavevector
           kz     = kzarray(ikz)   ! Toroidal wavevector
           kperp2 = kr**2 + kz**2  ! perpendicular wavevector
           be_2   = kperp2 * sigmae2_taue_o2 ! Kernel argument
 
           !! Compute moments and mixing terms
           ! term propto N_e^{p,j}
           TNapj = xNapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
           ! term propto N_e^{p+1,j}
           IF (ip+1 .LE. pmaxe+1) THEN ! OoB check
             TNapp1j = xNapp1j * moments_e(ip+1,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapp1j = 0._dp
           ENDIF
           ! term propto N_e^{p-1,j}
           IF (ip-1 .GE. 1) THEN ! OoB check
             TNapm1j = xNapm1j * moments_e(ip-1,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapm1j = 0._dp
           ENDIF
           ! term propto N_e^{p+2,j}
           IF (ip+2 .LE. pmaxe+1) THEN ! OoB check
             TNapp2j = xNapp2j * moments_e(ip+2,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapp2j = 0._dp
           ENDIF
           ! term propto N_e^{p-2,j}
           IF (ip-2 .GE. 1) THEN ! OoB check
             TNapm2j = xNapm2j * moments_e(ip-2,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapm2j = 0._dp
           ENDIF
           ! xterm propto N_e^{p,j+1}
           IF (ij+1 .LE. jmaxe+1) THEN ! OoB check
             TNapjp1 = xNapjp1 * moments_e(ip,ij+1,ikr,ikz,updatetlevel)
           ELSE
             TNapjp1 = 0._dp
           ENDIF
           ! term propto N_e^{p,j-1}
           IF (ij-1 .GE. 1) THEN ! OoB check
             TNapjm1 = xNapjm1 * moments_e(ip,ij-1,ikr,ikz,updatetlevel)
           ELSE
             TNapjm1 = 0._dp
           ENDIF
 
           !! Collision
           IF (CO .EQ. -2) THEN ! Dougherty Collision terms
             IF ( (pmaxe .GE. 2) ) THEN ! OoB check
               TColl20 = xCa20 * moments_e(3,1,ikr,ikz,updatetlevel)
             ELSE
               TColl20 = 0._dp
             ENDIF
             IF ( (jmaxe .GE. 1) ) THEN ! OoB check
               TColl01 = xCa01 * moments_e(1,2,ikr,ikz,updatetlevel)
             ELSE
               TColl01 = 0._dp
             ENDIF
             IF ( (pmaxe .GE. 1) ) THEN ! OoB check
               TColl10 = xCa10 * moments_e(2,1,ikr,ikz,updatetlevel)
             ELSE
               TColl10 = 0._dp
             ENDIF
 
             ! Total collisional term
             TColl =  xCapj* moments_e(ip,ij,ikr,ikz,updatetlevel)&
                    + TColl20 + TColl01 + TColl10
 
           ELSEIF (CO .EQ. -1) THEN ! Full Coulomb for electrons (COSOlver matrix)
 
             TColl = 0._dp ! Initialization
 
             ploopee: DO ip2 = 1,pmaxe+1 ! sum the electron-self and electron-ion test terms
               jloopee: DO ij2 = 1,jmaxe+1
                 TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
                    *( nu_e  * CeipjT(bare(ip-1,ij-1), bare(ip2-1,ij2-1)) &
                      +nu_ee * Ceepj (bare(ip-1,ij-1), bare(ip2-1,ij2-1)))
               ENDDO jloopee
             ENDDO ploopee
 
             ploopei: DO ip2 = 1,pmaxi+1 ! sum the electron-ion field terms
               jloopei: DO ij2 = 1,jmaxi+1
                 TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
                   *(nu_e * CeipjF(bare(ip-1,ij-1), bari(ip2-1,ij2-1)))
               END DO jloopei
             ENDDO ploopei
 
           ELSEIF (CO .EQ. 0) THEN ! Lenard Bernstein
             TColl = xCapj * moments_e(ip,ij,ikr,ikz,updatetlevel)
           ENDIF
 
           !! Electrical potential term
           IF ( (ip .LE. 3) ) THEN ! kronecker p0 p1 p2
             kernelj    = be_2**(ij-1) * exp(-be_2)/factj
             kerneljp1  = kernelj * be_2  /(ij_dp + 1._dp)
             IF ( be_2 .NE. 0 ) THEN
               kerneljm1  = kernelj * ij_dp / be_2
             ELSE
               kerneljm1 = 0._dp
             ENDIF
             Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
           ELSE
             Tphi = 0._dp
           ENDIF
 
           ! Sum of all linear terms
           moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
               -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
               -imagu * kpar*(TNapp1j + TNapm1j + xphijpar*kernelj*phi(ikr,ikz)) &
               - mu*kperp2**2 * moments_e(ip,ij,ikr,ikz,updatetlevel) &
               + TColl
 
           ! Adding non linearity
           IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN
             moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
-              moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz)*Nkr*Nkz
+              moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz)
           ENDIF
 
         END DO kzloope
       END DO krloope
 
     END DO jloope
   END DO ploope
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!! Ions moments RHS !!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ploopi : DO ip = ips_i, ipe_i  ! This loop is from 1 to pmaxi+1
     ip_dp = REAL(ip-1,dp) ! REAL index is one minus the loop index (0 to pmaxi)
 
     ! N_i^{p+1,j} coeff
     xNapp1j = sqrtTaui_qi * SQRT(ip_dp + 1)
     ! N_i^{p-1,j} coeff
     xNapm1j = sqrtTaui_qi * SQRT(ip_dp)
 
     ! x N_i^{p+2,j} coeff
     xNapp2j = taui_qi_etaB * SQRT((ip_dp + 1._dp) * (ip_dp + 2._dp))
     ! x N_i^{p-2,j} coeff
     xNapm2j = taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1._dp))
 
     factj = 1._dp ! Start of the recursive factorial
 
     jloopi : DO ij = ijs_i, ije_i  ! This loop is from 1 to jmaxi+1
       ij_dp = REAL(ij-1,dp) ! REAL index is one minus the loop index (0 to jmaxi)
 
       ! x N_i^{p,j+1} coeff
       xNapjp1 = -taui_qi_etaB * (ij_dp + 1._dp)
       ! x N_i^{p,j-1} coeff
       xNapjm1 = -taui_qi_etaB * ij_dp
       ! x N_i^{pj} coeff
       xNapj   = taui_qi_etaB * 2._dp*(ip_dp + ij_dp + 1._dp)
 
       !! Collision operator pj terms
       xCapj = -nu_i*(ip_dp + 2._dp*ij_dp) !DK Lenard-Bernstein basis
       ! Dougherty part
       IF ( CO .EQ. -2) THEN
         IF     ((ip .EQ. 3) .AND. (ij .EQ. 1)) THEN ! kronecker pj20
           xCa20 = nu_i * 2._dp/3._dp
           xCa01 = -SQRT2 * xCa20
           xCa10 = 0._dp
         ELSEIF ((ip .EQ. 1) .AND. (ij .EQ. 2)) THEN ! kronecker pj01
           xCa20 = -nu_i * SQRT2 * 2._dp/3._dp
           xCa01 = -SQRT2 * xCa20
           xCa10 = 0._dp
         ELSEIF ((ip .EQ. 2) .AND. (ij .EQ. 1)) THEN
           xCa20 = 0._dp
           xCa01 = 0._dp
           xCa10 = nu_i
         ELSE
           xCa20 = 0._dp; xCa01 = 0._dp; xCa10 = 0._dp
         ENDIF
       ENDIF
 
       !! Electrostatic potential pj terms
       IF (ip .EQ. 1) THEN ! krokecker p0
         xphij   =  (eta_n + 2._dp*ij_dp*eta_T - 2._dp*eta_B*(ij_dp+1._dp))
         xphijp1 = -(eta_T - eta_B)*(ij_dp+1._dp)
         xphijm1 = -(eta_T - eta_B)* ij_dp
         xphijpar =  0._dp
       ELSE IF (ip .EQ. 2) THEN ! kronecker p1
         xphijpar =  qi_sigmai_sqrtTaui
         xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp
       ELSE IF (ip .EQ. 3) THEN !krokecker p2
         xphij   =  (eta_T/SQRT2 - SQRT2*eta_B)
         xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar =  0._dp
       ELSE
         xphij = 0._dp; xphijp1 = 0._dp; xphijm1 = 0._dp; xphijpar =  0._dp
       ENDIF
 
       ! Recursive factorial
       IF (ij_dp .GT. 0) THEN
         factj = factj * ij_dp
       ELSE
         factj = 1._dp
       ENDIF
 
       ! Loop on kspace
       krloopi : DO ikr = ikrs,ikre
         kzloopi : DO ikz = ikzs,ikze
           kr     = krarray(ikr)   ! Poloidal wavevector
           kz     = kzarray(ikz)   ! Toroidal wavevector
           kperp2 = kr**2 + kz**2  ! perpendicular wavevector
           bi_2   = kperp2 * sigmai2_taui_o2 ! Kernel argument
 
           !! Compute moments and mixing terms
           ! term propto N_i^{p,j}
           TNapj = xNapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
           ! term propto N_i^{p+1,j}
           IF (ip+1 .LE. pmaxi+1) THEN ! OoB check
             TNapp1j = xNapp1j * moments_i(ip+1,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapp1j = 0._dp
           ENDIF
           ! term propto N_i^{p-1,j}
           IF (ip-1 .GE. 1) THEN ! OoB check
             TNapm1j = xNapm1j * moments_i(ip-1,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapm1j = 0._dp
           ENDIF
           ! term propto N_i^{p+2,j}
           IF (ip+2 .LE. pmaxi+1) THEN ! OoB check
             TNapp2j = xNapp2j * moments_i(ip+2,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapp2j = 0._dp
           ENDIF
           ! term propto N_i^{p-2,j}
           IF (ip-2 .GE. 1) THEN ! OoB check
             TNapm2j = xNapm2j * moments_i(ip-2,ij,ikr,ikz,updatetlevel)
           ELSE
             TNapm2j = 0._dp
           ENDIF
           ! xterm propto N_i^{p,j+1}
           IF (ij+1 .LE. jmaxi+1) THEN ! OoB check
             TNapjp1 = xNapjp1 * moments_i(ip,ij+1,ikr,ikz,updatetlevel)
           ELSE
             TNapjp1 = 0._dp
           ENDIF
           ! term propto N_i^{p,j-1}
           IF (ij-1 .GE. 1) THEN ! OoB check
             TNapjm1 = xNapjm1 * moments_i(ip,ij-1,ikr,ikz,updatetlevel)
           ELSE
             TNapjm1 = 0._dp
           ENDIF
 
           !! Collision
           IF (CO .EQ. -2) THEN  ! Dougherty Collision terms
             IF ( (pmaxi .GE. 2) ) THEN ! OoB check
               TColl20 = xCa20 * moments_i(3,1,ikr,ikz,updatetlevel)
             ELSE
               TColl20 = 0._dp
             ENDIF
             IF ( (jmaxi .GE. 1) ) THEN ! OoB check
               TColl01 = xCa01 * moments_i(1,2,ikr,ikz,updatetlevel)
             ELSE
               TColl01 = 0._dp
             ENDIF
             IF ( (pmaxi .GE. 1) ) THEN ! OoB check
               TColl10 = xCa10 * moments_i(2,1,ikr,ikz,updatetlevel)
             ELSE
               TColl10 = 0._dp
             ENDIF
 
             ! Total collisional term
             TColl =  xCapj* moments_i(ip,ij,ikr,ikz,updatetlevel)&
                    + TColl20 + TColl01 + TColl10
 
           ELSEIF (CO .EQ. -1) THEN !!! Full Coulomb for ions (COSOlver matrix) !!!
 
             TColl = 0._dp ! Initialization
 
             ploopii: DO ip2 = 1,pmaxi+1 ! sum the ion-self and ion-electron test terms
               jloopii: DO ij2 = 1,jmaxi+1
                 TColl = TColl + moments_i(ip2,ij2,ikr,ikz,updatetlevel) &
                     *( nu_ie * CiepjT(bari(ip-1,ij-1), bari(ip2-1,ij2-1)) &
                       +nu_i  * Ciipj (bari(ip-1,ij-1), bari(ip2-1,ij2-1)))
               ENDDO jloopii
             ENDDO ploopii
 
             ploopie: DO ip2 = 1,pmaxe+1 ! sum the ion-electron field terms
               jloopie: DO ij2 = 1,jmaxe+1
                 TColl = TColl + moments_e(ip2,ij2,ikr,ikz,updatetlevel) &
                   *(nu_ie * CiepjF(bari(ip-1,ij-1), bare(ip2-1,ij2-1)))
               ENDDO jloopie
             ENDDO ploopie
 
           ELSEIF (CO .EQ. 0) THEN! Lenhard Bernstein
             TColl = xCapj * moments_i(ip,ij,ikr,ikz,updatetlevel)
           ENDIF
 
           !! Electrical potential term
           IF ( (ip .LE. 3) ) THEN ! kronecker p0 p1 p2
             kernelj    = bi_2**(ij-1) * exp(-bi_2)/factj
             kerneljp1  = kernelj * bi_2  /(ij_dp + 1._dp)
             IF ( bi_2 .NE. 0 ) THEN
               kerneljm1  = kernelj * ij_dp / bi_2
             ELSE
               kerneljm1 = 0._dp
             ENDIF
             Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
           ELSE
             Tphi = 0._dp
           ENDIF
 
           ! Sum of linear terms
           moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
               -imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
               -imagu * kpar*(TNapp1j + TNapm1j + xphijpar*kernelj*phi(ikr,ikz)) &
               - mu*kperp2**2 * moments_i(ip,ij,ikr,ikz,updatetlevel) &
                + TColl
-               
+
           ! Adding non linearity
           IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN
            moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
-             moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz)*Nkr*Nkz
+             moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz)
           ENDIF
 
         END DO kzloopi
       END DO krloopi
 
     END DO jloopi
   END DO ploopi
 
 END SUBROUTINE moments_eq_rhs
diff --git a/src/poisson.F90 b/src/poisson.F90
index 760a9fe..c0fdbcf 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -1,94 +1,94 @@
 SUBROUTINE poisson
   ! Solve poisson equation to get phi
 
   USE time_integration, ONLY: updatetlevel
   USE array
   USE fields
   USE grid
   use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, DK
 
   USE prec_const
   IMPLICIT NONE
 
   INTEGER     :: ini,ine, i0j
   REAL(dp)    :: ini_dp, ine_dp
   REAL(dp)    :: kr, kz, kperp2
   REAL(dp)    :: Kne, Kni ! sub kernel factor for recursive build
   REAL(dp)    :: be_2, bi_2, alphaD
   REAL(dp)    :: sigmae2_taue_o2, sigmai2_taui_o2, qe2_taue, qi2_taui ! To avoid redondant computation
   REAL(dp)    :: sum_kernel2_e,    sum_kernel2_i    ! Store sum Kn^2
   COMPLEX(dp) :: sum_kernel_mom_e, sum_kernel_mom_i ! Store sum Kn*Napn
   REAL(dp)    :: gammaD
   COMPLEX(dp) :: gammaD_phi
 
   !Precompute species dependant factors
   sigmae2_taue_o2 = sigma_e**2 * tau_e/2._dp ! factor of the Kernel argument
   sigmai2_taui_o2 = sigma_i**2 * tau_i/2._dp ! (b_a/2)^2 = (kperp sqrt(2 tau_a) sigma_a/2)^2
                                              !           = kperp^2 tau_a sigma_a^2/2
   qe2_taue        = (q_e**2)/tau_e ! factor of the gammaD sum
   qi2_taui        = (q_i**2)/tau_i
 
   DO ikr=ikrs,ikre
     DO ikz=ikzs,ikze
       kr     = krarray(ikr)
       kz     = kzarray(ikz)
       kperp2 = kr**2 + kz**2
       ! non dim kernel argument (kperp2 sigma_a sqrt(2 tau_a))
       be_2  =  kperp2 * sigmae2_taue_o2
       bi_2  =  kperp2 * sigmai2_taui_o2
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       !!!!!!!!!!!!! Electrons sum(Kernel * Ne0n)  (skm) and sum(Kernel**2) (sk2)
       !! Sum(kernel * Moments_0n)
       ! Initialization for n = 0 (ine = 1)
       Kne  = exp(-be_2)
       sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel)
       sum_kernel2_e = Kne**2
       ! loop over n only if the max polynomial degree is 1 or more
       if (jmaxe .GT. 0) then
         DO ine=2,jmaxe+1 ! ine = n+1
           ine_dp = REAL(ine-1,dp)
           ! We update iteratively the kernel functions (to spare factorial computations)
           Kne  = Kne  * be_2/ine_dp
           sum_kernel_mom_e  = sum_kernel_mom_e  + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
           sum_kernel2_e     = sum_kernel2_e     + Kne**2 ! ... sum recursively ...
         END DO
       ENDIF
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       !!!!!!!!!!!!!!!!! Ions sum(Kernel * Ni0n)  (skm) and sum(Kernel**2) (sk2)
       !! Sum(kernel * Moments_0n)
       ! Initialization for n = 0 (ini = 1)
         Kni  = exp(-bi_2)
         sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel)
         sum_kernel2_i = Kni**2
         ! loop over n only if the max polynomial degree is 1 or more
         if (jmaxi .GT. 0) then
           DO ini=2,jmaxi+1
             ini_dp = REAL(ini-1,dp) ! Real index (0 to jmax)
             ! We update iteratively to spare factorial computations
             Kni  = Kni  * bi_2/ini_dp
             sum_kernel_mom_i  = sum_kernel_mom_i  + Kni * moments_i(1, ini, ikr, ikz, updatetlevel)
             sum_kernel2_i     = sum_kernel2_i     + Kni**2 ! ... sum recursively ...
           END DO
         ENDIF
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       !!!!!!!!!!!!!!! Assembling the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
       alphaD   = kperp2 * lambdaD**2
       gammaD   = alphaD + qe2_taue * (1._dp - sum_kernel2_e) & ! Called Poisson_ in MOLI
                         + qi2_taui * (1._dp - sum_kernel2_i)
 
       gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
 
       IF ( (gammaD .EQ. 0) .AND. (abs(kr)+abs(kz) .NE. 0._dp) ) THEN
-        write(*,*) 'Warning gammaD = 0', sum_kernel2_i
+        WRITE(*,*) 'Warning gammaD = 0', sum_kernel2_i
       ENDIF
 
       phi(ikr, ikz) =  gammaD_phi/gammaD
 
     END DO
   END DO
 
 ! Cancel origin singularity
 phi(ikr_0,ikz_0) = 0
 
 END SUBROUTINE poisson
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index a7257cc..7282015 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -1,69 +1,69 @@
 MODULE prec_const
 
   use mpi
 
   use, intrinsic :: iso_fortran_env, only: REAL32, REAL64, &
                                            stdin=>input_unit, &
                                            stdout=>output_unit, &
                                            stderr=>error_unit
 
   ! Precision for real and complex
   INTEGER, PARAMETER :: sp = REAL32 !Single precision, should not be used
   INTEGER, PARAMETER :: dp = REAL64 !real(dp), enforced through the code
 
   INTEGER, private :: dp_r, dp_p !Range and Aprecision of doubles
   INTEGER, private :: sp_r, sp_p !Range and precision of singles
 
 
   INTEGER, private :: MPI_SP !Single precision for MPI
   INTEGER, private :: MPI_DP !Double precision in MPI
   INTEGER, private :: MPI_SUM_DP !Sum reduction operation for DP datatype
   INTEGER, private :: MPI_MAX_DP !Max reduction operation for DP datatype
   INTEGER, private :: MPI_MIN_DP !Min reduction operation for DP datatype
 
 
   ! Some useful constants, to avoid recomputing then too often
   REAL(dp), PARAMETER :: PI=3.141592653589793238462643383279502884197_dp
   REAL(dp), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_dp
   REAL(dp), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_dp
   REAL(dp), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_dp
   REAL(dp), PARAMETER :: INVSQRT2=0.7071067811865475244008443621048490392848359377_dp
   REAL(dp), PARAMETER :: SQRT3=1.73205080756887729352744634150587236694281_dp
   REAL(dp), PARAMETER :: onetwelfth=0.08333333333333333333333333333333333333333333333_dp
   REAL(dp), PARAMETER :: onetwentyfourth=0.04166666666666666666666666666666666666666666667_dp
   REAL(dp), PARAMETER :: onethird=0.33333333333333333333333333333333333333333333333_dp
   REAL(dp), PARAMETER :: onesixth=0.1666666666666666666666666666666666666666666667_dp
   REAL(dp), PARAMETER :: fivesixths=0.8333333333333333333333333333333333333333333333_dp
   REAL(dp), PARAMETER :: sevensixths=1.1666666666666666666666666666666666666666666667_dp
   REAL(dp), PARAMETER :: elevensixths=1.833333333333333333333333333333333333333333333_dp
   REAL(dp), PARAMETER :: nineeighths=1.125_dp
   REAL(dp), PARAMETER :: onesixteenth=0.0625_dp
   REAL(dp), PARAMETER :: ninesixteenths=0.5625_dp
   REAL(dp), PARAMETER :: thirteentwelfths = 1.083333333333333333333333333333333333333333333_dp
   COMPLEX(dp),  PARAMETER :: imagu = (0,1)
   CONTAINS
 
     SUBROUTINE INIT_PREC_CONST
 
       IMPLICIT NONE
       integer :: ierr,me
       
       REAL(sp) :: a = 1_sp
       REAL(dp) :: b = 1_dp
 
       !Get range and precision of ISO FORTRAN sizes
       sp_r = range(a)
       sp_p = precision(a)
 
       dp_r = range(b)
       dp_p = precision(b)
 
-      call mpi_comm_rank(MPI_COMM_WORLD,me,ierr)
+      CALL mpi_comm_rank(MPI_COMM_WORLD,me,ierr)
 
       !Create MPI datatypes that support the specific size
-      call MPI_Type_create_f90_real(sp_p,sp_r,MPI_sp,ierr)
-      call MPI_Type_create_f90_real(dp_p,dp_r,MPI_dp,ierr)
+      CALL MPI_Type_create_f90_real(sp_p,sp_r,MPI_sp,ierr)
+      CALL MPI_Type_create_f90_real(dp_p,dp_r,MPI_dp,ierr)
 
     END SUBROUTINE INIT_PREC_CONST
 
 END MODULE prec_const
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index 44bfcf8..776f900 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -1,32 +1,37 @@
 SUBROUTINE readinputs
   ! Additional data specific for a new run
 
   USE grid,             ONLY: grid_readinputs
   USE diagnostics_par,  ONLY: output_par_readinputs
   USE model,            ONLY: model_readinputs
   USE initial_par,      ONLY: initial_readinputs
   USE time_integration, ONLY: time_integration_readinputs
 
   USE prec_const
   IMPLICIT NONE
 
   !WRITE(*,'(a/)') '=== Define additional data for a new run ==='
 
   ! The input file must be in the same order as the following read routines commands (eg spacegrid then output then model)
 
   ! Load grid data from input file
   CALL grid_readinputs
+  WRITE(*,*) 'check'
 
   ! Load diagnostic options from input file
   CALL output_par_readinputs
+  WRITE(*,*) 'check'
 
   ! Load model parameters from input file
   CALL model_readinputs
+  WRITE(*,*) 'check'
 
   ! Load initial condition parameters from input file
   CALL initial_readinputs
+  WRITE(*,*) 'check'
 
   ! Load parameters for time integration from input file
   CALL time_integration_readinputs
+  WRITE(*,*) 'check'
 
 END SUBROUTINE readinputs
diff --git a/src/stepon.F90 b/src/stepon.F90
index a900af8..f23c5ff 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -1,88 +1,86 @@
 SUBROUTINE stepon
   !   Advance one time step, (num_step=4 for Runge Kutta 4 scheme)
 
   USE basic
   USE time_integration
   USE fields, ONLY: moments_e, moments_i, phi
   USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj
   USE grid
   USE advance_field_routine, ONLY: advance_time_level, advance_field
   USE model
   USE utility, ONLY: checkfield
   use prec_const
   IMPLICIT NONE
 
   INTEGER :: num_step
 
-
    DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
-
       ! Compute right hand side of moments hierarchy equation
       CALL moments_eq_rhs
       ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
       CALL advance_time_level
       ! Update the moments with the hierarchy RHS (step by step)
       DO ip=ips_e,ipe_e
         DO ij=ijs_e,ije_e
           CALL advance_field(moments_e(ip,ij,:,:,:),moments_rhs_e(ip,ij,:,:,:))
         ENDDO
       ENDDO
       DO ip=ips_i,ipe_i
         DO ij=ijs_i,ije_i
           CALL advance_field(moments_i(ip,ij,:,:,:),moments_rhs_i(ip,ij,:,:,:))
         ENDDO
       ENDDO
       ! Update electrostatic potential
       CALL poisson
       ! Update nonlinear term
       IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN
         CALL compute_Sapj
       ENDIF
       !(The two routines above are called in inital for t=0)
 
       CALL enforce_symetry() ! Enforcing symmetry on kr = 0
 
       CALL checkfield_all()
    END DO
 
    CONTAINS
 
       SUBROUTINE checkfield_all ! Check all the fields for inf or nan
         IF(.NOT.nlend) THEN
            nlend=nlend .or. checkfield(phi,' phi')
            DO ip=ips_e,ipe_e
              DO ij=ijs_e,ije_e
               nlend=nlend .or. checkfield(moments_e(ip,ij,:,:,updatetlevel),' moments_e')
              ENDDO
            ENDDO
 
            DO ip=ips_i,ipe_i
              DO ij=ijs_i,ije_i
               nlend=nlend .or. checkfield(moments_i(ip,ij,:,:,updatetlevel),' moments_i')
              ENDDO
            ENDDO
         ENDIF
       END SUBROUTINE checkfield_all
 
       SUBROUTINE enforce_symetry
 
         DO ip=ips_e,ipe_e
           DO ij=ijs_e,ije_e
             DO ikz=2,Nkz/2 !symmetry at kr = 0
               moments_e( ip,ij,1,ikz, :) = CONJG(moments_e( ip,ij,1,Nkz+2-ikz, :))
             END DO
           END DO
         END DO
         DO ip=ips_i,ipe_i
           DO ij=ijs_i,ije_i
             DO ikz=2,Nkz/2 !symmetry at kr = 0
               moments_i( ip,ij,1,ikz, :) = CONJG(moments_i( ip,ij,1,Nkz+2-ikz, :))
             END DO
           END DO
         END DO
         DO ikz=2,Nkz/2 !symmetry at kr = 0
           phi(1,ikz) = phi(1,Nkz+2-ikz)
         END DO
       END SUBROUTINE enforce_symetry
 
 END SUBROUTINE stepon
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 3f47370..03e7098 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -1,227 +1,227 @@
 MODULE time_integration
 
   USE prec_const
   IMPLICIT NONE
   PRIVATE
 
   INTEGER, PUBLIC, PROTECTED :: ntimelevel=4 ! Total number of time levels required by the numerical scheme
   INTEGER, PUBLIC, PROTECTED :: updatetlevel ! Current time level to be updated
 
   real(dp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: A_E,A_I
   real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: b_E,b_Es,b_I
   real(dp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E,c_I !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17
 
   character(len=10),PUBLIC,PROTECTED :: numerical_scheme='RK4'
 
   PUBLIC :: set_updatetlevel, time_integration_readinputs, time_integration_outputinputs
 
 CONTAINS
 
   SUBROUTINE set_updatetlevel(new_updatetlevel)
     INTEGER, INTENT(in) :: new_updatetlevel
     updatetlevel = new_updatetlevel
   END SUBROUTINE set_updatetlevel
 
   SUBROUTINE time_integration_readinputs
     ! Read the input parameters
 
     USE prec_const
     USE basic, ONLY : lu_in
     IMPLICIT NONE
 
     NAMELIST /TIME_INTEGRATION_PAR/ numerical_scheme
 
     READ(lu_in,time_integration_par)
     !WRITE(*,time_integration_par)
 
     CALL set_numerical_scheme
 
   END SUBROUTINE time_integration_readinputs
 
 
   SUBROUTINE time_integration_outputinputs(fidres, str)
     ! Write the input parameters to the results_xx.h5 file
 
     USE prec_const
     USE futils, ONLY: attach
     IMPLICIT NONE
     INTEGER, INTENT(in) :: fidres
     CHARACTER(len=256), INTENT(in) :: str
 
     CALL attach(fidres, TRIM(str), "numerical_scheme", numerical_scheme)
 
   END SUBROUTINE time_integration_outputinputs
 
 
 
   SUBROUTINE set_numerical_scheme
     ! Initialize Butcher coefficient of numerical_scheme
 
     use basic
     IMPLICIT NONE
 
     SELECT CASE (numerical_scheme)
     CASE ('RK4')
       CALL RK4
     CASE ('DOPRI5')
        CALL DOPRI5
     CASE ('DOPRI5_ADAPT')
       CALL DOPRI5_ADAPT
     CASE DEFAULT
        WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
     END SELECT
 
-    write(*,*) " Time integration with ", numerical_scheme
+    WRITE(*,*) " Time integration with ", numerical_scheme
 
   END SUBROUTINE set_numerical_scheme
 
   SUBROUTINE RK4
     ! Butcher coeff for RK4 (default)
 
     USE basic
     USE prec_const
     IMPLICIT NONE
     INTEGER,PARAMETER :: nbstep = 4
 
     CALL allocate_array_dp1(c_E,1,nbstep)
     CALL allocate_array_dp1(b_E,1,nbstep)
     CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
 
     ntimelevel = 4
 
     c_E(1) = 0.0_dp
     c_E(2) = 1.0_dp/2.0_dp
     c_E(3) = 1.0_dp/2.0_dp
     c_E(4) = 1.0_dp
 
     b_E(1) = 1.0_dp/6.0_dp
     b_E(2) = 1.0_dp/3.0_dp
     b_E(3) = 1.0_dp/3.0_dp
     b_E(4) = 1.0_dp/6.0_dp
 
     A_E(2,1) = 1.0_dp/2.0_dp
     A_E(3,2) = 1.0_dp/2.0_dp
     A_E(4,3) = 1.0_dp
 
   END SUBROUTINE RK4
 
 
 
   SUBROUTINE DOPRI5
     ! Butcher coeff for DOPRI5 --> Stiffness detection
     ! DOPRI5 used for stiffness detection.
     ! 5 order method/7 stages
 
     USE basic
     IMPLICIT NONE
     INTEGER,PARAMETER :: nbstep =7
 
     CALL allocate_array_dp1(c_E,1,nbstep)
     CALL allocate_array_dp1(b_E,1,nbstep)
     CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
 
     ntimelevel = 7
     !
     c_E(1) = 0._dp
     c_E(2) = 1.0_dp/5.0_dp
     c_E(3) = 3.0_dp /10.0_dp
     c_E(4) = 4.0_dp/5.0_dp
     c_E(5) = 8.0_dp/9.0_dp
     c_E(6) = 1.0_dp
     c_E(7) = 1.0_dp
     !
     A_E(2,1) = 1.0_dp/5.0_dp
     A_E(3,1) = 3.0_dp/40.0_dp
     A_E(3,2) = 9.0_dp/40.0_dp
     A_E(4,1) = 	44.0_dp/45.0_dp
     A_E(4,2) = -56.0_dp/15.0_dp
     A_E(4,3) = 32.0_dp/9.0_dp
     A_E(5,1 ) = 19372.0_dp/6561.0_dp
     A_E(5,2) = -25360.0_dp/2187.0_dp
     A_E(5,3) = 64448.0_dp/6561.0_dp
     A_E(5,4) = -212.0_dp/729.0_dp
     A_E(6,1) = 9017.0_dp/3168.0_dp
     A_E(6,2)= -355.0_dp/33.0_dp
     A_E(6,3) = 46732.0_dp/5247.0_dp
     A_E(6,4) = 49.0_dp/176.0_dp
     A_E(6,5) = -5103.0_dp/18656.0_dp
     A_E(7,1) = 35.0_dp/384.0_dp
     A_E(7,3) = 500.0_dp/1113.0_dp
     A_E(7,4) = 125.0_dp/192.0_dp
     A_E(7,5) = -2187.0_dp/6784.0_dp
     A_E(7,6) = 11.0_dp/84.0_dp
     !
     b_E(1) = 35.0_dp/384.0_dp
     b_E(2) = 0._dp
     b_E(3) = 500.0_dp/1113.0_dp
     b_E(4) = 125.0_dp/192.0_dp
     b_E(5) = -2187.0_dp/6784.0_dp
     b_E(6) = 11.0_dp/84.0_dp
     b_E(7) = 0._dp
     !
   END SUBROUTINE DOPRI5
 
   SUBROUTINE DOPRI5_ADAPT
     ! Butcher coeff for DOPRI5 --> Stiffness detection
     ! DOPRI5 used for stiffness detection.
     ! 5 order method/7 stages
 
     USE basic
     IMPLICIT NONE
     INTEGER,PARAMETER :: nbstep =7
 
     CALL allocate_array_dp1(c_E,1,nbstep)
     CALL allocate_array_dp1(b_E,1,nbstep)
     CALL allocate_array_dp1(b_Es,1,nbstep)
     CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
 
     ntimelevel = 7
     !
     c_E(1) = 0._dp
     c_E(2) = 1.0_dp/5.0_dp
     c_E(3) = 3.0_dp /10.0_dp
     c_E(4) = 4.0_dp/5.0_dp
     c_E(5) = 8.0_dp/9.0_dp
     c_E(6) = 1.0_dp
     c_E(7) = 1.0_dp
     !
     A_E(2,1) = 1.0_dp/5.0_dp
     A_E(3,1) = 3.0_dp/40.0_dp
     A_E(3,2) = 9.0_dp/40.0_dp
     A_E(4,1) = 	44.0_dp/45.0_dp
     A_E(4,2) = -56.0_dp/15.0_dp
     A_E(4,3) = 32.0_dp/9.0_dp
     A_E(5,1 ) = 19372.0_dp/6561.0_dp
     A_E(5,2) = -25360.0_dp/2187.0_dp
     A_E(5,3) = 64448.0_dp/6561.0_dp
     A_E(5,4) = -212.0_dp/729.0_dp
     A_E(6,1) = 9017.0_dp/3168.0_dp
     A_E(6,2)= -355.0_dp/33.0_dp
     A_E(6,3) = 46732.0_dp/5247.0_dp
     A_E(6,4) = 49.0_dp/176.0_dp
     A_E(6,5) = -5103.0_dp/18656.0_dp
     A_E(7,1) = 35.0_dp/384.0_dp
     A_E(7,3) = 500.0_dp/1113.0_dp
     A_E(7,4) = 125.0_dp/192.0_dp
     A_E(7,5) = -2187.0_dp/6784.0_dp
     A_E(7,6) = 11.0_dp/84.0_dp
     !
     b_E(1) = 35.0_dp/384.0_dp
     b_E(2) = 0._dp
     b_E(3) = 500.0_dp/1113.0_dp
     b_E(4) = 125.0_dp/192.0_dp
     b_E(5) = -2187.0_dp/6784.0_dp
     b_E(6) = 11.0_dp/84.0_dp
     b_E(7) = 0._dp
     !
     b_Es(1) = 5179.0_dp/57600.0_dp
     b_Es(2) = 0._dp
     b_Es(3) = 7571.0_dp/16695.0_dp
     b_Es(4) = 393.0_dp/640.0_dp
     b_Es(5) = -92097.0_dp/339200.0_dp
     b_Es(6) = 187.0_dp/2100.0_dp
     b_Es(7) = 1._dp/40._dp
     !
   END SUBROUTINE DOPRI5_ADAPT
 
 END MODULE time_integration
diff --git a/wk/fort.90 b/wk/fort.90
index ca35060..c849db7 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,63 +1,63 @@
 &BASIC
   nrun   = 100000000
   dt     = 0.01
-  tmax   = 150
+  tmax   = 200
   RESTART = .false.
 /
 &GRID
   pmaxe =2
   jmaxe = 1
   pmaxi = 2
   jmaxi = 1
   Nr   = 256
-  Lr = 40
+  Lr = 66
   Nz   = 256
-  Lz = 40
+  Lz = 66
   kpar = 0
 /
 &OUTPUT_PAR
   nsave_0d = 0
   nsave_1d = 0
-  nsave_2d = 20
-  nsave_5d = 200
+  nsave_2d = 50
+  nsave_5d = 1000
   write_Na00    = .true.
   write_moments = .true.
   write_phi     = .true.
   write_non_lin     = .true.
   write_doubleprecision = .true.
-  resfile0      = 'ZP_forced_sym_256x128_L_40_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-03__mu_1e-04_'
-  rstfile0      = '../checkpoint/cp_ZP_forced_sym_256x128_L_40_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-03__mu_1e-04_'
-  job2load      = 0
+  resfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/outputs'
+  rstfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/checkpoint'
+  job2load      = 1
 /
 &MODEL_PAR
   ! Collisionality
-  CO      = 0
+  CO      = -2
   DK      = .false.
   NON_LIN = .true.
-  mu      = 0.0001
-  nu      = 0.001
+  mu      = 0.0005
+  nu      = 0.01
   tau_e   = 1
   tau_i   = 1
   sigma_e = 0.023338
   sigma_i = 1
   q_e     = -1
   q_i     = 1
   eta_n   = 1
   eta_T   = 0
   eta_B   = 0.5
   lambdaD = 0
   kr0KH   = 0
   A0KH    = 0
 /
 &INITIAL_CON
   only_Na00         =.false.
   initback_moments  =0
-  initnoise_moments =5e-05
+  initnoise_moments =1e-05
   iseed             =42
   selfmat_file      ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10.h5'
   eimat_file        ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5'
   iemat_file        ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5'
 /
 &TIME_INTEGRATION_PAR
   numerical_scheme='RK4'
 /
\ No newline at end of file
diff --git a/wk/parameters_ZP.m b/wk/parameters_ZP.m
index 45a749a..4033fbd 100644
--- a/wk/parameters_ZP.m
+++ b/wk/parameters_ZP.m
@@ -1,44 +1,43 @@
 clear all;
 addpath(genpath('../matlab')) % ... add
 default_plots_options
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1e-3;   % Collision frequency
+NU      = 1e-2;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 ETAB    = 0.5;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
-MU      = 1e-4;   % Hyper diffusivity coefficient
-LAMBDAD = 0.0; 
-NOISE0  = 5.0e-5;
+MU      = 5e-4;   % Hyper diffusivity coefficient
+NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
 N       = 256;     % Frequency gridpoints (Nkr = N/2)
-L       = 40;     % Size of the squared frequency domain
-KREQ0   = 0;      % put kr = 0
-PMAXE   = 02;     % Highest electron Hermite polynomial degree
-JMAXE   = 01;     % Highest ''       Laguerre ''
-PMAXI   = 02;     % Highest ion      Hermite polynomial degree
-JMAXI   = 01;     % Highest ''       Laguerre ''
-KPAR    = 0.0;    % Parellel wave vector component
+L       = 66;     % Size of the squared frequency domain
+PMAXE   = 2;     % Highest electron Hermite polynomial degree
+JMAXE   = 1;     % Highest ''       Laguerre ''
+PMAXI   = 2;     % Highest ion      Hermite polynomial degree
+JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS 
-TMAX    = 150;  % Maximal time unit
-DT    = 1e-2;   % Time step
-SPS2D   = 5;      % Sampling per time unit for 2D arrays
-SPS5D   = 0.5;    % Sampling per time unit for 5D arrays
+TMAX    = 200;  % Maximal time unit
+DT      = 1e-2;   % Time step
+SPS2D   = 2;      % Sampling per time unit for 2D arrays
+SPS5D   = 0.1;    % Sampling per time unit for 5D arrays
 RESTART = 0;      % To restart from last checkpoint
-JOB2LOAD= 00;
+JOB2LOAD= 01;
 %% OPTIONS
-SIMID   = 'ZP_forced_sym';  % Name of the simulation
-NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
-CO      = 0;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+SIMID   = 'ZP';  % Name of the simulation
+CO      = -2;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% unused
+%% unused
 KR0KH   = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
 NO_E    = 0;  % Remove electrons dynamic
 % DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
-
+KREQ0   = 0;      % put kr = 0
+KPAR    = 0.0;    % Parellel wave vector component
+LAMBDAD = 0.0; 
+NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 
 setup