From b5b92619b00fec4f0e798bcbbfbbb0258a82d10f Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Thu, 6 May 2021 14:24:33 +0200
Subject: [PATCH] Philippe 06/05/2021: budgets: add mask information in
 tbudiachrometadata + transmit this information

---
 src/MNH/modd_budget.f90            |   1 +
 src/MNH/mode_les_diachro.f90       | 109 +++++--
 src/MNH/write_aircraft_balloon.f90 |   2 +
 src/MNH/write_budget.f90           |  10 +-
 src/MNH/write_diachro.f90          |   8 +-
 src/MNH/write_les_budgetn.f90      |   8 +-
 src/MNH/write_les_rt_budgetn.f90   |   8 +-
 src/MNH/write_les_sv_budgetn.f90   |   6 +-
 src/MNH/write_lesn.f90             | 445 ++++++++++++++---------------
 src/MNH/write_profilern.f90        |   1 +
 src/MNH/write_seriesn.f90          |  37 ++-
 src/MNH/write_stationn.f90         |   1 +
 12 files changed, 361 insertions(+), 275 deletions(-)

diff --git a/src/MNH/modd_budget.f90 b/src/MNH/modd_budget.f90
index d0242d93f..4aed374f8 100644
--- a/src/MNH/modd_budget.f90
+++ b/src/MNH/modd_budget.f90
@@ -111,6 +111,7 @@ type :: tbudiachrometadata
   character(len=NCOMMENTLGTMAX) :: ccomment    = CNOTSET
   character(len=NBUNAMELGTMAX)  :: ccategory   = CNOTSET !budget, LES, aircraft, balloon, series, station, profiler
   character(len=NBUNAMELGTMAX)  :: cshape      = CNOTSET !Shape of the domain (mask, cartesian, vertical profile, point)
+  character(len=NBUNAMELGTMAX)  :: cmask       = CNOTSET !Mask defining where data is of meaning
   character(len=1)              :: cdirection  = ''        !Used for 2pt correlation and spectrum
   logical :: lmobile    = .false.                          !Is the domain moving? (ie for aircrafts and balloons)
   logical :: licompress = .false.
diff --git a/src/MNH/mode_les_diachro.f90 b/src/MNH/mode_les_diachro.f90
index ee778e0ec..5747d2340 100644
--- a/src/MNH/mode_les_diachro.f90
+++ b/src/MNH/mode_les_diachro.f90
@@ -748,9 +748,9 @@ end if
 
 end subroutine Les_diachro_2D
 
-!##################################################################################
-subroutine Les_diachro_3D( tpdiafile, tpfield, odoavg, odonorm, pfield, hsuffixes )
-!##################################################################################
+!##########################################################################################
+subroutine Les_diachro_3D( tpdiafile, tpfield, odoavg, odonorm, pfield, hsuffixes, 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,        &
@@ -763,6 +763,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) :: hsuffixes
+character(len=*),        dimension(:),     optional, intent(in) :: hmasks
 
 type(tfield_metadata_base) :: tzfield
 
@@ -779,9 +780,23 @@ if ( Any( tzfield%ndimlist(4:) /= NMNHDIM_UNUSED ) ) then
   tzfield%ndimlist(4:) = NMNHDIM_UNUSED
 end if
 
-if (         tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL                                                 &
-     .and.   tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME                                                  &
-     .and. ( tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_MASK .or.tzfield%ndimlist(3) == NMNHDIM_BUDGET_TERM ) ) then
+if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
+     .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME  &
+     .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_MASK  ) then
+  if ( .not. Present( hmasks ) ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
+                    'optional dummy argument hmasks is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
+
+  if ( Size( hmasks ) /= Size( pfield, 3) ) &
+    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_3D', 'wrong size for hmasks (' // Trim( tzfield%cmnhname ) // ')' )
+
+  tzfield%ndimlist(4) = NMNHDIM_UNUSED
+  call Les_diachro_common( tpdiafile, tzfield,                                                                &
+                           reshape( pfield, [ size( pfield, 1 ), size( pfield, 2 ), size( pfield, 3 ), 1 ] ), &
+                           odoavg, odonorm, hmasks = hmasks )
+else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
+          .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME  &
+          .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_TERM      ) then
   if ( .not. Present( hsuffixes ) ) &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
                     'optional dummy argument hsuffixes is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
@@ -792,7 +807,7 @@ if (         tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL
   tzfield%ndimlist(4) = NMNHDIM_UNUSED
   call Les_diachro_common( tpdiafile, tzfield,                                                                &
                            reshape( pfield, [ size( pfield, 1 ), size( pfield, 2 ), size( pfield, 3 ), 1 ] ), &
-                           odoavg, odonorm, hsuffixes )
+                           odoavg, odonorm, hsuffixes = hsuffixes )
 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
@@ -800,6 +815,10 @@ else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_3D', &
                     'optional dummy argument hsuffixes 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 ) // ')' )
+
   tzfield%ndimlist(4) = tzfield%ndimlist(3)
   tzfield%ndimlist(3) = NMNHDIM_UNUSED
   call Les_diachro_common( tpdiafile, tzfield,                                                                &
@@ -812,9 +831,9 @@ end if
 
 end subroutine Les_diachro_3D
 
-!##################################################################################
-subroutine Les_diachro_4D( tpdiafile, tpfield, odoavg, odonorm, pfield, hsuffixes )
-!##################################################################################
+!##########################################################################################
+subroutine Les_diachro_4D( tpdiafile, tpfield, odoavg, odonorm, pfield, hsuffixes, 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,        &
@@ -827,6 +846,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) :: hsuffixes
+character(len=*),        dimension(:),     optional, intent(in) :: hmasks
 
 type(tfield_metadata_base) :: tzfield
 
@@ -843,10 +863,22 @@ if ( Any( tzfield%ndimlist(5:) /= NMNHDIM_UNUSED ) ) then
   tzfield%ndimlist(5:) = NMNHDIM_UNUSED
 end if
 
-if (         tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL                                                 &
-     .and.   tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME                                                  &
-     .and. ( tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_MASK .or.tzfield%ndimlist(3) == NMNHDIM_BUDGET_TERM ) &
-     .and.   tzfield%ndimlist(4) == NMNHDIM_BUDGET_LES_SV                                                    ) then
+if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL&
+     .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME &
+     .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_LES_MASK &
+     .and. tzfield%ndimlist(4) == NMNHDIM_BUDGET_LES_SV   ) then
+  if ( .not. Present( hmasks ) ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
+                    'optional dummy argument hmasks is 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 ) // ')' )
+
+  call Les_diachro_common( tpdiafile, tzfield, pfield, odoavg, odonorm, hmasks = hmasks )
+else if (       tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL &
+          .and. tzfield%ndimlist(2) == NMNHDIM_BUDGET_LES_TIME  &
+          .and. tzfield%ndimlist(3) == NMNHDIM_BUDGET_TERM      &
+          .and. tzfield%ndimlist(4) == NMNHDIM_BUDGET_LES_SV    ) then
   if ( .not. Present( hsuffixes ) ) &
     call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
                     'optional dummy argument hsuffixes is needed for tpfield (' // Trim( tzfield%cmnhname ) // ')' )
