diff --git a/src/LIB/SURCOUCHE/src/modd_field.f90 b/src/LIB/SURCOUCHE/src/modd_field.f90
index b81b59f1f0354c2041f9908b0d3c7fb7230e0f3d..cf3e014068caca76b99fbcba63b96680894a39d6 100644
--- a/src/LIB/SURCOUCHE/src/modd_field.f90
+++ b/src/LIB/SURCOUCHE/src/modd_field.f90
@@ -12,6 +12,7 @@
 !  P. Wautelet 27/01/2020: create the tfield_metadata_base abstract datatype
 !  P. Wautelet 14/09/2020: add ndimlist field to tfield_metadata_base
 !  P. Wautelet 10/11/2020: new data structures for netCDF dimensions
+!  P. Wautelet 08/10/2021: add 2 new dimensions: LW_bands (NMNHDIM_NLWB) and SW_bands (NMNHDIM_NSWB)
 !-----------------------------------------------------------------
 module modd_field
 
@@ -42,49 +43,52 @@ integer, parameter :: NMNHDIM_TIME                = 9
 
 integer, parameter :: NMNHDIM_ONE                 = 10
 
-integer, parameter :: NMNHDIM_LASTDIM_NODIACHRO   = 10  ! Index of the last defined dimension for non-diachronic files
+integer, parameter :: NMNHDIM_NSWB                = 11
+integer, parameter :: NMNHDIM_NLWB                = 12
 
-integer, parameter :: NMNHDIM_COMPLEX             = 11
+integer, parameter :: NMNHDIM_LASTDIM_NODIACHRO   = 12 ! Index of the last defined dimension for non-diachronic files
 
-integer, parameter :: NMNHDIM_BUDGET_CART_NI      = 12
-integer, parameter :: NMNHDIM_BUDGET_CART_NJ      = 13
-integer, parameter :: NMNHDIM_BUDGET_CART_NI_U    = 14
-integer, parameter :: NMNHDIM_BUDGET_CART_NJ_U    = 15
-integer, parameter :: NMNHDIM_BUDGET_CART_NI_V    = 16
-integer, parameter :: NMNHDIM_BUDGET_CART_NJ_V    = 17
-integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL   = 18
-integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL_W = 19
+integer, parameter :: NMNHDIM_COMPLEX             = 13
 
-integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL   = 20
-integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL_W = 21
-integer, parameter :: NMNHDIM_BUDGET_MASK_NBUMASK = 22
+integer, parameter :: NMNHDIM_BUDGET_CART_NI      = 14
+integer, parameter :: NMNHDIM_BUDGET_CART_NJ      = 15
+integer, parameter :: NMNHDIM_BUDGET_CART_NI_U    = 16
+integer, parameter :: NMNHDIM_BUDGET_CART_NJ_U    = 17
+integer, parameter :: NMNHDIM_BUDGET_CART_NI_V    = 18
+integer, parameter :: NMNHDIM_BUDGET_CART_NJ_V    = 19
+integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL   = 20
+integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL_W = 21
 
-integer, parameter :: NMNHDIM_BUDGET_TIME         = 23
+integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL   = 22
+integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL_W = 23
+integer, parameter :: NMNHDIM_BUDGET_MASK_NBUMASK = 24
 
-integer, parameter :: NMNHDIM_BUDGET_LES_TIME     = 24
-integer, parameter :: NMNHDIM_BUDGET_LES_AVG_TIME = 25
-integer, parameter :: NMNHDIM_BUDGET_LES_LEVEL    = 26
-integer, parameter :: NMNHDIM_BUDGET_LES_SV       = 27
-integer, parameter :: NMNHDIM_BUDGET_LES_PDF      = 28
+integer, parameter :: NMNHDIM_BUDGET_TIME         = 25
+
+integer, parameter :: NMNHDIM_BUDGET_LES_TIME     = 26
+integer, parameter :: NMNHDIM_BUDGET_LES_AVG_TIME = 27
+integer, parameter :: NMNHDIM_BUDGET_LES_LEVEL    = 28
+integer, parameter :: NMNHDIM_BUDGET_LES_SV       = 29
+integer, parameter :: NMNHDIM_BUDGET_LES_PDF      = 30
 integer, parameter :: NMNHDIM_BUDGET_LES_MASK     = 100 ! This is not a true dimension
 
-integer, parameter :: NMNHDIM_SPECTRA_2PTS_NI     = 29
-integer, parameter :: NMNHDIM_SPECTRA_2PTS_NJ     = 30
-integer, parameter :: NMNHDIM_SPECTRA_SPEC_NI     = 31
-integer, parameter :: NMNHDIM_SPECTRA_SPEC_NJ     = 32
-integer, parameter :: NMNHDIM_SPECTRA_LEVEL       = 33
+integer, parameter :: NMNHDIM_SPECTRA_2PTS_NI     = 31
+integer, parameter :: NMNHDIM_SPECTRA_2PTS_NJ     = 32
+integer, parameter :: NMNHDIM_SPECTRA_SPEC_NI     = 33
+integer, parameter :: NMNHDIM_SPECTRA_SPEC_NJ     = 34
+integer, parameter :: NMNHDIM_SPECTRA_LEVEL       = 35
 
-integer, parameter :: NMNHDIM_SERIES_LEVEL        = 34
-integer, parameter :: NMNHDIM_SERIES_LEVEL_W      = 35
-integer, parameter :: NMNHDIM_SERIES_TIME         = 36  ! Time dimension for time series
+integer, parameter :: NMNHDIM_SERIES_LEVEL        = 36
+integer, parameter :: NMNHDIM_SERIES_LEVEL_W      = 37
+integer, parameter :: NMNHDIM_SERIES_TIME         = 38  ! Time dimension for time series
 
-integer, parameter :: NMNHDIM_FLYER_TIME          = 37  ! Time dimension for aircraft/balloon (dimension local to each flyer)
-integer, parameter :: NMNHDIM_PROFILER_TIME       = 38  ! Time dimension for profilers
-integer, parameter :: NMNHDIM_STATION_TIME        = 39  ! Time dimension for stations
+integer, parameter :: NMNHDIM_FLYER_TIME          = 39  ! Time dimension for aircraft/balloon (dimension local to each flyer)
+integer, parameter :: NMNHDIM_PROFILER_TIME       = 40  ! Time dimension for profilers
+integer, parameter :: NMNHDIM_STATION_TIME        = 41  ! Time dimension for stations
 
-integer, parameter :: NMNHDIM_PAIR                = 40  ! For values coming by pair (ie boundaries)
+integer, parameter :: NMNHDIM_PAIR                = 42  ! For values coming by pair (ie boundaries)
 
-integer, parameter :: NMNHDIM_LASTDIM_DIACHRO     = 40  ! Index of the last defined dimension for diachronic files
+integer, parameter :: NMNHDIM_LASTDIM_DIACHRO     = 42  ! Index of the last defined dimension for diachronic files
 
 integer, parameter :: NMNHDIM_BUDGET_NGROUPS      = 101 ! This is not a true dimension
 integer, parameter :: NMNHDIM_FLYER_PROC          = 102 ! This is not a true dimension
diff --git a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90 b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
index b48a37b2649d4d8ae06964e24a736db43439cbaf..9e95c6e518de2c22d3c373e9d8994cfe72e12a11 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -1079,6 +1079,8 @@ end subroutine IO_Ndimlist_reduce
     iresp     = 0
     iresp_lfi = 0
     iresp_nc4 = 0
+    iresp_tmp_lfi = 0
+    iresp_tmp_nc4 = 0
     GALLOC    = .FALSE.
     GALLOC_ll = .FALSE.
     IHEXTOT = 2*JPHEXT+1
