diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index ff48c1a260af3a821729e0d76fe6650d6752a40b..20f65d553291e415de6db3593c4b5023dfb19330 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -36,7 +36,7 @@ CONTAINS
         ! In GENE, there is a minus sign here...
         gamma_E = ExBrate*C_y*abs(Cyq0_x0/C_y)
         IF(abs(gamma_E) .GT. EPSILON(gamma_E)) THEN
-            CALL speak('-ExB background flow detected-')
+            CALL speak('-ExB background flow detected-',2)
             ExB    = .TRUE.
             t0     = deltakx/deltaky/gamma_E
             inv_t0 = 1._xp/t0
diff --git a/src/auxval.F90 b/src/auxval.F90
index c05fe97e56be2f4e3948e16ed20477ed1d25e86d..051228ff11574664e779294f9513898022fdfdc3 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -1,7 +1,7 @@
 subroutine auxval
   !   Set auxiliary values, at beginning of simulation
   USE, intrinsic :: iso_fortran_env, ONLY: OUTPUT_UNIT
-  USE basic,          ONLY: str, speak
+  USE basic,          ONLY: str, speak, VERBOSE_LVL
   USE grid,           ONLY: local_np, local_np_offset, total_np, local_nj, local_nj_offset,&
                             local_nky, local_nky_offset, total_nky, local_nkx, local_nkx_offset,&
                             local_nz, local_nz_offset, init_grids_data, set_grids
@@ -22,7 +22,6 @@ subroutine auxval
 #endif
   IMPLICIT NONE
   INTEGER :: i_, ierr
-  CALL speak('=== Set auxiliary values ===')
   ! Init the grids
   CALL init_grids_data(Na,EM,LINEARITY) 
   CALL set_grids(shear,ExBrate,Npol,LINEARITY,N_HD,N_HDz,EM,Na) 
@@ -54,38 +53,40 @@ subroutine auxval
   CALL init_CLA(local_nky,local_np*local_nj)
 #endif
   !! Display parallel settings
-  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-  DO i_ = 0,num_procs-1
+  IF(VERBOSE_LVL .GE. 1) THEN
     CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-    IF (my_id .EQ. i_) THEN
-      IF (my_id .EQ. 0) WRITE(*,*) ''
-      IF (my_id .EQ. 0) WRITE(*,*) '--------- Parallel environment ----------'
-      IF (my_id .EQ. 0) WRITE(*,'(A12,I3)') 'n_procs ', num_procs
-      IF (my_id .EQ. 0) WRITE(*,'(A12,I3,A14,I3,A14,I3)') 'num_procs_p   = ', num_procs_p, ', num_procs_ky   = ', num_procs_ky, ', num_procs_z   = ', num_procs_z
-      IF (my_id .EQ. 0) WRITE(*,*) ''
-      WRITE(*,'(A9,I3,A10,I3,A10,I3,A9,I3)')&
-       'my_id  = ', my_id, ', rank_p  = ', rank_p, ', rank_ky  = ', rank_ky,', rank_z  = ', rank_z
-       WRITE(*,'(A22,I3,A11,I3)')&
-       '         local_np   = ', local_np  , ', offset = ', local_np_offset
-       WRITE(*,'(A22,I3,A11,I3)')&
-       '         local_nj   = ', local_nj  , ', offset = ', local_nj_offset
-       WRITE(*,'(A22,I3,A11,I3)')&
-       '         local_nkx  = ', local_nkx , ', offset = ', local_nkx_offset
-       WRITE(*,'(A22,I3,A11,I3)')&
-       '         local_nky  = ', local_nky , ', offset = ', local_nky_offset
-       WRITE(*,'(A22,I3,A11,I3)')&
-       '         local_nz   = ', local_nz  , ', offset = ', local_nz_offset
-      IF (my_id .NE. num_procs-1) WRITE (*,*) ''
-      IF (my_id .EQ. num_procs-1) WRITE(*,*) '------------------------------------------'
-      IF (my_id .EQ. num_procs-1) CALL FLUSH(OUTPUT_UNIT)
-    ENDIF
-  ENDDO
-  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+    DO i_ = 0,num_procs-1
+      CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+      IF (my_id .EQ. i_) THEN
+        IF (my_id .EQ. 0) WRITE(*,*) ''
+        IF (my_id .EQ. 0) WRITE(*,*) '--------- Parallel environment ----------'
+        IF (my_id .EQ. 0) WRITE(*,'(A12,I3)') 'n_procs ', num_procs
+        IF (my_id .EQ. 0) WRITE(*,'(A12,I3,A14,I3,A14,I3)') 'num_procs_p   = ', num_procs_p, ', num_procs_ky   = ', num_procs_ky, ', num_procs_z   = ', num_procs_z
+        IF (my_id .EQ. 0) WRITE(*,*) ''
+        WRITE(*,'(A9,I3,A10,I3,A10,I3,A9,I3)')&
+        'my_id  = ', my_id, ', rank_p  = ', rank_p, ', rank_ky  = ', rank_ky,', rank_z  = ', rank_z
+        WRITE(*,'(A22,I3,A11,I3)')&
+        '         local_np   = ', local_np  , ', offset = ', local_np_offset
+        WRITE(*,'(A22,I3,A11,I3)')&
+        '         local_nj   = ', local_nj  , ', offset = ', local_nj_offset
+        WRITE(*,'(A22,I3,A11,I3)')&
+        '         local_nkx  = ', local_nkx , ', offset = ', local_nkx_offset
+        WRITE(*,'(A22,I3,A11,I3)')&
+        '         local_nky  = ', local_nky , ', offset = ', local_nky_offset
+        WRITE(*,'(A22,I3,A11,I3)')&
+        '         local_nz   = ', local_nz  , ', offset = ', local_nz_offset
+        IF (my_id .NE. num_procs-1) WRITE (*,*) ''
+        IF (my_id .EQ. num_procs-1) WRITE(*,*) '------------------------------------------'
+        IF (my_id .EQ. num_procs-1) CALL FLUSH(OUTPUT_UNIT)
+      ENDIF
+    ENDDO
+    CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+  ENDIF
   SELECT CASE(hierarchy_closure)
   CASE('truncation')