@@ -854,7 +886,7 @@ if (         tzfield%ndimlist(1) == NMNHDIM_BUDGET_LES_LEVEL
   if ( Size( hsuffixes ) /= Size( pfield, 3) ) &
     call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_4D', 'wrong size for hsuffixes (' // Trim( tzfield%cmnhname ) // ')' )
 
-  call Les_diachro_common( tpdiafile, tzfield, pfield, odoavg, odonorm, hsuffixes )
+  call Les_diachro_common( tpdiafile, tzfield, pfield, odoavg, odonorm, hsuffixes= hsuffixes )
 else
   call Print_msg( NVERB_ERROR, 'IO', 'Les_diachro_4D', &
                   'ndimlist configuration not yet implemented for ' // Trim( tzfield%cmnhname ) )
@@ -862,9 +894,9 @@ end if
 
 end subroutine Les_diachro_4D
 
-!#######################################################################################
-subroutine Les_diachro_common( tpdiafile, tpfield, pfield, odoavg, odonorm, hsuffixes )
-!#######################################################################################
+!##############################################################################################
+subroutine Les_diachro_common( tpdiafile, tpfield, pfield, odoavg, odonorm, hsuffixes, hmasks )
+!##############################################################################################
 
 use modd_field,         only: tfield_metadata_base
 use modd_io,            only: tfiledata
@@ -881,14 +913,14 @@ 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) :: hsuffixes
+character(len=*),           dimension(:),       optional, intent(in) :: hmasks
 
 character(len=100),         dimension(:),     allocatable :: ycomment                      ! Comment string
 character(len=100),         dimension(:),     allocatable :: ytitle                        ! Title
-integer                                                   :: iavg
 integer                                                   :: iles_k                        ! Number of vertical levels
 integer                                                   :: iil, iih, ijl, ijh, ikl, ikh  ! Cartesian area relatively to the
                                                                                            ! entire domain
-integer                                                   :: jk                            ! Vertical loop counter
+integer                                                   :: jp                            ! Process loop counter
 real,                       dimension(:,:,:), allocatable :: ztrajz                        ! x and y are not used for LES
 type(tfield_metadata_base), dimension(:),     allocatable :: tzfields
 !------------------------------------------------------------------------------
@@ -909,9 +941,18 @@ ikl = nles_levels(1)
 ikh = nles_levels(iles_k)
 
 if ( Present( hsuffixes ) ) then
+  if ( Present( hmasks ) ) &
+    call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'hsuffixes and hmasks optional arguments may not be present ' // &
+                    'at the same time (' // Trim( tpfield%cmnhname ) // ')' )
   if ( Size( hsuffixes ) /= Size( pfield, 3) ) &
     call Print_msg( NVERB_FATAL, 'IO', 'Les_diachro_common', 'wrong size for hsuffixes (' // Trim( tpfield%cmnhname ) // ')' )
   ycomment(:) = Trim( tpfield%ccomment(:) ) // hsuffixes(:)
+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 ) // ')' )
+  do jp = 1, Size( ycomment )
+    ycomment(jp) = Trim( tpfield%ccomment(:) ) // ' (' // Trim( hmasks(jp) ) // ')'
+  end do
 else
   ycomment(:) = tpfield%ccomment(:)
 end if
@@ -927,7 +968,9 @@ contains
 subroutine Les_diachro_common_intern( oavg, onorm )
 !#######################################################################################
 
-use modd_field,         only: NMNHDIM_BUDGET_LES_TIME, NMNHDIM_BUDGET_LES_AVG_TIME, NMNHDIM_BUDGET_LES_SV, NMNHDIM_UNUSED
+use modd_budget,        only: CNOTSET
+use modd_field,         only: NMNHDIM_BUDGET_LES_TIME, NMNHDIM_BUDGET_LES_AVG_TIME, NMNHDIM_BUDGET_LES_MASK, &
+                              NMNHDIM_BUDGET_LES_SV, NMNHDIM_UNUSED
 use modd_les,           only: nles_current_times
 
 use mode_write_diachro, only: Write_diachro
@@ -938,6 +981,7 @@ logical, intent(in) :: onorm
 character(len=10)                                    :: ygroup  ! Group title
 integer                                              :: iresp   ! Return code
 integer                                              :: ji
+integer                                              :: jk      ! Vertical loop counter
 integer                                              :: jp      ! Process loop counter
 integer                                              :: jsv     ! Scalar loop counter
 logical                                              :: gsv
@@ -1042,7 +1086,25 @@ if ( iresp == 0 .and. any( zfield /= XUNDEF ) ) then
   tzbudiachro%nkl        = ikl
   tzbudiachro%nkh        = ikh
 
-  call Write_diachro( tpdiafile, tzbudiachro, tzfields, tzdates, zwork6 )
+  if ( tzfields(1)%ndimlist(6) == NMNHDIM_BUDGET_LES_MASK ) then
+    tzfields(:)%ndimlist(6) = NMNHDIM_UNUSED
+
+    ! Loop on the different masks
+    ! Do not provide all tzfields once because they can be stored in different HDF groups (based on masks)
+    do jp = 1, Size( hmasks )
+      tzfields(jp)%clongname = Trim( ytitle(jp) ) // ' (' // Trim( hmasks(jp) ) // ')'
+      tzfields(jp)%ndims     = tzfields(jp)%ndims - 1
+
+      tzbudiachro%cmask      = hmasks(jp)
+
+      call Write_diachro( tpdiafile, tzbudiachro, [ tzfields(jp) ], tzdates, zwork6(:,:,:,:,:,jp:jp) )
+    end do
+  else
+    tzbudiachro%cmask      = 'cart'
+
+    call Write_diachro( tpdiafile, tzbudiachro, tzfields, tzdates, zwork6 )
+  end if
+
 end if
 
 !-------------------------------------------------------------------------------
@@ -1206,6 +1268,7 @@ else
 end if
 tzbudiachro%ccategory  = 'LES'
 tzbudiachro%cshape     = '2-point correlation'
+! tzbudiachro%cmask      = CNOTSET
 if ( tzfield%ndimlist(1) == NMNHDIM_SPECTRA_2PTS_NI ) then
   tzbudiachro%cdirection = 'I'
 else
@@ -1363,6 +1426,7 @@ tzbudiachro%cname      = ygroup
 tzbudiachro%ccomment   = tzfield%ccomment
 tzbudiachro%ccategory  = 'LES'
 tzbudiachro%cshape     = 'spectrum'
+! tzbudiachro%cmask      = CNOTSET
 if ( tzfield%ndimlist(1) == NMNHDIM_SPECTRA_SPEC_NI ) then
   tzbudiachro%cdirection = 'I'
 else
@@ -1400,6 +1464,7 @@ tzbudiachro%cname      = ygroup
 tzbudiachro%ccomment   = Trim( tzfield%ccomment ) // ' (time averaged)'
 tzbudiachro%ccategory  = 'LES'
 tzbudiachro%cshape     = 'spectrum'
+! tzbudiachro%cmask      = CNOTSET
 if ( tzfield%ndimlist(1) == NMNHDIM_SPECTRA_SPEC_NI ) then
   tzbudiachro%cdirection = 'I'
 else
diff --git a/src/MNH/write_aircraft_balloon.f90 b/src/MNH/write_aircraft_balloon.f90
index 15eefcb2c..60f5db69d 100644
--- a/src/MNH/write_aircraft_balloon.f90
+++ b/src/MNH/write_aircraft_balloon.f90
@@ -868,6 +868,7 @@ else
   tzbudiachro%ccategory  = 'unknown'
 end if
 tzbudiachro%cshape     = 'point'
+! tzbudiachro%cmask      = NOT SET (default values)
 tzbudiachro%lmobile    = .true.
 ! tzbudiachro%licompress = NOT SET (default values)
 ! tzbudiachro%ljcompress = NOT SET (default values)
@@ -908,6 +909,7 @@ tzbudiachro%cname      = ygroupz
 tzbudiachro%ccomment   = 'Vertical profiles at position of flyer ' // Trim( tpflyer%title )
 ! tzbudiachro%ccategory  =  !unchanged
 tzbudiachro%cshape     = 'vertical profile'
+! tzbudiachro%cmask      = NOT SET (default values)
 tzbudiachro%lmobile    = .true.
 tzbudiachro%licompress = .true.
 tzbudiachro%ljcompress = .true.
diff --git a/src/MNH/write_budget.f90 b/src/MNH/write_budget.f90
index 1baa0fe7e..99af4b02e 100644
--- a/src/MNH/write_budget.f90
+++ b/src/MNH/write_budget.f90
@@ -46,6 +46,8 @@ private
 
 public :: Write_budget
 
+character(len=*), parameter :: CMASK_VARNAME = 'MASKS'
+
 contains
 
 !#########################################################
@@ -272,7 +274,7 @@ subroutine Write_budget( tpdiafile, tpdtcur, ptstep, ksv )
         tzfile = tpdiafile
         tzfile%cformat = 'NETCDF4'
 
-        tzfield%cmnhname   = 'MASKS'
+        tzfield%cmnhname   = CMASK_VARNAME
         tzfield%cstdname   = ''
         tzfield%clongname  = Trim( tzfield%cmnhname )
         tzfield%cunits     = '1'
@@ -294,7 +296,7 @@ subroutine Write_budget( tpdiafile, tpdtcur, ptstep, ksv )
         !Write the data (partial write of the field with the given offset)
         call IO_Field_write( tzfile, tzfield, nbusurf(:,:,:,:), koffset= [ 0, 0, 0, ( nbutshift - 1 ) * nbusubwrite ] )
 
-        if ( nbutshift == 1 ) call Menu_diachro( tzfile, 'MASKS' )
+        if ( nbutshift == 1 ) call Menu_diachro( tzfile, CMASK_VARNAME )
       end if
   !
   END SELECT
@@ -561,9 +563,11 @@ subroutine Store_one_budget_rho( tpdiafile, tpdates, tprhodj, kp, knocompress, p
   tzbudiachro%ccategory  = 'budget'
   if ( ybutype == 'CART' ) then
     tzbudiachro%cshape   = 'cartesian'
+    ! tzbudiachro%cmask    = NOT SET (default values)
     tzbudiachro%lmobile  = .false.
   else
     tzbudiachro%cshape   = 'mask'
+    tzbudiachro%cmask    = CMASK_VARNAME
     !Masks are updated at each timestep (therefore the studied domains change during execution)
     tzbudiachro%lmobile  = .true.
   end if
@@ -817,9 +821,11 @@ subroutine Store_one_budget( tpdiafile, tpdates, tpbudget, prhodjn, knocompress,
   tzbudiachro%ccategory  = 'budget'
   if ( ybutype == 'CART' ) then
     tzbudiachro%cshape   = 'cartesian'
+    ! tzbudiachro%cmask    = NOT SET (default values)
     tzbudiachro%lmobile  = .false.
   else
     tzbudiachro%cshape   = 'mask'
+    tzbudiachro%cmask    = CMASK_VARNAME
     !Masks are updated at each timestep (therefore the studied domains change during execution)
     tzbudiachro%lmobile  = .true.
   end if
diff --git a/src/MNH/write_diachro.f90 b/src/MNH/write_diachro.f90
index b0e784da5..29a3f92be 100644
--- a/src/MNH/write_diachro.f90
+++ b/src/MNH/write_diachro.f90
@@ -1087,7 +1087,7 @@ select case ( idims )
 
   case (1)
 
-    if ( tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) then
+    if ( Any ( tpfields(1)%ndimlist(4) == [ NMNHDIM_BUDGET_LES_TIME, NMNHDIM_BUDGET_LES_AVG_TIME, NMNHDIM_SERIES_TIME ] ) ) then
       if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                    'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
 
@@ -1200,6 +1200,12 @@ select case ( idims )
       do ji = 1, Size( pvar, 6 )
         call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 4 ], gsplit, gdistributed )
       end do
+    else if (  ( tpfields(1)%ndimlist(3) == NMNHDIM_SERIES_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_SERIES_LEVEL_W ) &
+         .and. tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME ) then
+      !Correspond to WRITE_SERIES_n
+      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
+                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
+      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar(:,:,:,:,:,:), [ 3, 4 ], gsplit, gdistributed )
     else if (  tpfields(1)%ndimlist(4) == NMNHDIM_STATION_TIME &
          .and. tpfields(1)%ndimlist(6) == NMNHDIM_STATION_PROC ) then
       !Correspond to WRITE_STATION_n
diff --git a/src/MNH/write_les_budgetn.f90 b/src/MNH/write_les_budgetn.f90
index da600ed63..2897bd97f 100644
--- a/src/MNH/write_les_budgetn.f90
+++ b/src/MNH/write_les_budgetn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2000-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2000-2021 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.
@@ -436,7 +436,7 @@ tzfield%clongname = ygroup
 tzfield%ccomment  = 'resolved KE budget'
 tzfield%cunits    = 'm2 s-3'
 
-call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -688,7 +688,7 @@ tzfield%clongname = ygroup
 tzfield%ccomment  = 'thetal variance budget'
 tzfield%cunits    = 'K2 s-1'
 
-call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -1003,7 +1003,7 @@ tzfield%clongname = ygroup
 tzfield%ccomment  = 'thetal flux budget'
 tzfield%cunits    = 'm K s-2'
 