diff --git a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90 b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
index f7954781fd46754679b6f25aeb2c6fc018a19b6b..60cf44de04d57ca5056621b06c3fda899d42bdb5 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2016-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2016-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -19,6 +19,8 @@
 !  P. Wautelet 11/02/2020: bugfix: TDADFILE was wrongly constructed for output files
 !  S. Donnier  28/02/2020: type STREAM needed for use of ECOCLIMAP SG
 !  P. Wautelet 08/01/2021: allow output files with empty variable list (useful if IO_Field_user_write is not empty)
+!  P. Wautelet 18/03/2022: minor bugfix in ISTEP_MAX computation + adapt diagnostics messages
+                           (change verbosity level and remove some unnecessary warnings)
 !-----------------------------------------------------------------
 MODULE MODE_IO_MANAGE_STRUCT
 !
@@ -86,8 +88,8 @@ END IF
 DO IMI = 1, NMODEL
   IBAK_NUMB = 0
   IOUT_NUMB = 0
-  ISTEP_MAX = NINT(XSEGLEN/DYN_MODEL(IMI)%XTSTEP)+1
-  IF (IMI == 1) ISTEP_MAX = ISTEP_MAX - KSUP
+  !Reduce XSEGLEN by time added to XSEGLEN for 1st domain (see set_grid subroutine)
+  ISTEP_MAX = NINT( ( XSEGLEN - KSUP * DYN_MODEL(1)%XTSTEP ) / DYN_MODEL(IMI)%XTSTEP ) + 1
   !
   !* Insert regular backups/outputs into XBAK_TIME/XOUT_TIME arrays
   !
@@ -468,7 +470,7 @@ SUBROUTINE FIND_REMOVE_DUPLICATES(KNUMB,KSTEPS)
   DO JOUT = 1,KNUMB
     DO JKLOOP = JOUT+1,KNUMB
       IF ( KSTEPS(JKLOOP) == KSTEPS(JOUT) .AND. KSTEPS(JKLOOP) > 0 ) THEN
-        CALL PRINT_MSG(NVERB_INFO,'IO','FIND_REMOVE_DUPLICATES','found duplicated backup/output step (removed extra one)')
+        CALL PRINT_MSG(NVERB_DEBUG,'IO','FIND_REMOVE_DUPLICATES','found duplicated backup/output step (removed extra one)')
         KSTEPS(JKLOOP) = NNEGUNDEF
       END IF
     END DO
@@ -484,7 +486,8 @@ SUBROUTINE FIND_REMOVE_OUTOFTIMERANGE(KNUMB,KSTEPS)
   !
   DO JOUT = 1,KNUMB
     IF ( KSTEPS(JOUT) < 1 .OR. KSTEPS(JOUT) > ISTEP_MAX ) THEN
-      CALL PRINT_MSG(NVERB_WARNING,'IO','FIND_REMOVE_OUTOFTIMERANGE','found backup/output step outside of time range')
+      IF ( KSTEPS(JOUT) /= NNEGUNDEF ) &
+        CALL PRINT_MSG(NVERB_WARNING,'IO','FIND_REMOVE_OUTOFTIMERANGE','found backup/output step outside of time range')
       KSTEPS(JOUT) = NNEGUNDEF
     END IF
   END DO
diff --git a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
index 477c389dbf6ffe6691029ff11a627e431b2e0a84..b436b0d18cf7763e6d5cdfc213050d36ec6cfec5 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -19,6 +19,7 @@
 !  P. Wautelet 18/03/2021: workaround for an intel compiler bug
 !  P. Wautelet 04/05/2021: improve IO_Vdims_fill_nc4 if l2d and lpack
 !  P. Wautelet 27/05/2021: improve IO_Mnhname_clean to autocorrect names to be CF compliant
+!  P. Wautelet 08/10/2021: add 2 new dimensions: LW_bands and SW_bands
 !-----------------------------------------------------------------
 #ifdef MNH_IOCDF4
 module mode_io_tools_nc4
@@ -251,10 +252,10 @@ USE MODD_CONF,          ONLY: CPROGRAM, l2d, lpack
 USE MODD_CONF_n,        ONLY: CSTORAGE_TYPE
 USE MODD_DIM_n,         ONLY: NIMAX_ll, NJMAX_ll, NKMAX
 use modd_dyn,           only: xseglen
-use modd_dyn_n,         only: xtstep
+use modd_dyn_n,         only: dyn_model
 use modd_field,         only: NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_NI_U, NMNHDIM_NJ_U, NMNHDIM_NI_V, NMNHDIM_NJ_V,   &
                               NMNHDIM_LEVEL, NMNHDIM_LEVEL_W, NMNHDIM_TIME,                                     &
-                              NMNHDIM_ONE,  NMNHDIM_COMPLEX,                                                    &
+                              NMNHDIM_ONE,  NMNHDIM_NSWB, NMNHDIM_NLWB, NMNHDIM_COMPLEX,                        &
                               NMNHDIM_BUDGET_CART_NI, NMNHDIM_BUDGET_CART_NJ, NMNHDIM_BUDGET_CART_NI_U,         &
                               NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NI_V, NMNHDIM_BUDGET_CART_NJ_V,     &
                               NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W,                           &
@@ -274,7 +275,9 @@ use modd_les,           only: lles_pdf, nles_k, npdf, nspectra_k, xles_temp_mean
 use modd_les_n,         only: nles_times, nspectra_ni, nspectra_nj
 use modd_nsv,           only: nsv
 USE MODD_PARAMETERS_ll, ONLY: JPHEXT, JPVEXT
+use modd_param_n,       only: crad
 use modd_profiler_n,    only: numbprofiler, tprofiler
+use modd_radiations_n,  only: nlwb_mnh, nswb_mnh
 use modd_series,        only: lseries
 use modd_series_n,      only: nsnbstept
 use modd_station_n,     only: numbstat, tstation
@@ -336,6 +339,11 @@ if ( tpfile%ctype == 'MNHDIACHRONIC' .or. ( lpack .and. l2d ) ) then
   call IO_Add_dim_nc4( tpfile, NMNHDIM_ONE, 'one', 1 )
 end if
 
+if ( tpfile%ctype /= 'MNHDIACHRONIC' .and. crad /= 'NONE' ) then
+  call IO_Add_dim_nc4( tpfile, NMNHDIM_NSWB, 'SW_bands', nswb_mnh ) !number of SW bands practically used
+  call IO_Add_dim_nc4( tpfile, NMNHDIM_NLWB, 'LW_bands', nlwb_mnh ) !number of LW bands practically used
+end if
+
 !Write dimensions used in diachronic files
 if ( tpfile%ctype == 'MNHDIACHRONIC' ) then
   !Dimension of size 2 used for NMNHDIM_COMPLEX
@@ -413,13 +421,13 @@ if ( tpfile%ctype == 'MNHDIACHRONIC' ) then
 
   !Dimension for the number of profiler times
   if ( numbprofiler > 0 ) then
-    iprof = Int ( ( xseglen - xtstep ) / tprofiler%step ) + 1
+    iprof = Nint ( ( xseglen - dyn_model(1)%xtstep ) / tprofiler%step ) + 1
     call IO_Add_dim_nc4( tpfile, NMNHDIM_PROFILER_TIME, 'time_profiler', iprof )
   end if
 
   !Dimension for the number of station times
   if ( numbstat > 0 ) then