-    CALL speak('Truncation closure')
+    CALL speak('Truncation closure',2)
   CASE('max_degree')
-    CALL speak('Max degree closure -> Maximal Napj degree is D = '// str(dmax))
+    CALL speak('Max degree closure -> Maximal Napj degree is D = '// str(dmax),2)
   END SELECT
 
 END SUBROUTINE auxval
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index a6cf8a8b33c70e5e68a62031def3a1093680b5bd..3909738950c853c5ac660b696d4e5f132342ca62 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -10,6 +10,7 @@ MODULE basic
   real(xp), PUBLIC, PROTECTED :: dt         = 1.0      ! Time step
   real(xp), PUBLIC, PROTECTED :: maxruntime = 1e9      ! Maximum simulation CPU time
   INTEGER,  PUBLIC, PROTECTED :: job2load   = 99       ! jobnum of the checkpoint to load
+  INTEGER,  PUBLIC, PROTECTED :: VERBOSE_LVL= 1        ! tune the amount of std out (0: only time and fluxes, 1: MPI distribution, 2: all stages of init)
   ! Auxiliary variables
   real(xp), PUBLIC, PROTECTED :: time   = 0            ! Current simulation time (Init from restart file)
 
@@ -23,9 +24,11 @@ MODULE basic
   INTEGER, PUBLIC :: iframe2d ! counting the number of times 2d datasets are outputed (for diagnose)
   INTEGER, PUBLIC :: iframe3d ! counting the number of times 3d datasets are outputed (for diagnose)
   INTEGER, PUBLIC :: iframe5d ! counting the number of times 5d datasets are outputed (for diagnose)
+  ! input file name
+  CHARACTER(len=32), PUBLIC, PROTECTED :: input_file = "params.in" 
   !  List of logical file units
-  INTEGER, PUBLIC, PROTECTED  :: lu_in   = 90              ! File duplicated from STDIN
-  INTEGER, PUBLIC, PROTECTED  :: lu_stop = 91              ! stop file, see subroutine TESEND
+  INTEGER, PUBLIC, PROTECTED :: lu_in   = 90 ! File duplicated from STDIN
+  INTEGER, PUBLIC, PROTECTED :: lu_stop = 91 ! stop file, see subroutine TESEND
   ! To measure computation time
   type :: chrono
     real(xp) :: tstart !start of the chrono
@@ -51,9 +54,9 @@ MODULE basic
   CHARACTER(len=*), PUBLIC, PARAMETER :: maindir = ""
 #endif
   ! Routines interfaces
-  PUBLIC :: allocate_array, basic_outputinputs,basic_data,&
+  PUBLIC :: allocate_array, basic_outputinputs,basic_data, show_title,&
             speak, str, increase_step, increase_cstep, increase_time, display_h_min_s,&
-            set_basic_cp, daytim, start_chrono, stop_chrono, change_dt
+            set_basic_cp, daytim, start_chrono, stop_chrono, change_dt, day_and_time_str
   ! Interface for allocating arrays, these routines allocate and initialize directly to zero
   INTERFACE allocate_array
     MODULE PROCEDURE allocate_array_xp1,allocate_array_xp2,allocate_array_xp3, &
@@ -76,7 +79,7 @@ CONTAINS
     use prec_const
     IMPLICIT NONE
 
-    NAMELIST /BASIC/  nrun, dt, tmax, maxruntime, job2load
+    NAMELIST /BASIC/  nrun, dt, tmax, maxruntime, job2load, VERBOSE_LVL
 
     CALL find_input_file
 
@@ -162,35 +165,78 @@ CONTAINS
     timer%ttot = timer%ttot + (timer%tstop-timer%tstart)
   END SUBROUTINE
   !================================================================================
-  ! routine to speak in the terminal
-  SUBROUTINE speak(message)
+  ! routine to speak in the terminal according to the verbose level
+  ! The larger the verbose level, the more std output
+  SUBROUTINE speak(message,priority)
     USE parallel, ONLY: my_id
     IMPLICIT NONE
-    CHARACTER(len=*), INTENT(in) :: message
-    IF(my_id .EQ. 0) write(*,*) message
+    ! input variables
+    CHARACTER(len=*),  INTENT(IN) :: message
+    INTEGER,           INTENT(IN) :: priority
+    ! local variables
+    CHARACTER(len=16) :: indent
+    INTEGER :: i_
+    indent = ''
+    DO i_ = 1,priority
+      indent = trim(indent)//'-'
+    ENDDO
+    IF (priority .LE. VERBOSE_LVL) THEN
+      IF(my_id .EQ. 0) write(*,*) trim(indent)//message
+    ENDIF
   END SUBROUTINE
   !================================================================================
   SUBROUTINE find_input_file
     USE parallel, ONLY: my_id
     IMPLICIT NONE
-    CHARACTER(len=32) :: str_, input_file
+    CHARACTER(len=32) :: str_
     INTEGER :: nargs, fileid, l, ierr
-    LOGICAL :: mlexist
+    LOGICAL :: mlexist, file_found
     nargs = COMMAND_ARGUMENT_COUNT()
-    IF((nargs .EQ. 1) .OR. (nargs .EQ. 4)) THEN
+    file_found = .false. ! to check if the input file was correctly set
+
+    ! Input param, option 1: input file is defined with an index as command argument
+    IF((nargs .EQ. 1) .OR. (nargs .EQ. 4)) THEN 
       CALL GET_COMMAND_ARGUMENT(nargs, str_, l, ierr)
       READ(str_(1:l),'(i3)')  fileid
+      ! Old verson with "fort_XX.90" format (for retro compatibility)
       WRITE(input_file,'(a,a1,i2.2,a3)') 'fort','_',fileid,'.90'
-
       INQUIRE(file=input_file, exist=mlexist)
       IF( mlexist ) THEN
-        IF(my_id.EQ.0) write(*,*) 'Reading input ', input_file,'...'
+        CALL speak('Reading input '// input_file,1)
         OPEN(lu_in, file=input_file)
+        file_found = .true.
       ELSE
-        IF(my_id.EQ.0) write(*,*) 'Reading input fort.90...'
+        ! new verson with "params_XX.in" format
+        WRITE(input_file,'(a,a1,i2.2,a3)') 'params','_',fileid,'.in'
+        INQUIRE(file=input_file, exist=mlexist)
+        IF( mlexist ) THEN
+          CALL speak('Reading input '// input_file,1)
+          OPEN(lu_in, file=input_file)
+          file_found = .true.
+        ENDIF
+      ENDIF
+      IF (.NOT. file_found) THEN
+        WRITE(*,*)  'Error stop: ' // input_file // ' not found.'
+        ERROR STOP
+      ENDIF
+    ELSE
+      ! Input param, option 2: a fort.90 or parameter.in input file is present
+      INQUIRE(file='fort.90', exist=mlexist)
+      IF ( mlexist ) THEN
+        CALL speak('Reading default input file (fort.90)',1)
         OPEN(lu_in, file='fort.90')
+        file_found = .true.
       ENDIF
+      ! Test a parameters.in input file
+      INQUIRE(file='params.in', exist=mlexist)
+      IF ( mlexist ) THEN
+        CALL speak('Reading default input file (params.in)',1)
+        OPEN(lu_in, file='params.in')
+        file_found = .true.
+      ENDIF        
     ENDIF
+    IF(.NOT. file_found) ERROR STOP "Error stop: basic input file not found (a file named params.in or fort.90 must be present)"
+
   END SUBROUTINE find_input_file
   !================================================================================
   SUBROUTINE daytim(str)
@@ -210,6 +256,21 @@ CONTAINS
       WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), time(1:12)
     !
   END SUBROUTINE daytim
+    !================================================================================
+  CHARACTER(64) FUNCTION day_and_time_str
+    !   Print date and time
+    USE parallel, ONLY: my_id
+    use prec_const
+    IMPLICIT NONE
+    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)
+    day_and_time_str = trim(dat(1:10)//' '//trim(time(1:12)))
+    RETURN
+  END
   !================================================================================
   SUBROUTINE display_h_min_s(time)
     USE parallel, ONLY: my_id
@@ -443,4 +504,23 @@ CONTAINS
     a=.false.
   END SUBROUTINE allocate_array_l4
 
+  SUBROUTINE show_title
+  USE parallel, ONLY: my_id
+  IMPLICIT NONE
+  INCLUDE 'srcinfo.h' ! for the code version
+  IF(my_id .EQ. 0) THEN
+    write(*,*) "=============================================="
+    write(*,*) "   ______                                     "
+    write(*,*) "  / ____/_  ______ __________  ____ ___  ____ "
+    write(*,*) " / / __/ / / / __ `/ ___/ __ \/ __ `__ \/ __ \"
+    write(*,*) "/ /_/ / /_/ / /_/ / /__/ /_/ / / / / / / /_/ /"
+    write(*,*) "\____/\__, /\__,_/\___/\____/_/ /_/ /_/\____/ "
+    write(*,*) "     /____/                                   "
+    ! Write the git version of the code
+    write(*,*) 'This is GYACOMO code version:'
+    write(*,*)  VERSION
+    write(*,*) "=============================================="
+  ENDIF
+  END SUBROUTINE show_title
+
 END MODULE basic
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 820c393b9b109e93027dcd023463fac9a7172b96..7816ba24775a4d493b87202d4574531cc1910b2d 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -84,7 +84,7 @@ SUBROUTINE set_closure_model
   SELECT CASE(nonlinear_closure)
   CASE('truncation','monomial')
     IF(nmax .LT. 0) THEN
-      CALL speak("Set nonlinear truncation to anti Laguerre aliasing")
+      CALL speak("Set nonlinear truncation to anti Laguerre aliasing",2)
       DO ij = 1,local_nj
         nmaxarray(ij) = jmax - jarray(ij+ngj_o2)
       ENDDO
diff --git a/src/control.F90 b/src/control.F90
index 0095c2eb2aa84a045234a6f68008d360c0c6f0df..4614cba723e5ff6e7f56caf33cde3d6f39a8817a 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -4,7 +4,7 @@ SUBROUTINE control
   use basic,          ONLY: str,daytim,speak,basic_data,&
                         nlend,step,increase_step,increase_time,increase_cstep,&
                         chrono_runt,chrono_step, chrono_diag, chrono_ExBs,&
-                        start_chrono, stop_chrono
+                        start_chrono, stop_chrono, day_and_time_str, show_title
   use prec_const,     ONLY: xp, stdout
   USE parallel,       ONLY: ppinit
   USE initial,        ONLY: initialize
@@ -17,52 +17,54 @@ SUBROUTINE control
   ! start the chronometer for the total runtime
   CALL start_chrono(chrono_runt)
   !________________________________________________________________________________
-  !              1.   Prologue
+  !           ]   1.   Prologue
+  ! Print out title
+  CALL show_title
+
+  ! Write the current time
+  CALL speak('Start time: '//day_and_time_str(),0)
+
   !                   1.1     Initialize the parallel environment
+  CALL speak('Init MPI [',2)
   CALL ppinit
-  CALL speak('MPI initialized')
+  CALL speak('] MPI initialized',2)
   CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-
-  CALL daytim('Start at ')
   
   !                   1.2     Define data specific to run
-  CALL speak( 'Load basic data...')
+  CALL speak('Load basic data [',2)
   CALL basic_data
-  ! CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-  CALL speak('...basic data loaded.')
+  CALL speak('] basic data loaded.',2)
   
   !                   1.3   Read input parameters from input file
-  CALL speak('Read input parameters...')
+  CALL speak('Read input parameters [',2)
   CALL readinputs
-  ! CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-  CALL speak('...input parameters read')
+  CALL speak('] input parameters read',2)
 
   !                   1.4     Set auxiliary values (allocate arrays, set grid, ...)
-  CALL speak('Calculate auxval...')
+  CALL speak('Setting auxiliary values [',2)
   CALL auxval
-  ! CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-  CALL speak('...auxval calculated')
+  CALL speak('] auxval set ',2)
   
   !                   1.5     Initial conditions
-  CALL speak( 'Create initial state...')
+  CALL speak( 'Create initial state [',2)
   CALL initialize
   ! CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-  CALL speak('...initial state created')
+  CALL speak('] initial state created',2)
   
   !                   1.6     Initial diagnostics
-  CALL speak( 'Initial diagnostics...')
+  CALL speak( 'Initial diagnostics [',2)
   CALL cpu_time(t_init_diag_0) ! Measure the time of the init diag
   CALL diagnose(0)
   CALL cpu_time(t_init_diag_1)
   ! CALL mpi_barrier(MPI_COMM_WORLD, ierr)
-  CALL speak('...initial diagnostics done')
-  CALL speak('('//str(t_init_diag_1-t_init_diag_0)//'[s])')
+  CALL speak('] initial diagnostics done',2)
+  CALL speak('('//str(t_init_diag_1-t_init_diag_0)//'[s])',2)
 
   CALL FLUSH(stdout)
   CALL mpi_barrier(MPI_COMM_WORLD, ierr)
   !________________________________________________________________________________
 
-  CALL speak( 'Time integration loop..')
+  CALL speak( 'Time integration loop [',2)
   !________________________________________________________________________________
   !              2.   Main loop
   DO
@@ -101,7 +103,7 @@ SUBROUTINE control
 
   END DO
 
-  CALL speak('...time integration done')
+  CALL speak('] time integration done',2)
   !________________________________________________________________________________
   !              9.   Epilogue
   ! Stop total run chronometer (must be done before the last diagnostic)
@@ -111,7 +113,8 @@ SUBROUTINE control
   ! end the run
   CALL endrun
   ! display final time
-  CALL daytim('Done at ')
+  ! CALL daytim('Run terminated at ')
+  CALL speak('Finish time: '//day_and_time_str(),0)
   ! close mpi environement
   CALL ppexit
 
diff --git a/src/cosolver_interface_mod.F90 b/src/cosolver_interface_mod.F90
index f59906757f58aab7147bce1256273a18d6e559b3..01a53dd675a59f3f17b5c2561d46b9c7ee6c5f39 100644
--- a/src/cosolver_interface_mod.F90
+++ b/src/cosolver_interface_mod.F90
@@ -286,7 +286,7 @@ CONTAINS
       DEALLOCATE (Caa__kp); DEALLOCATE (CabT_kp); DEALLOCATE (CabF_kp)
 
       IF( .NOT. INTERSPECIES ) THEN
-        CALL speak("--Like Species operator--")
+        CALL speak("--Like Species operator--",2)
         Cab_F = 0._xp;
         Cab_T = 0._xp;
       ENDIF
@@ -301,7 +301,7 @@ CONTAINS
           ENDIF
         ENDDO
       ENDDO
-      CALL speak('============DONE===========')
+      CALL speak('============DONE===========',2)
 
     END SUBROUTINE load_COSOlver_mat
     !******************************************************************************!
diff --git a/src/diagnostics_mod.F90 b/src/diagnostics_mod.F90
index 26e07218e7ec39965296677e39b532a3ba156683..ec0ed7cd66f4d7db03165458f95524bdf68df66b 100644
--- a/src/diagnostics_mod.F90
+++ b/src/diagnostics_mod.F90
@@ -90,7 +90,7 @@ CONTAINS
     ELSE
       CALL creatf(file, fid, mpicomm=comm)
     END IF
-    CALL speak(TRIM(file)//' created')
+    CALL speak(TRIM(file)//' created',2)
     !  basic data group
     CALL creatg(fid, "/data", "data")
     !  File group
diff --git a/src/endrun.F90 b/src/endrun.F90
index 80413fe64732e7554f1bd9e2e0453733972fd203..935a28af792ddbad136ba577e0326daa4a1543a4 100644
--- a/src/endrun.F90
+++ b/src/endrun.F90
@@ -9,12 +9,12 @@ SUBROUTINE endrun
   IF( nlend ) THEN
      !----------------------------------------------------------------------
      !              1.   Normal end of run
-     CALL speak('   Normal exit')
+     CALL speak('Normal exit',0)
 
      !----------------------------------------------------------------------
      !              2.   Abnormal exit
   ELSE
-     WRITE(*,'(/a)') '   Abnormal exit'
+     CALL speak('Abnormal exit',0)
   END IF
 
 
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 8a40f441b2387a9f5c4c8db395c3e758da668038..54e15d62e78cee7b18358cc9db710e3881ce22b4 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -598,7 +598,7 @@ END SUBROUTINE fft1D_plans
     SUBROUTINE finalize_plans
         USE basic, ONLY: speak
         IMPLICIT NONE
-        CALL speak('..plan Destruction.')
+        CALL speak('FFTW3 plan destruction]',2)
         call fftw_destroy_plan(planb)
         call fftw_destroy_plan(planf)
         call fftw_mpi_cleanup()
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 9bf73ed26dcff2e65950924ef20cd24c91351c47..4b239fd756e55adccff28beb7f7ace6f9df92fef 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -101,7 +101,7 @@ CONTAINS
       CASE DEFAULT
         ERROR STOP '>> ERROR << Parallel BC not recognized'
     END SELECT
-    CALL speak('Parallel BC : '//parallel_bc)
+    CALL speak('Parallel BC : '//parallel_bc,2)
   END SUBROUTINE geometry_readinputs
 
   subroutine eval_magnetic_geometry
@@ -125,19 +125,19 @@ CONTAINS
 
     !
     IF( (total_nky .EQ. 1) .AND. (total_nz .EQ. 1)) THEN !1D perp linear run
-      CALL speak('1D perpendicular geometry')
+      CALL speak('1D perpendicular geometry',2)
       call eval_1D_geometry
     ELSE
       SELECT CASE(geom)
         CASE('s-alpha','salpha')
-          CALL speak('s-alpha geometry')
+          CALL speak('s-alpha geometry',2)
           IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for s-alpha geometry"
           IF(MODULO(FLOOR(Npol),2) .EQ. 0)   ERROR STOP "Npol must be odd for s-alpha"
           call eval_salpha_geometry
           C_y     = 1._xp
           Cyq0_x0 = 1._xp
         CASE('z-pinch','Z-pinch','zpinch','Zpinch')
-          CALL speak('Z-pinch geometry')
+          CALL speak('Z-pinch geometry',2)
           call eval_zpinch_geometry
           SHEARED = .FALSE.
           shear   = 0._xp
@@ -147,7 +147,7 @@ CONTAINS
           C_y     = 1._xp
           Cyq0_x0 = 1._xp
         CASE('miller','Miller','Miller_global','miller_global')
-          CALL speak('Miller geometry')
+          CALL speak('Miller geometry',2)
           IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for Miller geometry"
           IF(MODULO(FLOOR(Npol),2) .EQ. 0)   THEN
             write(*,*) "Npol must be odd for Miller, (Npol = ",Npol,")"
@@ -158,7 +158,7 @@ CONTAINS
                           C_y,C_xy,Cyq0_x0,xpdx_pm_geom,gxx,gxy,gxz,gyy,gyz,gzz,&
                           dBdx,dBdy,dBdz,hatB,jacobian,hatR,hatZ,dxdR,dxdZ)
         CASE('circular','circ')
-          CALL speak('Circular geometry')
+          CALL speak('Circular geometry',2)
           IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for circular geometry"
           IF(MODULO(FLOOR(Npol),2) .EQ. 0)   THEN
             write(*,*) "Npol must be odd for circular, (Npol = ",Npol,")"
@@ -174,9 +174,8 @@ CONTAINS
     ! inv squared of the magnetic gradient bckg amplitude (for fast kperp update)
     inv_hatB2 = 1/hatB/hatB
     ! Reset kx grid (to account for Cyq0_x0 factor)
-    CALL speak('--Reset kx grid according to Cyq0_x0 factor--')
+    CALL speak('Reset kx grid according to Cyq0_x0 factor',1)
     CALL set_kxgrid(shear,ExBrate,Npol,Cyq0_x0,LINEARITY,N_HD) 
-    CALL speak('-- done --')
     !
     ! Evaluate perpendicular wavenumber
     !  k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 7b7fc5b9963f8dde07d70139855e1674a841f290..4b0b2243b97aafc4d2e1361fc4ed262cc794a7f3 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -199,10 +199,10 @@ CONTAINS
     !!----------------- SPATIAL GRIDS (parallelized)
     !! Parallel distribution of kx ky grid
     IF (LINEARITY .NE. 'linear') THEN ! we let FFTW distribute if we use it
-      IF (my_id .EQ. 0) write(*,*) 'FFTW3 y-grid distribution'
+      CALL speak('FFTW3 y-grid distribution',2)
       CALL init_grid_distr_and_plans(Nx,Ny,comm_ky,local_nx_ptr,local_nx_ptr_offset,local_nky_ptr,local_nky_ptr_offset)
     ELSE ! otherwise we distribute equally
-      IF (my_id .EQ. 0) write(*,*) 'Manual y-grid distribution'
+      CALL speak('Manual y-grid distribution',2)
       ! balanced distribution among the processes
       CALL  decomp1D( Ny/2+1, num_procs_ky, rank_ky, ikys, ikye )
       local_nky_ptr        = ikye - ikys + 1
@@ -239,7 +239,7 @@ CONTAINS
     !!---------------- PARALLEL Z GRID (parallelized)
     total_nz = Nz
     IF (SG) THEN
-      CALL speak('--2 staggered z grids--')
+      CALL speak('--2 staggered z grids--',2)
       ! indices for even p and odd p grids (used in kernel, jacobian, gij etc.)
       ieven  = 1
       iodd   = 2
@@ -499,17 +499,17 @@ CONTAINS
     INTEGER :: ikx, iky
     REAL(xp):: Lx_adapted
     IF(shear .GT. 0) THEN
-      CALL speak('Magnetic shear detected: set up sheared kx grid..')
+      CALL speak('Magnetic shear detected: set up sheared kx grid..',1)
       ! mininal size of box in x to respect dkx = 2pi shear dky
       Lx_adapted = Ly/(2._xp*pi*shear*Npol*Cyq0_x0)
       ! Put Nexc to 0 so that it is computed from a target value Lx
       IF(Nexc .EQ. 0) THEN
         Nexc = CEILING(0.9 * Lx/Lx_adapted)
-        CALL speak('Adapted Nexc ='//str(Nexc))
+        CALL speak('Adapted Nexc ='//str(Nexc),1)
       ENDIF
       ! x length is adapted
       Lx = Lx_adapted*Nexc
-      CALL speak('Simulation Lx ='//str(Lx))
+      CALL speak('Simulation Lx ='//str(Lx),1)
     ENDIF
     deltakx = 2._xp*PI/Lx
     inv_dkx = 1._xp/deltakx 
diff --git a/src/initial_mod.F90 b/src/initial_mod.F90
index c5668dfac54b7bf69ff3030b76c950d5982fecfa..355d68b0e0f3c08bf808448a196a5739a5fc2894 100644
--- a/src/initial_mod.F90
+++ b/src/initial_mod.F90
@@ -80,7 +80,7 @@ CONTAINS
     !!!!!! Set the moments arrays Nepj, Nipj and phi!!!!!!
     ! through loading a previous state
     IF ( job2load .GE. 0 ) THEN
-      CALL speak('Load moments')
+      CALL speak('Load moments',2)
       CALL load_moments ! get N_0
       CALL apply_closure_model
       CALL update_ghosts_moments
@@ -90,44 +90,44 @@ CONTAINS
       SELECT CASE (INIT_OPT)
       ! set phi with noise and set moments to 0
       CASE ('phi')
-        CALL speak('Init noisy phi')
+        CALL speak('Init noisy phi',2)
         CALL init_phi
       CASE ('phi_ppj')
-        CALL speak('Init noisy phi')
+        CALL speak('Init noisy phi',2)
         CALL init_phi_ppj
       ! set moments_00 (GC density) with noise and compute phi afterwards
       CASE('mom00')
-        CALL speak('Init noisy gyrocenter density')
+        CALL speak('Init noisy gyrocenter density',2)
         CALL init_gyrodens ! init only gyrocenter density
         CALL update_ghosts_moments
         CALL solve_EM_fields
       ! set moments_00 (GC density) with a single kx0,ky0 mode
       ! kx0 is setup but by the noise value and ky0 by the background value
       CASE('mom00_modes','mom00_mode')
-        CALL speak('Init modes in the gyrocenter density')
+        CALL speak('Init modes in the gyrocenter density',2)
         CALL init_modes ! init single mode in gyrocenter density
         CALL update_ghosts_moments
         CALL solve_EM_fields
       ! init all moments randomly (unadvised)
       CASE('allmom')
-        CALL speak('Init noisy moments')
+        CALL speak('Init noisy moments',2)
         CALL init_moments ! init all moments
         CALL update_ghosts_moments
         CALL solve_EM_fields
       ! init a gaussian blob in gyrodens
       CASE('blob','blob_nokx')
-        CALL speak('--init a blob')
+        CALL speak('Blob initialization',2)
         CALL initialize_blob
         CALL update_ghosts_moments
         CALL solve_EM_fields
       ! init moments 00 with a power law similarly to GENE
       CASE('ppj')
-        CALL speak('ppj init ~ GENE')
+        CALL speak('ppj init ~ GENE',2)
         call init_ppj
         CALL update_ghosts_moments
         CALL solve_EM_fields
       CASE('ricci')
-        CALL speak('Init Ricci')
+        CALL speak('Init Ricci',2)
         CALL init_ricci ! init only gyrocenter density
         CALL update_ghosts_moments
         CALL solve_EM_fields
@@ -136,10 +136,10 @@ CONTAINS
       END SELECT
     ENDIF
     ! closure of j>J, p>P and j<0, p<0 moments
-    CALL speak('Apply closure')
+    CALL speak('Apply closure',2)
     CALL apply_closure_model
     ! ghosts for p parallelization
-    CALL speak('Ghosts communication')
+    CALL speak('Ghosts communication',2)
     CALL update_ghosts_moments
     CALL update_ghosts_EM
     !! End of phi and moments initialization
@@ -147,7 +147,7 @@ CONTAINS
     CALL init_collision
     !! Preparing auxiliary arrays at initial state
     ! particle density, fluid velocity and temperature (used in diagnose)
-    CALL speak('Computing fluid moments and transport')
+    CALL speak('Computing fluid moments and transport',2)
     CALL compute_fluid_moments
     CALL compute_radial_transport
     CALL compute_radial_heatflux
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index 5df1d153fa1c8f9a2b8c79bea474f06f5950c6b7..7855e8cc29da4fecd9b8ce12d96ac5fb849c5dac 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -330,18 +330,15 @@ CONTAINS
 
     !definition of radial coordinate! dx_drho=1 --> x = r
     dx_drho=1. !drPsi/Psi0*Lnorm*q0
-    if (my_id==0) write(*,"(A,ES12.4)") 'Using radial coordinate with dx/dr = ',dx_drho
+    CALL speak('Using radial coordinate with dx/dr = '//str(dx_drho),2)
     dxPsi  = drPsi/dx_drho
     C_y    = dxPsi*sign_Ip_CW
     C_xy   = abs(B0*dxPsi/C_y)
 
-    if (my_id==0) then
-       write(*,"(A,ES12.4,A,ES12.4,A,ES12.4)") &
-            "Setting C_xy = ",C_xy,' C_y = ', C_y," C_x' = ", 1./dxPsi
-       write(*,'(A,ES12.4)') "B_unit/Bref conversion factor = ", q0/rho*drPsi
-       write(*,'(A,ES12.4)') "dPsi/dr = ", drPsi
-       if (thetaShift.ne.0._xp) write(*,'(A,ES12.4)') "thetaShift = ", thetaShift
-    endif
+    CALL speak("Setting C_xy = "//str(C_xy)//' C_y = '//str(C_y)//" C_x' = "//str(1./dxPsi),2)
+    CALL speak("B_unit/Bref conversion factor = "//str(q0/rho*drPsi),2)
+    CALL speak("dPsi/dr = "//str(drPsi),2)
+    IF (thetaShift.ne.0._xp) CALL speak("thetaShift = "//str(thetaShift),2)
 
 
     !--------shear is expected to be defined as rho/q*dq/drho--------!
@@ -382,9 +379,8 @@ CONTAINS
 
     ffprime=-(sign_Ip_CW*dq_dpsi*2._xp*pi*Npol_ext+D0_full+D2_full*pprime)/D1_full
 
-    if (my_id==0) then
-       write(*,'(A,ES12.4)') "ffprime = ", ffprime
-    endif
+    CALL speak("ffprime = "//str(ffprime),2)
+
     D0=D0-D0_mid
     D1=D1-D1_mid
     D2=D2-D2_mid
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 529dd217cb4ca08896eb6cbfc4a147be106654d2..7e84045ba8578baf3404d280393cce6e2d0d4101 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -77,15 +77,15 @@ CONTAINS
     IF(Na .EQ. 1) THEN
       IF((.NOT. ADIAB_E) .AND. (.NOT. ADIAB_I)) ERROR STOP "With one species, ADIAB_E or ADIAB_I must be set to .true. STOP"
       IF(ADIAB_E) THEN
-        CALL speak('Adiabatic electron model -> beta = 0')
+        CALL speak('Adiabatic electron model -> beta = 0',1)
         beta = 0._xp
         EM   = .false.
       ENDIF
-      IF(ADIAB_I) CALL speak('Adiabatic ion model')
+      IF(ADIAB_I) CALL speak('Adiabatic ion model',1)
     ENDIF
 
     IF(beta .GT. 0) THEN
-      CALL speak('Electromagnetic effects are included')
+      CALL speak('Electromagnetic effects are included',1)
       EM = .TRUE.
     ELSE
       EM = .FALSE.
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 579ca824b6a9259327d46a5f1d4c3d0aa7e865de..4ab23a6396d28dcca8245e1cb3f7006d4c66e23c 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -43,7 +43,7 @@ CONTAINS
     CALL cpu_time(timer_tot_1)
     ! Checkpoint filename
     WRITE(rstfile,'(a,i2.2,a3)') 'outputs_',job2load,'.h5'
-    CALL speak("Resume from "//rstfile)
+    CALL speak("Resume from "//rstfile,0)
     ! Open file
 #ifdef SINGLE_PRECISION
     CALL openf(rstfile, fidrst, mode='r', mpicomm=comm0)
@@ -58,21 +58,21 @@ CONTAINS
     CALL getatt(fidrst,"/data/input/grid" , "Nkx", Nkx_cp)
     CALL getatt(fidrst,"/data/input/grid" ,  "Nz", Nz_cp)
     !! Check dimensional compatibility with the planned simulation
-    IF(Na_cp  .NE. total_Na ) CALL speak("NOTE: Na  has changed")
-    IF(Np_cp  .NE. total_Np ) CALL speak("NOTE: Np  has changed")
-    IF(Nj_cp  .NE. total_Nj ) CALL speak("NOTE: Nj  has changed")
-    IF(Nky_cp .NE. total_Nky) CALL speak("NOTE: Nky has changed")
-    IF(Nkx_cp .NE. total_Nkx) CALL speak("NOTE: Nkx has changed")
-    IF(Nz_cp  .NE. total_Nz)  CALL speak("NOTE: Nz  has changed")
+    IF(Na_cp  .NE. total_Na ) CALL speak("NOTE: Na  has changed",1)
+    IF(Np_cp  .NE. total_Np ) CALL speak("NOTE: Np  has changed",1)
+    IF(Nj_cp  .NE. total_Nj ) CALL speak("NOTE: Nj  has changed",1)
+    IF(Nky_cp .NE. total_Nky) CALL speak("NOTE: Nky has changed",1)
+    IF(Nkx_cp .NE. total_Nkx) CALL speak("NOTE: Nkx has changed",1)
+    IF(Nz_cp  .NE. total_Nz)  CALL speak("NOTE: Nz  has changed",1)
     ! Get the x,y fourier modes and the z space and check grids
     ALLOCATE(ky_cp(Nky_cp),kx_cp(Nkx_cp),z_cp(Nz_cp))
     CALL getarr(fidrst,"/data/grid/coordky",ky_cp); Ly_cp = 2._xp*PI/ky_cp(2)
     CALL getarr(fidrst,"/data/grid/coordkx",kx_cp); Lx_cp = 2._xp*PI/kx_cp(2)
     CALL getarr(fidrst,"/data/grid/coordz" , z_cp); Lz_cp =-2._xp*z_cp(1)
     ! check changes
-    IF(ABS(ky_cp(2)-kyarray_full(2)) .GT. 1e-3) CALL speak("NOTE: Ly  has changed")
-    IF(ABS(kx_cp(2)-kxarray_full(2)) .GT. 1e-3) CALL speak("NOTE: Lx  has changed")
-    IF(ABS(z_cp(1) - zarray_full(1)) .GT. 1e-3) CALL speak("NOTE: Lz  has changed")
+    IF(ABS(ky_cp(2)-kyarray_full(2)) .GT. 1e-3) CALL speak("NOTE: Ly  has changed",1)
+    IF(ABS(kx_cp(2)-kxarray_full(2)) .GT. 1e-3) CALL speak("NOTE: Lx  has changed",1)
+    IF(ABS(z_cp(1) - zarray_full(1)) .GT. 1e-3) CALL speak("NOTE: Lz  has changed",1)
     dz_cp = (z_cp(2)- z_cp(1)) !! check deltaz changes
     ! IF(ABS(deltaz - dz_cp) .GT. 1e-3) ERROR STOP "ERROR STOP: change of deltaz is not implemented."
     ! compute the mapping for z grid adaptation
@@ -102,7 +102,7 @@ CONTAINS
     CALL getatt(fidrst, dset_name, 'jobnum', jobnum_cp)
     ! Set the cstep, step, time and iframnd variables in basic from the checkpoint
     CALL set_basic_cp(cstep_cp,time_cp,jobnum_cp)
-    CALL speak('.. restart from t = '//str(time))
+    CALL speak('.. restart from t = '//str(time),1)
     ! Read state of system from checkpoint file
     ! Brute force loading: load the full moments and take what is needed (RAM dangerous...)
     ! (other possibility would be to load all on process 0 and send the slices)
@@ -150,11 +150,11 @@ CONTAINS
     !! deallocate the full moment variable
     DEALLOCATE(moments_cp)
     CALL closef(fidrst)
-    CALL speak("Reading from restart file "//TRIM(rstfile)//" completed!")
+    CALL speak("Reading from restart file "//TRIM(rstfile)//" completed!",1)
     CALL mpi_barrier(MPI_COMM_WORLD, ierr)
     ! stop time measurement
     CALL cpu_time(timer_tot_2)
-    CALL speak('** Total load time : '// str(timer_tot_2 - timer_tot_1)//' **')
+    CALL speak('** Total load time : '// str(timer_tot_2 - timer_tot_1)//' **',1)
   END SUBROUTINE load_moments
   !******************************************************************************!
   ! Auxiliary routine
diff --git a/src/tesend.F90 b/src/tesend.F90
index 0dedde74007e58b393d4e9b1694eb001eea04440..ef1e6bb747e37c26329e12747fe92aae7e354532 100644
--- a/src/tesend.F90
+++ b/src/tesend.F90
@@ -16,8 +16,8 @@ SUBROUTINE tesend
   IF( mlend ) THEN
     nlend   = .TRUE.
     crashed = .TRUE.
-    CALL speak('rhs are NaN/Inf')
-    CALL speak('Run terminated at cstep='//str(cstep))
+    CALL speak('rhs is NaN/Inf',0)
+    CALL speak('Run terminated at cstep='//str(cstep),0)
     RETURN
   END IF
 
@@ -25,7 +25,7 @@ SUBROUTINE tesend
   !                   2.  Test on NRUN
   nlend = step .GE. nrun
   IF ( nlend ) THEN
-     CALL speak('NRUN steps done')
+     CALL speak('NRUN steps done',0)
      RETURN
   END IF
 
@@ -34,7 +34,7 @@ SUBROUTINE tesend
   !                   3.  Test on TMAX
   nlend = time .GE. tmax
   IF ( nlend ) THEN
-     CALL speak('TMAX reached')
+     CALL speak('TMAX reached',0)
      RETURN
   END IF
   !
@@ -48,7 +48,7 @@ SUBROUTINE tesend
 
   CALL mpi_allreduce(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr)
   IF ( nlend ) THEN
-     CALL speak('Max run time reached')
+     CALL speak('Max run time reached',0)
      RETURN
   END IF
   !________________________________________________________________________________
@@ -59,7 +59,7 @@ SUBROUTINE tesend
      IF( mlexist ) THEN
         OPEN(lu_stop, file=stop_file)
         mlend = mlexist ! Send stop status asa the file exists
-        CALL speak('Stop file found -> finishing..')
+        CALL speak('Stop file found -> finishing..',0)
         CLOSE(lu_stop, status='delete')
      END IF
   END IF
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index ce550814a7491af6f5b3eeff59599b58826dd5ee..b76fd8f1119fb07b3ba80daea06c56745b4932ca 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -68,6 +68,7 @@ CONTAINS
   SUBROUTINE set_numerical_scheme
     ! Initialize Butcher coefficient of set_numerical_scheme
     use parallel, ONLY: my_id
+    use basic,    ONLY: speak
     IMPLICIT NONE
     SELECT CASE (numerical_scheme)
     ! Order 1 method
@@ -93,9 +94,9 @@ CONTAINS
       time_scheme_is_adaptive = .true.
       CALL DOPRI5
     CASE DEFAULT
-       IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
+       ERROR STOP 'Cannot initialize time integration scheme. Name invalid.'
     END SELECT
-    IF (my_id .EQ. 0) WRITE(*,*) " Time integration with ", numerical_scheme
+    CALL speak("-Time integration with "// numerical_scheme,2)
   END SUBROUTINE set_numerical_scheme
 
   !!! first order time schemes
@@ -326,8 +327,8 @@ CONTAINS
 
     if(should_discard_step) then
         ! Note this relates to the previous step!
-        if(dt/=old_dt) CALL speak("step discarded. dt update "//str(old_dt)//" => "//str(dt))
-        if(dt/=old_dt) CALL speak("errscale was "//str(adaptive_accuracy_ratio))
+        if(dt/=old_dt) CALL speak("step discarded. dt update "//str(old_dt)//" => "//str(dt),2)
+        if(dt/=old_dt) CALL speak("errscale was "//str(adaptive_accuracy_ratio),2)
     end if
 
     ! Setting `adaptive_accuracy_ratio` is equivalent to saying that the step