diff --git a/src/MNH/mode_les_diachro.f90 b/src/MNH/mode_les_diachro.f90
index b3c518177fc78530bfe1736e101976db2f861e36..f56b2ce0f607015dde207609bd257cbdcf1050da 100644
--- a/src/MNH/mode_les_diachro.f90
+++ b/src/MNH/mode_les_diachro.f90
@@ -753,9 +753,10 @@ end if
 
 end subroutine Les_diachro_2D
 
-!###################################################################################################################
-subroutine Les_diachro_3D( tpdiafile, tpfield, hgroup, hgroupcomment, odoavg, odonorm, pfield, hfieldnames, hmasks )
-!###################################################################################################################
+!###############################################################################################
+subroutine Les_diachro_3D( tpdiafile, tpfield, hgroup, hgroupcomment, odoavg, odonorm, pfield, &
+                           hfieldnames, hfieldcomments, hmasks )
+!###############################################################################################
 
 use modd_field, only: NMNHDIM_BUDGET_LES_LEVEL, NMNHDIM_BUDGET_LES_MASK, NMNHDIM_BUDGET_LES_SV, &
                       NMNHDIM_BUDGET_LES_TIME,  NMNHDIM_BUDGET_TERM,     NMNHDIM_UNUSED,        &
@@ -770,6 +771,7 @@ logical,                                   intent(in) :: odoavg     ! Compute an
 logical,                                   intent(in) :: odonorm    ! Compute and store normalized field
 real,                    dimension(:,:,:), intent(in) :: pfield     ! Data array
 character(len=*),        dimension(:),     optional, intent(in) :: hfieldnames
+character(len=*),        dimension(:),     optional, intent(in) :: hfieldcomments
 character(len=*),        dimension(:),     optional, intent(in) :: hmasks
 
 type(tfield_metadata_base) :: tzfield
@@ -808,13 +810,20 @@ else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
                     'optional dummy argument hfieldnames is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
 
+  if ( .not. Present( hfieldcomments ) ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
+                    'optional dummy argument hfieldcomments is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
   if ( Size( hfieldnames ) /= Size( pfield, 3) ) &
     call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_3D', 'wrong size for hfieldnames (' // Trim( tzfield%cmnhname ) // ')' )
 
+  if ( Size( hfieldcomments ) /= Size( pfield, 3) ) &
+    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_3D', 'wrong size for hfieldcomments (' // Trim( tzfield%cmnhname ) // ')' )
+
   tzfield%ndimlist(4) = NMNHDIM_UNUSED
   call Les_diachro_common( tpdiafile, tzfield, hgroup, hgroupcomment,                                         &
                            reshape( pfield, [ size( pfield, 1 ), size( pfield, 2 ), size( pfield, 3 ), 1 ] ), &
-                           odoavg, odonorm, hfieldnames = hfieldnames )
+                           odoavg, odonorm, hfieldnames = hfieldnames, hfieldcomments = hfieldcomments )
 else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
           .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME  &
           .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_SV    ) then
@@ -822,6 +831,10 @@ else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
                     'optional dummy argument hfieldnames is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
 