-call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/write_les_rt_budgetn.f90 b/src/MNH/write_les_rt_budgetn.f90
index 595556a56..3bda050f1 100644
--- a/src/MNH/write_les_rt_budgetn.f90
+++ b/src/MNH/write_les_rt_budgetn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2002-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-2021 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.
@@ -363,7 +363,7 @@ tzfield%clongname = ygroup
 tzfield%ccomment  = 'Rt variance budget'
 tzfield%cunits    = 'kg2 kg-2 s-1'
 
-call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -673,7 +673,7 @@ tzfield%clongname = ygroup
 tzfield%ccomment  = 'Rt flux budget'
 tzfield%cunits    = 'm kg kg-1 s-2'
 
-call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -923,7 +923,7 @@ tzfield%clongname = ygroup
 tzfield%ccomment  = 'Thl-Rt covariance budget'
 tzfield%cunits    = 'K kg kg-1 s-1'
 
-call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/write_les_sv_budgetn.f90 b/src/MNH/write_les_sv_budgetn.f90
index e2434073f..377a5ea84 100644
--- a/src/MNH/write_les_sv_budgetn.f90
+++ b/src/MNH/write_les_sv_budgetn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2002-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-2021 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.
@@ -399,7 +399,7 @@ 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, gdoavg, gdonorm, zles_budget(:, :, :iles, :), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles, :), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
@@ -781,7 +781,7 @@ 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, gdoavg, gdonorm, zles_budget(:, :, :iles, :), ysubtitle(:iles) )
+call Les_diachro( tpdiafile, tzfield, gdoavg, gdonorm, zles_budget(:, :, :iles, :), hsuffixes = ysubtitle(:iles) )
 
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/write_lesn.f90 b/src/MNH/write_lesn.f90
index 09724ca6a..b2a3891f6 100644
--- a/src/MNH/write_lesn.f90
+++ b/src/MNH/write_lesn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2000-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2000-2021 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.
@@ -107,8 +107,9 @@ INTEGER :: IMASK
 INTEGER :: JSV       ! scalar loop counter
 INTEGER :: JI        ! loop counter
 !
-CHARACTER(len=9), DIMENSION(NLES_MASKS) :: YSUBTITLE
+character(len=3)                        :: ynum
 CHARACTER(len=5)                        :: YGROUP
+character(len=7), dimension(nles_masks) :: ymasks
 !
 logical :: gdoavg    ! Compute and store time average
 logical :: gdonorm   ! Compute and store normalized field
@@ -234,32 +235,33 @@ tfield%ntype = TYPEREAL
 !*      2.   (z,t) profiles (all masks)
 !            --------------
 IMASK = 1
-YSUBTITLE(IMASK) = " (cart)"
+ymasks(imask) = 'cart'
 IF (LLES_NEB_MASK) THEN
   IMASK=IMASK+1
-  YSUBTITLE(IMASK) = " (neb)"
+  ymasks(imask) = 'neb'
   IMASK=IMASK+1
-  YSUBTITLE(IMASK) = " (clear)"
+  ymasks(imask) = 'clear'
 END IF
 IF (LLES_CORE_MASK) THEN
   IMASK=IMASK+1
-  YSUBTITLE(IMASK) = " (core)"
+  ymasks(imask) = 'core'
   IMASK=IMASK+1
-  YSUBTITLE(IMASK) = " (env)"
+  ymasks(imask) = 'env'
 END IF
 IF (LLES_MY_MASK) THEN
    DO JI=1,NLES_MASKS_USER
         IMASK=IMASK+1
-        YSUBTITLE(IMASK) = " (user)"
+        Write( ynum, '( i3.3 )' ) ji
+        ymasks(imask) = 'user' // ynum
    END DO
 END IF
 IF (LLES_CS_MASK) THEN
   IMASK=IMASK+1
-  YSUBTITLE(IMASK) = " (cs1)"
+  ymasks(imask) = 'cs1'
   IMASK=IMASK+1
-  YSUBTITLE(IMASK) = " (cs2)"
+  ymasks(imask) = 'cs2'
   IMASK=IMASK+1
-  YSUBTITLE(IMASK) = " (cs3)"
+  ymasks(imask) = 'cs3'
 END IF
 !
 !*      2.0  averaging diagnostics
@@ -281,10 +283,10 @@ tfield%ndimlist(4:) = NMNHDIM_UNUSED
 ldoavg  = xles_temp_mean_start /= XUNDEF .and. xles_temp_mean_end /= XUNDEF
 ldonorm = .false.
 
-call Les_diachro_write( tpdiafile, zavg_pts_ll,                'AVG_PTS',  'number of points used for averaging',   '1', ysubtitle)
-call Les_diachro_write( tpdiafile, zavg_pts_ll / zcart_pts_ll, 'AVG_PTSF', 'fraction of points used for averaging', '1', ysubtitle)
-call Les_diachro_write( tpdiafile, zund_pts_ll,                'UND_PTS',  'number of points below orography',      '1', ysubtitle)
-call Les_diachro_write( tpdiafile, zund_pts_ll / zcart_pts_ll, 'UND_PTSF', 'fraction of points below orography',    '1', ysubtitle)
+call Les_diachro_write( tpdiafile, zavg_pts_ll,                'AVG_PTS',  'number of points used for averaging',   '1', ymasks )
+call Les_diachro_write( tpdiafile, zavg_pts_ll / zcart_pts_ll, 'AVG_PTSF', 'fraction of points used for averaging', '1', ymasks )
+call Les_diachro_write( tpdiafile, zund_pts_ll,                'UND_PTS',  'number of points below orography',      '1', ymasks )
+call Les_diachro_write( tpdiafile, zund_pts_ll / zcart_pts_ll, 'UND_PTSF', 'fraction of points below orography',    '1', ymasks )
 
 DEALLOCATE(ZAVG_PTS_ll)
 DEALLOCATE(ZUND_PTS_ll)
@@ -301,54 +303,54 @@ tfield%ndimlist(4:) = NMNHDIM_UNUSED
 ldoavg  = xles_temp_mean_start /= XUNDEF .and. xles_temp_mean_end /= XUNDEF
 ldonorm = trim(cles_norm_type) /= 'NONE'
 