-    istation = Int ( ( xseglen - xtstep ) / tstation%step ) + 1
+    istation = Nint ( ( xseglen - dyn_model(1)%xtstep ) / tstation%step ) + 1
     call IO_Add_dim_nc4( tpfile, NMNHDIM_STATION_TIME, 'time_station', istation )
   end if
 
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index 53053069c936e71c6df416594babe000686782a7..4d30222db2026b5f5a85fe3ae84242ae36771db7 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -28,6 +28,7 @@
 !  P. Wautelet 14/01/2021: add IO_Field_write_nc4_N4, IO_Field_partial_write_nc4_N2,
 !                          IO_Field_partial_write_nc4_N3 and IO_Field_partial_write_nc4_N4 subroutines
 !  P. Wautelet 30/03/2021: budgets: LES cartesian subdomain limits are defined in the physical domain
+!  P. Wautelet 22/03/2022: correct time_les_avg and time_les_avg_bounds coordinates
 !-----------------------------------------------------------------
 #ifdef MNH_IOCDF4
 module mode_io_write_nc4
@@ -86,6 +87,8 @@ use modd_parameters, only: jphext
 
 use mode_tools_ll,   only: Get_globaldims_ll
 
+use NETCDF,          only: NF90_FLOAT
+
 type(tfiledata),       intent(in) :: tpfile
 type(tfielddata),      intent(in) :: tpfield
 integer,               intent(in) :: knblocks
@@ -112,7 +115,11 @@ call IO_Mnhname_clean( tpfield%cmnhname, yvarname )
 istatus = NF90_INQ_VARID( incid, yvarname, ivarid )
 if ( istatus /= NF90_NOERR ) then
 
-  istatus = NF90_DEF_VAR( incid, yvarname, MNHREAL_NF90, ivarid)
+  if ( tpfile%lncreduce_float_precision ) then
+    istatus = NF90_DEF_VAR( incid, yvarname, NF90_FLOAT, ivarid )
+  else
+    istatus = NF90_DEF_VAR( incid, yvarname, MNHREAL_NF90, ivarid )
+  end if
 
   if ( tpfield%ndims /= 3 ) call Print_msg( NVERB_FATAL, 'IO', 'IO_Field_header_split_write_nc4', &
                   trim( tpfile%cname )//': '//trim( yvarname )//': NDIMS should be 3' )
@@ -1444,9 +1451,9 @@ use modd_grid,       only: xlatori, xlonori
 use modd_grid_n,     only: lsleve, xxhat, xyhat, xzhat
 use modd_les,        only: cles_level_type, cspectra_level_type, nlesn_iinf, nlesn_isup, nlesn_jinf, nlesn_jsup, &
                            nles_k, nles_levels, nspectra_k, nspectra_levels,                                     &
-                           xles_altitudes, xles_temp_mean_start, xles_temp_mean_end, xles_temp_mean_step,        &
-                           xspectra_altitudes
-use modd_les_n,      only: nles_times, nspectra_ni, nspectra_nj, tles_dates
+                           xles_altitudes, xspectra_altitudes
+use modd_les_n,      only: nles_dtcount, nles_mean_end, nles_mean_start, nles_mean_step, nles_mean_times, &
+                           nles_times, nspectra_ni, nspectra_nj, tles_dates, xles_times
 use modd_netcdf,     only: tdimnc
 use modd_parameters, only: jphext, JPVEXT
 use modd_profiler_n, only: numbprofiler, tprofiler
@@ -1470,14 +1477,13 @@ character(len=:),                         allocatable :: yprogram
 integer                                               :: iiu, iju, iku
 integer                                               :: id, iid, iresp
 integer                                               :: imi
-integer                                               :: iavg
 integer                                               :: ji
 integer                                               :: jt
+integer                                               :: jtb, jte
 integer(kind=cdfint)                                  :: incid
 logical                                               :: gchangemodel
 logical                                               :: gdealloc
 logical,                         pointer              :: gsleve
-real                                                  :: zles_temp_mean_start, zles_temp_mean_end
 real,            dimension(:),   pointer              :: zxhat, zyhat, zzhat
 real,            dimension(:),            allocatable :: zxhatm, zyhatm, zzhatm !Coordinates at mass points in the transformed space
 real,            dimension(:),            allocatable :: zles_levels
@@ -1740,31 +1746,34 @@ if ( tpfile%lmaster ) then
       call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_LES_TIME), 'time axis for LES budgets', tles_dates )
 
     !Coordinates for the number of LES budget time averages
-    iavg = int( xles_temp_mean_end - 1.e-10 - xles_temp_mean_start ) / xles_temp_mean_step + 1
     !Condition also on nles_times to not create this coordinate when not used (no time average if nles_times=0)
-    if ( nles_times > 0 .and. iavg > 0 ) then
-      Allocate( tzdates(iavg) )
-      Allocate( tzdates_bound(2, iavg) )
-
-      do jt = 1, iavg
-        zles_temp_mean_start = xles_temp_mean_start + ( jt - 1 ) * xles_temp_mean_step
-        zles_temp_mean_end   = Min( xles_temp_mean_end, zles_temp_mean_start + xles_temp_mean_step )
+    if ( nles_times > 0 .and. nles_mean_times > 0 ) then
+      Allocate( tzdates(nles_mean_times) )
+      Allocate( tzdates_bound(2, nles_mean_times) )
+
+      do jt = 1, nles_mean_times
+        jtb = ( nles_mean_start + ( jt - 1 ) * nles_mean_step ) / nles_dtcount
+        jte = MIN( jtb + nles_mean_step / nles_dtcount, nles_mean_end / nles_dtcount, nles_times )
+        ! jtb could be 0 if nles_mean_start is smaller than the first LES measurement
+        ! For example, it occurs if xles_temp_mean_start is smaller than xles_temp_sampling (if xles_temp_mean_start=0.)
+        ! Do this correction only after computation of jte
+        if ( jtb < 1 ) jtb = 1
 
         tzdates(jt)%nyear  = tdtseg%nyear
         tzdates(jt)%nmonth = tdtseg%nmonth
         tzdates(jt)%nday   = tdtseg%nday
-        tzdates(jt)%xtime        = tdtseg%xtime + ( zles_temp_mean_start + zles_temp_mean_end ) / 2.
+        tzdates(jt)%xtime        = tdtseg%xtime + ( xles_times(jtb) + xles_times(jte) ) / 2.
         !Not necessary:  call Datetime_correctdate( tzdates(jt ) )
 
         tzdates_bound(1, jt)%nyear  = tdtseg%nyear
         tzdates_bound(1, jt)%nmonth = tdtseg%nmonth
         tzdates_bound(1, jt)%nday   = tdtseg%nday
-        tzdates_bound(1, jt)%xtime        = tdtseg%xtime + zles_temp_mean_start
+        tzdates_bound(1, jt)%xtime  = tdtseg%xtime + xles_times(jtb)
 
         tzdates_bound(2, jt)%nyear  = tdtseg%nyear
         tzdates_bound(2, jt)%nmonth = tdtseg%nmonth
         tzdates_bound(2, jt)%nday   = tdtseg%nday
-        tzdates_bound(2, jt)%xtime        = tdtseg%xtime + zles_temp_mean_end
+        tzdates_bound(2, jt)%xtime  = tdtseg%xtime + xles_times(jte)
       end do
       call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_LES_AVG_TIME), 'time axis for LES budget time averages', &
                              tzdates, tzdates_bound )
diff --git a/src/MNH/contrav.f90 b/src/MNH/contrav.f90
index 50e81468802efef8d210ae1736c96011e3bbbede..c890435f523160d94c0da938e847bf36c21b1da6 100644
--- a/src/MNH/contrav.f90
+++ b/src/MNH/contrav.f90
@@ -201,8 +201,20 @@ END IF
 !              ------------------------------------
 !
 IF (LFLAT) THEN