+  if ( Present( hfieldcomments ) ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
+                    'optional dummy argument hfieldcomments is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
   if ( Present( hmasks ) ) &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
                     'optional dummy argument hmasks is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
@@ -838,9 +851,10 @@ end if
 
 end subroutine Les_diachro_3D
 
-!###################################################################################################################
-subroutine Les_diachro_4D( tpdiafile, tpfield, hgroup, hgroupcomment, odoavg, odonorm, pfield, hfieldnames, hmasks )
-!###################################################################################################################
+!###############################################################################################
+subroutine Les_diachro_4D( tpdiafile, tpfield, hgroup, hgroupcomment, odoavg, odonorm, pfield, &
+                           hfieldnames, hfieldcomments, hmasks )
+!###############################################################################################
 
 use modd_field, only: NMNHDIM_BUDGET_LES_LEVEL, NMNHDIM_BUDGET_LES_MASK, NMNHDIM_BUDGET_LES_PDF,  NMNHDIM_BUDGET_LES_SV, &
                       NMNHDIM_BUDGET_LES_TIME,  NMNHDIM_BUDGET_TERM,     NMNHDIM_UNUSED,                                 &
@@ -855,6 +869,7 @@ logical,                                     intent(in) :: odoavg     ! Compute
 logical,                                     intent(in) :: odonorm    ! Compute and store normalized field
 real,                    dimension(:,:,:,:), intent(in) :: pfield     ! Data array
 character(len=*),        dimension(:),     optional, intent(in) :: hfieldnames
+character(len=*),        dimension(:),     optional, intent(in) :: hfieldcomments
 character(len=*),        dimension(:),     optional, intent(in) :: hmasks
 
 type(tfield_metadata_base) :: tzfield
@@ -884,6 +899,10 @@ if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL&
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
                     'optional dummy argument hfieldnames is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
 
+  if ( Present( hfieldcomments ) ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
+                    'optional dummy argument hfieldcomments is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
   if ( Size( hmasks ) /= Size( pfield, 3) ) &
     call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_4D', 'wrong size for hmasks (' // Trim( tzfield%cmnhname ) // ')' )
 
@@ -895,6 +914,10 @@ else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
                     'optional dummy argument hfieldnames is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
 
+  if ( .not. Present( hfieldcomments ) ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
+                    'optional dummy argument hfieldcomments is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
   if ( Present( hmasks ) ) &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
                     'optional dummy argument hmasks is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
@@ -902,6 +925,9 @@ else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
   if ( Size( hfieldnames ) /= Size( pfield, 3) ) &
     call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_4D', 'wrong size for hfieldnames (' // Trim( tzfield%cmnhname ) // ')' )
 
+  if ( Size( hfieldcomments ) /= Size( pfield, 3) ) &
+    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_4D', 'wrong size for hfieldcomments (' // Trim( tzfield%cmnhname ) // ')' )
+
 else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
           .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME  &
           .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_MASK      &
@@ -914,6 +940,10 @@ else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
                     'optional dummy argument hfieldnames is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
 
+  if ( Present( hfieldcomments ) ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
+                    'optional dummy argument hfieldcomments is not needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
   if ( Size( hmasks ) /= Size( pfield, 3) ) &
     call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_4D', 'wrong size for hmasks (' // Trim( tzfield%cmnhname ) // ')' )
 
@@ -922,14 +952,15 @@ else
                   'ndimlist configuration not yet implemented for ' // Trim( tzfield%cmnhname ) )
 end if
 
-call Les_diachro_common( tpdiafile, tzfield, hgroup, hgroupcomment, pfield, odoavg, odonorm, &
-                         hfieldnames = hfieldnames, hmasks = hmasks )
+call Les_diachro_common( tpdiafile, tzfield, hgroup, hgroupcomment, pfield, odoavg, odonorm,         &
+                         hfieldnames = hfieldnames, hfieldcomments = hfieldcomments, hmasks = hmasks )
 
 end subroutine Les_diachro_4D
 
-!#######################################################################################################################
-subroutine Les_diachro_common( tpdiafile, tpfield, hgroup, hgroupcomment, pfield, odoavg, odonorm, hfieldnames, hmasks )
-!#######################################################################################################################
+!###################################################################################################
+subroutine Les_diachro_common( tpdiafile, tpfield, hgroup, hgroupcomment, pfield, odoavg, odonorm, &
+                               hfieldnames, hfieldcomments, hmasks )
+!###################################################################################################
 
 use modd_field,         only: tfield_metadata_base
 use modd_io,            only: tfiledata
@@ -948,6 +979,7 @@ real,                       dimension(:,:,:,:),           intent(in) :: pfield
 logical,                                                  intent(in) :: odoavg    ! Compute and store time average
 logical,                                                  intent(in) :: odonorm   ! Compute and store normalized field
 character(len=*),           dimension(:),       optional, intent(in) :: hfieldnames
+character(len=*),           dimension(:),       optional, intent(in) :: hfieldcomments
 character(len=*),           dimension(:),       optional, intent(in) :: hmasks
 
 character(len=100),         dimension(:),     allocatable :: ycomment                      ! Comment string
@@ -975,13 +1007,13 @@ ijh = nles_current_jsup
 ikl = nles_levels(1)
 ikh = nles_levels(iles_k)
 
-if ( Present( hfieldnames ) ) then
+if ( Present( hfieldcomments ) ) then
   if ( Present( hmasks ) ) &
-    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'hfieldnames and hmasks optional arguments may not be present ' // &
+    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'hfieldcomments and hmasks optional arguments may not be present ' // &
                     'at the same time (' // Trim( tpfield%cmnhname ) // ')' )
-  if ( Size( hfieldnames ) /= Size( pfield, 3) ) &
-    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'wrong size for hfieldnames (' // Trim( tpfield%cmnhname ) // ')' )
-  ycomment(:) = Trim( tpfield%ccomment(:) ) // ': ' // hfieldnames(:)
+  if ( Size( hfieldcomments ) /= Size( pfield, 3) ) &
+    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'wrong size for hfieldcomments (' // Trim( tpfield%cmnhname ) // ')' )
+  ycomment(:) = Trim( tpfield%ccomment(:) ) // ': ' // hfieldcomments(:)
 else if ( Present( hmasks ) ) then
   if ( Size( hmasks ) /= Size( pfield, 3) ) &
     call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'wrong size for hmasks (' // Trim( tpfield%cmnhname ) // ')' )
diff --git a/src/MNH/write_les_budgetn.f90 b/src/MNH/write_les_budgetn.f90
index 9745903d66dc6893fe21f4badd918576f7b54e46..7827ba184b18d710839d63155d807c186b67eb01 100644
--- a/src/MNH/write_les_budgetn.f90
+++ b/src/MNH/write_les_budgetn.f90
@@ -98,7 +98,8 @@ INTEGER :: ILES_P1, ILES_P2
 INTEGER :: JK ! vertical loop counter
 INTEGER :: JT ! temporal loop counter
 !
-CHARACTER(len=9), DIMENSION(NMAX_ILES) :: YFIELDNAMES
+CHARACTER(len=9),   DIMENSION(NMAX_ILES) :: YFIELDNAMES
+CHARACTER(len=100), DIMENSION(NMAX_ILES) :: YFIELDCOMMENTS
 character(len=:), allocatable          :: ygroup
 character(len=:), allocatable          :: ygroupcomment
 !
@@ -141,16 +142,18 @@ ILES_STA=ILES
 !     --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TEND'
+YFIELDNAMES(ILES)    = 'SBG_TEND'
+YFIELDCOMMENTS(ILES) = 'subgrid tendency'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_TEND)
 !
 !
-!* 1.2 production by mean wind gradient
+!* 1.2 dynamic production by mean wind gradient
 !     ---------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 !
 ZLES_BUDGET(:,:,ILES)= - XLES_SUBGRID_WU (:,:,1) * XLES_MEAN_DUDZ(:,:,1)  &
                        - XLES_SUBGRID_WV (:,:,1) * XLES_MEAN_DVDZ(:,:,1)  &
@@ -160,7 +163,8 @@ ZLES_BUDGET(:,:,ILES)= - XLES_SUBGRID_WU (:,:,1) * XLES_MEAN_DUDZ(:,:,1)  &
 !     ------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_DP) - ZLES_BUDGET(:,:,2)
 !
@@ -171,7 +175,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_DP) - ZLES_BUDGET(:,:,2)
 !
 IF ( ANY(XLES_BU_SBG_TKE(:,:,NLES_ADVM)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_ADVM'
+YFIELDNAMES(ILES)    = 'SBG_ADVM'
+YFIELDCOMMENTS(ILES) = 'subgrid advection by mean flow'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_ADVM)
 END IF
@@ -182,7 +187,8 @@ END IF
 !
 IF ( ANY(XLES_BU_SBG_TKE(:,:,NLES_FORC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_FORC'
+YFIELDNAMES(ILES)    = 'SBG_FORC'
+YFIELDCOMMENTS(ILES) = 'subgrid advection by large-scale W forcing'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_FORC)
 END IF
@@ -192,7 +198,8 @@ END IF
 !      -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_TR)
 !
@@ -201,7 +208,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_TR)
 !      -----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_ADVR'
+YFIELDNAMES(ILES)    = 'SBG_ADVR'
+YFIELDCOMMENTS(ILES) = 'subgrid advection by resolved flow'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_ADVR)
 !
@@ -211,7 +219,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_ADVR)
 !
 IF ( ANY(XLES_BU_SBG_TKE(:,:,NLES_PRES)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_PRES'
+YFIELDNAMES(ILES)    = 'SBG_PRES'
+YFIELDCOMMENTS(ILES) = 'subgrid pressure-correlation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_PRES)
 END IF
@@ -221,7 +230,8 @@ END IF
 !      ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TP'
+YFIELDNAMES(ILES)    = 'SBG_TP'
+YFIELDCOMMENTS(ILES) = 'subgrid thermal production'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_TP)
 !
@@ -230,7 +240,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_TP)
 !       -----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DISS'
+YFIELDNAMES(ILES)    = 'SBG_DISS'
+YFIELDCOMMENTS(ILES) = 'subgrid dissipation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_DISS)
 !
@@ -240,7 +251,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_DISS)
 !
 IF ( ANY(XLES_BU_SBG_TKE(:,:,NLES_DIFF)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_NUMD'
+YFIELDNAMES(ILES)    = 'SBG_NUMD'
+YFIELDCOMMENTS(ILES) = 'subgrid numerical diffusion'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_DIFF)
 END IF
@@ -250,7 +262,8 @@ END IF
 !
 IF ( ANY(XLES_BU_SBG_TKE(:,:,NLES_RELA)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RELA'
+YFIELDNAMES(ILES)    = 'SBG_RELA'
+YFIELDCOMMENTS(ILES) = 'subgrid sponge layer relaxation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_RELA)
 END IF
@@ -260,7 +273,8 @@ END IF
 !
 IF ( ANY(XLES_BU_SBG_TKE(:,:,NLES_NEST)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_NEST'
+YFIELDNAMES(ILES)    = 'SBG_NEST'
+YFIELDCOMMENTS(ILES) = 'subgrid average from smaller nested models'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_NEST)
 END IF
@@ -271,7 +285,8 @@ END IF
 !
 IF ( ANY(XLES_BU_SBG_TKE(:,:,NLES_MISC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_MISC'
+YFIELDNAMES(ILES)    = 'SBG_MISC'
+YFIELDCOMMENTS(ILES) = 'subgrid: other effects'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_SBG_TKE(:,:,NLES_MISC)
 END IF
@@ -281,7 +296,8 @@ END IF
 !       --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -294,7 +310,8 @@ ILES_STA=ILES
 !       --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_TEND)
 !
@@ -304,7 +321,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_TEND)
 !
 IF ( ANY(XLES_BU_RES_KE(:,:,NLES_ADVM)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_ADV'
+YFIELDNAMES(ILES)    = 'RES_ADV'
+YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_Ke(:,:,NLES_ADVM)
 END IF
@@ -315,17 +333,19 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_KE(:,:,NLES_FORC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_FORC'
+YFIELDNAMES(ILES)    = 'RES_FORC'
+YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_Ke(:,:,NLES_FORC)
 END IF
 !
 !
-!* 1.19 production by mean wind gradient
+!* 1.19 dynamic production by mean wind gradient
 !       --------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Ke(:,:,NLES_DP)
 !
@@ -334,7 +354,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Ke(:,:,NLES_DP)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_TR)
 !
@@ -344,7 +365,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_TR)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_PRES'
+YFIELDNAMES(ILES)    = 'RES_PRES'
+YFIELDCOMMENTS(ILES) = 'resolved pressure-correlation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_PRES)
 !
@@ -353,7 +375,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_PRES)
 !       ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TP'
+YFIELDNAMES(ILES)    = 'RES_TP'
+YFIELDCOMMENTS(ILES) = 'resolved thermal production'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_GRAV)
 !
@@ -362,7 +385,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_GRAV)
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_VTURB) + XLES_BU_RES_KE(:,:,NLES_HTURB)
 !
@@ -371,7 +395,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_VTURB) + XLES_BU_RES_KE(:,:,NLES
 !
 IF ( ANY(XLES_BU_RES_KE(:,:,NLES_COR)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_CORI'
+YFIELDNAMES(ILES)    = 'RES_CORI'
+YFIELDCOMMENTS(ILES) = 'resolved Coriolis effect'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_COR)
 END IF
@@ -381,7 +406,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_KE(:,:,NLES_DIFF)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NUMD'
+YFIELDNAMES(ILES)    = 'RES_NUMD'
+YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_DIFF)
 END IF
@@ -391,7 +417,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_KE(:,:,NLES_RELA)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RELA'
+YFIELDNAMES(ILES)    = 'RES_RELA'
+YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_RELA)
 END IF
@@ -401,7 +428,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_KE(:,:,NLES_NEST)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NEST'
+YFIELDNAMES(ILES)    = 'RES_NEST'
+YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_NEST)
 END IF