-call Les_diachro_write( tpdiafile, XLES_MEAN_U,      'MEAN_U',      'Mean U Profile',                        'm s-1',  ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_V,      'MEAN_V',      'Mean V Profile',                        'm s-1',  ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_W,      'MEAN_W',      'Mean W Profile',                        'm s-1',  ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_P,      'MEAN_PRE',    'Mean pressure Profile',                 'Pa',     ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_DP,     'MEAN_DP',     'Mean Dyn production TKE Profile',       'm2 s-3', ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_TP,     'MEAN_TP',     'Mean Thermal production TKE Profile',   'm2 s-3', ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_TR,     'MEAN_TR',     'Mean transport production TKE Profile', 'm2 s-3', ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_DISS,   'MEAN_DISS',   'Mean Dissipation TKE Profile',          'm2 s-3', ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_LM,     'MEAN_LM',     'Mean mixing length Profile',            'm',      ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_RHO,    'MEAN_RHO',    'Mean density Profile',                  'kg m-3', ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_Th,     'MEAN_TH',     'Mean potential temperature Profile',    'K',      ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_MEAN_Mf,     'MEAN_MF',     'Mass-flux Profile',                     'm s-1',  ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_U,      'MEAN_U',      'Mean U Profile',                        'm s-1',  ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_V,      'MEAN_V',      'Mean V Profile',                        'm s-1',  ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_W,      'MEAN_W',      'Mean W Profile',                        'm s-1',  ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_P,      'MEAN_PRE',    'Mean pressure Profile',                 'Pa',     ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_DP,     'MEAN_DP',     'Mean Dyn production TKE Profile',       'm2 s-3', ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_TP,     'MEAN_TP',     'Mean Thermal production TKE Profile',   'm2 s-3', ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_TR,     'MEAN_TR',     'Mean transport production TKE Profile', 'm2 s-3', ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_DISS,   'MEAN_DISS',   'Mean Dissipation TKE Profile',          'm2 s-3', ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_LM,     'MEAN_LM',     'Mean mixing length Profile',            'm',      ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_RHO,    'MEAN_RHO',    'Mean density Profile',                  'kg m-3', ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Th,     'MEAN_TH',     'Mean potential temperature Profile',    'K',      ymasks )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Mf,     'MEAN_MF',     'Mass-flux Profile',                     'm s-1',  ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Thl,    'MEAN_THL',    'Mean liquid potential temperature Profile',  'K', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Thl,    'MEAN_THL',    'Mean liquid potential temperature Profile',  'K', ymasks )
 if ( luserv ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Thv,    'MEAN_THV',    'Mean virtual potential temperature Profile', 'K', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Thv,    'MEAN_THV',    'Mean virtual potential temperature Profile', 'K', ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rt,     'MEAN_RT',     'Mean Rt Profile', 'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rt,     'MEAN_RT',     'Mean Rt Profile', 'kg kg-1', ymasks )
 if ( luserv ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rv,     'MEAN_RV',     'Mean Rv Profile', 'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rv,     'MEAN_RV',     'Mean Rv Profile', 'kg kg-1', ymasks )
 if ( luserv ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rehu,   'MEAN_REHU',   'Mean Rh Profile', 'percent', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rehu,   'MEAN_REHU',   'Mean Rh Profile', 'percent', ymasks )
 if ( luserv ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Qs,     'MEAN_QS',     'Mean Qs Profile', 'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Qs,     'MEAN_QS',     'Mean Qs Profile', 'kg kg-1', ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_KHt,    'MEAN_KHT',    'Eddy-diffusivity (temperature) Profile', 'm2 s-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_KHt,    'MEAN_KHT',    'Eddy-diffusivity (temperature) Profile', 'm2 s-1', ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_KHr,    'MEAN_KHR',    'Eddy-diffusivity (wvapor) Profile',      'm2 s-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_KHr,    'MEAN_KHR',    'Eddy-diffusivity (wvapor) Profile',      'm2 s-1', ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rc,     'MEAN_RC',     'Mean Rc Profile', 'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rc,     'MEAN_RC',     'Mean Rc Profile',              'kg kg-1', ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Cf,     'MEAN_CF',     'Mean Cf Profile',              '1',       ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Cf,     'MEAN_CF',     'Mean Cf Profile',              '1',       ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_INDCf,  'MEAN_INDCF',  'Mean Cf>1-6 Profile (0 or 1)', '1',       ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_INDCf,  'MEAN_INDCF',  'Mean Cf>1-6 Profile (0 or 1)', '1',       ymasks )
 if ( luserc ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_INDCf2, 'MEAN_INDCF2', 'Mean Cf>1-5 Profile (0 or 1)', '1',       ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_INDCf2, 'MEAN_INDCF2', 'Mean Cf>1-5 Profile (0 or 1)', '1',       ymasks )
 if ( luserr ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rr,     'MEAN_RR',     'Mean Rr Profile',              'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rr,     'MEAN_RR',     'Mean Rr Profile',              'kg kg-1', ymasks )
 if ( luserr ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_RF,     'MEAN_RF',     'Mean RF Profile',              '1',       ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_RF,     'MEAN_RF',     'Mean RF Profile',              '1',       ymasks )
 if ( luseri ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Ri,     'MEAN_RI',     'Mean Ri Profile',              'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Ri,     'MEAN_RI',     'Mean Ri Profile',              'kg kg-1', ymasks )
 if ( lusers ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rs,     'MEAN_RS',     'Mean Rs Profile',              'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rs,     'MEAN_RS',     'Mean Rs Profile',              'kg kg-1', ymasks )
 if ( luserg ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rg,     'MEAN_RG',     'Mean Rg Profile',              'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rg,     'MEAN_RG',     'Mean Rg Profile',              'kg kg-1', ymasks )
 if ( luserh ) &
-call Les_diachro_write( tpdiafile, XLES_MEAN_Rh,     'MEAN_RH',     'Mean Rh Profile',              'kg kg-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_Rh,     'MEAN_RH',     'Mean Rh Profile',              'kg kg-1', ymasks )
 
 if ( nsv > 0 ) then
   tfield%ndims = 4
@@ -358,7 +360,7 @@ if ( nsv > 0 ) then
   tfield%ndimlist(4)  = NMNHDIM_BUDGET_LES_SV
   tfield%ndimlist(5:) = NMNHDIM_UNUSED
 
-  call Les_diachro_write( tpdiafile, XLES_MEAN_Sv, 'MEAN_SV', 'Mean Sv Profiles', 'kg kg-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_MEAN_Sv, 'MEAN_SV', 'Mean Sv Profiles', 'kg kg-1', ymasks )
 
   tfield%ndims = 3
   !tfield%ndimlist(1)  = NMNHDIM_BUDGET_LES_LEVEL
@@ -368,28 +370,28 @@ if ( nsv > 0 ) then
   !tfield%ndimlist(5:) = NMNHDIM_UNUSED
 end if
 
-call Les_diachro_write( tpdiafile, XLES_MEAN_WIND, 'MEANWIND',       'Profile of Mean Modulus of Wind', 'm s-1',      ysubtitle )
-call Les_diachro_write( tpdiafile, XLES_RESOLVED_MASSFX, 'MEANMSFX', 'Total updraft mass flux',         'kg m-2 s-1', ysubtitle )
+call Les_diachro_write( tpdiafile, XLES_MEAN_WIND, 'MEANWIND',       'Profile of Mean Modulus of Wind', 'm s-1',      ymasks )
+call Les_diachro_write( tpdiafile, XLES_RESOLVED_MASSFX, 'MEANMSFX', 'Total updraft mass flux',         'kg m-2 s-1', ymasks )
 
 if ( lles_pdf ) then
-  call Les_diachro_write( tpdiafile,   XLES_PDF_TH,  'PDF_TH',  'Pdf potential temperature Profiles', '1', ysubtitle )
-  call Les_diachro_write( tpdiafile,   XLES_PDF_W,   'PDF_W',   'Pdf vertical velocity Profiles',     '1', ysubtitle )
-  call Les_diachro_write( tpdiafile,   XLES_PDF_THV, 'PDF_THV', 'Pdf virtual pot. temp. Profiles',    '1', ysubtitle )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_TH,  'PDF_TH',  'Pdf potential temperature Profiles', '1', ymasks )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_W,   'PDF_W',   'Pdf vertical velocity Profiles',     '1', ymasks )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_THV, 'PDF_THV', 'Pdf virtual pot. temp. Profiles',    '1', ymasks )
   if ( luserv ) &
-  call Les_diachro_write( tpdiafile,   XLES_PDF_RV,  'PDF_RV',  'Pdf Rv Profiles',                    '1', ysubtitle )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_RV,  'PDF_RV',  'Pdf Rv Profiles',                    '1', ymasks )
   if ( luserc ) then
-    call Les_diachro_write( tpdiafile, XLES_PDF_RC,  'PDF_RC',  'Pdf Rc Profiles',                    '1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_PDF_RT,  'PDF_RT',  'Pdf Rt Profiles',                    '1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_PDF_THL, 'PDF_THL', 'Pdf Thl Profiles',                   '1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_PDF_RC,  'PDF_RC',  'Pdf Rc Profiles',                    '1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_PDF_RT,  'PDF_RT',  'Pdf Rt Profiles',                    '1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_PDF_THL, 'PDF_THL', 'Pdf Thl Profiles',                   '1', ymasks )
   end if
   if ( luserr ) &
-  call Les_diachro_write( tpdiafile,   XLES_PDF_RR,  'PDF_RR',  'Pdf Rr Profiles',                    '1', ysubtitle )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_RR,  'PDF_RR',  'Pdf Rr Profiles',                    '1', ymasks )
   if ( luseri ) &
-  call Les_diachro_write( tpdiafile,   XLES_PDF_RI,  'PDF_RI',  'Pdf Ri Profiles',                    '1', ysubtitle )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_RI,  'PDF_RI',  'Pdf Ri Profiles',                    '1', ymasks )
   if ( lusers ) &
-  call Les_diachro_write( tpdiafile,   XLES_PDF_RS,  'PDF_RS',  'Pdf Rs Profiles',                    '1', ysubtitle )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_RS,  'PDF_RS',  'Pdf Rs Profiles',                    '1', ymasks )
   if ( luserg ) &
-  call Les_diachro_write( tpdiafile,   XLES_PDF_RG,  'PDF_RG',  'Pdf Rg Profiles',                    '1', ysubtitle )
+  call Les_diachro_write( tpdiafile,   XLES_PDF_RG,  'PDF_RG',  'Pdf Rg Profiles',                    '1', ymasks )
 end if
 !
 !*      2.2  resolved quantities
@@ -406,89 +408,79 @@ if ( lles_resolved ) then
   tfield%ndimlist(3)  = NMNHDIM_BUDGET_LES_MASK
   tfield%ndimlist(4:) = NMNHDIM_UNUSED
 
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_U2, 'RES_U2',  'Resolved <u2> variance',        'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_V2, 'RES_V2',  'Resolved <v2> variance',        'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2, 'RES_W2',  'Resolved <w2> variance',        'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UV, 'RES_UV',  'Resolved <uv> Flux',            'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WU, 'RES_WU',  'Resolved <wu> Flux',            'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WV, 'RES_WV',  'Resolved <wv> Flux',            'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_Ke, 'RES_KE',  'Resolved TKE Profile',          'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_P2, 'RES_P2',  'Resolved pressure variance',    'Pa2',    ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UP, 'RES_UPZ', 'Resolved <up> horizontal Flux', 'Pa s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_VP, 'RES_VPZ', 'Resolved <vp> horizontal Flux', 'Pa s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WP, 'RES_WPZ', 'Resolved <wp> vertical Flux',   'Pa s-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_U2, 'RES_U2',  'Resolved <u2> variance',        'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_V2, 'RES_V2',  'Resolved <v2> variance',        'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2, 'RES_W2',  'Resolved <w2> variance',        'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UV, 'RES_UV',  'Resolved <uv> Flux',            'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WU, 'RES_WU',  'Resolved <wu> Flux',            'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WV, 'RES_WV',  'Resolved <wv> Flux',            'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_Ke, 'RES_KE',  'Resolved TKE Profile',          'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_P2, 'RES_P2',  'Resolved pressure variance',    'Pa2',    ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UP, 'RES_UPZ', 'Resolved <up> horizontal Flux', 'Pa s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_VP, 'RES_VPZ', 'Resolved <vp> horizontal Flux', 'Pa s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WP, 'RES_WPZ', 'Resolved <wp> vertical Flux',   'Pa s-1', ymasks )
 
   if ( luserv ) &
   call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThThv, 'RES_THTV', &
-                          'Resolved potential temperature - virtual potential temperature covariance',        'K2', ysubtitle )
+                          'Resolved potential temperature - virtual potential temperature covariance',        'K2', ymasks )
   if ( luserc ) &
   call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlThv, 'RES_TLTV', &
