From 27689d9b20e3d5e0d94929817f5d1405a8a14d47 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 28 Oct 2020 11:25:46 +0100
Subject: [PATCH] Philippe 28/10/2020: IO: add Write_diachro_nc4 subroutine

---
 src/MNH/write_diachro.f90 | 1144 ++++++++++++++++++++++++++++++++++++-
 1 file changed, 1143 insertions(+), 1 deletion(-)

diff --git a/src/MNH/write_diachro.f90 b/src/MNH/write_diachro.f90
index a2fe6cec6..15cfb6000 100644
--- a/src/MNH/write_diachro.f90
+++ b/src/MNH/write_diachro.f90
@@ -549,10 +549,1152 @@ call IO_Field_write( tpdiafile, tzfield, zdatime )
 
 deallocate( zdatime )
 !
+#ifdef MNH_IOCDF4
+if ( tpdiafile%cformat == 'NETCDF4' .or. tpdiafile%cformat == 'LFICDF4' ) &
+  call Write_diachro_nc4( tpdiafile, tpfields, hgroup, htype, tpdates, pvar, gicp, gjcp, gkcp, kil, kih, kjl, kjh, kkl, kkh )
+#endif
+
 CALL MENU_DIACHRO(TPDIAFILE,HGROUP)
 LPACK=GPACK
+
+end subroutine Write_diachro
+
 !-----------------------------------------------------------------------------