@@ -413,7 +441,8 @@ END IF
 IF ( ANY( XLES_BU_RES_KE(:,:,NLES_MISC) &
          +XLES_BU_RES_KE(:,:,NLES_CURV) /= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_MISC'
+YFIELDNAMES(ILES)    = 'RES_MISC'
+YFIELDCOMMENTS(ILES) = 'resolved: other effects'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_KE(:,:,NLES_MISC)  &
                       + XLES_BU_RES_KE(:,:,NLES_CURV)
@@ -423,7 +452,8 @@ END IF
 !       ------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -440,8 +470,8 @@ tzfield%clongname = ygroup !clongname will be overwritten by yfieldnames(:) in L
 tzfield%ccomment  = 'resolved KE budget' !ccomment will be completed with yfieldnames(:) in Les_diachro
 tzfield%cunits    = 'm2 s-3'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -455,11 +485,12 @@ ILES=0
 !
 ILES_STA=ILES
 !
-!* 2.1 production by mean gradients
+!* 2.1 dynamic production by mean gradients
 !      ----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 ILES_P1=ILES
 !
 ZLES_BUDGET(:,:,ILES)= - 2. * XLES_SUBGRID_WThl(:,:,1) * XLES_MEAN_dThldz(:,:,1)
@@ -470,7 +501,8 @@ ZLES_BUDGET(:,:,ILES)= - 2. * XLES_SUBGRID_WThl(:,:,1) * XLES_MEAN_dThldz(:,:,1)
 !
 IF ( ANY(XLES_SUBGRID_WThl2(:,:,1)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES) = - ( XLES_SUBGRID_WThl2 (JK+1,:,1)      &
@@ -487,7 +519,8 @@ END IF
 !      --------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 ILES_P2=ILES
 !
 ZLES_BUDGET(:,:,ILES)= - 2. * XLES_RES_ddxa_Thl_SBG_UaThl(:,:,1)  &
@@ -498,7 +531,8 @@ ZLES_BUDGET(:,:,ILES)= - 2. * XLES_RES_ddxa_Thl_SBG_UaThl(:,:,1)  &
 !      -----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DISS'
+YFIELDNAMES(ILES)    = 'SBG_DISS'
+YFIELDCOMMENTS(ILES) = 'subgrid dissipation'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_DISS_Thl2(:,:,1)
 !
@@ -507,7 +541,8 @@ ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_DISS_Thl2(:,:,1)
 !      --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -520,7 +555,8 @@ ILES_STA=ILES
 !      --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_TEND)
 !
@@ -530,7 +566,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_TEND)
 !
 IF ( ANY(XLES_BU_RES_Thl2(:,:,NLES_ADVM)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_ADV'
+YFIELDNAMES(ILES)    = 'RES_ADV'
+YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_ADVM)
 END IF
@@ -541,7 +578,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Thl2(:,:,NLES_FORC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_FORC'
+YFIELDNAMES(ILES)    = 'RES_FORC'
+YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_FORC)
 END IF
@@ -551,7 +589,8 @@ END IF
 !      ----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_DP)
 
@@ -560,7 +599,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_DP)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_TR)
 !
@@ -569,7 +609,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_TR)
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_VTURB) + XLES_BU_RES_Thl2(:,:,NLES_HTURB)
 !
@@ -578,7 +619,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_VTURB) + XLES_BU_RES_Thl2(:,:,
 !
 IF ( ANY(XLES_BU_RES_Thl2(:,:,NLES_DIFF)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NUMD'
+YFIELDNAMES(ILES)    = 'RES_NUMD'
+YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_DIFF)
 END IF
@@ -588,7 +630,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Thl2(:,:,NLES_RELA)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RELA'
+YFIELDNAMES(ILES)    = 'RES_RELA'
+YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_RELA)
 END IF
@@ -598,7 +641,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Thl2(:,:,NLES_NEST)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NEST'
+YFIELDNAMES(ILES)    = 'RES_NEST'
+YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_NEST)
 END IF