-                          'Resolved liquid potential temperature - virtual potential temperature covariance', 'K2',  ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_Th2, 'RES_TH2', 'Resolved potential temperature variance', 'K2', ysubtitle )
+                          'Resolved liquid potential temperature - virtual potential temperature covariance', 'K2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_Th2, 'RES_TH2', 'Resolved potential temperature variance', 'K2', ymasks )
   if ( luserc ) &
   call Les_diachro_write( tpdiafile, XLES_RESOLVED_Thl2, 'RES_THL2', 'Resolved liquid potential temperature variance', 'K2',&
-                          ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UTh, 'RES_UTH', 'Resolved <uth> horizontal Flux', 'm K s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_VTh, 'RES_VTH', 'Resolved <vth> horizontal Flux', 'm K s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WTh, 'RES_WTH', 'Resolved <wth> vertical Flux',   'm K s-1', ysubtitle )
+                          ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UTh, 'RES_UTH', 'Resolved <uth> horizontal Flux', 'm K s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_VTh, 'RES_VTH', 'Resolved <vth> horizontal Flux', 'm K s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WTh, 'RES_WTH', 'Resolved <wth> vertical Flux',   'm K s-1', ymasks )
 
   if ( luserc ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_UThl, 'RES_UTHL', 'Resolved <uthl> horizontal Flux', 'm K s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VThl, 'RES_VTHL', 'Resolved <vthl> horizontal Flux', 'm K s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThl, 'RES_WTHL', 'Resolved <wthl> vertical Flux',   'm K s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Rt2, 'RES_RT2',  'Resolved total water variance',    'kg2 kg-2', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRt, 'RES_WRT',  'Resolved <wrt> vertical Flux', 'm kg kg-1 s-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_UThl, 'RES_UTHL', 'Resolved <uthl> horizontal Flux', 'm K s-1',       ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VThl, 'RES_VTHL', 'Resolved <vthl> horizontal Flux', 'm K s-1',       ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThl, 'RES_WTHL', 'Resolved <wthl> vertical Flux',   'm K s-1',       ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Rt2,  'RES_RT2',  'Resolved total water variance',   'kg2 kg-2',      ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRt,  'RES_WRT',  'Resolved <wrt> vertical Flux',    'm kg kg-1 s-1', ymasks )
   end if
 
   if ( luserv ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_UThv,  'RES_UTHV', 'Resolved <uthv> horizontal Flux', 'm K s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VThv,  'RES_VTHV', 'Resolved <vthv> horizontal Flux', 'm K s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThv,  'RES_WTHV', 'Resolved <wthv> vertical Flux',   'm K s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Rv2,   'RES_RV2',  'Resolved water vapor variance',   'kg2 kg-2', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThRv,  'RES_THRV', 'Resolved <thrv> covariance',      'K kg kg-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_UThv,  'RES_UTHV', 'Resolved <uthv> horizontal Flux', 'm K s-1',       ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VThv,  'RES_VTHV', 'Resolved <vthv> horizontal Flux', 'm K s-1',       ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThv,  'RES_WTHV', 'Resolved <wthv> vertical Flux',   'm K s-1',       ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Rv2,   'RES_RV2',  'Resolved water vapor variance',   'kg2 kg-2',      ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThRv,  'RES_THRV', 'Resolved <thrv> covariance',      'K kg kg-1',     ymasks )
     if ( luserc ) &
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlRv, 'RES_TLRV', 'Resolved <thlrv> covariance',     'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvRv, 'RES_TVRV', 'Resolved <thvrv> covariance',     'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_URv,   'RES_URV',  'Resolved <urv> horizontal flux',  'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VRv,   'RES_VRV',  'Resolved <vrv> horizontal flux',  'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRv,   'RES_WRV',  'Resolved <wrv> vertical flux',    'm kg kg-1 s-1', &
-                            ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlRv, 'RES_TLRV', 'Resolved <thlrv> covariance',     'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvRv, 'RES_TVRV', 'Resolved <thvrv> covariance',     'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_URv,   'RES_URV',  'Resolved <urv> horizontal flux',  'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VRv,   'RES_VRV',  'Resolved <vrv> horizontal flux',  'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRv,   'RES_WRV',  'Resolved <wrv> vertical flux',    'm kg kg-1 s-1', ymasks )
   end if
 
   if ( luserc ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Rc2,   'RES_RC2',  'Resolved cloud water variance',   'kg2 kg-2', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThRc,  'RES_THRC', 'Resolved <thrc> covariance',      'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlRc, 'RES_TLRC', 'Resolved <thlrc> covariance',     'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvRc, 'RES_TVRC', 'Resolved <thvrc> covariance',     'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_URc,   'RES_URC',  'Resolved <urc> horizontal flux',  'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VRc,   'RES_VRC',  'Resolved <vrc> horizontal flux',  'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRc,   'RES_WRC',  'Resolved <wrc> vertical flux',    'm kg kg-1 s-1', &
-                            ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Rc2,   'RES_RC2',  'Resolved cloud water variance',  'kg2 kg-2',      ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThRc,  'RES_THRC', 'Resolved <thrc> covariance',     'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlRc, 'RES_TLRC', 'Resolved <thlrc> covariance',    'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvRc, 'RES_TVRC', 'Resolved <thvrc> covariance',    'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_URc,   'RES_URC',  'Resolved <urc> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VRc,   'RES_VRC',  'Resolved <vrc> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRc,   'RES_WRC',  'Resolved <wrc> vertical flux',   'm kg kg-1 s-1', ymasks )
   end if
 
   if ( luseri ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Ri2,   'RES_RI2',  'Resolved cloud ice variance',     'kg2 kg-2', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThRi,  'RES_THRI', 'Resolved <thri> covariance',      'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlRi, 'RES_TLRI', 'Resolved <thlri> covariance',     'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvRi, 'RES_TVRI', 'Resolved <thvri> covariance',     'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_URi,   'RES_URI',  'Resolved <uri> horizontal flux',  'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VRi,   'RES_VRI',  'Resolved <vri> horizontal flux',  'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRi,   'RES_WRI',  'Resolved <wri> vertical flux',    'm kg kg-1 s-1', &
-                            ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Ri2,   'RES_RI2',  'Resolved cloud ice variance',    'kg2 kg-2',      ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThRi,  'RES_THRI', 'Resolved <thri> covariance',     'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlRi, 'RES_TLRI', 'Resolved <thlri> covariance',    'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvRi, 'RES_TVRI', 'Resolved <thvri> covariance',    'K kg kg-1',     ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_URi,   'RES_URI',  'Resolved <uri> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VRi,   'RES_VRI',  'Resolved <vri> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRi,   'RES_WRI',  'Resolved <wri> vertical flux',   'm kg kg-1 s-1', ymasks )
   end if
 
   if ( luserr ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRr,   'RES_WRR',  'Resolved <wrr> vertical flux',    'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_INPRR3D,     'INPRR3D',   'Precipitation flux',  'm s-1',       ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_MAX_INPRR3D, 'MAXINPR3D', 'Max Precip flux',     'm s-1',       ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_EVAP3D,      'EVAP3D',    'Evaporation profile', 'kg kg-1 s-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRr,   'RES_WRR',   'Resolved <wrr> vertical flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_INPRR3D,        'INPRR3D',   'Precipitation flux',           'm s-1',         ymasks )
+    call Les_diachro_write( tpdiafile, XLES_MAX_INPRR3D,    'MAXINPR3D', 'Max Precip flux',              'm s-1',         ymasks )
+    call Les_diachro_write( tpdiafile, XLES_EVAP3D,         'EVAP3D',    'Evaporation profile',          'kg kg-1 s-1',   ymasks )
   end if
 
   if ( nsv > 0 ) then
@@ -499,19 +491,15 @@ if ( lles_resolved ) then
     tfield%ndimlist(4)  = NMNHDIM_BUDGET_LES_SV
     tfield%ndimlist(5:) = NMNHDIM_UNUSED
 
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Sv2,   'RES_SV2',  'Resolved scalar variables variances', 'kg2 kg-2', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThSv,  'RES_THSV', 'Resolved <ThSv> variance',  'K kg kg-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_Sv2,   'RES_SV2',  'Resolved scalar variables variances', 'kg2 kg-2', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThSv,  'RES_THSV', 'Resolved <ThSv> variance',  'K kg kg-1',          ymasks )
     if ( luserc ) &
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlSv, 'RES_TLSV', 'Resolved <ThlSv> variance', 'K kg kg-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlSv, 'RES_TLSV', 'Resolved <ThlSv> variance', 'K kg kg-1',          ymasks )
     if ( luserv ) &
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvSv, 'RES_TVSV', 'Resolved <ThvSv> variance', 'K kg kg-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_USv,   'RES_USV',  'Resolved <uSv> horizontal flux', 'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VSv,   'RES_VSV',  'Resolved <vSv> horizontal flux', 'm kg kg-1 s-1', &
-                            ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WSv,   'RES_WSV',  'Resolved <wSv> vertical flux', 'm kg kg-1 s-1',   &
-                            ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThvSv, 'RES_TVSV', 'Resolved <ThvSv> variance', 'K kg kg-1',          ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_USv,   'RES_USV',  'Resolved <uSv> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_VSv,   'RES_VSV',  'Resolved <vSv> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WSv,   'RES_WSV',  'Resolved <wSv> vertical flux', 'm kg kg-1 s-1',   ymasks )
 
     tfield%ndims = 3
     !tfield%ndimlist(1)  = NMNHDIM_BUDGET_LES_LEVEL
@@ -521,37 +509,37 @@ if ( lles_resolved ) then
     !tfield%ndimlist(5:) = NMNHDIM_UNUSED
   end if
 
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_U3, 'RES_U3', 'Resolved <u3>', 'm3 s-3', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_V3, 'RES_V3', 'Resolved <v3>', 'm3 s-3', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W3, 'RES_W3', 'Resolved <w3>', 'm3 s-3', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_U4, 'RES_U4', 'Resolved <u4>', 'm4 s-4', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_V4, 'RES_V4', 'Resolved <v4>', 'm4 s-4', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W4, 'RES_W4', 'Resolved <>w4', 'm4 s-4', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_U3, 'RES_U3', 'Resolved <u3>', 'm3 s-3', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_V3, 'RES_V3', 'Resolved <v3>', 'm3 s-3', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W3, 'RES_W3', 'Resolved <w3>', 'm3 s-3', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_U4, 'RES_U4', 'Resolved <u4>', 'm4 s-4', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_V4, 'RES_V4', 'Resolved <v4>', 'm4 s-4', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W4, 'RES_W4', 'Resolved <>w4', 'm4 s-4', ymasks )
 
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThl2, 'RES_WTL2', 'Resolved <wThl2',  'm K2 s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Thl, 'RES_W2TL', 'Resolved <w2Thl>', 'm2 K s-2', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThl2, 'RES_WTL2', 'Resolved <wThl2',  'm K2 s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Thl, 'RES_W2TL', 'Resolved <w2Thl>', 'm2 K s-2', ymasks )
 
   if ( luserv ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRv2,   'RES_WRV2', 'Resolved <wRv2>',   'm kg2 kg-2 s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Rv,   'RES_W2RV', 'Resolved <w2Rv>',   'm2 kg kg-1 s-2',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRt2,   'RES_WRT2', 'Resolved <wRt2>',   'm kg2 kg-2 s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Rt,   'RES_W2RT', 'Resolved <w2Rt>',   'm2 kg kg-1 s-2',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRv, 'RE_WTLRV', 'Resolved <wThlRv>', 'm K kg kg-1 s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRt, 'RE_WTLRT', 'Resolved <wThlRt>', 'm K kg kg-1 s-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRv2,   'RES_WRV2', 'Resolved <wRv2>',   'm kg2 kg-2 s-1',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Rv,   'RES_W2RV', 'Resolved <w2Rv>',   'm2 kg kg-1 s-2',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRt2,   'RES_WRT2', 'Resolved <wRt2>',   'm kg2 kg-2 s-1',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Rt,   'RES_W2RT', 'Resolved <w2Rt>',   'm2 kg kg-1 s-2',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRv, 'RE_WTLRV', 'Resolved <wThlRv>', 'm K kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRt, 'RE_WTLRT', 'Resolved <wThlRt>', 'm K kg kg-1 s-1', ymasks )
   end if
 
   if ( luserc ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRc2,   'RES_WRC2', 'Resolved <wRc2>',   'm kg2 kg-2 s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Rc,   'RES_W2RC', 'Resolved <w2Rc>',   'm2 kg kg-1 s-2',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRc, 'RE_WTLRC', 'Resolved <wThlRc>', 'm K kg kg-1 s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRvRc,  'RE_WRVRC', 'Resolved <wRvRc>',  'm kg2 kg-2 s-1',  ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRc2,   'RES_WRC2', 'Resolved <wRc2>',   'm kg2 kg-2 s-1',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Rc,   'RES_W2RC', 'Resolved <w2Rc>',   'm2 kg kg-1 s-2',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRc, 'RE_WTLRC', 'Resolved <wThlRc>', 'm K kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRvRc,  'RE_WRVRC', 'Resolved <wRvRc>',  'm kg2 kg-2 s-1',  ymasks )
   end if
 
   if ( luseri ) then
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRi2,   'RES_WRI2', 'Resolved <wRi2>',   'm kg2 kg-2 s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Ri,   'RES_W2RI', 'Resolved <w2Ri>',   'm2 kg kg-1 s-2',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRi, 'RE_WTLRI', 'Resolved <wThlRi>', 'm K kg kg-1 s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRvRi,  'RE_WRVRI', 'Resolved <wRvRi>',  'm kg2 kg-2 s-1',  ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRi2,   'RES_WRI2', 'Resolved <wRi2>',   'm kg2 kg-2 s-1',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Ri,   'RES_W2RI', 'Resolved <w2Ri>',   'm2 kg kg-1 s-2',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlRi, 'RE_WTLRI', 'Resolved <wThlRi>', 'm K kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRvRi,  'RE_WRVRI', 'Resolved <wRvRi>',  'm kg2 kg-2 s-1',  ymasks )
   end if
 
   if ( nsv > 0 ) then
@@ -562,11 +550,11 @@ if ( lles_resolved ) then
     tfield%ndimlist(4)  = NMNHDIM_BUDGET_LES_SV
     tfield%ndimlist(5:) = NMNHDIM_UNUSED
 
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WSv2,   'RES_WSV2', 'Resolved <wSv2>',   'm kg2 kg-2 s-1',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Sv,   'RES_W2SV', 'Resolved <w2Sv>',   'm2 kg kg-1 s-2',  ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlSv, 'RE_WTLSV', 'Resolved <wThlSv>', 'm K kg kg-1 s-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WSv2,   'RES_WSV2', 'Resolved <wSv2>',   'm kg2 kg-2 s-1',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_W2Sv,   'RES_W2SV', 'Resolved <w2Sv>',   'm2 kg kg-1 s-2',  ymasks )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WThlSv, 'RE_WTLSV', 'Resolved <wThlSv>', 'm K kg kg-1 s-1', ymasks )
     if ( luserv ) &
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRvSv,  'RE_WRVSV', 'Resolved <wRvSv>',  'm kg2 kg-2 s-1',  ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_WRvSv,  'RE_WRVSV', 'Resolved <wRvSv>',  'm kg2 kg-2 s-1',  ymasks )
 
     tfield%ndims = 3
     !tfield%ndimlist(1)  = NMNHDIM_BUDGET_LES_LEVEL
@@ -576,15 +564,15 @@ if ( lles_resolved ) then
     !tfield%ndimlist(5:) = NMNHDIM_UNUSED
   end if
 
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlPz, 'RES_TLPZ', 'Resolved <Thldp/dz>', 'K Pa m-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_ThlPz, 'RES_TLPZ', 'Resolved <Thldp/dz>', 'K Pa m-1',        ymasks )
   if ( luserc ) &
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RtPz, 'RES_RTPZ', 'Resolved <Rtdp/dz>', 'kg2 kg-2 Pa m-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RtPz,  'RES_RTPZ', 'Resolved <Rtdp/dz>',  'kg2 kg-2 Pa m-1', ymasks )
   if ( luserv ) &
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RvPz, 'RES_RVPZ', 'Resolved <Rvdp/dz>', 'kg2 kg-2 Pa m-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RvPz,  'RES_RVPZ', 'Resolved <Rvdp/dz>',  'kg2 kg-2 Pa m-1', ymasks )
   if ( luserc ) &
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RcPz, 'RES_RCPZ', 'Resolved <Rcdp/dz>', 'kg2 kg-2 Pa m-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RcPz,  'RES_RCPZ', 'Resolved <Rcdp/dz>',  'kg2 kg-2 Pa m-1', ymasks )
   if ( luseri ) &
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RiPz, 'RES_RIPZ', 'Resolved <Ridp/dz>', 'kg2 kg-2 Pa m-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_RiPz,  'RES_RIPZ', 'Resolved <Ridp/dz>',  'kg2 kg-2 Pa m-1', ymasks )
 
   if ( nsv > 0 ) then
     tfield%ndims = 4
@@ -594,7 +582,7 @@ if ( lles_resolved ) then
     tfield%ndimlist(4)  = NMNHDIM_BUDGET_LES_SV
     tfield%ndimlist(5:) = NMNHDIM_UNUSED
 
-    call Les_diachro_write( tpdiafile, XLES_RESOLVED_SvPz, 'RES_SVPZ', 'Resolved <Svdp/dz>', 'kg2 kg-2 Pa m-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_RESOLVED_SvPz, 'RES_SVPZ', 'Resolved <Svdp/dz>', 'kg2 kg-2 Pa m-1', ymasks )
 
     tfield%ndims = 3
     !tfield%ndimlist(1)  = NMNHDIM_BUDGET_LES_LEVEL
@@ -604,9 +592,9 @@ if ( lles_resolved ) then
     !tfield%ndimlist(5:) = NMNHDIM_UNUSED
   end if
 
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UKe, 'RES_UKE', 'Resolved flux of resolved kinetic energy', 'm3 s-3', ysubtitle)
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_VKe, 'RES_VKE', 'Resolved flux of resolved kinetic energy', 'm3 s-3', ysubtitle)
-  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WKe, 'RES_WKE', 'Resolved flux of resolved kinetic energy', 'm3 s-3', ysubtitle)
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_UKe, 'RES_UKE', 'Resolved flux of resolved kinetic energy', 'm3 s-3', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_VKe, 'RES_VKE', 'Resolved flux of resolved kinetic energy', 'm3 s-3', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_RESOLVED_WKe, 'RES_WKE', 'Resolved flux of resolved kinetic energy', 'm3 s-3', ymasks )
 end if
 !
 !