-  PRWCT(:,:,:) = PRWT(:,:,:) / PDZZ(:,:,:)
-  RETURN
+   !
+   PRWCT(:,:,:) = PRWT(:,:,:) / PDZZ(:,:,:)
+   !
+   IF (KADV_ORDER == 4 ) THEN
+      NULLIFY(TZFIELD_U)
+      NULLIFY(TZFIELD_V)
+      CALL ADD3DFIELD_ll( TZFIELD_U, PRUCT, 'CONTRAV::PRUCT' )
+      CALL ADD3DFIELD_ll( TZFIELD_V, PRVCT, 'CONTRAV::PRVCT' )
+      CALL UPDATE_HALO_ll(TZFIELD_U,IINFO_ll)
+      CALL UPDATE_HALO_ll(TZFIELD_V,IINFO_ll)
+   ENDIF
+   !
+   RETURN
+   !
 END IF
 !
 !*       3.    Compute the vertical contravariant components (general case)
diff --git a/src/MNH/ground_paramn.f90 b/src/MNH/ground_paramn.f90
index b504d4b1beeac8a1347455eac33184fbd2a63e7f..68e38a1c36babb073228eb2a0a08dc11ada394e0 100644
--- a/src/MNH/ground_paramn.f90
+++ b/src/MNH/ground_paramn.f90
@@ -111,6 +111,7 @@ END MODULE MODI_GROUND_PARAM_n
 !!     (V. Vionnet)           18/07/2017 add coupling for blowing snow module 
 !!     (Bielli S.) 02/2019  Sea salt : significant sea wave height influences salt emission; 5 salt modes
 !  P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!  P. Wautelet 09/02/2022: bugfix: add missing XCURRENT_LEI computation
 !-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -122,7 +123,7 @@ USE MODI_GET_HALO
 USE MODI_MNH_OASIS_RECV
 USE MODI_MNH_OASIS_SEND
 USE MODD_SFX_OASIS, ONLY : LOASIS
-USE MODD_DYN, ONLY : XSEGLEN
+USE MODD_DYN,        ONLY: DYN_MODEL, XSEGLEN
 #endif
 !
 USE MODD_LUNIT_n, ONLY: TLUOUT
@@ -335,7 +336,8 @@ REAL, DIMENSION(:),   ALLOCATABLE :: ZP_PET_B_COEF
 REAL, DIMENSION(:),   ALLOCATABLE :: ZP_PEQ_B_COEF
 REAL, DIMENSION(:),   ALLOCATABLE :: ZP_RN        ! net radiation           (W/m2)
 REAL, DIMENSION(:),   ALLOCATABLE :: ZP_H         ! sensible heat flux      (W/m2)
-REAL, DIMENSION(:),   ALLOCATABLE :: ZP_LE        ! latent heat flux        (W/m2)
+REAL, DIMENSION(:),   ALLOCATABLE :: ZP_LE        ! Total latent heat flux  (W/m2)
+REAL, DIMENSION(:),   ALLOCATABLE :: ZP_LEI       ! Solid Latent heat flux  (W/m2)
 REAL, DIMENSION(:),   ALLOCATABLE :: ZP_GFLUX     ! ground flux             (W/m2)
 REAL, DIMENSION(:),   ALLOCATABLE :: ZP_T2M       ! Air temperature at 2 meters (K)
 REAL, DIMENSION(:),   ALLOCATABLE :: ZP_Q2M       ! Air humidity at 2 meters    (kg/kg)
@@ -573,7 +575,7 @@ CALL DATETIME_DISTANCE(TDTSEG,TDTCUR,ZTIMEC)
 #ifdef CPLOASIS
 IF (LOASIS) THEN
   IF ( MOD(ZTIMEC,1.0) .LE. 1E-2 .OR. (1.0 - MOD(ZTIMEC,1.0)) .LE. 1E-2 ) THEN
-    IF ( NINT(ZTIMEC-(XSEGLEN-XTSTEP)) .LT. 0 ) THEN
+    IF ( NINT(ZTIMEC-(XSEGLEN-DYN_MODEL(1)%XTSTEP)) .LT. 0 ) THEN
       WRITE(ILUOUT,*) '----------------------------'
       WRITE(ILUOUT,*) ' Reception des champs avec OASIS'
       WRITE(ILUOUT,*) 'NINT(ZTIMEC)=', NINT(ZTIMEC)