@@ -611,7 +655,8 @@ IF ( ANY( XLES_BU_RES_Thl2(:,:,NLES_MISC) &
          +XLES_BU_RES_Thl2(:,:,NLES_MICR) &
          + XLES_BU_RES_Thl2(:,:,NLES_PREF) /= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_MISC'
+YFIELDNAMES(ILES)    = 'RES_MISC'
+YFIELDCOMMENTS(ILES) = 'resolved: other effects'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Thl2(:,:,NLES_MISC) &
                       + XLES_BU_RES_Thl2(:,:,NLES_RAD ) &
@@ -624,7 +669,8 @@ END IF
 !       ---------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -635,7 +681,8 @@ END DO
 !       ------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_TEND'
+YFIELDNAMES(ILES)    = 'NSG_TEND'
+YFIELDCOMMENTS(ILES) = 'neglected tendency'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 IF (NLES_TIMES>2) THEN
@@ -654,7 +701,8 @@ END IF
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVM'
+YFIELDNAMES(ILES)    = 'NSG_ADVM'
+YFIELDCOMMENTS(ILES) = 'neglected advection by mean flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)=  -XLES_MEAN_W(JK,:,1)                 &
@@ -671,7 +719,8 @@ ZLES_BUDGET(NLES_K,:,ILES) = ZLES_BUDGET(NLES_K-1,:,ILES)
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVR'
+YFIELDNAMES(ILES)    = 'NSG_ADVR'
+YFIELDCOMMENTS(ILES) = 'neglected advection by resolved flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)= - ( XLES_RES_W_SBG_Thl2   (JK+1,:,1)    &
@@ -694,8 +743,8 @@ tzfield%clongname = ygroup !clongname will be overwritten by yfieldnames(:) in L
 tzfield%ccomment  = 'thetal variance budget' !ccomment will be completed with yfieldnames(:) in Les_diachro
 tzfield%cunits    = 'K2 s-1'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -708,11 +757,12 @@ ILES=0
 !
 ILES_STA=ILES
 !
-!* 3.1 production by mean gradients
+!* 3.1 dynamic production by mean gradients
 !     -----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 ILES_P1=ILES
 !
 ZLES_BUDGET(:,:,ILES) =  - XLES_SUBGRID_W2(:,:,1) * XLES_MEAN_dThldz(:,:,1)
@@ -722,7 +772,8 @@ ZLES_BUDGET(:,:,ILES) =  - XLES_SUBGRID_W2(:,:,1) * XLES_MEAN_dThldz(:,:,1)
 !     -------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 !
 ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddz_Thl_SBG_W2(:,:,1) &
                       - ZLES_BUDGET(:,:,ILES_P1)
@@ -733,7 +784,8 @@ ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddz_Thl_SBG_W2(:,:,1) &
 !
 IF ( ANY(XLES_SUBGRID_W2Thl(:,:,1)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES) = - ( XLES_SUBGRID_W2Thl (JK+1,:,1)       &
@@ -752,7 +804,8 @@ END IF
 !
 IF ( ANY(XLES_SUBGRID_ThlPz(:,:,1)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_PRES'
+YFIELDNAMES(ILES)    = 'SBG_PRES'
+YFIELDCOMMENTS(ILES) = 'subgrid pressure-correlation'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_ThlPz(:,:,1)
 END IF
@@ -762,7 +815,8 @@ END IF
 !      ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TP'
+YFIELDNAMES(ILES)    = 'SBG_TP'
+YFIELDCOMMENTS(ILES) = 'subgrid thermal production'
 !
 IF (LUSERV) THEN
   ZLES_BUDGET(:,:,ILEs) =  XG * XLES_SUBGRID_ThlThv(:,:,1)   &
@@ -779,7 +833,8 @@ END IF
 !      --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -792,7 +847,8 @@ ILES_STA=ILES
 !      --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_TEND)
 !
@@ -801,7 +857,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_TEND)
 !
 IF ( ANY(XLES_BU_RES_WThl(:,:,NLES_ADVM)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_ADV'
+YFIELDNAMES(ILES)    = 'RES_ADV'
+YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_WThl(:,:,NLES_ADVM)
 END IF
@@ -811,7 +868,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WThl(:,:,NLES_FORC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_FORC'
+YFIELDNAMES(ILES)    = 'RES_FORC'
+YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_WThl(:,:,NLES_FORC)
 END IF
@@ -821,7 +879,8 @@ END IF
 !       ----------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_WThl(:,:,NLES_DP)
 !
@@ -829,7 +888,8 @@ ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_WThl(:,:,NLES_DP)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_WThl(:,:,NLES_TR)
 !
@@ -838,7 +898,8 @@ ZLES_BUDGET(:,:,ILES) =  XLES_BU_RES_WThl(:,:,NLES_TR)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_PRES'
+YFIELDNAMES(ILES)    = 'RES_PRES'
+YFIELDCOMMENTS(ILES) = 'resolved pressure-correlation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_PRES)
 !
@@ -847,7 +908,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_PRES)
 !       ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TP'
+YFIELDNAMES(ILES)    = 'RES_TP'
+YFIELDCOMMENTS(ILES) = 'resolved thermal production'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_GRAV)
 !
@@ -856,7 +918,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_GRAV)
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_VTURB) + XLES_BU_RES_WThl(:,:,NLES_HTURB)
 !
@@ -865,7 +928,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_VTURB) + XLES_BU_RES_WThl(:,:,
 !
 IF ( ANY(XLES_BU_RES_WThl(:,:,NLES_COR)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_CORI'
+YFIELDNAMES(ILES)    = 'RES_CORI'
+YFIELDCOMMENTS(ILES) = 'resolved Coriolis effect'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_COR)
 END IF
@@ -875,7 +939,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WThl(:,:,NLES_DIFF)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NUMD'
+YFIELDNAMES(ILES)    = 'RES_NUMD'
+YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_DIFF)
 END IF
@@ -885,7 +950,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WThl(:,:,NLES_RELA)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RELA'
+YFIELDNAMES(ILES)    = 'RES_RELA'
+YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_RELA)
 END IF
@@ -895,7 +961,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WThl(:,:,NLES_NEST)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NEST'
+YFIELDNAMES(ILES)    = 'RES_NEST'
+YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_NEST)
 END IF
@@ -909,7 +976,8 @@ IF ( ANY( XLES_BU_RES_WThl(:,:,NLES_MISC) &
          +XLES_BU_RES_WThl(:,:,NLES_PREF) &
          +XLES_BU_RES_WThl(:,:,NLES_CURV) /= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_MISC'
+YFIELDNAMES(ILES)    = 'RES_MISC'
+YFIELDCOMMENTS(ILES) = 'resolved: other effects'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WThl(:,:,NLES_MISC) &
                       + XLES_BU_RES_WThl(:,:,NLES_RAD ) &
@@ -923,7 +991,8 @@ END IF
 !       --------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -934,7 +1003,8 @@ END DO
 !       ------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_TEND'
+YFIELDNAMES(ILES)    = 'NSG_TEND'
+YFIELDCOMMENTS(ILES) = 'neglected tendency'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 IF (NLES_TIMES>2) THEN
@@ -955,7 +1025,8 @@ END IF
 !       ------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVM'
+YFIELDNAMES(ILES)    = 'NSG_ADVM'
+YFIELDCOMMENTS(ILES) = 'neglected advection by mean flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)= - XLES_MEAN_W(JK,:,1)                               &
@@ -972,7 +1043,8 @@ ZLES_BUDGET(NLES_K,:,ILES) = ZLES_BUDGET(NLES_K-1,:,ILES)
 !       ------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVR'
+YFIELDNAMES(ILES)    = 'NSG_ADVR'
+YFIELDCOMMENTS(ILES) = 'neglected advection by resolved flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)= - ( XLES_RES_W_SBG_WThl(JK+1,:,1)   &
@@ -988,7 +1060,8 @@ ZLES_BUDGET(NLES_K,:,ILES) = ZLES_BUDGET(NLES_K-1,:,ILES)
 !       ----------------------------------------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_DPGW'
+YFIELDNAMES(ILES)    = 'NSG_DPGW'
+YFIELDCOMMENTS(ILES) = 'neglected production by gradient of vertical velocity for subgrid quantity'
 !
 ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddxa_W_SBG_UaThl(:,:,1)
 !
@@ -997,7 +1070,8 @@ ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddxa_W_SBG_UaThl(:,:,1)
 !       -------------------------------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_DPGT'
+YFIELDNAMES(ILES)    = 'NSG_DPGT'
+YFIELDCOMMENTS(ILES) = 'neglected production by horizontal gradient of Thl for subgrid quantity'
 !
 ZLES_BUDGET(:,:,ILES)=-XLES_RES_ddxa_Thl_SBG_UaW(:,:,1)       &
                       -ZLES_BUDGET(:,:,ILES_P1) -ZLES_BUDGET(:,:,ILES_P2)
@@ -1014,8 +1088,8 @@ tzfield%clongname = ygroup !clongname will be overwritten by yfieldnames(:) in L
 tzfield%ccomment  = 'thetal flux budget' !ccomment will be completed with yfieldnames(:) in Les_diachro
 tzfield%cunits    = 'm K s-2'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/write_les_rt_budgetn.f90 b/src/MNH/write_les_rt_budgetn.f90
index cb6e13e3e5813226c6464f7929279ec37f17b6e6..114d39cda70bc47d713bcfbddc3f550c9a27d491 100644
--- a/src/MNH/write_les_rt_budgetn.f90
+++ b/src/MNH/write_les_rt_budgetn.f90
@@ -91,9 +91,10 @@ INTEGER :: ILES_P1, ILES_P2
 INTEGER :: JK ! vertical loop counter
 INTEGER :: JT ! temporal loop counter
 !
-CHARACTER(len=9), DIMENSION(NMAX_ILES) :: YFIELDNAMES
-character(len=:), allocatable          :: ygroup
-character(len=:), allocatable          :: ygroupcomment
+CHARACTER(len=9),   DIMENSION(NMAX_ILES) :: YFIELDNAMES
+CHARACTER(len=100), DIMENSION(NMAX_ILES) :: YFIELDCOMMENTS
+character(len=:),   allocatable          :: ygroup
+character(len=:),   allocatable          :: ygroupcomment
 !
 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZLES_BUDGET
 !
@@ -132,11 +133,12 @@ ILES=0
 ILES_STA=ILES
 !
 !
-!* 2.1 production by mean gradients
+!* 2.1 dynamic production by mean gradients
 !      ----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 ILES_P1=ILES
 !
 ZLES_BUDGET(:,:,ILES)= - 2. * XLES_SUBGRID_WRt(:,:,1) * XLES_MEAN_dRtdz(:,:,1)
@@ -146,7 +148,8 @@ ZLES_BUDGET(:,:,ILES)= - 2. * XLES_SUBGRID_WRt(:,:,1) * XLES_MEAN_dRtdz(:,:,1)
 !      --------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 !
 ZLES_BUDGET(:,:,ILES)= - 2. * XLES_RES_ddxa_Rt_SBG_UaRt(:,:,1)  &
                           - ZLES_BUDGET(:,:,ILES_P1)
@@ -157,7 +160,8 @@ ZLES_BUDGET(:,:,ILES)= - 2. * XLES_RES_ddxa_Rt_SBG_UaRt(:,:,1)  &
 !
 IF ( ANY(XLES_SUBGRID_WRt2(:,:,1)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES) = - ( XLES_SUBGRID_WRt2 (JK+1,:,1)      &
@@ -174,7 +178,8 @@ END IF
 !      -----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DISS'
+YFIELDNAMES(ILES)    = 'SBG_DISS'
+YFIELDCOMMENTS(ILES) = 'subgrid dissipation'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_DISS_Rt2(:,:,1)
 !
@@ -183,7 +188,8 @@ ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_DISS_Rt2(:,:,1)
 !      --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -196,7 +202,8 @@ ILES_STA=ILES
 !      --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_TEND)
 !
@@ -206,7 +213,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_TEND)
 !
 IF ( ANY(XLES_BU_RES_Rt2(:,:,NLES_ADVM)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_ADV'
+YFIELDNAMES(ILES)    = 'RES_ADV'
+YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_ADVM)
 END IF
@@ -216,7 +224,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Rt2(:,:,NLES_FORC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_FORC'
+YFIELDNAMES(ILES)    = 'RES_FORC'
+YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_FORC)
 END IF
@@ -226,7 +235,8 @@ END IF
 !      ----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_DP)
 
@@ -235,7 +245,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_DP)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_TR)
 !
@@ -244,7 +255,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_TR)
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_VTURB) + XLES_BU_RES_Rt2(:,:,NLES_HTURB)
 !