@@ -624,73 +612,73 @@ if ( lles_subgrid ) then
   tfield%ndimlist(3)  = NMNHDIM_BUDGET_LES_MASK
   tfield%ndimlist(4:) = NMNHDIM_UNUSED
 
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_Tke,  'SBG_TKE',  'Subgrid TKE',           'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_U2,   'SBG_U2',   'Subgrid <u2> variance', 'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_V2,   'SBG_V2',   'Subgrid <v2> variance', 'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_W2,   'SBG_W2',   'Subgrid <w2> variance', 'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_UV,   'SBG_UV',   'Subgrid <uv> flux',     'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WU,   'SBG_WU',   'Subgrid <wu> flux',     'm2 s-2', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WV,   'SBG_WV',   'Subgrid <wv> flux',     'm2 s-2', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_Tke,  'SBG_TKE',  'Subgrid TKE',           'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_U2,   'SBG_U2',   'Subgrid <u2> variance', 'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_V2,   'SBG_V2',   'Subgrid <v2> variance', 'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_W2,   'SBG_W2',   'Subgrid <w2> variance', 'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_UV,   'SBG_UV',   'Subgrid <uv> flux',     'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WU,   'SBG_WU',   'Subgrid <wu> flux',     'm2 s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WV,   'SBG_WV',   'Subgrid <wv> flux',     'm2 s-2', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_Thl2, 'SBG_THL2', 'Subgrid liquid potential temperature variance', &
-                          'K2', ysubtITLE )
+                          'K2', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_UThl, 'SBG_UTHL', 'Subgrid hor. flux of liquid potential temperature',  &
-                          'm K s-1', YSUBTITLE )
+                          'm K s-1', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_VThl, 'SBG_VTHL', 'Subgrid hor. flux of liquid potential temperature',  &
-                          'm K s-1', YSUBTITLE )
+                          'm K s-1', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_WThl, 'SBG_WTHL', 'Subgrid vert. flux of liquid potential temperature', &
-                          'm K s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WP,   'SBG_WP',   'Subgrid <wp> vertical Flux', 'm Pa s-1', ysubtitle )
-
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_THLUP_MF, 'THLUP_MF', 'Subgrid <thl> of updraft',    'K',          ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RTUP_MF,  'RTUP_MF',  'Subgrid <rt> of updraft',     'kg kg-1',    ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RVUP_MF,  'RVUP_MF',  'Subgrid <rv> of updraft',     'kg kg-1',    ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RCUP_MF,  'RCUP_MF',  'Subgrid <rc> of updraft',     'kg kg-1',    ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RIUP_MF,  'RIUP_MF',  'Subgrid <ri> of updraft',     'kg kg-1',    ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WUP_MF,   'WUP_MF',   'Subgrid <w> of updraft',      'm s-1',      ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_MASSFLUX, 'MAFLX_MF', 'Subgrid <MF> of updraft',     'kg m-2 s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_DETR,     'DETR_MF',  'Subgrid <detr> of updraft',   'kg m-3 s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_ENTR,     'ENTR_MF',  'Subgrid <entr> of updraft',   'kg m-3 s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_FRACUP,   'FRCUP_MF', 'Subgrid <FracUp> of updraft', '1',          ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_THVUP_MF, 'THVUP_MF', 'Subgrid <thv> of updraft',    'K',          ysubtitle )
+                          'm K s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WP,   'SBG_WP',   'Subgrid <wp> vertical Flux', 'm Pa s-1', ymasks )
+
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_THLUP_MF, 'THLUP_MF', 'Subgrid <thl> of updraft',    'K',          ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RTUP_MF,  'RTUP_MF',  'Subgrid <rt> of updraft',     'kg kg-1',    ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RVUP_MF,  'RVUP_MF',  'Subgrid <rv> of updraft',     'kg kg-1',    ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RCUP_MF,  'RCUP_MF',  'Subgrid <rc> of updraft',     'kg kg-1',    ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_RIUP_MF,  'RIUP_MF',  'Subgrid <ri> of updraft',     'kg kg-1',    ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WUP_MF,   'WUP_MF',   'Subgrid <w> of updraft',      'm s-1',      ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_MASSFLUX, 'MAFLX_MF', 'Subgrid <MF> of updraft',     'kg m-2 s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_DETR,     'DETR_MF',  'Subgrid <detr> of updraft',   'kg m-3 s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_ENTR,     'ENTR_MF',  'Subgrid <entr> of updraft',   'kg m-3 s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_FRACUP,   'FRCUP_MF', 'Subgrid <FracUp> of updraft', '1',          ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_THVUP_MF, 'THVUP_MF', 'Subgrid <thv> of updraft',    'K',          ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_WTHLMF,   'WTHL_MF',  'Subgrid <wthl> of mass flux convection scheme', &
-                          'm K s-1', ysubtitle )
+                          'm K s-1', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_WRTMF,    'WRT_MF',   'Subgrid <wrt> of mass flux convection scheme',  &
-                          'm kg kg-1 s-1', ysubtitle )
+                          'm kg kg-1 s-1', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_WTHVMF,   'WTHV_MF',  'Subgrid <wthv> of mass flux convection scheme', &
-                          'm K s-1', ysubtitle )
+                          'm K s-1', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_WUMF,     'WU_MF',    'Subgrid <wu> of mass flux convection scheme',   &
-                          'm2 s-2', ysubtitle )
+                          'm2 s-2', ymasks )
   call Les_diachro_write( tpdiafile, XLES_SUBGRID_WVMF,     'WV_MF',    'Subgrid <wv> of mass flux convection scheme',   &
-                          'm2 s-2', ysubtitle )
+                          'm2 s-2', ymasks )
 
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_PHI3,  'SBG_PHI3', 'Subgrid Phi3 function',         '1',      ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_LMix,  'SBG_LMIX', 'Subgrid Mixing Length',         '1',      ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_LDiss, 'SBG_LDIS', 'Subgrid Dissipation Length',    '1',      ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_Km,    'SBG_KM',   'Eddy diffusivity for momentum', 'm2 s-1', ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_Kh,    'SBG_KH',   'Eddy diffusivity for heat',     'm2 s-1', ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_PHI3,  'SBG_PHI3', 'Subgrid Phi3 function',         '1',      ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_LMix,  'SBG_LMIX', 'Subgrid Mixing Length',         '1',      ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_LDiss, 'SBG_LDIS', 'Subgrid Dissipation Length',    '1',      ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_Km,    'SBG_KM',   'Eddy diffusivity for momentum', 'm2 s-1', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_Kh,    'SBG_KH',   'Eddy diffusivity for heat',     'm2 s-1', ymasks )
 
   if ( luserv ) then
     call Les_diachro_write( tpdiafile, XLES_SUBGRID_WThv,  'SBG_WTHV', 'Subgrid vert. flux of liquid potential temperature', &
-                            'm K s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_SUBGRID_Rt2,   'SBG_RT2',  'Subgrid total water variance', 'kg2 kg-2', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_SUBGRID_ThlRt, 'SBG_TLRT', 'Subgrid <thlrt> covariance',  'K kg kg-1', ysubtitle )
+                            'm K s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_SUBGRID_Rt2,   'SBG_RT2',  'Subgrid total water variance', 'kg2 kg-2', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_SUBGRID_ThlRt, 'SBG_TLRT', 'Subgrid <thlrt> covariance',  'K kg kg-1', ymasks )
     call Les_diachro_write( tpdiafile, XLES_SUBGRID_URt,   'SBG_URT',  'Subgrid total water horizontal flux', &
-                            'm kg kg-1 s-1', ysubtitle )
+                            'm kg kg-1 s-1', ymasks )
     call Les_diachro_write( tpdiafile, XLES_SUBGRID_VRt,   'SBG_VRT',  'Subgrid total water horizontal flux', &
-                            'm kg kg-1 s-1', ysubtitle )
+                            'm kg kg-1 s-1', ymasks )
     call Les_diachro_write( tpdiafile, XLES_SUBGRID_WRt,   'SBG_WRT',  'Subgrid total water vertical flux',   &
-                            'm kg kg-1 s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_SUBGRID_PSI3,  'SBG_PSI3', 'Subgrid Psi3 function', '1', ysubtitle )
+                            'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_SUBGRID_PSI3,  'SBG_PSI3', 'Subgrid Psi3 function', '1', ymasks )
   end if
 
   if ( luserc ) then
-    call Les_diachro_write( tpdiafile, XLES_SUBGRID_Rc2, 'SBG_RC2', 'Subgrid cloud water variance',        'kg2 kg-2', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_SUBGRID_Rc2, 'SBG_RC2', 'Subgrid cloud water variance',        'kg2 kg-2', ymasks )
     call Les_diachro_write( tpdiafile, XLES_SUBGRID_URc, 'SBG_URC', 'Subgrid cloud water horizontal flux', 'm kg kg-1 s-1', &
-                            ysubtitle )
+                            ymasks )
     call Les_diachro_write( tpdiafile, XLES_SUBGRID_VRc, 'SBG_VRC', 'Subgrid cloud water horizontal flux', 'm kg kg-1 s-1', &
-                            ysubtitle )
+                            ymasks )
     call Les_diachro_write( tpdiafile, XLES_SUBGRID_WRc, 'SBG_WRC', 'Subgrid cloud water vertical flux',   'm kg kg-1 s-1', &
-                            ysubtitle )
+                            ymasks )
   end if
 
   if ( nsv > 0 ) then