+#ifdef MNH_IOCDF4
+subroutine Write_diachro_nc4( tpdiafile, tpfields, hgroup, htype, tpdates, pvar, oicp, ojcp, okcp, kil, kih, kjl, kjh, kkl, kkh )
+
+use NETCDF,            only: NF90_DEF_DIM, NF90_DEF_GRP, NF90_DEF_VAR, NF90_INQ_NCID, NF90_PUT_ATT, NF90_PUT_VAR, &
+                             NF90_GLOBAL, NF90_NOERR, NF90_STRERROR
+
+use modd_field
+use modd_io,           only: isp, tfiledata
+use modd_les,          only: nles_masks
+use modd_parameters,   only: jphext
+use modd_precision,    only: CDFINT, MNHREAL_NF90
+use modd_type_date,    only: date_time
+
+use mode_io_field_write, only: IO_Field_write, IO_Field_write_box
+use mode_io_tools_nc4,   only: IO_Err_handle_nc4
+use mode_msg
+
+type(tfiledata),                                     intent(in)           :: tpdiafile        ! File to write
+class(tfield_metadata_base), dimension(:),           intent(in)           :: tpfields
+character(len=*),                                    intent(in)           :: hgroup, htype
+type(date_time),             dimension(:),           intent(in)           :: tpdates
+real,                        dimension(:,:,:,:,:,:), intent(in)           :: pvar
+logical,                                             intent(in)           :: oicp, ojcp, okcp
+integer,                                             intent(in), optional :: kil, kih
+integer,                                             intent(in), optional :: kjl, kjh
+integer,                                             intent(in), optional :: kkl, kkh
+
+character(len=3)                          :: ynum
+integer                                   :: icompx, icompy, icompz
+integer                                   :: idims
+integer                                   :: ji
+integer(kind=CDFINT)                      :: isavencid
+integer(kind=CDFINT)                      :: idimid
+integer(kind=CDFINT)                      :: igrpid
+integer(kind=CDFINT)                      :: istatus
+integer(kind=CDFINT)                      :: idimtimeid
+real                                      :: zdata0d
+real, dimension(:),           allocatable :: zdata1d
+real, dimension(:,:),         allocatable :: zdata2d
+real, dimension(:,:,:),       allocatable :: zdata3d
+real, dimension(:,:,:,:),     allocatable :: zdata4d
+real, dimension(:,:,:,:,:),   allocatable :: zdata5d
+real, dimension(:,:,:,:,:,:), allocatable :: zdata6d
+type(tfielddata)                          :: tzfield
+type(tfiledata)                           :: tzfile
+
+
+if ( trim ( htype ) == 'CART' .or. trim ( htype ) == 'MASK' .or. trim ( htype ) == 'SPXY') then
+  if (        .not. Present( kil ) .or. .not. Present( kih ) &
+         .or. .not. Present( kjl ) .or. .not. Present( kjh ) &
+         .or. .not. Present( kkl ) .or. .not. Present( kkh ) ) then
+    call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
+                      'kil, kih, kjl, kjh, kkl or kkh not provided for variable ' // Trim( tpfields(1)%cmnhname ) )
+  end if
+end if
+
+tzfile = tpdiafile
+
+if ( oicp ) then
+  icompx = 1
+else
+  icompx = 0
+endif
+if ( ojcp ) then
+  icompy = 1
+else
+  icompy = 0
+endif
+if ( okcp ) then
+  icompz = 1
+else
+  icompz = 0
+endif
+
+MASTER: if ( isp == tzfile%nmaster_rank) then
+  istatus = NF90_INQ_NCID( tzfile%nncid, trim( hgroup ), igrpid )
+  if ( istatus == NF90_NOERR ) then
+    call Print_msg( NVERB_WARNING, 'IO', 'Write_diachro_nc4', trim(tzfile%cname)//': group '//trim(hgroup)//' already defined' )
+  else
+    istatus = NF90_DEF_GRP( tzfile%nncid, trim( hgroup ), igrpid )
+    if ( istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_DEF_GRP', 'for '//trim(hgroup)//' group' )
+  end if
+
+  !Save id of the file root group ('/' group)
+  isavencid = tzfile%nncid
+  tzfile%nncid = igrpid
+
+  !Write only in netCDF files
+  tzfile%cformat = 'NETCDF4'
+
+  istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'type', trim( htype ) )
+  if (istatus /= NF90_NOERR ) &
+    call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'type for '//trim(hgroup)//' group' )
+
+  if ( trim ( htype ) == 'CART' .or. trim ( htype ) == 'MASK' .or. trim ( htype ) == 'SPXY') then
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'min x index', kil )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'min x index for '//trim(hgroup)//' group' )
+
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'max x index', kih )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'max x index for '//trim(hgroup)//' group' )
+
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'min y index', kjl )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'min y index for '//trim(hgroup)//' group' )
+
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'max y index', kjh )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'max y index for '//trim(hgroup)//' group' )
+
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'min z index', kkl )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'min z index for '//trim(hgroup)//' group' )
+
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'max z index', kkh )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'max z index for '//trim(hgroup)//' group' )
+  end if
+
+  if ( trim ( htype ) == 'CART' ) then
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'averaged on x dimension', icompx )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'averaged on x dimension '//trim(hgroup)//' group' )
+
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'averaged on y dimension', icompy )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'averaged on y dimension '//trim(hgroup)//' group' )
+  end if
+
+  if ( trim ( htype ) == 'CART' .or. trim ( htype ) == 'MASK' .or. trim ( htype ) == 'SPXY') then
+    istatus = NF90_PUT_ATT( igrpid, NF90_GLOBAL, 'averaged on z dimension', icompz )
+    if (istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', 'averaged on z dimension '//trim(hgroup)//' group' )
+  end if
+
+end if MASTER
+
+if ( Count( tpfields(1)%ndimlist(:) == NMNHDIM_UNKNOWN ) /= 0 ) &
+  call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
+                  'some dimensions are unknown for variable '//trim(tpfields(1)%cmnhname) )
+
+idims = Count( tpfields(1)%ndimlist(:) /= NMNHDIM_UNUSED )
+
+if ( tpfields(1)%ndims /= idims ) &
+  call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
+                  'ndims is not coherent with ndimlist for variable '//trim(tpfields(1)%cmnhname) )
+
+do ji = 2, NMNHMAXDIMS
+  if ( tpfields(1)%ndimlist(ji) == NMNHDIM_UNKNOWN ) then
+    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
+                    'NMNHDIM_UNKNOWN for some dimensions of variable '//trim(tpfields(1)%cmnhname) )
+    idims = 6
+  end if
+end do
+
+!Check that if 'CART' and no horizontal compression, parameters are as expected
+if ( htype == 'CART' .and. .not. oicp .and. .not. ojcp ) then
+  if ( idims /= 3 .and. idims /= 4 ) then
+    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4',                                                 &
+                    'unexpected number of dimensions for CART without horizontal compression for variable ' &
+                    // Trim( tpfields(1)%cmnhname ) )
+  end if
+
+  if (      (       tpfields(1)%ndimlist(1) /= NMNHDIM_BUDGET_CART_NI              &
+              .and. tpfields(1)%ndimlist(1) /= NMNHDIM_BUDGET_CART_NI_U            &
+              .and. tpfields(1)%ndimlist(1) /= NMNHDIM_BUDGET_CART_NI_V          ) &
+       .or. (       tpfields(1)%ndimlist(2) /= NMNHDIM_BUDGET_CART_NJ              &
+              .and. tpfields(1)%ndimlist(2) /= NMNHDIM_BUDGET_CART_NJ_U            &
+              .and. tpfields(1)%ndimlist(2) /= NMNHDIM_BUDGET_CART_NJ_V          ) &
+       .or. (       tpfields(1)%ndimlist(3) /= NMNHDIM_BUDGET_CART_LEVEL           &
+              .and. tpfields(1)%ndimlist(3) /= NMNHDIM_BUDGET_CART_LEVEL_W       ) &
+       .or. ( idims == 4 .and. tpfields(1)%ndimlist(6) /= NMNHDIM_BUDGET_NGROUPS ) ) then
+    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4',                                       &
+                    'unexpected dimensions for CART without horizontal compression for variable ' &
+                    // Trim( tpfields(1)%cmnhname ) )
+  end if
+
+  if ( .not. Present( kil ) .or. .not. Present( kih ) .or. .not. Present( kjl ) .or. .not. Present( kjh ) ) then
+    call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4',                                                &
+                    'kil or kih or kjl or kjh not provided for variable ' // Trim( tpfields(1)%cmnhname ) )
+  end if
+end if
+
+select case ( idims )
+  case (0)
+    zdata0d = pvar(1, 1, 1, 1, 1, 1)
+
+    TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+    TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+    TZFIELD%CLONGNAME  = tpfields(1)%clongname
+    TZFIELD%CUNITS     = tpfields(1)%cunits
+    TZFIELD%CDIR       = '--'
+    TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+    TZFIELD%NGRID      = tpfields(1)%ngrid
+    TZFIELD%NTYPE      = tpfields(1)%ntype
+    TZFIELD%NDIMS      = 1
+    TZFIELD%LTIMEDEP   = .FALSE.
+    CALL IO_Field_write( tzfile, tzfield, zdata0d )
+
+  case (1)
+    if ( tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+
+      allocate( zdata1d( size(pvar,4) ) )
+
+      zdata1d(:) = pvar(1, 1, 1, :, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 1
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata1d )
+
+      deallocate( zdata1d )
+    else if (      tpfields(1)%ndimlist(1) == NMNHDIM_NI   &
+              .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_U   &
+              .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_V   &
+              .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI   &
+              .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_U &
+              .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_V ) then
+      tzfield%ndimlist(1) = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+
+      allocate( zdata1d( size(pvar,1) ) )
+
+      zdata1d(:) = pvar(:, 1, 1, 1, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 1
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata1d )
+
+      deallocate( zdata1d )
+    else if (      tpfields(1)%ndimlist(2) == NMNHDIM_NJ   &
+              .or. tpfields(1)%ndimlist(2) == NMNHDIM_NJ_U   &
+              .or. tpfields(1)%ndimlist(2) == NMNHDIM_NJ_V   &
+              .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ   &
+              .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_U &
+              .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_V ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+
+      allocate( zdata1d( size(pvar,2) ) )
+
+      zdata1d(:) = pvar(1, :, 1, 1, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 1
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata1d )
+
+      deallocate( zdata1d )
+    else if (      tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL   &
+              .or. tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL_W   &
+              .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL   &
+              .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL_W ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+
+      allocate( zdata1d( size(pvar,3) ) )
+
+      zdata1d(:) = pvar(1, 1, :, 1, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 1
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata1d )
+
+      deallocate( zdata1d )
+    else
+      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
+                      'case not yet implemented (variable '//trim(tpfields(1)%cmnhname)//')' )
+    end if
+
+
+  case (2)
+    if (       ( tpfields(1)%ndimlist(1) == NMNHDIM_NI .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_U .or. &
+                 tpfields(1)%ndimlist(1) == NMNHDIM_NI_V .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI &
+                 .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_U &
+                 .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_V )     &
+         .and. ( tpfields(1)%ndimlist(2) == NMNHDIM_NJ .or. tpfields(1)%ndimlist(2) == NMNHDIM_NJ_U .or. &
+                  tpfields(1)%ndimlist(2) == NMNHDIM_NJ_V .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ .or. &
+                  tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_U &
+                 .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_V ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,1), size(pvar,2) ) )
+
+      zdata2d(:,:) = pvar(:, :, 1, 1, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 2
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata2d )
+
+      deallocate( zdata2d )
+
+    else if (       ( tpfields(1)%ndimlist(1) == NMNHDIM_NI .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_U &
+                      .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_V .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI &
+                      .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_U &
+                 .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_V )     &
+       .and.  ( tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL_W &
+           .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL &
+           .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL_W ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,1), size(pvar,3) ) )
+
+      zdata2d(:,:) = pvar(:, 1, :, 1, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 2
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata2d )
+
+      deallocate( zdata2d )
+
+    else if (       ( tpfields(1)%ndimlist(2) == NMNHDIM_NJ .or. tpfields(1)%ndimlist(2) == NMNHDIM_NJ_U &
+                     .or. tpfields(1)%ndimlist(2) == NMNHDIM_NJ_V .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ &
+                     .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_U &
+                 .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_V )     &
+       .and.  ( tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL_W &
+           .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL &
+           .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL_W ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,2), size(pvar,3) ) )
+
+      zdata2d(:,:) = pvar(1, :, :, 1, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 2
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata2d )
+
+      deallocate( zdata2d )
+
+    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,3), size(pvar,4) ) )
+
+      zdata2d(:,:) = pvar(1, 1, :, :, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 2
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata2d )
+
+      deallocate( zdata2d )
+
+    else if (  (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+         .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV       ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,4), size(pvar,5) ) )
+
+      zdata2d(:,:) = pvar(1, 1, 1, :, :, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 2
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata2d )
+
+      deallocate( zdata2d )
+    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_FLYER_TIME &
+         .and. tpfields(1)%ndimlist(6) == NMNHDIM_FLYER_PROC ) then
+      !Correspond to FLYER_DIACHRO
+
+      !Create local time dimension
+      if ( isp == tzfile%nmaster_rank) then
+        istatus = NF90_DEF_DIM( igrpid, 'time_flyer', Size( pvar, 4), idimid )
+        if ( istatus /= NF90_NOERR ) &
+          call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_DEF_DIM', Trim( tpfields(1)%cmnhname ) )
+      end if
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+
+      allocate( zdata1d( size(pvar,4) ) )
+
+      ! Loop on the processes (1 written variable per process)
+      do ji = 1, size(pvar,6)
+        zdata1d(:) = pvar(1, 1, 1, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 1
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata1d )
+      end do
+
+      deallocate( zdata1d )
+    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_STATION_TIME &
+         .and. tpfields(1)%ndimlist(6) == NMNHDIM_STATION_PROC ) then
+      !Correspond to STATION_DIACHRO_n
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+
+      allocate( zdata1d( size(pvar,4) ) )
+
+      ! Loop on the processes (1 written variable per process)
+      do ji = 1, size(pvar,6)
+        zdata1d(:) = pvar(1, 1, 1, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 1
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata1d )
+      end do
+
+      deallocate( zdata1d )
+    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME &
+         .and. tpfields(1)%ndimlist(6) == NMNHDIM_SERIES_PROC ) then
+      !Correspond to WRITE_SERIES_n
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+
+      allocate( zdata1d( size(pvar,4) ) )
+
+      ! Loop on the processes (1 written variable per process)
+      do ji = 1, size(pvar,6)
+        zdata1d(:) = pvar(1, 1, 1, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 1
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata1d )
+      end do
+
+      deallocate( zdata1d )
+    else
+      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
+                      'case not yet implemented (variable '//trim(tpfields(1)%cmnhname)//')' )
+    end if
+
+
+  case (3)
+    if (       ( tpfields(1)%ndimlist(1) == NMNHDIM_NI .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_U .or. &
+                 tpfields(1)%ndimlist(1) == NMNHDIM_NI_V .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI &
+                 .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_U &
+                 .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_V )     &
+         .and. ( tpfields(1)%ndimlist(2) == NMNHDIM_NJ .or. tpfields(1)%ndimlist(2) == NMNHDIM_NJ_U .or. &
+                  tpfields(1)%ndimlist(2) == NMNHDIM_NJ_V .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ .or. &
+                  tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_U &
+                 .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_V ) &
+       .and.  ( tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL_W &
+           .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL &
+           .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL_W ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,1), size(pvar,2), size(pvar,3) ) )
+
+      zdata3d(:,:,:) = pvar(:, :, :, 1, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 3
+      TZFIELD%LTIMEDEP   = .FALSE.
+      if ( htype == 'CART' .and. .not. oicp .and. .not. ojcp ) then
+        !Data is distributed between all the processes
+        TZFIELD%CDIR     = 'XY'
+        CALL IO_Field_write_BOX( tzfile, tzfield, 'BUDGET', zdata3d, &
+                                 KIL + JPHEXT, KIH + JPHEXT, KJL + JPHEXT, KJH + JPHEXT )
+      else
+        !Data is already collected on the master process
+        TZFIELD%CDIR     = '--'
+        CALL IO_Field_write( tzfile, tzfield, zdata3d )
+      end if
+
+      deallocate( zdata3d )
+    else if (  ( tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL &
+                  .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL_W ) &
+       .and.    tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_MASK_TIME &
+       .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_MASK_NBUMASK      ) then
+      !Correspond to Store_one_budget_rho (MASK)
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,3), size(pvar,4), size(pvar,5) ) )
+
+      zdata3d(:,:,:) = pvar(1, 1, :, :, :, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 3
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata3d )
+
+      deallocate( zdata3d )
+    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+       .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_LES_MASK      ) then
+      if ( nles_masks /= Size( pvar, 6 ) ) &
+        call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4',                             &
+                        'last dimension size of pvar is not equal to nles_masks (variable ' &
+                        // Trim( tpfields(1)%cmnhname ) // ')' )
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,3), size(pvar,4) ) )
+
+      ! Loop on the masks (1 written variable per mask)
+      do ji = 1, size(pvar,6)
+        zdata2d(:,:) = pvar(1, 1, :, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 2
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata2d )
+      end do
+
+      deallocate( zdata2d )
+    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+       .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_TERM      ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,3), size(pvar,4) ) )
+
+      ! Loop on the masks (1 written variable per mask)
+      do ji = 1, size(pvar,6)
+        zdata2d(:,:) = pvar(1, 1, :, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 2
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata2d )
+      end do
+
+      deallocate( zdata2d )
+    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+       .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV      ) then
+      !Correspond to Les_diachro_sv_new
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,3), size(pvar,4), size(pvar,5) ) )
+
+      zdata3d(:,:,:) = pvar(1, 1, :, :, :, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 3
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata3d )
+
+      deallocate( zdata3d )
+    else if (       tpfields(1)%ndimlist(1) == NMNHDIM_SPECTRA_2PTS_NI                   &
+              .and. tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL                &
+              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
+                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,1), size(pvar,3), size(pvar,4) ) )
+
+      zdata3d(:,:,:) = pvar(:, 1, :, :, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 3
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata3d )
+
+      deallocate( zdata3d )
+    else if (       tpfields(1)%ndimlist(2) == NMNHDIM_SPECTRA_2PTS_NJ                   &
+              .and. tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL                &
+              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
+                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,2), size(pvar,3), size(pvar,4) ) )
+
+      zdata3d(:,:,:) = pvar(1, :, :, :, 1, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 3
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata3d )
+
+      deallocate( zdata3d )
+    else if (  tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL      &
+         .and. tpfields(1)%ndimlist(4) == NMNHDIM_FLYER_TIME &
+         .and. tpfields(1)%ndimlist(6) == NMNHDIM_FLYER_PROC ) then
+      !Correspond to FLYER_DIACHRO
+
+      !Create local time dimension
+      if ( isp == tzfile%nmaster_rank) then
+        istatus = NF90_DEF_DIM( igrpid, 'time_flyer', Size( pvar, 4), idimid )
+        if ( istatus /= NF90_NOERR ) &
+          call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_DEF_DIM', Trim( tpfields(1)%cmnhname ) )
+      end if
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,3), size(pvar,4) ) )
+
+      ! Loop on the processes (1 written variable per process)
+      do ji = 1, size(pvar,6)
+        zdata2d(:, :) = pvar(1, 1, :, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 2
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata2d )
+      end do
+
+      deallocate( zdata2d )
+    else if (  tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL      &
+         .and. tpfields(1)%ndimlist(4) == NMNHDIM_PROFILER_TIME &
+         .and. tpfields(1)%ndimlist(6) == NMNHDIM_PROFILER_PROC ) then
+      !Correspond to PROFILER_DIACHRO_n
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,3), size(pvar,4) ) )
+
+      ! Loop on the processes (1 written variable per process)
+      do ji = 1, size(pvar,6)
+        zdata2d(:, :) = pvar(1, 1, :, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 2
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata2d )
+      end do
+
+      deallocate( zdata2d )
+    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 &
+         .and. tpfields(1)%ndimlist(6) == NMNHDIM_SERIES_PROC ) then
+      !Correspond to PROFILER_DIACHRO_n
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,3), size(pvar,4) ) )
+
+      ! Loop on the processes (1 written variable per process)
+      do ji = 1, size(pvar,6)
+        zdata2d(:, :) = pvar(1, 1, :, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 2
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata2d )
+      end do
+
+      deallocate( zdata2d )
+    else if (  ( tpfields(1)%ndimlist(1) == NMNHDIM_NI .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_U )      &
+         .and. tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME &
+         .and. tpfields(1)%ndimlist(6) == NMNHDIM_SERIES_PROC ) then
+      !Correspond to PROFILER_DIACHRO_n
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3:) = NMNHDIM_UNUSED
+
+      allocate( zdata2d( size(pvar,1), size(pvar,4) ) )
+
+      ! Loop on the processes (1 written variable per process)
+      do ji = 1, size(pvar,6)
+        zdata2d(:, :) = pvar(:, 1, 1, :, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 2
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata2d )
+      end do
+
+      deallocate( zdata2d )
+    else
+      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
+                      'case not yet implemented (variable '//trim(tpfields(1)%cmnhname)//')' )
+    end if
+
+  case (4)
+    if (       ( tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_U &
+                 .or. tpfields(1)%ndimlist(1) == NMNHDIM_BUDGET_CART_NI_V )     &
+         .and. ( tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_U &
+                 .or. tpfields(1)%ndimlist(2) == NMNHDIM_BUDGET_CART_NJ_V ) &
+       .and.  ( tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL &
+                  .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_CART_LEVEL_W ) &
+       .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS      ) then
+      !Correspond to Store_one_budget (CART)
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,1), size(pvar,2), size(pvar,3) ) )
+
+      ! Loop on the processes
+      do ji = 1, size(pvar,6)
+        zdata3d(:,:,:) = pvar(:, :, :, 1, 1, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 3
+        TZFIELD%LTIMEDEP   = .FALSE.
+        if ( htype == 'CART' .and. .not. oicp .and. .not. ojcp ) then
+          !Data is distributed between all the processes
+          TZFIELD%CDIR     = 'XY'
+          CALL IO_Field_write_BOX( tzfile, tzfield, 'BUDGET', zdata3d, &
+                                   KIL + JPHEXT, KIH + JPHEXT, KJL + JPHEXT, KJH + JPHEXT )
+        else
+          !Data is already collected on the master process
+          TZFIELD%CDIR     = '--'
+          CALL IO_Field_write( tzfile, tzfield, zdata3d )
+        end if
+      end do
+
+      deallocate( zdata3d )
+
+    elseif (  ( tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL &
+                  .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL_W ) &
+       .and.    tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_MASK_TIME &
+       .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_MASK_NBUMASK &
+       .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS      ) then
+      !Correspond to Store_one_budget (MASK)
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,3), size(pvar,4), size(pvar,5) ) )
+
+      ! Loop on the processes
+      do ji = 1, size(pvar,6)
+        zdata3d(:,:,:) = pvar(1, 1, :, :, :, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 3
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata3d )
+      end do
+
+      deallocate( zdata3d )
+    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+       .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV      &
+       .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_LES_MASK      ) then
+      if ( nles_masks /= Size( pvar, 6 ) ) &
+        call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4',                             &
+                        'last dimension size of pvar is not equal to nles_masks (variable ' &
+                        // Trim( tpfields(1)%cmnhname ) // ')' )
+
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,3), size(pvar,4), size(pvar,5) ) )
+
+      ! Loop on the masks (1 written variable per mask)
+      do ji = 1, size(pvar,6)
+        zdata3d(:,:,:) = pvar(1, 1, :, :, :, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 3
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata3d )
+      end do
+
+      deallocate( zdata3d )
+    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+       .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV      &
+       .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_TERM      ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(4:) = NMNHDIM_UNUSED
+
+      allocate( zdata3d( size(pvar,3), size(pvar,4), size(pvar,5) ) )
+
+      ! Loop on the masks (1 written variable per mask)
+      do ji = 1, size(pvar,6)
+        zdata3d(:,:,:) = pvar(1, 1, :, :, :, ji)
+
+        TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+        TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+        TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+        TZFIELD%CUNITS     = tpfields(ji)%cunits
+        TZFIELD%CDIR       = '--'
+        TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+        TZFIELD%NGRID      = tpfields(ji)%ngrid
+        TZFIELD%NTYPE      = tpfields(ji)%ntype
+        TZFIELD%NDIMS      = 3
+        TZFIELD%LTIMEDEP   = .FALSE.
+        CALL IO_Field_write( tzfile, tzfield, zdata3d )
+      end do
+
+      deallocate( zdata3d )
+    else if (       tpfields(1)%ndimlist(1) == NMNHDIM_SPECTRA_SPEC_NI     &
+       .and.        tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+       .and. tpfields(1)%ndimlist(5) == NMNHDIM_COMPLEX      ) then
+      !Correspond to LES_DIACHRO_SPEC
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(4)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(5:) = NMNHDIM_UNUSED
+
+      allocate( zdata4d( size(pvar,1), size(pvar,3), size(pvar,4), size(pvar,5) ) )
+
+      zdata4d(:,:,:,:) = pvar(:, 1, :, :, :, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 4
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata4d )
+
+      deallocate( zdata4d )
+    else if (       tpfields(1)%ndimlist(2) == NMNHDIM_SPECTRA_SPEC_NJ     &
+       .and.        tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL &
+       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
+               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
+       .and. tpfields(1)%ndimlist(5) == NMNHDIM_COMPLEX      ) then
+      !Correspond to LES_DIACHRO_SPEC
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(4)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(5:) = NMNHDIM_UNUSED
+
+      allocate( zdata4d( size(pvar,2), size(pvar,3), size(pvar,4), size(pvar,5) ) )
+
+      zdata4d(:,:,:,:) = pvar(1, :, :, :, :, 1)
+
+      TZFIELD%CMNHNAME   = tpfields(1)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(1)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(1)%clongname
+      TZFIELD%CUNITS     = tpfields(1)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(1)%ccomment
+      TZFIELD%NGRID      = tpfields(1)%ngrid
+      TZFIELD%NTYPE      = tpfields(1)%ntype
+      TZFIELD%NDIMS      = 4
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, zdata4d )
+
+      deallocate( zdata4d )
+    else
+      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
+                      'case not yet implemented (variable '//trim(tpfields(1)%cmnhname)//')' )
+    end if
+
+!   case (5)
+
+!   case (6)
+
+  case default
+    if ( All( tpfields(1)%ndimlist(:) /= NMNHDIM_UNKNOWN ) ) then
+      tzfield%ndimlist(1)  = tpfields(1)%ndimlist(1)
+      tzfield%ndimlist(2)  = tpfields(1)%ndimlist(2)
+      tzfield%ndimlist(3)  = tpfields(1)%ndimlist(3)
+      tzfield%ndimlist(4)  = tpfields(1)%ndimlist(4)
+      tzfield%ndimlist(5)  = tpfields(1)%ndimlist(5)
+      tzfield%ndimlist(6:) = NMNHDIM_UNUSED
+    end if
+
+    do ji = 1, Size( pvar, 6 )
+      TZFIELD%CMNHNAME   = tpfields(ji)%cmnhname
+      TZFIELD%CSTDNAME   = tpfields(ji)%cstdname
+      TZFIELD%CLONGNAME  = tpfields(ji)%clongname
+      TZFIELD%CUNITS     = tpfields(ji)%cunits
+      TZFIELD%CDIR       = '--'
+      TZFIELD%CCOMMENT   = tpfields(ji)%ccomment
+      TZFIELD%NGRID      = tpfields(ji)%ngrid
+      TZFIELD%NTYPE      = tpfields(ji)%ntype
+      TZFIELD%NDIMS      = 5
+      TZFIELD%LTIMEDEP   = .FALSE.
+      CALL IO_Field_write( tzfile, tzfield, pvar(:, :, :, :, :, ji ) )
+    end do
+end select
+
+!Reset ndimlist (to prevent problems later)
+tzfield%ndimlist(:) = NMNHDIM_UNKNOWN
+
+TZFIELD%CMNHNAME   = 'dates'
+TZFIELD%CSTDNAME   = ''
+TZFIELD%CLONGNAME  = 'dates'
+TZFIELD%CUNITS     = 'seconds since YYYY-MM-DD HH:MM:SS.S'
+TZFIELD%CDIR       = '--'
+TZFIELD%CCOMMENT   = 'Dates at the middle of the budget timesteps'
+TZFIELD%NGRID      = 0
+TZFIELD%NTYPE      = TYPEDATE
+TZFIELD%NDIMS      = 1
+TZFIELD%LTIMEDEP   = .FALSE.
+
+if ( tpfields(1)%ndimlist(4) /= NMNHDIM_UNKNOWN .and. tpfields(1)%ndimlist(4) /= NMNHDIM_UNUSED) then
+  tzfield%ndimlist(1) = tpfields(1)%ndimlist(4)
+  tzfield%ndimlist(2:) = NMNHDIM_UNUSED
+end if
+
+CALL IO_Field_write( tzfile, tzfield, tpdates(:) )
+
+!Reset ndimlist
+tzfield%ndimlist(:) = NMNHDIM_UNKNOWN
+
 !
-END SUBROUTINE Write_diachro
+
+!Restore id of the file root group ('/' group)
+tzfile%nncid = isavencid
+
+end  subroutine Write_diachro_nc4
+#endif
 
 end module mode_write_diachro
-- 
GitLab