@@ -254,7 +266,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_VTURB) + XLES_BU_RES_Rt2(:,:,NL
 !
 IF ( ANY(XLES_BU_RES_Rt2(:,:,NLES_DIFF)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NUMD'
+YFIELDNAMES(ILES)    = 'RES_NUMD'
+YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_DIFF)
 END IF
@@ -264,7 +277,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Rt2(:,:,NLES_RELA)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RELA'
+YFIELDNAMES(ILES)    = 'RES_RELA'
+YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_RELA)
 END IF
@@ -274,7 +288,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Rt2(:,:,NLES_NEST)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NEST'
+YFIELDNAMES(ILES)    = 'RES_NEST'
+YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_NEST)
 END IF
@@ -285,7 +300,8 @@ END IF
 IF ( ANY( XLES_BU_RES_Rt2(:,:,NLES_MISC) &
          +XLES_BU_RES_Rt2(:,:,NLES_MICR)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_MISC'
+YFIELDNAMES(ILES)    = 'RES_MISC'
+YFIELDCOMMENTS(ILES) = 'resolved: other effects'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_Rt2(:,:,NLES_MISC) &
                       + XLES_BU_RES_Rt2(:,:,NLES_MICR)
@@ -296,7 +312,8 @@ END IF
 !       ---------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -307,7 +324,8 @@ END DO
 !       ------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_TEND'
+YFIELDNAMES(ILES)    = 'NSG_TEND'
+YFIELDCOMMENTS(ILES) = 'neglected tendency'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 IF (NLES_TIMES>2) THEN
@@ -326,7 +344,8 @@ END IF
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVM'
+YFIELDNAMES(ILES)    = 'NSG_ADVM'
+YFIELDCOMMENTS(ILES) = 'neglected advection by mean flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)=  -XLES_MEAN_W(JK,:,1)                &
@@ -343,7 +362,8 @@ ZLES_BUDGET(NLES_K,:,ILES) = ZLES_BUDGET(NLES_K-1,:,ILES)
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVR'
+YFIELDNAMES(ILES)    = 'NSG_ADVR'
+YFIELDCOMMENTS(ILES) = 'neglected advection by resolved flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)= - ( XLES_RES_W_SBG_Rt2   (JK+1,:,1)    &
@@ -365,8 +385,8 @@ tzfield%clongname = ygroup !clongname will be overwritten by yfieldnames(:) in L
 tzfield%ccomment  = 'Rt variance budget' !ccomment will be completed with yfieldnames(:) in Les_diachro
 tzfield%cunits    = 'kg2 kg-2 s-1'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -379,11 +399,12 @@ ygroupcomment = 'Total water flux budget'
 ILES=0
 ILES_STA=ILES
 !
-!* 3.1 production by mean gradients
+!* 3.1 dynamic production by mean gradients
 !     -----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 ILES_P1=ILES
 !
 ZLES_BUDGET(:,:,ILES) =  - XLES_SUBGRID_W2(:,:,1) * XLES_MEAN_dRtdz(:,:,1)
@@ -393,7 +414,8 @@ ZLES_BUDGET(:,:,ILES) =  - XLES_SUBGRID_W2(:,:,1) * XLES_MEAN_dRtdz(:,:,1)
 !     -------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(2) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 ILES_P2=ILES
 !
 ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddz_Rt_SBG_W2(:,:,1) &
@@ -406,7 +428,8 @@ ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddz_Rt_SBG_W2(:,:,1) &
 !
 IF ( ANY(XLES_SUBGRID_W2Rt(:,:,1)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES) = - ( XLES_SUBGRID_W2Rt (JK+1,:,1)       &
@@ -424,7 +447,8 @@ END IF
 !      -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_PRES'
+YFIELDNAMES(ILES)    = 'SBG_PRES'
+YFIELDCOMMENTS(ILES) = 'subgrid pressure-correlation'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_RtPz(:,:,1)
 !
@@ -433,7 +457,8 @@ ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_RtPz(:,:,1)
 !      ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TP'
+YFIELDNAMES(ILES)    = 'SBG_TP'
+YFIELDCOMMENTS(ILES) = 'subgrid thermal production'
 !
 ZLES_BUDGET(:,:,ILES) =  XG * XLES_SUBGRID_RtThv(:,:,1)   &
                             / XLES_MEAN_Thv     (:,:,1)
@@ -444,7 +469,8 @@ ZLES_BUDGET(:,:,ILES) =  XG * XLES_SUBGRID_RtThv(:,:,1)   &
 !
 !PW: not in the documentation, but set to 0 anyway
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DISS'
+YFIELDNAMES(ILES)    = 'SBG_DISS'
+YFIELDCOMMENTS(ILES) = 'subgrid dissipation'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 !
@@ -453,7 +479,8 @@ ZLES_BUDGET(:,:,ILES) = 0.
 !      --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -466,7 +493,8 @@ ILES_STA=ILES
 !      --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_TEND)
 !
@@ -475,7 +503,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_TEND)
 !
 IF ( ANY(XLES_BU_RES_WRt(:,:,NLES_ADVM)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_ADV'
+YFIELDNAMES(ILES)    = 'RES_ADV'
+YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_ADVM)
 END IF
@@ -485,7 +514,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WRt(:,:,NLES_FORC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_FORC'
+YFIELDNAMES(ILES)    = 'RES_FORC'
+YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_FORC)
 END IF
@@ -494,7 +524,8 @@ END IF
 !       ----------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_DP)
 !
@@ -502,7 +533,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_DP)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_TR)
 !
@@ -511,7 +543,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_TR)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_PRES'
+YFIELDNAMES(ILES)    = 'RES_PRES'
+YFIELDCOMMENTS(ILES) = 'resolved pressure-correlation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_PRES)
 !
@@ -520,7 +553,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_PRES)
 !       ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TP'
+YFIELDNAMES(ILES)    = 'RES_TP'
+YFIELDCOMMENTS(ILES) = 'resolved thermal production'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_GRAV)
 !