@@ -701,9 +689,9 @@ if ( lles_subgrid ) then
     tfield%ndimlist(4)  = NMNHDIM_BUDGET_LES_SV
     tfield%ndimlist(5:) = NMNHDIM_UNUSED
 
-    call Les_diachro_write( tpdiafile, XLES_SUBGRID_USv, 'SBG_USV', 'Subgrid <uSv> horizontal flux', 'm kg kg-1 s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_SUBGRID_VSv, 'SBG_VSV', 'Subgrid <vSv> horizontal flux', 'm kg kg-1 s-1', ysubtitle )
-    call Les_diachro_write( tpdiafile, XLES_SUBGRID_WSv, 'SBG_WSV', 'Subgrid <wSv> vertical flux',   'm kg kg-1 s-1', ysubtitle )
+    call Les_diachro_write( tpdiafile, XLES_SUBGRID_USv, 'SBG_USV', 'Subgrid <uSv> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_SUBGRID_VSv, 'SBG_VSV', 'Subgrid <vSv> horizontal flux', 'm kg kg-1 s-1', ymasks )
+    call Les_diachro_write( tpdiafile, XLES_SUBGRID_WSv, 'SBG_WSV', 'Subgrid <wSv> vertical flux',   'm kg kg-1 s-1', ymasks )
 
     tfield%ndims = 3
     !tfield%ndimlist(1)  = NMNHDIM_BUDGET_LES_LEVEL
@@ -715,16 +703,11 @@ if ( lles_subgrid ) then
 
   end if
 
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_UTke,  'SBG_UTKE', 'Subgrid flux of subgrid kinetic energy', 'm3 s-3',   &
-                          ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_VTke,  'SBG_VTKE', 'Subgrid flux of subgrid kinetic energy', 'm3 s-3',   &
-                          ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WTke,  'SBG_WTKE', 'Subgrid flux of subgrid kinetic energy', 'm3 s-3',   &
-                          ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_W2Thl, 'SBG_W2TL', 'Subgrid flux of subgrid kinetic energy', 'm2 K s-2', &
-                          ysubtitle )
-  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WThl2, 'SBG_WTL2', 'Subgrid flux of subgrid kinetic energy', 'm K2 s-1', &
-                          ysubtitle )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_UTke,  'SBG_UTKE', 'Subgrid flux of subgrid kinetic energy', 'm3 s-3',   ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_VTke,  'SBG_VTKE', 'Subgrid flux of subgrid kinetic energy', 'm3 s-3',   ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WTke,  'SBG_WTKE', 'Subgrid flux of subgrid kinetic energy', 'm3 s-3',   ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_W2Thl, 'SBG_W2TL', 'Subgrid flux of subgrid kinetic energy', 'm2 K s-2', ymasks )
+  call Les_diachro_write( tpdiafile, XLES_SUBGRID_WThl2, 'SBG_WTL2', 'Subgrid flux of subgrid kinetic energy', 'm K2 s-1', ymasks )
 end if
 
 
@@ -1118,12 +1101,14 @@ if ( nspectra_k > 0 ) then
     call Les_diachro_2pt_write( tpdiafile, XCORRi_WRi,   XCORRj_WRi,   'WRI',  'W*ri    2 points correlations', 'm kg s-1 kg-1' )
   end if
 