@@ -603,7 +605,7 @@ CALL COUPLING_SURF_ATM_n(YSURF_CUR,'MESONH', 'E',ZTIMEC,
 #ifdef CPLOASIS
 IF (LOASIS) THEN
   IF ( MOD(ZTIMEC,1.0) .LE. 1E-2 .OR. (1.0 - MOD(ZTIMEC,1.0)) .LE. 1E-2 ) THEN
-    IF (NINT(ZTIMEC-(XSEGLEN-XTSTEP)) .LT. 0) THEN
+    IF (NINT(ZTIMEC-(XSEGLEN-DYN_MODEL(1)%XTSTEP)) .LT. 0) THEN
       WRITE(ILUOUT,*) '----------------------------'
       WRITE(ILUOUT,*) ' Envoi des champs avec OASIS'
       WRITE(ILUOUT,*) 'NINT(ZTIMEC)=', NINT(ZTIMEC)
@@ -616,9 +618,9 @@ END IF
 !
 IF (CPROGRAM=='DIAG  ' .OR. LDIAG_IN_RUN) THEN
   CALL DIAG_SURF_ATM_n(YSURF_CUR,'MESONH')
-  CALL  MNHGET_SURF_PARAM_n(PRN=ZP_RN,PH=ZP_H,PLE=ZP_LE,PGFLUX=ZP_GFLUX, &
-                           PT2M=ZP_T2M,PQ2M=ZP_Q2M,PHU2M=ZP_HU2M,        &
-                           PZON10M=ZP_ZON10M,PMER10M=ZP_MER10M           )
+  CALL  MNHGET_SURF_PARAM_n( PRN = ZP_RN,         PH = ZP_H,           PLE = ZP_LE,   PLEI = ZP_LEI,   &
+                             PGFLUX = ZP_GFLUX,   PT2M = ZP_T2M,       PQ2M = ZP_Q2M, PHU2M = ZP_HU2M, &
+                             PZON10M = ZP_ZON10M, PMER10M = ZP_MER10M                                  )
 END IF
 !
 ! Transform 1D output fields into 2D:
@@ -841,6 +843,7 @@ ALLOCATE(ZP_QSURF   (KDIM1D))
 ALLOCATE(ZP_RN      (KDIM1D))
 ALLOCATE(ZP_H       (KDIM1D))
 ALLOCATE(ZP_LE      (KDIM1D))
+ALLOCATE(ZP_LEI     (KDIM1D))
 ALLOCATE(ZP_GFLUX   (KDIM1D))
 ALLOCATE(ZP_T2M     (KDIM1D))
 ALLOCATE(ZP_Q2M     (KDIM1D))
@@ -964,6 +967,7 @@ IF (LDIAG_IN_RUN) THEN
   XCURRENT_RN      (IIB:IIE,IJB:IJE)  = RESHAPE(ZP_RN(:),     ISHAPE_2)
   XCURRENT_H       (IIB:IIE,IJB:IJE)  = RESHAPE(ZP_H (:),     ISHAPE_2)
   XCURRENT_LE      (IIB:IIE,IJB:IJE)  = RESHAPE(ZP_LE(:),     ISHAPE_2)
+  XCURRENT_LEI     (IIB:IIE,IJB:IJE)  = RESHAPE(ZP_LEI(:),    ISHAPE_2)
   XCURRENT_GFLUX   (IIB:IIE,IJB:IJE)  = RESHAPE(ZP_GFLUX(:),  ISHAPE_2)
   XCURRENT_T2M     (IIB:IIE,IJB:IJE)  = RESHAPE(ZP_T2M(:),    ISHAPE_2)
   XCURRENT_Q2M     (IIB:IIE,IJB:IJE)  = RESHAPE(ZP_Q2M(:),    ISHAPE_2)
@@ -1012,6 +1016,7 @@ DEALLOCATE(ZP_EMIS    )
 DEALLOCATE(ZP_RN      )
 DEALLOCATE(ZP_H       )
 DEALLOCATE(ZP_LE      )
+DEALLOCATE(ZP_LEI     )
 DEALLOCATE(ZP_GFLUX   )
 DEALLOCATE(ZP_T2M     )
 DEALLOCATE(ZP_Q2M     )
diff --git a/src/MNH/ini_aircraft_balloon.f90 b/src/MNH/ini_aircraft_balloon.f90
index 0b22d34031958913fa591c450709ede432d98301..cdb990bb73ae06e3440e85858b68161bc9806c92 100644
--- a/src/MNH/ini_aircraft_balloon.f90
+++ b/src/MNH/ini_aircraft_balloon.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2000-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2000-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -360,7 +360,7 @@ IF (IMI /= TPFLYER%NMODEL .AND. .NOT. (IMI==1 .AND. TPFLYER%NMODEL==0) ) RETURN
 IF ( CPROGRAM == 'DIAG  ' ) THEN
   ISTORE = INT ( NTIME_AIRCRAFT_BALLOON / TPFLYER%STEP ) + 1
 ELSE
-  ISTORE = INT ( (PSEGLEN-XTSTEP) / TPFLYER%STEP ) + 1
+  ISTORE = NINT ( ( PSEGLEN - DYN_MODEL(1)%XTSTEP ) / TPFLYER%STEP ) + 1
 ENDIF
 !
 IF (TPFLYER%NMODEL == 0) ISTORE=0
diff --git a/src/MNH/ini_lesn.f90 b/src/MNH/ini_lesn.f90
index 45d7b9f9ac85e947fda34a09e56fe335c79d1d1c..84f90dba5a61e12523a215ab1e7f07ae3627a66d 100644
--- a/src/MNH/ini_lesn.f90
+++ b/src/MNH/ini_lesn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2000-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2000-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -40,6 +40,7 @@
 !  P. Wautelet 04/01/2021: bugfix: nles_k was used instead of nspectra_k for a loop index
 !  P. Wautelet 30/03/2021: budgets: LES cartesian subdomain limits are defined in the physical domain
 !  P. Wautelet 09/07/2021: bugfix: altitude levels are on the correct grid position (mass point)
+!  P. Wautelet 22/03/2022: LES averaging periods are more reliable (compute with integers instead of reals)
 ! --------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -321,7 +322,7 @@ XLES_TEMP_SAMPLING = XTSTEP * NLES_DTCOUNT
 !            ----------------------------------------
 !
 !
-NLES_TIMES = ( INT( (XSEGLEN-XTSTEP+1.E-6) / XTSTEP ) ) / NLES_DTCOUNT
+NLES_TIMES = ( NINT( ( XSEGLEN - DYN_MODEL(1)%XTSTEP ) / XTSTEP ) ) / NLES_DTCOUNT
 !
 !*      3.5  current LES time counter
 !            ------------------------
@@ -342,6 +343,46 @@ IF (NLES_TIMES==0) THEN
   RETURN
 END IF
 !
+!*     3.8  Averaging
+!           ---------
+IF (     XLES_TEMP_MEAN_END   == XUNDEF &
+    .OR. XLES_TEMP_MEAN_START == XUNDEF &
+    .OR. XLES_TEMP_MEAN_STEP  == XUNDEF ) THEN
+  !No LES temporal averaging
+  NLES_MEAN_TIMES = 0
+  NLES_MEAN_STEP  = NNEGUNDEF
+  NLES_MEAN_START = NNEGUNDEF
+  NLES_MEAN_END   = NNEGUNDEF
+ELSE
+  !LES temporal averaging is enabled
+  !Ensure that XLES_TEMP_MEAN_END is not after segment end
+  XLES_TEMP_MEAN_END = MIN( XLES_TEMP_MEAN_END, XSEGLEN - DYN_MODEL(1)%XTSTEP )
+
+  NLES_MEAN_START = NINT( XLES_TEMP_MEAN_START / XTSTEP )
+
+  IF ( MODULO( NLES_MEAN_START, NLES_DTCOUNT ) /= 0 ) THEN
+    CMNHMSG(1) = 'XLES_TEMP_MEAN_START is not a multiple of XLES_TEMP_SAMPLING'
+    CMNHMSG(2) = 'LES averaging periods could be wrong'
+    CALL Print_msg( NVERB_WARNING, 'IO', 'INI_LES_n' )
+  END IF
+
+  NLES_MEAN_END = NINT( XLES_TEMP_MEAN_END / XTSTEP )
+
+  NLES_MEAN_STEP = NINT( XLES_TEMP_MEAN_STEP / XTSTEP )
+
+  IF ( NLES_MEAN_STEP < NLES_DTCOUNT ) &
+    CALL Print_msg( NVERB_ERROR, 'IO', 'INI_LES_n', 'XLES_TEMP_MEAN_STEP < XLES_TEMP_SAMPLING not allowed' )
+
+  IF ( MODULO( NLES_MEAN_STEP, NLES_DTCOUNT ) /= 0 ) THEN
+    CMNHMSG(1) = 'XLES_TEMP_MEAN_STEP is not a multiple of XLES_TEMP_SAMPLING'
+    CMNHMSG(2) = 'LES averaging periods could be wrong'
+    CALL Print_msg( NVERB_WARNING, 'IO', 'INI_LES_n' )
+  END IF
+
+  NLES_MEAN_TIMES = ( NLES_MEAN_END - NLES_MEAN_START ) / NLES_MEAN_STEP
+  !Add 1 averaging period if the last one is incomplete (for example: start=0., end=10., step=3.)
+  IF ( MODULO( NLES_MEAN_END - NLES_MEAN_START, NLES_MEAN_STEP ) > 0 ) NLES_MEAN_TIMES = NLES_MEAN_TIMES + 1
+END IF
 !-------------------------------------------------------------------------------
 !
 !*      4.   Number of vertical levels for local diagnostics
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 718f11ce615bb12340a0c91c7d2b514630be805a..d5a6ffa05b054a5ebe8d8ef848ec4cf44b3cef0e 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -1798,7 +1798,7 @@ gles = lles_mean .or. lles_resolved  .or. lles_subgrid .or. lles_updraft &
                  .or. lles_downdraft .or. lles_spectra
 !Called if budgets are enabled via NAM_BUDGET
 !or if LES budgets are enabled via NAM_LES (condition on kmi==1 to call it max once)
-if ( ( cbutype /= "NONE" .and. nbumod == kmi ) .or. ( gles .and. kmi == 1 ) .or. LCHECK ) THEN
+if ( ( cbutype /= "NONE" .and. nbumod == kmi ) .or. ( ( gles .or. lcheck ) .and. kmi == 1 ) ) THEN
   call Budget_preallocate()
 end if
 
diff --git a/src/MNH/ini_posprofilern.f90 b/src/MNH/ini_posprofilern.f90
index 08b25fba13d2dfdc78cdaf669b31ee81a1bdd62f..00279c8a81488360ebcc1e2ce4bb677094e55361 100644
--- a/src/MNH/ini_posprofilern.f90
+++ b/src/MNH/ini_posprofilern.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -156,7 +156,7 @@ SUBROUTINE ALLOCATE_PROFILER_n(TPROFILER)
 !
 TYPE(PROFILER), INTENT(INOUT) :: TPROFILER
 !
-ISTORE = INT ( (PSEGLEN-XTSTEP) / TPROFILER%STEP ) + 1
+ISTORE = NINT( ( PSEGLEN - DYN_MODEL(1)%XTSTEP ) / TPROFILER%STEP ) + 1
 !
 allocate( tprofiler%tpdates( istore ) )
 ALLOCATE(TPROFILER%ERROR (NUMBPROFILER))
diff --git a/src/MNH/ini_surfstationn.f90 b/src/MNH/ini_surfstationn.f90
index f53ee35d11e73fd6873e8e7abf40d90a227579af..5b79bb289026b25a3b04bdcc223ec9434345cefd 100644
--- a/src/MNH/ini_surfstationn.f90
+++ b/src/MNH/ini_surfstationn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -171,11 +171,7 @@ if ( tstation%step < xtstep ) then
   tstation%step = xtstep
 end if
 
-IF (KMI==1) THEN
- ISTORE = NINT ( (PSEGLEN-XTSTEP) / TSTATION%STEP ) + 1
-ELSE 
- ISTORE = NINT ( (PSEGLEN-XTSTEP * NDTRATIO(KMI)) / TSTATION%STEP ) + 1
-END IF
+ISTORE = NINT ( ( PSEGLEN - DYN_MODEL(1)%XTSTEP ) / TSTATION%STEP ) + 1
 
 allocate( tstation%tpdates( istore ) )
 ALLOCATE(TSTATION%ERROR (NUMBSTAT))
diff --git a/src/MNH/initial_guess.f90 b/src/MNH/initial_guess.f90
index 004d6be13168a80d05113e16940d8d9283283f72..cba1058bc1b554620e674215843be01b2218f9ae 100644
--- a/src/MNH/initial_guess.f90
+++ b/src/MNH/initial_guess.f90
@@ -153,6 +153,7 @@ use modd_budget,     only: lbudget_u,  lbudget_v,  lbudget_w,  lbudget_th, lbudg
                            NBUDGET_RR, NBUDGET_RI, NBUDGET_RS, NBUDGET_RG, NBUDGET_RH,  NBUDGET_SV1,             &
                            lbu_beg, lbu_enable, tbudgets
 USE MODD_CONF
+use modd_conf_n, only: luserv
 USE MODD_GRID_n
 
 use mode_budget,     only: Budget_store_init, Budget_store_end
@@ -213,11 +214,14 @@ IF (SIZE(PTKET,1) /= 0) THEN
   PRTKES(:,:,:) = PTKET(:,:,:) * ZINVTSTEP * PRHODJ(:,:,:)
 END IF
 !
-! Case with KRR moist variables 
-DO JRR=1,KRR
-  PRRS(:,:,:,JRR) = PRT(:,:,:,JRR) * ZINVTSTEP * PRHODJ(:,:,:) 
-END DO
-CALL MPPDB_CHECK3DM("initial_guess:PRRS/RT/RHO",PRECISION,PRRS(:,:,:,1) , PRT(:,:,:,1) , PRHODJ)
+! Case with KRR moist variables
+!
+IF (LUSERV) THEN
+   DO JRR=1,KRR
+      PRRS(:,:,:,JRR) = PRT(:,:,:,JRR) * ZINVTSTEP * PRHODJ(:,:,:)
+   END DO
+   CALL MPPDB_CHECK3DM("initial_guess:PRRS/RT/RHO",PRECISION,PRRS(:,:,:,1) , PRT(:,:,:,1) , PRHODJ)
+ENDIF
 !
 ! *** passive tracers
 !
diff --git a/src/MNH/isocom.f b/src/MNH/isocom.f
index d92460fe1ac19f237be0ac2bd6c7a9fa2c79deab..42481f83a3f489686c70fa8a696f9ffe8966ff04 100644
--- a/src/MNH/isocom.f
+++ b/src/MNH/isocom.f
@@ -1,4 +1,4 @@
-CMNH_LIC Copyright 1996-2021 CNRS, Meteo-France and Universite Paul Sabatier
+CMNH_LIC Copyright 1996-2022 CNRS, Meteo-France and Universite Paul Sabatier
 CMNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 CMNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 CMNH_LIC for details. version 1.
diff --git a/src/MNH/modd_lesn.f90 b/src/MNH/modd_lesn.f90
index 28db43c4d13ba4f85754ace59ad988145b5f1c75..c7b54dd8d4922ec77317e5abaa5711ca49f32e11 100644
--- a/src/MNH/modd_lesn.f90
+++ b/src/MNH/modd_lesn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -42,6 +42,7 @@
 !  P. Wautelet 08/02/2019: add missing NULL association for pointers
 !  C. Lac         02/2019: add rain fraction as a LES diagnostic
 !  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
+!  P. Wautelet 22/03/2022: LES averaging periods are more reliable (compute with integers instead of reals)
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -61,6 +62,10 @@ TYPE LES_t
   INTEGER :: NLES_TIMES         ! number of LES computations in time
   INTEGER :: NLES_DTCOUNT       ! number of time steps between two LES comp.
   INTEGER :: NLES_TCOUNT        ! current time counter for LES comp.
+  INTEGER :: NLES_MEAN_TIMES    ! Number of LES averaging periods
+  INTEGER :: NLES_MEAN_STEP     ! number of time steps between two LES average comp.
+  INTEGER :: NLES_MEAN_START    ! First time step number taken into account for LES averaging
+  INTEGER :: NLES_MEAN_END      ! Last  time step number taken into account for LES averaging
   INTEGER :: NSPECTRA_NI        ! number of wave lengths in I direction
   INTEGER :: NSPECTRA_NJ        ! number of wave lengths in J direction
 !
@@ -660,6 +665,10 @@ TYPE(LES_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: LES_MODEL
 INTEGER, POINTER :: NLES_TIMES=>NULL()
 INTEGER, POINTER :: NLES_DTCOUNT=>NULL()
 INTEGER, POINTER :: NLES_TCOUNT=>NULL()
+INTEGER, POINTER :: NLES_MEAN_TIMES => NULL()
+INTEGER, POINTER :: NLES_MEAN_STEP  => NULL()
+INTEGER, POINTER :: NLES_MEAN_START => NULL()
+INTEGER, POINTER :: NLES_MEAN_END   => NULL()
 INTEGER, POINTER :: NSPECTRA_NI=>NULL()
 INTEGER, POINTER :: NSPECTRA_NJ=>NULL()
 type(date_time), dimension(:), pointer :: tles_dates => null()
@@ -1509,6 +1518,10 @@ LES_MODEL(KFROM)%XLES_RADEFF=>XLES_RADEFF
 NLES_TIMES=>LES_MODEL(KTO)%NLES_TIMES
 NLES_DTCOUNT=>LES_MODEL(KTO)%NLES_DTCOUNT
 NLES_TCOUNT=>LES_MODEL(KTO)%NLES_TCOUNT
+NLES_MEAN_TIMES => LES_MODEL(KTO)%NLES_MEAN_TIMES
+NLES_MEAN_STEP  => LES_MODEL(KTO)%NLES_MEAN_STEP
+NLES_MEAN_START => LES_MODEL(KTO)%NLES_MEAN_START
+NLES_MEAN_END   => LES_MODEL(KTO)%NLES_MEAN_END
 NSPECTRA_NI=>LES_MODEL(KTO)%NSPECTRA_NI
 NSPECTRA_NJ=>LES_MODEL(KTO)%NSPECTRA_NJ
 tles_dates=>les_model(kto)%tles_dates
diff --git a/src/MNH/mode_les_diachro.f90 b/src/MNH/mode_les_diachro.f90
index f56b2ce0f607015dde207609bd257cbdcf1050da..7d8d698722653f0bb5255fcbe99e88a33992bf0e 100644
--- a/src/MNH/mode_les_diachro.f90
+++ b/src/MNH/mode_les_diachro.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -12,6 +12,7 @@
 !  P. Wautelet    10/2020: restructure subroutines to use tfield_metadata_base type
 !  P. Wautelet 03/03/2021: budgets: add tbudiachrometadata type (useful to pass more information to Write_diachro)
 !  P. Wautelet 11/03/2021: budgets: remove ptrajx/y/z optional dummy arguments of Write_diachro
+!  P. Wautelet 22/03/2022: LES averaging periods are more reliable (compute with integers instead of reals)
 !-----------------------------------------------------------------
 !#######################
 MODULE MODE_LES_DIACHRO
@@ -524,8 +525,7 @@ END SUBROUTINE LES_TIME_AVG
 subroutine Les_time_avg_4D( pwork4, tpdates, kresp )
 !########################################################
 
-use modd_les,        only: nles_current_times, xles_temp_mean_start, xles_temp_mean_end, xles_temp_mean_step
-use modd_parameters, only: XUNDEF
+use modd_les_n,      only: nles_dtcount, nles_mean_start, nles_mean_end, nles_mean_step, nles_mean_times, nles_times
 use modd_time,       only: tdtseg
 use modd_type_date,  only: date_time
 
@@ -537,60 +537,33 @@ real,            dimension(:,:,:,:), allocatable, intent(inout) :: pwork4 ! cont
 type(date_time), dimension(:),       allocatable, intent(inout) :: tpdates
 integer,                                          intent(out)   :: kresp  ! return code (0 is ok)
 !------------------------------------------------------------------------------
-integer                               :: jt                   ! time counter
-integer                               :: itime                ! nb of avg. points
-integer                               :: iavg                 ! nb of avg. periods
 integer                               :: javg                 ! loop counter on avg. periods
 integer                               :: jk                   ! vertical loop counter
 integer                               :: jp                   ! process loop counter
 integer                               :: jsv                  ! scalar loop counter
-integer                               :: jx                   ! first  spatial or spectral coordinate loop counter
-integer                               :: jy                   ! second spatial or spectral coordinate loop counter
 integer                               :: jtb, jte
-real                                  :: zles_temp_mean_start ! initial and end times
-real                                  :: zles_temp_mean_end   ! of one avergaing preiod
 real, dimension(:,:,:,:), allocatable :: zwork4               ! contains averaged physical field
 !------------------------------------------------------------------------------
 
-if (     xles_temp_mean_end   == XUNDEF &
-    .or. xles_temp_mean_start == XUNDEF &
-    .or. xles_temp_mean_step  == XUNDEF ) then
-  kresp = -1
-  return
-end if
-
-iavg = Int( xles_temp_mean_end - 1.e-10 - xles_temp_mean_start ) / xles_temp_mean_step + 1
-if ( iavg <= 0 ) then
+if ( nles_mean_times == 0 ) then
   kresp = -1
   return
 end if
 
 Deallocate( tpdates )
 
-Allocate( tpdates(iavg) )
-Allocate( zwork4(Size( pwork4, 1 ), iavg, Size( pwork4, 3 ), Size( pwork4, 4 )) )
+Allocate( tpdates(nles_mean_times) )
+Allocate( zwork4(Size( pwork4, 1 ), nles_mean_times, Size( pwork4, 3 ), Size( pwork4, 4 )) )
 
 zwork4(:, :, :, :) = 0.
 
-do javg = 1, iavg
-  zles_temp_mean_start = xles_temp_mean_start + (javg - 1) * xles_temp_mean_step
-  zles_temp_mean_end   = Min( xles_temp_mean_end, zles_temp_mean_start + xles_temp_mean_step )
-
-  jtb = -1
-  jte = -2
-  do jt = 1, nles_current_times
-    if ( xles_times(jt) >= zles_temp_mean_start ) then
-      jtb = jt
-      exit
-    end if
-  end do
-  do jt = jtb, nles_current_times
-    if ( xles_times(jt) <= zles_temp_mean_end ) then
-      jte = jt
-    else
-      exit
-    end if
-  end do
+do javg = 1, nles_mean_times
+  jtb = ( nles_mean_start + ( javg - 1 ) * nles_mean_step ) / nles_dtcount
+  jte = MIN( jtb + nles_mean_step / nles_dtcount, nles_mean_end / nles_dtcount, nles_times )
+  ! jtb could be 0 if nles_mean_start is smaller than the first LES measurement
+  ! For example, it occurs if xles_temp_mean_start is smaller than xles_temp_sampling (if xles_temp_mean_start=0.)
+  ! Do this correction only after computation of jte
+  if ( jtb < 1 ) jtb = 1
 
   do jp = 1, Size( pwork4, 4 )
     do jsv = 1, Size( pwork4, 3 )
@@ -603,12 +576,12 @@ do javg = 1, iavg
   tpdates(javg)%nyear  = tdtseg%nyear
   tpdates(javg)%nmonth = tdtseg%nmonth
   tpdates(javg)%nday   = tdtseg%nday
-  tpdates(javg)%xtime  = tdtseg%xtime + ( zles_temp_mean_start + zles_temp_mean_end ) / 2.
+  tpdates(javg)%xtime  = tdtseg%xtime + ( xles_times(jtb) + xles_times(jte) ) / 2.
   call Datetime_correctdate( tpdates(javg) )
 end do
 
 Deallocate( pwork4 )
-Allocate( pwork4(Size( zwork4, 1 ), iavg, Size( zwork4, 3 ), Size( zwork4, 4 )) )
+Allocate( pwork4(Size( zwork4, 1 ), nles_mean_times, Size( zwork4, 3 ), Size( zwork4, 4 )) )
 
 pwork4 = zwork4
 
diff --git a/src/MNH/profilern.f90 b/src/MNH/profilern.f90
index c1eabd5e17de3dff4cbce78e0f26fbb09c915329..cd276826f40564aab3c831b571789555869614fc 100644
--- a/src/MNH/profilern.f90
+++ b/src/MNH/profilern.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2002-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -85,6 +85,9 @@ END MODULE MODI_PROFILER_n
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
 !  M. Taufour  05/07/2021: modify RARE for hydrometeors containing ice and add bright band calculation for RARE
+!  P. Wautelet 09/02/2022: add message when some variables not computed
+!                          + bugfix: put values in variables in this case
+!                          + move some operations outside a do loop
 ! --------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
@@ -103,6 +106,7 @@ USE MODD_TIME,           only: tdtexp
 USE MODD_TIME_n,         only: tdtcur
 !
 USE MODE_ll
+USE MODE_MSG
 !
 USE MODI_GPS_ZENITH_GRID
 USE MODI_LIDAR
@@ -515,6 +519,14 @@ IF (GSTORE) THEN
         TPROFILER%ZTD(IN,I)= ZZTD_PROFILER
         TPROFILER%ZWD(IN,I)= ZZWD_PROFILER
         TPROFILER%ZHD(IN,I)= ZZHD_PROFILER
+      ELSE
+        CMNHMSG(1) = 'altitude of profiler ' // TRIM( TPROFILER%NAME(I) ) // ' is too far from orography'
+        CMNHMSG(2) = 'some variables are therefore not computed (IWV, ZTD, ZWD, ZHD)'
+        CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'PROFILER_n' )
+        TPROFILER%IWV(IN,I)= XUNDEF
+        TPROFILER%ZTD(IN,I)= XUNDEF
+        TPROFILER%ZWD(IN,I)= XUNDEF
+        TPROFILER%ZHD(IN,I)= XUNDEF
       END IF
       TPROFILER%ZON (IN,:,I) = ZU_PROFILER(:) * COS(ZGAM) + ZV_PROFILER(:) * SIN(ZGAM)
       TPROFILER%MER (IN,:,I) = - ZU_PROFILER(:) * SIN(ZGAM) + ZV_PROFILER(:) * COS(ZGAM)
@@ -886,10 +898,7 @@ IF (GSTORE) THEN
   CALL DISTRIBUTE_PROFILER(TPROFILER%CRARE(IN,JK,I))
   CALL DISTRIBUTE_PROFILER(TPROFILER%CRARE_ATT(IN,JK,I))
   CALL DISTRIBUTE_PROFILER(TPROFILER%CIZ(IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%IWV(IN,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZTD(IN,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZHD(IN,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZWD(IN,I))
+
   !
   IF (LDIAG_IN_RUN) CALL DISTRIBUTE_PROFILER(TPROFILER%TKE_DISS(IN,JK,I))
   !
@@ -901,6 +910,11 @@ IF (GSTORE) THEN
   END DO
   IF (SIZE(PTKE)>0) CALL DISTRIBUTE_PROFILER(TPROFILER%TKE  (IN,JK,I))
  ENDDO
+
+  CALL DISTRIBUTE_PROFILER(TPROFILER%IWV(IN,I))
+  CALL DISTRIBUTE_PROFILER(TPROFILER%ZTD(IN,I))
+  CALL DISTRIBUTE_PROFILER(TPROFILER%ZHD(IN,I))
+  CALL DISTRIBUTE_PROFILER(TPROFILER%ZWD(IN,I))
 ENDDO
 !
 END IF
diff --git a/src/MNH/read_chem_data_netcdf_case.f90 b/src/MNH/read_chem_data_netcdf_case.f90
index 547a55b535791a53149ef6e20f232706ddfbff71..3ecc3f594c36347e046515c7e9a5995490a3ad02 100644
--- a/src/MNH/read_chem_data_netcdf_case.f90
+++ b/src/MNH/read_chem_data_netcdf_case.f90
@@ -87,6 +87,7 @@ END MODULE MODI_READ_CHEM_DATA_NETCDF_CASE
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 18/09/2019: correct support of 64bit integers (MNH_INT=8)
 !  P. Wautelet 09/03/2021: move some chemistry initializations to ini_nsv
+!  P. Wautelet 01/04/2022: add error if CDUMMY1 not set correctly
 !-------------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -395,6 +396,9 @@ ELSEIF (CDUMMY1=="18") THEN
        itimeindex=3
 ELSEIF ((CDUMMY1=="24").OR.(CDUMMY1=="00")) THEN
        itimeindex=4
+ELSE
+  call Print_msg( NVERB_ERROR, 'GEN', 'READ_CHEM_DATA_NETCDF_CASE', 'CDUMMY1 is not set correctly (or not set at all)' )
+  itimeindex=1
 ENDIF
  istart3d(4) = itimeindex
 !
@@ -423,7 +427,7 @@ enddo
   istatus = nf90_get_var(incid, ips_varid, ZPSMOZ(:,:), start=istart2d, count=icount2d)
   if (istatus /= nf90_noerr) call handle_err(istatus)
 
-  
+
 !------------------------------------------------------------------------
 !* 3 Interpolation of MOZART variable
 !---------------------------------------------------------------------
@@ -580,7 +584,7 @@ DO JI = 1,IMOZ               !for every MNH species existing in MOZ1.nam
     ENDIF      
 
   ENDDO ! JNCHEM
-  DO JNAER = NSV_AERBEG, NSV_AEREND 
+  DO JNAER = NSV_AERBEG, NSV_AEREND
     IF (trim(CAERONAMES(JNAER-NSV_AERBEG+1))==trim(YSPCMNH(JI))) THEN !MNH mechanism species
        IF (ISPCMOZ(JI)==1) THEN
          istatus = nf90_inq_varid(incid, trim(YCHANGE(JI,1)), ind_netcdf)
@@ -661,7 +665,7 @@ DO JI = 1,IMOZ               !for every MNH species existing in MOZ1.nam
     ENDIF         
   ENDDO ! JNAER
 ENDDO  ! JIDO JNCHEM = NSV_CHEMBEG, NSV_CHEMEND  !loop on all MNH species
-DEALLOCATE(YSPCMNH) 
+DEALLOCATE(YSPCMNH)
 DEALLOCATE(TZSTOC)
 DEALLOCATE(ISPCMOZ) 
 DEALLOCATE(ZCOEFMOZART)
diff --git a/src/MNH/read_surf_mnh.f90 b/src/MNH/read_surf_mnh.f90
index f051ed7393a5eb3042200c278fb0c6af3ebcbf0a..30b2a93712df54151a1f2dd8927e3e2de626d8d9 100644
--- a/src/MNH/read_surf_mnh.f90
+++ b/src/MNH/read_surf_mnh.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2003-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2003-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -1023,7 +1023,7 @@ IF (.NOT. GCOVER_PACKED) THEN
   TZFIELD%CSTDNAME   = ''
   TZFIELD%CLONGNAME  = TRIM(YREC)
   TZFIELD%CUNITS     = ''
-  TZFIELD%CDIR       = YDIR1
+  TZFIELD%CDIR       = YDIR
   TZFIELD%CCOMMENT   = 'X_Y_'//TRIM(YREC)
   TZFIELD%NGRID      = 4
   TZFIELD%NTYPE      = TYPEREAL
diff --git a/src/MNH/stationn.f90 b/src/MNH/stationn.f90
index 4de047e716b3d60dda687a4abe41342d8f6dabcb..8774f316af1407a1f08b151fd469aa351461ae21 100644
--- a/src/MNH/stationn.f90
+++ b/src/MNH/stationn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2002-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-2022 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -80,6 +80,7 @@ END MODULE MODI_STATION_n
 !  P. Wautelet  05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet  13/09/2019: budget: simplify and modernize date/time management
 !  R. Schoetter    11/2019: use LCARTESIAN instead of LSTATLAT for multiproc in cartesian
+!  P. Wautelet  09/05/2022: bugfix: use correct indices for U and V interpolation
 !
 ! --------------------------------------------------------------------------
 !
@@ -491,7 +492,7 @@ ELSEIF (L1D) THEN
      JI=2
      JJ=2
 ELSE
-     JI=II(I)
+     JI=IU(I)
      JJ=IJ(I)
 END IF
 !
@@ -521,7 +522,7 @@ ELSEIF (L1D) THEN
      JJ=2  
 ELSE
   JI=II(I)
-  JJ=IJ(I)
+  JJ=IV(I)
 END IF
 !
 IF ((JI .GT. 0).AND. (JI .LT. SIZE(PA,1)) .AND. &
diff --git a/src/SURFEX/pgd_grid.F90 b/src/SURFEX/pgd_grid.F90
index 5a7158fc109dd922e185beef9c6b51261a77ee8e..9c17d153a3428715d0cb1593b01b253b90c94fb8 100644
--- a/src/SURFEX/pgd_grid.F90
+++ b/src/SURFEX/pgd_grid.F90
@@ -316,7 +316,7 @@ ALLOCATE(UG%XJPDIR       (NL))
  CALL MPPDB_CHECK_SURFEX2D(UG%G%XLAT,"PGD_GRID after LATLON_GRID:XLAT",PRECISION,ILUOUT)
  CALL MPPDB_CHECK_SURFEX2D(UG%G%XLON,"PGD_GRID after LATLON_GRID:XLON",PRECISION,ILUOUT)
  CALL MPPDB_CHECK_SURFEX2D(UG%G%XMESH_SIZE,"PGD_GRID after LATLON_GRID:XMESH_SIZE",PRECISION,ILUOUT)
-#endif!
+#endif
 !------------------------------------------------------------------------------
 !
 !*    7.      Average grid length (in degrees)