@@ -529,7 +563,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_GRAV)
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_VTURB) + XLES_BU_RES_WRt(:,:,NLES_HTURB)
 !
@@ -538,7 +573,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_VTURB) + XLES_BU_RES_WRt(:,:,NL
 !
 IF ( ANY(XLES_BU_RES_WRt(:,:,NLES_COR)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_CORI'
+YFIELDNAMES(ILES)    = 'RES_CORI'
+YFIELDCOMMENTS(ILES) = 'resolved Coriolis effect'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_COR)
 END IF
@@ -548,7 +584,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WRt(:,:,NLES_DIFF)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NUMD'
+YFIELDNAMES(ILES)    = 'RES_NUMD'
+YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_DIFF)
 END IF
@@ -558,7 +595,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WRt(:,:,NLES_RELA)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RELA'
+YFIELDNAMES(ILES)    = 'RES_RELA'
+YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_RELA)
 END IF
@@ -568,7 +606,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WRt(:,:,NLES_NEST)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NEST'
+YFIELDNAMES(ILES)    = 'RES_NEST'
+YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_NEST)
 END IF
@@ -579,7 +618,8 @@ END IF
 IF ( ANY( XLES_BU_RES_WRt(:,:,NLES_MISC) &
          +XLES_BU_RES_WRt(:,:,NLES_MICR) /= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_MISC'
+YFIELDNAMES(ILES)    = 'RES_MISC'
+YFIELDCOMMENTS(ILES) = 'resolved: other effects'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_WRt(:,:,NLES_MISC) &
                       + XLES_BU_RES_WRt(:,:,NLES_MICR)
@@ -589,7 +629,8 @@ END IF
 !       --------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -600,7 +641,8 @@ END DO
 !       ------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TEND'
+YFIELDNAMES(ILES)    = 'NSG_TEND'
+YFIELDCOMMENTS(ILES) = 'neglected tendency'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 IF (NLES_TIMES>2) THEN
@@ -621,7 +663,8 @@ END IF
 !       ------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVM'
+YFIELDNAMES(ILES)    = 'NSG_ADVM'
+YFIELDCOMMENTS(ILES) = 'neglected advection by mean flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)= - XLES_MEAN_W(JK,:,1)                              &
@@ -638,7 +681,8 @@ ZLES_BUDGET(NLES_K,:,ILES) = ZLES_BUDGET(NLES_K-1,:,ILES)
 !       ------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVR'
+YFIELDNAMES(ILES)    = 'NSG_ADVR'
+YFIELDCOMMENTS(ILES) = 'neglected advection by resolved flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)= - ( XLES_RES_W_SBG_WRt(JK+1,:,1)   &
@@ -654,7 +698,8 @@ ZLES_BUDGET(NLES_K,:,ILES) = ZLES_BUDGET(NLES_K-1,:,ILES)
 !       ----------------------------------------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_DPGW'
+YFIELDNAMES(ILES)    = 'NSG_DPGW'
+YFIELDCOMMENTS(ILES) = 'neglected production by gradient of vertical velocity for subgrid quantity'
 !
 ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddxa_W_SBG_UaRt(:,:,1)
 !
@@ -663,7 +708,8 @@ ZLES_BUDGET(:,:,ILES)=- XLES_RES_ddxa_W_SBG_UaRt(:,:,1)
 !       -------------------------------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_DPGT'
+YFIELDNAMES(ILES)    = 'NSG_DPGT'
+YFIELDCOMMENTS(ILES) = 'neglected production by horizontal gradient of Thl for subgrid quantity'
 !
 ZLES_BUDGET(:,:,ILES)=-XLES_RES_ddxa_Rt_SBG_UaW(:,:,1)       &
                       -ZLES_BUDGET(:,:,ILES_P1) -ZLES_BUDGET(:,:,ILES_P2)
@@ -678,8 +724,8 @@ tzfield%clongname = ygroup !clongname will be overwritten by yfieldnames(:) in L
 tzfield%ccomment  = 'Rt flux budget' !ccomment will be completed with yfieldnames(:) in Les_diachro
 tzfield%cunits    = 'm kg kg-1 s-2'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -694,11 +740,12 @@ ILES=0
 ILES_STA=ILES
 !
 !
-!* 2.1 production by mean gradients
+!* 2.1 dynamic production by mean gradients
 !      ----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 ILES_P1=ILES
 !
 ZLES_BUDGET(:,:,ILES)=-XLES_SUBGRID_WRt (:,:,1) * XLES_MEAN_dThldz(:,:,1) &
@@ -709,7 +756,8 @@ ZLES_BUDGET(:,:,ILES)=-XLES_SUBGRID_WRt (:,:,1) * XLES_MEAN_dThldz(:,:,1) &
 !      --------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 !
 ZLES_BUDGET(:,:,ILES)= - XLES_RES_ddxa_Rt_SBG_UaThl(:,:,1)  &
                        - XLES_RES_ddxa_Thl_SBG_UaRt(:,:,1)  &
@@ -721,7 +769,8 @@ ZLES_BUDGET(:,:,ILES)= - XLES_RES_ddxa_Rt_SBG_UaThl(:,:,1)  &
 !
 IF ( ANY(XLES_SUBGRID_WThlRt(:,:,1)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES) = - ( XLES_SUBGRID_WThlRt (JK+1,:,1)      &
@@ -738,7 +787,8 @@ END IF
 !      -----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DISS'
+YFIELDNAMES(ILES)    = 'SBG_DISS'
+YFIELDCOMMENTS(ILES) = 'subgrid dissipation'
 !
 ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_DISS_ThlRt(:,:,1)
 !
@@ -747,7 +797,8 @@ ZLES_BUDGET(:,:,ILES) =  XLES_SUBGRID_DISS_ThlRt(:,:,1)
 !      --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -760,7 +811,8 @@ ILES_STA=ILES
 !      --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_TEND)
 !
@@ -770,7 +822,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_TEND)
 !
 IF ( ANY(XLES_BU_RES_ThlRt(:,:,NLES_ADVM)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(7) = 'RES_ADV'
+YFIELDNAMES(ILES)    = 'RES_ADV'
+YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_ADVM)
 END IF
@@ -780,7 +833,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_ThlRt(:,:,NLES_FORC)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_FORC'
+YFIELDNAMES(ILES)    = 'RES_FORC'
+YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_FORC)
 END IF
@@ -789,7 +843,8 @@ END IF
 !      ----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_DP)
 !
@@ -797,7 +852,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_DP)
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_TR)
 !
@@ -806,7 +862,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_TR)
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_VTURB) + XLES_BU_RES_ThlRt(:,:,NLES_HTURB)
 !
@@ -816,7 +873,8 @@ ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_VTURB) + XLES_BU_RES_ThlRt(:,
 !
 IF ( ANY(XLES_BU_RES_ThlRt(:,:,NLES_DIFF)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NUMD'
+YFIELDNAMES(ILES)    = 'RES_NUMD'
+YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_DIFF)
 END IF
@@ -826,7 +884,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_ThlRt(:,:,NLES_RELA)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RELA'
+YFIELDNAMES(ILES)    = 'RES_RELA'
+YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_RELA)
 END IF
@@ -836,7 +895,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_ThlRt(:,:,NLES_NEST)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_NEST'
+YFIELDNAMES(ILES)    = 'RES_NEST'
+YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_NEST)
 END IF