+!PW: TODO: ameliorer le ygroup (tenir compte de ce qu'est la variable scalaire et pas juste son jsv!)
   do jsv = 1, nsv
     Write( ygroup, fmt = "( a2, i3.3 )" ) "SS", jsv
     call Les_diachro_2pt_write( tpdiafile, XCORRi_SvSv(:,:,:,JSV), XCORRj_SvSv(:,:,:,JSV), ygroup, &
                                 'Sv*Sv   2 points correlations','kg2 kg-2' )
   end do
 
+!PW: TODO: ameliorer le ygroup (tenir compte de ce qu'est la variable scalaire et pas juste son jsv!)
   do jsv = 1, nsv
     Write( ygroup, fmt = "( a2, i3.3 )" ) "WS", jsv
     call Les_diachro_2pt_write( tpdiafile, XCORRi_WSv(:,:,:,JSV), XCORRj_WSv(:,:,:,JSV), ygroup, &
@@ -1205,7 +1190,7 @@ end subroutine Les_diachro_write_2D
 
 !------------------------------------------------------------------------------
 
-subroutine Les_diachro_write_3D( tpdiafile, pdata, hmnhname, hcomment, hunits, hsuffixes )
+subroutine Les_diachro_write_3D( tpdiafile, pdata, hmnhname, hcomment, hunits, hmasks )
 
 use modd_io,          only: tfiledata
 
@@ -1216,20 +1201,20 @@ real,             dimension(:,:,:),       intent(in) :: pdata
 character(len=*),                         intent(in) :: hmnhname
 character(len=*),                         intent(in) :: hcomment
 character(len=*),                         intent(in) :: hunits
-character(len=*), dimension(:), optional, intent(in) :: hsuffixes
+character(len=*), dimension(:), optional, intent(in) :: hmasks
 
 tfield%cmnhname  = hmnhname
 tfield%clongname = hmnhname
 tfield%ccomment  = hcomment
 tfield%cunits    = hunits
 
-call Les_diachro( tpdiafile, tfield, ldoavg, ldonorm, pdata, hsuffixes )
+call Les_diachro( tpdiafile, tfield, ldoavg, ldonorm, pdata, hmasks = hmasks )
 
 end subroutine Les_diachro_write_3D
 
 !------------------------------------------------------------------------------
 
-subroutine Les_diachro_write_4D( tpdiafile, pdata, hmnhname, hcomment, hunits, hsuffixes )
+subroutine Les_diachro_write_4D( tpdiafile, pdata, hmnhname, hcomment, hunits, hmasks )
 
 use modd_io,          only: tfiledata
 
@@ -1240,14 +1225,14 @@ real,             dimension(:,:,:,:),       intent(in) :: pdata
 character(len=*),                           intent(in) :: hmnhname
 character(len=*),                           intent(in) :: hcomment
 character(len=*),                           intent(in) :: hunits
-character(len=*), dimension(:),   optional, intent(in) :: hsuffixes
+character(len=*), dimension(:),   optional, intent(in) :: hmasks
 
 tfield%cmnhname  = hmnhname
 tfield%clongname = hmnhname
 tfield%ccomment  = hcomment
 tfield%cunits    = hunits
 
-call Les_diachro( tpdiafile, tfield, ldoavg, ldonorm, pdata, hsuffixes )
+call Les_diachro( tpdiafile, tfield, ldoavg, ldonorm, pdata, hmasks = hmasks )
 
 end subroutine Les_diachro_write_4D
 
diff --git a/src/MNH/write_profilern.f90 b/src/MNH/write_profilern.f90
index 325c4d3a6..7bda486cc 100644
--- a/src/MNH/write_profilern.f90
+++ b/src/MNH/write_profilern.f90
@@ -655,6 +655,7 @@ tzbudiachro%cname      = ygroup
 tzbudiachro%ccomment   = 'Vertical profiles at position of profiler ' // Trim( ygroup )
 tzbudiachro%ccategory  = 'profiler'
 tzbudiachro%cshape     = 'vertical profile'
+! tzbudiachro%cmask     = NOT SET (default values)
 tzbudiachro%lmobile    = .false.
 tzbudiachro%licompress = .true.
 tzbudiachro%ljcompress = .true.
diff --git a/src/MNH/write_seriesn.f90 b/src/MNH/write_seriesn.f90
index 7cf5d5d2f..6b8b2bcfa 100644
--- a/src/MNH/write_seriesn.f90
+++ b/src/MNH/write_seriesn.f90
@@ -95,7 +95,7 @@ TYPE(TFILEDATA),    INTENT(IN) :: TPDIAFILE    ! file to write
 !
 INTEGER                              :: IIB,IJB,IKB     ! Begin of physical dimensions
 INTEGER                              :: IIE,IJE,IKE     ! End   of physical dimensions
-INTEGER                              :: JS,JT,JJ,JI      ! Loop indices 
+INTEGER                              :: JS, JT, JJ, JI, JP ! Loop indices
 INTEGER                              :: ISB1,ISB2
 CHARACTER (LEN=2)                    :: YS             ! String for y-slice
 CHARACTER (LEN=3)                    :: YSL,YSH        ! Strings for y-slice
@@ -251,19 +251,22 @@ tzfields(:)%cunits    = csunit1(:)
 tzfields(:)%ccomment  = cscomment1(:)
 tzfields(:)%ngrid     = nsgridd1(:)
 tzfields(:)%ntype     = TYPEREAL
-tzfields(:)%ndims     = 2
+! tzfields(:)%ndims     = 2
+tzfields(:)%ndims     = 1 !Set to 1 because write are done in a loop (1 write per "process")
 tzfields(:)%ndimlist(1) = NMNHDIM_UNUSED
 tzfields(:)%ndimlist(2) = NMNHDIM_UNUSED
 tzfields(:)%ndimlist(3) = NMNHDIM_UNUSED
 tzfields(:)%ndimlist(4) = NMNHDIM_SERIES_TIME
 tzfields(:)%ndimlist(5) = NMNHDIM_UNUSED
-tzfields(:)%ndimlist(6) = NMNHDIM_SERIES_PROC
+! tzfields(:)%ndimlist(6) = NMNHDIM_SERIES_PROC
+tzfields(:)%ndimlist(6) = NMNHDIM_UNUSED !Set to unused because write are done in a loop (1 write per "process")
 
 tzbudiachro%cgroupname = 'TSERIES'
 tzbudiachro%cname      = 'TSERIES'
 tzbudiachro%ccomment   = 'Time series of horizontally and vertically averaged fields'
 tzbudiachro%ccategory  = 'time series'
 tzbudiachro%cshape     = 'cartesian' !It is based on a cartesian domain (with compression in all directions)
+! tzbudiachro%cmask    =  set in the process loop
 tzbudiachro%lmobile    = .false.
 tzbudiachro%licompress = .true.
 tzbudiachro%ljcompress = .true.
@@ -277,8 +280,14 @@ tzbudiachro%njh        = njboxh
 tzbudiachro%nkl        = 1
 tzbudiachro%nkh        = ikmax
 
-call Write_diachro( tpdiafile, tzbudiachro, tzfields, tpsdates(1:nsnbstept), &
-                    xsseries1(1:1,1:1,1:1,1:nsnbstept,1:1,:)                 )
+! Loop on the different masks
+! Do not provide all tzfields once because they can be stored in different HDF groups (based on masks)
+do jp = 1 , nstemp_serie1
+  tzbudiachro%cmask    = csmask1(jp)
+
+  call Write_diachro( tpdiafile, tzbudiachro, [ tzfields(jp) ], tpsdates(1:nsnbstept), &
+                      xsseries1(1:1,1:1,1:1,1:nsnbstept,1:1,jp:jp)                 )
+end do
 
 deallocate( tzfields )
 !
@@ -338,7 +347,8 @@ tzfields(:)%cunits    = csunit2(:)
 tzfields(:)%ccomment  = cscomment2(:)
 tzfields(:)%ngrid     = nsgridd2(:)
 tzfields(:)%ntype     = TYPEREAL
-tzfields(:)%ndims     = 3
+! tzfields(:)%ndims     = 3
+tzfields(:)%ndims     = 2 !Set to 2 because write are done in a loop (1 write per "process")
 tzfields(:)%ndimlist(1) = NMNHDIM_UNUSED
 tzfields(:)%ndimlist(2) = NMNHDIM_UNUSED
 do ji = 1, nstemp_serie2
@@ -353,13 +363,15 @@ do ji = 1, nstemp_serie2
 end do
 tzfields(:)%ndimlist(4) = NMNHDIM_SERIES_TIME
 tzfields(:)%ndimlist(5) = NMNHDIM_UNUSED
-tzfields(:)%ndimlist(6) = NMNHDIM_SERIES_PROC
+! tzfields(:)%ndimlist(6) = NMNHDIM_SERIES_PROC
+tzfields(:)%ndimlist(6) = NMNHDIM_UNUSED !Set to unused because write are done in a loop (1 write per "process")
 
 tzbudiachro%cgroupname = 'ZTSERIES'
 tzbudiachro%cname      = 'ZTSERIES'
 tzbudiachro%ccomment   = 'Time series of horizontally averaged vertical profile'
 tzbudiachro%ccategory  = 'time series'
 tzbudiachro%cshape     = 'cartesian'  !It is based on a cartesian domain (with horizontal compression)
+! tzbudiachro%cmask    =  set in the process loop
 tzbudiachro%lmobile    = .false.
 tzbudiachro%licompress = .true.
 tzbudiachro%ljcompress = .true.
@@ -373,8 +385,14 @@ tzbudiachro%njh        = njboxh
 tzbudiachro%nkl        = 1
 tzbudiachro%nkh        = ikmax
 
-call Write_diachro( tpdiafile, tzbudiachro, tzfields, tpsdates(1:nsnbstept), &
-                    xsseries2(1:1,1:1,1:ikmax,1:nsnbstept,1:1,:)             )
+! Loop on the different masks
+! Do not provide all tzfields once because they can be stored in different HDF groups (based on masks)
+do jp = 1 , nstemp_serie2
+  tzbudiachro%cmask    = csmask2(jp)
+
+  call Write_diachro( tpdiafile, tzbudiachro, [ tzfields(jp) ], tpsdates(1:nsnbstept), &
+                      xsseries2(1:1,1:1,1:ikmax,1:nsnbstept,1:1,jp:jp)                 )
+end do
 
 deallocate( tzfields )
 !
@@ -460,6 +478,7 @@ DO JS=1,NBJSLICE
   tzbudiachro%ccomment   = 'Time series of y-horizontally averaged fields at one level or vertically averaged between 2 levels'
   tzbudiachro%ccategory  = 'time series'
   tzbudiachro%cshape     = 'cartesian' !It is based on a cartesian domain (with compression in 1 direction)
+!  tzbudiachro%cmask      = NOT SET (default values)
   tzbudiachro%lmobile    = .false.
   tzbudiachro%licompress = .false.
   tzbudiachro%ljcompress = .true.
diff --git a/src/MNH/write_stationn.f90 b/src/MNH/write_stationn.f90
index 36c118023..d20e5ef0c 100644
--- a/src/MNH/write_stationn.f90
+++ b/src/MNH/write_stationn.f90
@@ -740,6 +740,7 @@ tzbudiachro%cname      = ygroup
 tzbudiachro%ccomment   = 'Values at position of station ' // Trim( ygroup )
 tzbudiachro%ccategory  = 'station'
 tzbudiachro%cshape     = 'point'
+! tzbudiachro%cmask     = NOT SET (default values)
 tzbudiachro%lmobile    = .false.
 tzbudiachro%licompress = .true.
 tzbudiachro%ljcompress = .true.
-- 
GitLab