@@ -849,7 +909,8 @@ IF ( ANY( XLES_BU_RES_ThlRt(:,:,NLES_MISC) &
          +XLES_BU_RES_ThlRt(:,:,NLES_RAD ) &
          +XLES_BU_RES_ThlRt(:,:,NLES_MICR) /= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_MISC'
+YFIELDNAMES(ILES)    = 'RES_MISC'
+YFIELDCOMMENTS(ILES) = 'resolved: other effects'
 !
 ZLES_BUDGET(:,:,ILES) = XLES_BU_RES_ThlRt(:,:,NLES_MISC) &
                       + XLES_BU_RES_ThlRt(:,:,NLES_PREF) &
@@ -862,7 +923,8 @@ END IF
 !       ---------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 DO JLES=ILES_STA+1,ILES-1
@@ -873,7 +935,8 @@ END DO
 !       ------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_TEND'
+YFIELDNAMES(ILES)    = 'NSG_TEND'
+YFIELDCOMMENTS(ILES) = 'neglected tendency'
 !
 ZLES_BUDGET(:,:,ILES) = 0.
 IF (NLES_TIMES>2) THEN
@@ -892,7 +955,8 @@ END IF
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVM'
+YFIELDNAMES(ILES)    = 'NSG_ADVM'
+YFIELDCOMMENTS(ILES) = 'neglected advection by mean flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)=  -XLES_MEAN_W(JK,:,1)                &
@@ -909,7 +973,8 @@ ZLES_BUDGET(NLES_K,:,ILES) = ZLES_BUDGET(NLES_K-1,:,ILES)
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVR'
+YFIELDNAMES(ILES)    = 'NSG_ADVR'
+YFIELDCOMMENTS(ILES) = 'neglected advection by resolved flow'
 !
 DO JK=2,NLES_K-1
   ZLES_BUDGET(JK,:,ILES)= - ( XLES_RES_W_SBG_ThlRt (JK+1,:,1)    &
@@ -931,8 +996,8 @@ tzfield%clongname = ygroup !clongname will be overwritten by yfieldnames(:) in L
 tzfield%ccomment  = 'Thl-Rt covariance budget' !ccomment will be completed with yfieldnames(:) in Les_diachro
 tzfield%cunits    = 'K kg kg-1 s-1'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/write_les_sv_budgetn.f90 b/src/MNH/write_les_sv_budgetn.f90
index bb7a206016d11dd4c932441bcef711b91142408a..6a3997964fff08e950d3de7c18fa1603f5e392ad 100644
--- a/src/MNH/write_les_sv_budgetn.f90
+++ b/src/MNH/write_les_sv_budgetn.f90
@@ -92,9 +92,10 @@ INTEGER :: JT ! temporal loop counter
 INTEGER :: JSV! scalar loop counter
 INTEGER :: JP ! process loop counter
 !
-CHARACTER(len=9), DIMENSION(NMAX_ILES)      :: YFIELDNAMES
-character(len=:), allocatable               :: ygroup
-character(len=:), allocatable               :: ygroupcomment
+CHARACTER(len=9),   DIMENSION(NMAX_ILES) :: YFIELDNAMES
+CHARACTER(len=100), DIMENSION(NMAX_ILES) :: YFIELDCOMMENTS
+character(len=:),   allocatable          :: ygroup
+character(len=:),   allocatable          :: ygroupcomment
 !
 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZLES_BUDGET
 !
@@ -121,11 +122,12 @@ ygroupcomment = 'Total scalar variance budget'
 ILES=0
 ILES_STA=ILES
 !
-!* 2.1 production by mean gradients
+!* 2.1 dynamic production by mean gradients
 !      ----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 ILES_P1=ILES
 !
 DO JSV=1,NSV
@@ -137,7 +139,8 @@ END DO
 !      --------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV)= - 2. * XLES_RES_ddxa_Sv_SBG_UaSv(:,:,1,JSV)  &
@@ -150,7 +153,8 @@ END DO
 !
 IF ( ANY(XLES_SUBGRID_WSv2(:,:,1,:)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 DO JSV=1,NSV
   DO JK=2,NLES_K-1
@@ -169,7 +173,8 @@ END IF
 !      -----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DISS'
+YFIELDNAMES(ILES)    = 'SBG_DISS'
+YFIELDCOMMENTS(ILES) = 'subgrid dissipation'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) =  XLES_SUBGRID_DISS_Sv2(:,:,1,JSV)
@@ -180,7 +185,8 @@ END DO
 !      --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = 0.
@@ -195,7 +201,8 @@ ILES_STA=ILES
 !      --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_TEND,JSV)
@@ -206,7 +213,8 @@ END DO
 !
 IF ( ANY(XLES_BU_RES_Sv2(:,:,NLES_ADVM,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_ADV'
+  YFIELDNAMES(ILES)    = 'RES_ADV'
+  YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
 !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_ADVM,JSV)
@@ -218,7 +226,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Sv2(:,:,NLES_FORC,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_FORC'
+  YFIELDNAMES(ILES)    = 'RES_FORC'
+  YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
 !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_FORC,JSV)
@@ -230,7 +239,8 @@ END IF
 !      ----------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_DP,JSV)
@@ -240,7 +250,8 @@ END DO
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_TR,JSV)
@@ -251,7 +262,8 @@ END DO
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_VTURB,JSV) &
@@ -263,7 +275,8 @@ END DO
 !
 IF ( ANY(XLES_BU_RES_Sv2(:,:,NLES_DIFF,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_NUMD'
+  YFIELDNAMES(ILES)    = 'RES_NUMD'
+  YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
 !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_DIFF,JSV)
@@ -275,7 +288,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Sv2(:,:,NLES_RELA,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_RELA'
+  YFIELDNAMES(ILES)    = 'RES_RELA'
+  YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
 !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_RELA,JSV)
@@ -287,7 +301,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Sv2(:,:,NLES_NEST,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_NEST'
+  YFIELDNAMES(ILES)    = 'RES_NEST'
+  YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
 !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_NEST,JSV)
@@ -299,7 +314,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_Sv2(:,:,NLES_MISC,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_MISC'
+  YFIELDNAMES(ILES)    = 'RES_MISC'
+  YFIELDCOMMENTS(ILES) = 'resolved: other effects'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_Sv2(:,:,NLES_MISC,JSV)
@@ -310,7 +326,8 @@ END IF
 !       ---------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = 0.
@@ -323,7 +340,8 @@ END DO
 !       ------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_TEND'
+YFIELDNAMES(ILES)    = 'NSG_TEND'
+YFIELDCOMMENTS(ILES) = 'neglected tendency'
 !
 IF (NLES_TIMES>2) THEN
   DO JSV=1,NSV
@@ -344,7 +362,8 @@ END IF
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVM'
+YFIELDNAMES(ILES)    = 'NSG_ADVM'
+YFIELDCOMMENTS(ILES) = 'neglected advection by mean flow'
 !
 DO JSV=1,NSV
   DO JK=2,NLES_K-1
@@ -362,7 +381,8 @@ END DO
 !       ----------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVR'
+YFIELDNAMES(ILES)    = 'NSG_ADVR'
+YFIELDCOMMENTS(ILES) = 'neglected advection by resolved flow'
 !
 DO JSV=1,NSV
   DO JK=2,NLES_K-1
@@ -398,8 +418,8 @@ tzfield%ndimlist(5:) = NMNHDIM_UNUSED
 gdoavg  = xles_temp_mean_start /= XUNDEF .and. xles_temp_mean_end /= XUNDEF
 gdonorm = trim(cles_norm_type) /= 'NONE'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles, :), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles, :), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -414,11 +434,12 @@ ygroupcomment = 'Total water flux budget'
 ILES=0
 ILES_STA=ILES
 !
-!* 3.1 production by mean gradients
+!* 3.1 dynamic production by mean gradients
 !     -----------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_M'
+YFIELDNAMES(ILES)    = 'SBG_DP_M'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by mean gradient'
 ILES_P1=ILES
 !
 DO JSV=1,NSV
@@ -431,7 +452,8 @@ END DO
 !     -------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_DP_R'
+YFIELDNAMES(ILES)    = 'SBG_DP_R'
+YFIELDCOMMENTS(ILES) = 'subgrid dynamic production by resolved fluctuations'
 ILES_P2=ILES
 !
 DO JSV=1,NSV
@@ -446,7 +468,8 @@ END DO
 !
 IF ( ANY(XLES_SUBGRID_W2Sv(:,:,1,:)/= 0.) ) THEN
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TR'
+YFIELDNAMES(ILES)    = 'SBG_TR'
+YFIELDCOMMENTS(ILES) = 'subgrid turbulent transport'
 !
 DO JSV=1,NSV
   DO JK=2,NLES_K-1
@@ -465,7 +488,8 @@ END IF
 !      -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_PRES'
+YFIELDNAMES(ILES)    = 'SBG_PRES'
+YFIELDCOMMENTS(ILES) = 'subgrid pressure-correlation'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) =  XLES_SUBGRID_SvPz(:,:,1,JSV)
@@ -476,7 +500,8 @@ END DO
 !      ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_TP'
+YFIELDNAMES(ILES)    = 'SBG_TP'
+YFIELDCOMMENTS(ILES) = 'subgrid thermal production'
 !
 IF (LUSERV) THEN
   DO JSV=1,NSV
@@ -496,7 +521,8 @@ END IF
 !
 ILES=ILES+1
 !PW: not in documentation. Always set to 0
-YFIELDNAMES(ILES) = 'SBG_DISS'
+YFIELDNAMES(ILES)    = 'SBG_DISS'
+YFIELDCOMMENTS(ILES) = 'subgrid dissipation'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = 0.
@@ -507,7 +533,8 @@ END DO
 !      --------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'SBG_RESI'
+YFIELDNAMES(ILES)    = 'SBG_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of subgrid budget'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = 0.
@@ -522,7 +549,8 @@ ILES_STA=ILES
 !      --------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TEND'
+YFIELDNAMES(ILES)    = 'RES_TEND'
+YFIELDCOMMENTS(ILES) = 'resolved tendency'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_TEND,JSV)
@@ -533,7 +561,8 @@ END DO
 !
 IF ( ANY(XLES_BU_RES_WSv(:,:,NLES_ADVM,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_ADV'
+  YFIELDNAMES(ILES)    = 'RES_ADV'
+  YFIELDCOMMENTS(ILES) = 'resolved advection by mean flow'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_ADVM,JSV)
@@ -545,7 +574,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WSv(:,:,NLES_FORC,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_FORC'
+  YFIELDNAMES(ILES)    = 'RES_FORC'
+  YFIELDCOMMENTS(ILES) = 'resolved advection by large-scale W forcing'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_FORC,JSV)
@@ -556,7 +586,8 @@ END IF
 !       ----------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_DP'
+YFIELDNAMES(ILES)    = 'RES_DP'
+YFIELDCOMMENTS(ILES) = 'resolved dynamic production by mean gradient'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_DP,JSV)
@@ -566,7 +597,8 @@ END DO
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TR'
+YFIELDNAMES(ILES)    = 'RES_TR'
+YFIELDCOMMENTS(ILES) = 'turbulent transport of resolved flux by itself'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_TR,JSV)
@@ -577,7 +609,8 @@ END DO
 !       -------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_PRES'
+YFIELDNAMES(ILES)    = 'RES_PRES'
+YFIELDCOMMENTS(ILES) = 'resolved pressure-correlation'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_PRES,JSV)
@@ -588,7 +621,8 @@ END DO
 !       ------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_TP'
+YFIELDNAMES(ILES)    = 'RES_TP'
+YFIELDCOMMENTS(ILES) = 'resolved thermal production'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_GRAV,JSV)
@@ -599,7 +633,8 @@ END DO
 !       ----------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_SBGT'
+YFIELDNAMES(ILES)    = 'RES_SBGT'
+YFIELDCOMMENTS(ILES) = 'resolved sink due to subgrid turbulence'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_VTURB,JSV) + XLES_BU_RES_WSv(:,:,NLES_HTURB,JSV)
@@ -610,7 +645,8 @@ END DO
 !
 IF ( ANY(XLES_BU_RES_WSv(:,:,NLES_COR,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_CORI'
+  YFIELDNAMES(ILES)    = 'RES_CORI'
+  YFIELDCOMMENTS(ILES) = 'resolved Coriolis effect'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_COR,JSV)
@@ -622,7 +658,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WSv(:,:,NLES_DIFF,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_NUMD'
+  YFIELDNAMES(ILES)    = 'RES_NUMD'
+  YFIELDCOMMENTS(ILES) = 'resolved numerical diffusion'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_DIFF,JSV)
@@ -634,7 +671,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WSv(:,:,NLES_RELA,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_RELA'
+  YFIELDNAMES(ILES)    = 'RES_RELA'
+  YFIELDCOMMENTS(ILES) = 'resolved sponge layer relaxation'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_RELA,JSV)
@@ -646,7 +684,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WSv(:,:,NLES_NEST,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_NEST'
+  YFIELDNAMES(ILES)    = 'RES_NEST'
+  YFIELDCOMMENTS(ILES) = 'resolved average from smaller nested models'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_NEST,JSV)
@@ -658,7 +697,8 @@ END IF
 !
 IF ( ANY(XLES_BU_RES_WSv(:,:,NLES_MISC,:)/= 0.) ) THEN
   ILES=ILES+1
-  YFIELDNAMES(ILES) = 'RES_MISC'
+  YFIELDNAMES(ILES)    = 'RES_MISC'
+  YFIELDCOMMENTS(ILES) = 'resolved: other effects'
   !
   DO JSV=1,NSV
     ZLES_BUDGET(:,:,ILES,JSV) = XLES_BU_RES_WSv(:,:,NLES_MISC,JSV)
@@ -669,7 +709,8 @@ END IF
 !       -------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'RES_RESI'
+YFIELDNAMES(ILES)    = 'RES_RESI'
+YFIELDCOMMENTS(ILES) = 'residual of resolved budget'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = 0.
@@ -682,7 +723,8 @@ END DO
 !       ------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_TEND'
+YFIELDNAMES(ILES)    = 'NSG_TEND'
+YFIELDCOMMENTS(ILES) = 'neglected tendency'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV) = 0.
@@ -705,7 +747,8 @@ END DO
 !       ------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVM'
+YFIELDNAMES(ILES)    = 'NSG_ADVM'
+YFIELDCOMMENTS(ILES) = 'neglected advection by mean flow'
 !
 DO JSV=1,NSV
   DO JK=2,NLES_K-1
@@ -723,7 +766,8 @@ END DO
 !       ------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_ADVR'
+YFIELDNAMES(ILES)    = 'NSG_ADVR'
+YFIELDCOMMENTS(ILES) = 'neglected advection by resolved flow'
 !
 DO JSV=1,NSV
   DO JK=2,NLES_K-1
@@ -740,7 +784,8 @@ END DO
 !       ----------------------------------------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_DPGW'
+YFIELDNAMES(ILES)    = 'NSG_DPGW'
+YFIELDCOMMENTS(ILES) = 'neglected production by gradient of vertical velocity for subgrid quantity'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV)=- XLES_RES_ddxa_W_SBG_UaSv(:,:,1,JSV)
@@ -751,7 +796,8 @@ END DO
 !       -------------------------------------------------------------------------
 !
 ILES=ILES+1
-YFIELDNAMES(ILES) = 'NSG_DPGT'
+YFIELDNAMES(ILES)    = 'NSG_DPGT'
+YFIELDCOMMENTS(ILES) = 'neglected production by horizontal gradient of Thl for subgrid quantity'
 !
 DO JSV=1,NSV
   ZLES_BUDGET(:,:,ILES,JSV)=-XLES_RES_ddxa_Sv_SBG_UaW(:,:,1,JSV)       &
@@ -781,8 +827,8 @@ tzfield%ndimlist(5:) = NMNHDIM_UNUSED
 gdoavg  = xles_temp_mean_start /= XUNDEF .and. xles_temp_mean_end /= XUNDEF
 gdonorm = trim(cles_norm_type) /= 'NONE'
 
-call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, &
-                  zles_budget(:, :, :iles, :), hfieldnames = yfieldnames(:iles) )
+call Les_diachro( tpdiafile, tzfield, ygroup, ygroupcomment, gdoavg, gdonorm, zles_budget(:, :, :iles, :), &
+                  hfieldnames = yfieldnames(:iles), hfieldcomments = yfieldcomments(:iles) )
 
 !-------------------------------------------------------------------------------
 !