From a5de46e3c6f13db07f2e205896a636e6a39aa683 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 20 Mar 2024 16:32:17 +0100
Subject: [PATCH] Philippe 20/03/2024: add boxes for output files (not yet
 fully implemented)

---
 src/LIB/SURCOUCHE/src/modd_field.f90          |  16 +
 src/LIB/SURCOUCHE/src/modd_io.f90             |   4 +
 src/LIB/SURCOUCHE/src/mode_io_field_write.f90 | 950 ++++++++++--------
 .../SURCOUCHE/src/mode_io_manage_struct.f90   |  70 +-
 src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90   | 229 ++++-
 src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90   |  66 ++
 src/MNH/modd_bakout.f90                       |  15 +-
 src/MNH/modd_outn.f90                         |  19 +-
 src/MNH/modn_output.f90                       |  72 +-
 9 files changed, 989 insertions(+), 452 deletions(-)

diff --git a/src/LIB/SURCOUCHE/src/modd_field.f90 b/src/LIB/SURCOUCHE/src/modd_field.f90
index 9f34d92a2..aabe37203 100644
--- a/src/LIB/SURCOUCHE/src/modd_field.f90
+++ b/src/LIB/SURCOUCHE/src/modd_field.f90
@@ -17,6 +17,7 @@
 !  P. Wautelet 14/10/2021: dynamically allocate tfieldlist (+ reallocate if necessary)
 !  P. Wautelet 04/11/2021: add TFIELDMETADATA type
 !  P. Wautelet 15/02/2024: add time dimension for Lagrangian trajectories
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-----------------------------------------------------------------
 module modd_field
 
@@ -98,6 +99,21 @@ integer, parameter :: NMNHDIM_PAIR                = 43  ! For values coming by p
 
 integer, parameter :: NMNHDIM_LASTDIM_DIACHRO     = 43  ! Index of the last defined dimension for diachronic files
 
+! Box dimensions are after the NMNHDIM_LASTDIM_NODIACHRO/NMNHDIM_LASTDIM_DIACHRO indices
+! because they are allocated in separated structures (1 for each box instead of being in the global dimension list)
+integer, parameter :: NMNHDIM_BOX_FIRST_ENTRY     = 44
+
+integer, parameter :: NMNHDIM_BOX_NI              = NMNHDIM_BOX_FIRST_ENTRY
+integer, parameter :: NMNHDIM_BOX_NJ              = 45
+integer, parameter :: NMNHDIM_BOX_NI_U            = 46
+integer, parameter :: NMNHDIM_BOX_NJ_U            = 47
+integer, parameter :: NMNHDIM_BOX_NI_V            = 48
+integer, parameter :: NMNHDIM_BOX_NJ_V            = 49
+integer, parameter :: NMNHDIM_BOX_LEVEL           = 50
+integer, parameter :: NMNHDIM_BOX_LEVEL_W         = 51
+
+integer, parameter :: NMNHDIM_BOX_LAST_ENTRY      = NMNHDIM_BOX_LEVEL_W
+
 integer, parameter :: NMNHDIM_BUDGET_NGROUPS      = 101 ! This is not a true dimension
 integer, parameter :: NMNHDIM_FLYER_PROC          = 102 ! This is not a true dimension
 integer, parameter :: NMNHDIM_PROFILER_PROC       = 103 ! This is not a true dimension
diff --git a/src/LIB/SURCOUCHE/src/modd_io.f90 b/src/LIB/SURCOUCHE/src/modd_io.f90
index 4a2e25727..f4548dc12 100644
--- a/src/LIB/SURCOUCHE/src/modd_io.f90
+++ b/src/LIB/SURCOUCHE/src/modd_io.f90
@@ -13,6 +13,8 @@
 !  P. Wautelet 22/09/2020: add ldimreduced in tfiledata
 !  P. Wautelet 10/11/2020: new data structures for netCDF dimensions
 !  P. Wautelet 14/12/2023: add lossy compression for output files
+!  P. Wautelet 07/02/2024: add compression for all netCDF files
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-----------------------------------------------------------------
 
 #define MNH_REDUCE_DIMENSIONS_IN_FILES 1
@@ -133,6 +135,8 @@ TYPE TFILEDATA
   INTEGER(KIND=CDFINT)   :: NNCCOMPRESS_LOSSY_ALGO = NF90_QUANTIZE_GRANULARBR ! Lossy compression algorithm
   INTEGER(KIND=CDFINT)   :: NNCCOMPRESS_LOSSY_NSD  = 3                  ! Number of Significant Digits (or Bits)
   TYPE(TDIMSNC), POINTER :: TNCDIMS => NULL()     ! Dimensions of netCDF file
+  INTEGER(KIND=CDFINT), DIMENSION(:), ALLOCATABLE :: NBOXNCID   ! Box HDF group identifiers     (used for MNHOUTPUT files)
+  TYPE(TDIMSNC),        DIMENSION(:), ALLOCATABLE :: TBOXNCDIMS ! Box dimensions of netCDF file (used for MNHOUTPUT files)
 #endif
   !
   !Fields for other files
diff --git a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90 b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
index 25ab1ff7c..e32f00474 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
@@ -21,6 +21,7 @@
 !  P. Wautelet 14/01/2021: add IO_Field_write_byname_N4 and IO_Field_write_byfield_N4 subroutines
 !  P. Wautelet 07/04/2023: correct IO_Field_user_write examples
 !  P. Wautelet 08/01/2024: add zero-size check for monoprocess runs
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-----------------------------------------------------------------
 
 #define MNH_SCALARS_IN_SPLITFILES 0
@@ -3559,7 +3560,7 @@ end subroutine IO_Ndimlist_reduce
     CALL IO_Format_write_select(TPFILE,GLFI,GNC4)
     !
     if ( Present( koffset ) .and. glfi ) then
-      call Print_msg( NVERB_ERROR, 'IO', 'IO_Field_write_box_byfield_X3', Trim( tpfile%cname ) &
+      call Print_msg( NVERB_ERROR, 'IO', 'IO_Field_write_box_byfield_X2', Trim( tpfile%cname ) &
                       // ': LFI format not supported' )
       glfi = .false.
     end if
@@ -3612,7 +3613,8 @@ end subroutine IO_Ndimlist_reduce
   END SUBROUTINE IO_Field_write_box_byfield_X2
 
 
-  SUBROUTINE IO_Field_write_box_byfield_X3( TPFILE, TPFIELD, HBUDGET, PFIELD, KXOBOX, KXEBOX, KYOBOX, KYEBOX, KRESP, koffset )
+  SUBROUTINE IO_Field_write_box_byfield_X3( TPFILE, TPFIELD, HBUDGET, PFIELD, KXOBOX, KXEBOX, KYOBOX, KYEBOX, &
+                                            KZOBOX, KZEBOX, KRESP, koffset )
     !
     USE MODD_IO, ONLY: GSMONOPROC, ISP
     !
@@ -3629,12 +3631,15 @@ end subroutine IO_Ndimlist_reduce
     INTEGER,                         INTENT(IN)  :: KXEBOX   ! Global coordinates of the box
     INTEGER,                         INTENT(IN)  :: KYOBOX   !
     INTEGER,                         INTENT(IN)  :: KYEBOX   !
+    INTEGER,               OPTIONAL, INTENT(IN)  :: KZOBOX   !
+    INTEGER,               OPTIONAL, INTENT(IN)  :: KZEBOX   !
     INTEGER,               OPTIONAL, INTENT(OUT) :: KRESP    ! return-code
     integer, dimension(3), optional, intent(in)  :: koffset
     !
     !*      0.2   Declarations of local variables
     !
     integer                             :: iresp, iresp_lfi, iresp_nc4, iresp_glob
+    INTEGER                             :: IZOBOX, IZEBOX
     REAL, DIMENSION(:,:,:), POINTER     :: ZFIELDP
     LOGICAL                             :: GALLOC
     LOGICAL                             :: GLFI, GNC4
@@ -3659,10 +3664,22 @@ end subroutine IO_Ndimlist_reduce
     end if
 
     IF (IRESP==0) THEN
-       IF (GSMONOPROC) THEN ! sequential execution
+      IF ( PRESENT( KZOBOX ) ) THEN
+        IZOBOX = KZOBOX
+      ELSE
+        IZOBOX = LBOUND( PFIELD, 3 )
+      END IF
+
+      IF ( PRESENT( KZEBOX ) ) THEN
+        IZEBOX = KZEBOX
+      ELSE
+        IZEBOX = UBOUND( PFIELD, 3 )
+      END IF
+
+      IF (GSMONOPROC) THEN ! sequential execution
           IF (HBUDGET /= 'BUDGET') THEN
              ! take the sub-section of PFIELD defined by the box
-             ZFIELDP=>PFIELD(KXOBOX:KXEBOX,KYOBOX:KYEBOX,:)
+             ZFIELDP=>PFIELD(KXOBOX:KXEBOX,KYOBOX:KYEBOX,IZOBOX:IZEBOX)
           ELSE
              ! take the field as a budget
              ZFIELDP=>PFIELD
@@ -3677,14 +3694,14 @@ end subroutine IO_Ndimlist_reduce
        ELSE ! multiprocesses execution
           IF (ISP == TPFILE%NMASTER_RANK)  THEN
              ! Allocate the box
-             ALLOCATE(ZFIELDP(KXEBOX-KXOBOX+1,KYEBOX-KYOBOX+1,SIZE(PFIELD,3)))
+             ALLOCATE( ZFIELDP(KXEBOX-KXOBOX+1, KYEBOX-KYOBOX+1, IZEBOX-IZOBOX+1) )
              GALLOC = .TRUE.
           ELSE
              ALLOCATE(ZFIELDP(0,0,0))
              GALLOC = .TRUE.
           END IF
           !
-          CALL GATHER_XYFIELD(PFIELD,ZFIELDP,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM,&
+          CALL GATHER_XYFIELD(PFIELD(:,:,IZOBOX:IZEBOX),ZFIELDP,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM,&
                & KXOBOX,KXEBOX,KYOBOX,KYEBOX,HBUDGET)
           !
           IF (ISP == TPFILE%NMASTER_RANK)  THEN
@@ -3879,431 +3896,524 @@ end subroutine IO_Ndimlist_reduce
 
 
 SUBROUTINE IO_Fieldlist_write(TPOUTPUT)
-!
-USE MODD_OUT_n,          ONLY: NOUT_FIELDLIST
-!
+
+USE MODD_IO,             ONLY: ISP
+USE MODD_OUT_n,          ONLY: NOUT_FIELDLIST, NOUT_NBOXES, TOUT_BOXES
+USE MODD_PRECISION,      ONLY: CDFINT
+
 USE MODE_MODELN_HANDLER, ONLY: GET_CURRENT_MODEL_INDEX
-!
+
 IMPLICIT NONE
-!
+
 TYPE(TFILEDATA), INTENT(IN) :: TPOUTPUT !Output file
-!
-INTEGER :: IDX
-INTEGER :: IMI
-INTEGER :: JI
+
+INTEGER              :: IMI
+INTEGER              :: JBOX
+INTEGER              :: JI
+INTEGER(KIND=CDFINT) :: IGROUPID_ROOT
+TYPE(TFILEDATA)      :: TZOUTPUT
 !
 IMI = GET_CURRENT_MODEL_INDEX()
 !
 DO JI = 1, SIZE( NOUT_FIELDLIST )
-  IDX = NOUT_FIELDLIST(JI)
-  NDIMS: SELECT CASE (TFIELDLIST(IDX)%NDIMS)
-    !
-    !0D output
-    !
-    CASE (0)
-      NTYPE0D: SELECT CASE (TFIELDLIST(IDX)%NTYPE)
-        !
-        !0D real
-        !
-        CASE (TYPEREAL)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_X0D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X0D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_X0D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X0D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_X0D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 0D real fields' )
-          END IF
-        !
-        !0D integer
-        !
-        CASE (TYPEINT)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_N0D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N0D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_N0D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N0D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_N0D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 0D integer fields' )
-          END IF
-        !
-        !0D logical
-        !
-        CASE (TYPELOG)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_L0D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_L0D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_L0D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_L0D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_L0D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 0D logical fields' )
-          END IF
-        !
-        !0D string
-        !
-        CASE (TYPECHAR)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_C0D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_C0D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_C0D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_C0D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_C0D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 0D character fields' )
-          END IF
-        !
-        !0D date/time
-        !
-        CASE (TYPEDATE)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_T0D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_T0D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_T0D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_T0D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_T0D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 0D date/time fields' )
-          END IF
-        !
-        !0D other types
-        !
-        CASE DEFAULT
-          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': type not yet supported for 0D output' )
-      END SELECT NTYPE0D
-    !
-    !1D output
-    !
-    CASE (1)
-      NTYPE1D: SELECT CASE (TFIELDLIST(IDX)%NTYPE)
-        !
-        !1D real
-        !
-        CASE (TYPEREAL)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_X1D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X1D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_X1D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X1D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_X1D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 1D real fields' )
-          END IF
-        !
-        !1D integer
-        !
-        CASE (TYPEINT)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_N1D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N1D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_N1D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N1D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_N1D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 1D integer fields' )
-          END IF
-        !
-        !1D logical
-        !
-        CASE (TYPELOG)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_L1D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_L1D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_L1D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_L1D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_L1D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 1D logical fields' )
-          END IF
-        !
-        !1D string
-        !
-        CASE (TYPECHAR)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_C1D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_C1D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_C1D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_C1D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_C1D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 1D character fields' )
-          END IF
-        !
-        !1D date/time
-        !
-        CASE (TYPEDATE)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_T1D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_T1D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_T1D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_T1D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_T1D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 1D date/time fields' )
-          END IF
-        !
-        !1D other types
-        !
-        CASE DEFAULT
-          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': type not yet supported for 1D output' )
-      END SELECT NTYPE1D
-    !
-    !2D output
-    !
-    CASE (2)
-      NTYPE2D: SELECT CASE (TFIELDLIST(IDX)%NTYPE)
-        !
-        !2D real
-        !
-        CASE (TYPEREAL)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_X2D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X2D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_X2D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X2D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_X2D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 2D real fields' )
-          END IF
-        !
-        !2D integer
-        !
-        CASE (TYPEINT)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_N2D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N2D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_N2D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N2D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_N2D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not allowed for 2D integer fields' )
-          END IF
-        !
-        !2D other types
-        !
-        CASE DEFAULT
-          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': type not yet supported for 2D output' )
-      END SELECT NTYPE2D
-    !
-    !3D output
-    !
-    CASE (3)
-      NTYPE3D: SELECT CASE (TFIELDLIST(IDX)%NTYPE)
-        !
-        !3D real
-        !
-        CASE (TYPEREAL)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_X3D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X3D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_X3D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X3D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_X3D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not (yet) allowed for 3D real fields' )
-            !CALL IO_Field_write_lb(TPOUTPUT,TFIELDLIST(IDX),***,TFIELDLIST(IDX)%TFIELD_X3D(IMI)%DATA)
-          END IF
-        !
-        !3D integer
-        !
-        CASE (TYPEINT)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_N3D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N3D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_N3D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_N3D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_N3D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not (yet) allowed for 3D integer fields' )
-            !CALL IO_Field_write_lb(TPOUTPUT,TFIELDLIST(IDX),***,TFIELDLIST(IDX)%TFIELD_N3D(IMI)%DATA)
-          END IF
-        !
-        !3D other types
-        !
-        CASE DEFAULT
-          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': type not yet supported for 3D output' )
-      END SELECT NTYPE3D
-    !
-    !4D output
-    !
-    CASE (4)
-      NTYPE4D: SELECT CASE (TFIELDLIST(IDX)%NTYPE)
-        !
-        !4D real
-        !
-        CASE (TYPEREAL)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_X4D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X4D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_X4D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X4D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_X4D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not (yet) allowed for 4D real fields' )
-            !CALL IO_Field_write_lb(TPOUTPUT,TFIELDLIST(IDX),***,TFIELDLIST(IDX)%TFIELD_X4D(IMI)%DATA)
-          END IF
-        !
-        !4D other types
-        !
-        CASE DEFAULT
-          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': type not yet supported for 4D output' )
-      END SELECT NTYPE4D
-    !
-    !5D output
-    !
-    CASE (5)
-      NTYPE5D: SELECT CASE (TFIELDLIST(IDX)%NTYPE)
-        !
-        !5D real
-        !
-        CASE (TYPEREAL)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_X5D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X5D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_X5D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X5D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_X5D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not (yet) allowed for 5D real fields' )
-            !CALL IO_Field_write_lb(TPOUTPUT,TFIELDLIST(IDX),***,TFIELDLIST(IDX)%TFIELD_X5D(IMI)%DATA)
-          END IF
-        !
-        !5D other types
-        !
-        CASE DEFAULT
-          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': type not yet supported for 5D output' )
-      END SELECT NTYPE5D
-    !
-    !6D output
-    !
-    CASE (6)
-      NTYPE6D: SELECT CASE (TFIELDLIST(IDX)%NTYPE)
-        !
-        !6D real
-        !
-        CASE (TYPEREAL)
-          IF ( .NOT.ALLOCATED(TFIELDLIST(IDX)%TFIELD_X6D) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X6D is NOT allocated ' )
-          END IF
-          IF ( .NOT.ASSOCIATED(TFIELDLIST(IDX)%TFIELD_X6D(IMI)%DATA) ) THEN
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': TFIELD_X6D%DATA is NOT associated' )
-          END IF
-          IF ( TFIELDLIST(IDX)%CLBTYPE == 'NONE' ) THEN
-            CALL IO_Field_write(TPOUTPUT,TFIELDLIST(IDX),TFIELDLIST(IDX)%TFIELD_X6D(IMI)%DATA)
-          ELSE
-            call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                            ': CLBTYPE/=NONE not (yet) allowed for 6D real fields' )
-            !CALL IO_Field_write_lb(TPOUTPUT,TFIELDLIST(IDX),***,TFIELDLIST(IDX)%TFIELD_X6D(IMI)%DATA)
-          END IF
-        !
-        !6D other types
-        !
-        CASE DEFAULT
-          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': type not yet supported for 6D output' )
-      END SELECT NTYPE6D
-    !
-    !Other number of dimensions
-    !
-    CASE DEFAULT
-      call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', trim(tfieldlist(idx)%cmnhname)// &
-                          ': number of dimensions not yet supported' )
-  END SELECT NDIMS
+  CALL IO_Fieldlist_1field_write( TPOUTPUT, IMI, TFIELDLIST(NOUT_FIELDLIST(JI)) )
 END DO
-!
+
+! Treat boxes (subdomains)
+IF ( NOUT_NBOXES > 0 ) THEN
+  ! Only available for netCDF files
+  TZOUTPUT = TPOUTPUT
+  IF ( TZOUTPUT%CFORMAT == 'LFI') THEN
+    RETURN
+  ELSE IF ( TZOUTPUT%CFORMAT == 'LFICDF4') THEN
+    TZOUTPUT%CFORMAT = 'NETCDF4'
+  ELSE IF ( TZOUTPUT%CFORMAT /= 'NETCDF4' ) THEN
+    CALL PRINT_MSG( NVERB_ERROR, 'IO', 'IO_Fieldlist_write', &
+                    TRIM(TZOUTPUT%CNAME) // ': unknown fileformat: ' // TRIM(TZOUTPUT%CFORMAT) )
+    RETURN
+  END IF
+
+  IGROUPID_ROOT = TZOUTPUT%NNCID
+
+  DO JBOX = 1, NOUT_NBOXES
+    ! Go to the group
+    IF ( ISP == TZOUTPUT%NMASTER_RANK ) TZOUTPUT%NNCID = TZOUTPUT%NBOXNCID(JBOX)
+
+    ! Write fields common to all boxes
+    DO JI = 1, SIZE( NOUT_FIELDLIST )
+      CALL IO_Fieldlist_1field_write( TZOUTPUT, IMI, TFIELDLIST(NOUT_FIELDLIST(JI)), TOUT_BOXES(JBOX) )
+    END DO
+
+      ! Restore the root group (not really necessary but cleaner)
+    IF ( ISP == TZOUTPUT%NMASTER_RANK ) TZOUTPUT%NNCID = IGROUPID_ROOT
+  END DO
+END IF
+
 END SUBROUTINE IO_Fieldlist_write
 
 
+SUBROUTINE IO_Fieldlist_1field_write( TPOUTPUT, KMI, TPFIELD, TPBOX )
+
+USE MODD_FIELD,      ONLY: NMNHDIM_BOX_NI, NMNHDIM_BOX_NJ, NMNHDIM_BOX_NI_U, NMNHDIM_BOX_NJ_U, NMNHDIM_BOX_NI_V, &
+                           NMNHDIM_BOX_NJ_V, NMNHDIM_BOX_LEVEL, NMNHDIM_BOX_LEVEL_W, NMNHDIM_NOTLISTED,          &
+                           TFIELDDATA, TFIELDMETADATA
+USE MODD_OUT_n,      ONLY: TOUTBOXMETADATA
+USE MODD_PARAMETERS, ONLY: JPHEXT, JPVEXT
+
+TYPE(TFILEDATA),                 INTENT(IN) :: TPOUTPUT !Output file
+INTEGER,                         INTENT(IN) :: KMI
+TYPE(TFIELDDATA),                INTENT(IN) :: TPFIELD
+TYPE(TOUTBOXMETADATA), OPTIONAL, INTENT(IN) :: TPBOX
+
+TYPE(TFIELDMETADATA) :: TZFIELDMD
+
+IF ( PRESENT(TPBOX) ) THEN
+  IF ( TPFIELD%NDIMS /= 3 .AND. TPFIELD%NTYPE /= TYPEREAL ) THEN
+    CALL Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', 'TPBOX optional dummy argument not allowed for field ' &
+                    // TRIM(TPFIELD%CMNHNAME) )
+    RETURN
+  END IF
+END IF
+
+NDIMS: SELECT CASE (TPFIELD%NDIMS)
+  !
+  !0D output
+  !
+  CASE (0)
+    NTYPE0D: SELECT CASE (TPFIELD%NTYPE)
+      !
+      !0D real
+      !
+      CASE (TYPEREAL)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_X0D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X0D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_X0D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X0D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_X0D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 0D real fields' )
+        END IF
+      !
+      !0D integer
+      !
+      CASE (TYPEINT)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_N0D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N0D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_N0D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N0D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_N0D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 0D integer fields' )
+        END IF
+      !
+      !0D logical
+      !
+      CASE (TYPELOG)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_L0D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_L0D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_L0D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_L0D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_L0D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 0D logical fields' )
+        END IF
+      !
+      !0D string
+      !
+      CASE (TYPECHAR)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_C0D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_C0D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_C0D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_C0D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_C0D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 0D character fields' )
+        END IF
+      !
+      !0D date/time
+      !
+      CASE (TYPEDATE)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_T0D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_T0D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_T0D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_T0D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_T0D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 0D date/time fields' )
+        END IF
+      !
+      !0D other types
+      !
+      CASE DEFAULT
+        call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': type not yet supported for 0D output' )
+    END SELECT NTYPE0D
+  !
+  !1D output
+  !
+  CASE (1)
+    NTYPE1D: SELECT CASE (TPFIELD%NTYPE)
+      !
+      !1D real
+      !
+      CASE (TYPEREAL)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_X1D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X1D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_X1D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X1D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_X1D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 1D real fields' )
+        END IF
+      !
+      !1D integer
+      !
+      CASE (TYPEINT)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_N1D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N1D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_N1D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N1D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_N1D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 1D integer fields' )
+        END IF
+      !
+      !1D logical
+      !
+      CASE (TYPELOG)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_L1D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_L1D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_L1D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_L1D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_L1D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 1D logical fields' )
+        END IF
+      !
+      !1D string
+      !
+      CASE (TYPECHAR)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_C1D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_C1D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_C1D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_C1D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_C1D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 1D character fields' )
+        END IF
+      !
+      !1D date/time
+      !
+      CASE (TYPEDATE)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_T1D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_T1D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_T1D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_T1D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_T1D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 1D date/time fields' )
+        END IF
+      !
+      !1D other types
+      !
+      CASE DEFAULT
+        call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': type not yet supported for 1D output' )
+    END SELECT NTYPE1D
+  !
+  !2D output
+  !
+  CASE (2)
+    NTYPE2D: SELECT CASE (TPFIELD%NTYPE)
+      !
+      !2D real
+      !
+      CASE (TYPEREAL)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_X2D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X2D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_X2D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X2D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_X2D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 2D real fields' )
+        END IF
+      !
+      !2D integer
+      !
+      CASE (TYPEINT)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_N2D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N2D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_N2D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N2D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_N2D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not allowed for 2D integer fields' )
+        END IF
+      !
+      !2D other types
+      !
+      CASE DEFAULT
+        call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': type not yet supported for 2D output' )
+    END SELECT NTYPE2D
+  !
+  !3D output
+  !
+  CASE (3)
+    NTYPE3D: SELECT CASE (TPFIELD%NTYPE)
+      !
+      !3D real
+      !
+      CASE (TYPEREAL)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_X3D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X3D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_X3D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X3D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          IF ( .NOT.PRESENT( TPBOX ) ) THEN
+            CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_X3D(KMI)%DATA)
+          ELSE
+            TZFIELDMD = TFIELDMETADATA( TPFIELD ) !Copy only metadata (TZFIELD is of TYPE(TFIELDMETADATA))
+            IF ( TZFIELDMD%CDIR /= 'XY' ) &
+              CALL Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)//': must have CDIR="XY"' )
+            SELECT CASE ( TZFIELDMD%NGRID )
+              CASE( 1, 4 )
+                TZFIELDMD%NDIMLIST(1) = NMNHDIM_BOX_NI
+              CASE( 2 )
+                TZFIELDMD%NDIMLIST(1) = NMNHDIM_BOX_NI_U
+              CASE( 3 )
+                TZFIELDMD%NDIMLIST(1) = NMNHDIM_BOX_NI_V
+            END SELECT
+            SELECT CASE ( TZFIELDMD%NGRID )
+              CASE( 1, 4 )
+                TZFIELDMD%NDIMLIST(2) = NMNHDIM_BOX_NJ
+              CASE( 2 )
+                TZFIELDMD%NDIMLIST(2) = NMNHDIM_BOX_NJ_U
+              CASE( 3 )
+                TZFIELDMD%NDIMLIST(2) = NMNHDIM_BOX_NJ_V
+            END SELECT
+            SELECT CASE ( TZFIELDMD%NGRID )
+              CASE( 1, 2, 3 )
+                TZFIELDMD%NDIMLIST(3) = NMNHDIM_BOX_LEVEL
+              CASE( 4 )
+                TZFIELDMD%NDIMLIST(3) = NMNHDIM_BOX_LEVEL_W
+            END SELECT
+            CALL IO_Field_write_BOX( TPOUTPUT, TZFIELDMD, 'OTHER', TPFIELD%TFIELD_X3D(KMI)%DATA,   &
+                                     KXOBOX = JPHEXT + TPBOX%NIINF, KXEBOX = JPHEXT + TPBOX%NISUP, &
+                                     KYOBOX = JPHEXT + TPBOX%NJINF, KYEBOX = JPHEXT + TPBOX%NJSUP, &
+                                     KZOBOX = JPVEXT + TPBOX%NKINF, KZEBOX = JPVEXT + TPBOX%NKSUP  )
+          END IF
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not (yet) allowed for 3D real fields' )
+          !CALL IO_Field_write_lb(TPOUTPUT,TPFIELD,***,TPFIELD%TFIELD_X3D(KMI)%DATA)
+        END IF
+      !
+      !3D integer
+      !
+      CASE (TYPEINT)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_N3D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N3D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_N3D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_N3D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_N3D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not (yet) allowed for 3D integer fields' )
+          !CALL IO_Field_write_lb(TPOUTPUT,TPFIELD,***,TPFIELD%TFIELD_N3D(KMI)%DATA)
+        END IF
+      !
+      !3D other types
+      !
+      CASE DEFAULT
+        call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': type not yet supported for 3D output' )
+    END SELECT NTYPE3D
+  !
+  !4D output
+  !
+  CASE (4)
+    NTYPE4D: SELECT CASE (TPFIELD%NTYPE)
+      !
+      !4D real
+      !
+      CASE (TYPEREAL)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_X4D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X4D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_X4D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X4D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_X4D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not (yet) allowed for 4D real fields' )
+          !CALL IO_Field_write_lb(TPOUTPUT,TPFIELD,***,TPFIELD%TFIELD_X4D(KMI)%DATA)
+        END IF
+      !
+      !4D other types
+      !
+      CASE DEFAULT
+        call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': type not yet supported for 4D output' )
+    END SELECT NTYPE4D
+  !
+  !5D output
+  !
+  CASE (5)
+    NTYPE5D: SELECT CASE (TPFIELD%NTYPE)
+      !
+      !5D real
+      !
+      CASE (TYPEREAL)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_X5D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X5D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_X5D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X5D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_X5D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not (yet) allowed for 5D real fields' )
+          !CALL IO_Field_write_lb(TPOUTPUT,TPFIELD,***,TPFIELD%TFIELD_X5D(KMI)%DATA)
+        END IF
+      !
+      !5D other types
+      !
+      CASE DEFAULT
+        call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': type not yet supported for 5D output' )
+    END SELECT NTYPE5D
+  !
+  !6D output
+  !
+  CASE (6)
+    NTYPE6D: SELECT CASE (TPFIELD%NTYPE)
+      !
+      !6D real
+      !
+      CASE (TYPEREAL)
+        IF ( .NOT.ALLOCATED(TPFIELD%TFIELD_X6D) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X6D is NOT allocated ' )
+        END IF
+        IF ( .NOT.ASSOCIATED(TPFIELD%TFIELD_X6D(KMI)%DATA) ) THEN
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': TFIELD_X6D%DATA is NOT associated' )
+        END IF
+        IF ( TPFIELD%CLBTYPE == 'NONE' ) THEN
+          CALL IO_Field_write(TPOUTPUT,TPFIELD,TPFIELD%TFIELD_X6D(KMI)%DATA)
+        ELSE
+          call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                          ': CLBTYPE/=NONE not (yet) allowed for 6D real fields' )
+          !CALL IO_Field_write_lb(TPOUTPUT,TPFIELD,***,TPFIELD%TFIELD_X6D(KMI)%DATA)
+        END IF
+      !
+      !6D other types
+      !
+      CASE DEFAULT
+        call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': type not yet supported for 6D output' )
+    END SELECT NTYPE6D
+  !
+  !Other number of dimensions
+  !
+  CASE DEFAULT
+    call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
+                        ': number of dimensions not yet supported' )
+END SELECT NDIMS
+
+END SUBROUTINE IO_Fieldlist_1field_write
+
+
 SUBROUTINE IO_Field_user_write( TPOUTPUT )
 !
 #if 0
diff --git a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90 b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
index cf21d21c4..e4745932f 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
@@ -26,6 +26,7 @@
 !  P. Wautelet 17/01/2024: add IO_File_remove_from_list subroutine
 !  P. Wautelet 02/02/2024: restructure backups/outputs lists
 !  P. Wautelet 07/02/2024: add compression for all netCDF files
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-----------------------------------------------------------------
 MODULE MODE_IO_MANAGE_STRUCT
 !
@@ -61,7 +62,7 @@ USE MODD_DYN_n,       ONLY: DYN_MODEL
 USE MODD_FIELD,       ONLY: TFIELDLIST
 USE MODD_NESTING,     ONLY: NDAD
 USE MODD_SUB_MODEL_N, ONLY: NFILE_BACKUP_CURRENT, NFILE_OUTPUT_CURRENT, SUB_MODEL_MODEL
-USE MODD_OUT_n,       ONLY: OUT_GOTO_MODEL, OUT_MODEL
+USE MODD_OUT_n,       ONLY: OUT_GOTO_MODEL, OUT_MODEL, TOUT_BOXES
 USE MODD_VAR_ll,      ONLY: IP
 
 use mode_field, only: Find_field_id_from_mnhname
@@ -83,6 +84,7 @@ INTEGER           :: ISTEP
 INTEGER           :: ISTEP_MAX        ! Number of timesteps
 INTEGER           :: IPOS,IFIELD      ! Indices
 INTEGER           :: JOUT,IDX         ! Loop indices
+INTEGER           :: JI
 INTEGER           :: IRESP
 INTEGER           :: ISTEPDADFIRST
 INTEGER, DIMENSION(:,:), ALLOCATABLE :: IBAK_STEPLIST
@@ -523,6 +525,9 @@ DO IMI = 1, NMODEL
     END IF
   END IF
   !
+  ! Treat the boxes (sub-domains) for output files
+  IF ( OUT_MODEL(IMI)%NOUT_NUMB > 0 ) CALL IO_BOX_PREPARE( IMI )
+  !
   IF ( IP == 1 ) THEN
     ! Backup information
     WRITE( *, '( "-------------------------------------------------" )' )
@@ -578,6 +583,13 @@ DO IMI = 1, NMODEL
           WRITE(*, '( "  ", A )' ) TRIM(TFIELDLIST(IDX)%CMNHNAME)
         END DO
       END IF
+      WRITE( *, '( "Number of boxes:", I9 )' ) OUT_MODEL(IMI)%NOUT_NBOXES
+      IF ( OUT_MODEL(IMI)%NOUT_NBOXES > 0 ) THEN
+        WRITE( *, '( "  List of boxes:" )' )
+        DO JI = 1, OUT_MODEL(IMI)%NOUT_NBOXES
+          WRITE(*, '( "    ", A )' ) TRIM(OUT_MODEL(IMI)%TOUT_BOXES(JI)%CNAME)
+        END DO
+      END IF
     END IF
 
     WRITE( *, '( "-------------------------------------------------" )' )
@@ -722,6 +734,62 @@ SUBROUTINE SORT_ENTRIES(KNUMB,KSTEPS)
   END DO
 END SUBROUTINE SORT_ENTRIES
 !
+!#########################################################################
+SUBROUTINE IO_BOX_PREPARE( KMI )
+!#########################################################################
+
+  USE MODD_DIM_n, ONLY: NIMAX_ll, NJMAX_ll, NKMAX
+
+  INTEGER, INTENT(IN) :: KMI
+
+  OUT_MODEL(IMI)%NOUT_NBOXES = NOUT_BOXES(IMI)
+  ALLOCATE( OUT_MODEL(IMI)%TOUT_BOXES(NOUT_BOXES(IMI)) )
+  TOUT_BOXES => OUT_MODEL(IMI)%TOUT_BOXES
+  DO JI = 1, OUT_MODEL(IMI)%NOUT_NBOXES
+    IF ( LEN_TRIM(COUT_BOX_NAME(IMI,JI)) > 0 ) THEN
+      TOUT_BOXES(JI)%CNAME = COUT_BOX_NAME(IMI,JI)
+    ELSE
+      WRITE( TOUT_BOXES(JI)%CNAME, '( "Box_", I4.4 )' ) JI
+    END IF
+    TOUT_BOXES(JI)%NIINF = NOUT_BOX_IINF(IMI,JI)
+    TOUT_BOXES(JI)%NISUP = NOUT_BOX_ISUP(IMI,JI)
+    TOUT_BOXES(JI)%NJINF = NOUT_BOX_JINF(IMI,JI)
+    TOUT_BOXES(JI)%NJSUP = NOUT_BOX_JSUP(IMI,JI)
+    TOUT_BOXES(JI)%NKINF = NOUT_BOX_KINF(IMI,JI)
+    TOUT_BOXES(JI)%NKSUP = NOUT_BOX_KSUP(IMI,JI)
+
+    !Set default values to physical domain boundaries
+    IF ( TOUT_BOXES(JI)%NIINF == NNEGUNDEF ) TOUT_BOXES(JI)%NIINF = 1
+    IF ( TOUT_BOXES(JI)%NISUP == NNEGUNDEF ) TOUT_BOXES(JI)%NIINF = NIMAX_ll
+    IF ( TOUT_BOXES(JI)%NJINF == NNEGUNDEF ) TOUT_BOXES(JI)%NJINF = 1
+    IF ( TOUT_BOXES(JI)%NJSUP == NNEGUNDEF ) TOUT_BOXES(JI)%NJINF = NJMAX_ll
+    IF ( TOUT_BOXES(JI)%NKINF == NNEGUNDEF ) TOUT_BOXES(JI)%NKINF = 1
+    IF ( TOUT_BOXES(JI)%NKSUP == NNEGUNDEF ) TOUT_BOXES(JI)%NKINF = NKMAX
+
+    !Check that selected indices are in physical domain
+    IF ( TOUT_BOXES(JI)%NIINF < 1 )        CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_IINF too small (<1)' )
+    IF ( TOUT_BOXES(JI)%NIINF > NIMAX_ll ) CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_IINF too large (>NIMAX)')
+    IF ( TOUT_BOXES(JI)%NISUP < 1 )        CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_ISUP too small (<1)' )
+    IF ( TOUT_BOXES(JI)%NISUP > NIMAX_ll ) CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_ISUP too large (>NIMAX)')
+    IF ( TOUT_BOXES(JI)%NISUP < TOUT_BOXES(JI)%NIINF ) &
+                                           CALL Print_msg( NVERB_ERROR, 'GEN', 'INI_LES_n', 'NOUT_BOX_ISUP < NOUT_BOX_IINF' )
+
+    IF ( TOUT_BOXES(JI)%NJINF < 1 )        CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_JINF too small (<1)' )
+    IF ( TOUT_BOXES(JI)%NJINF > NJMAX_ll ) CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_JINF too large (>NJMAX)')
+    IF ( TOUT_BOXES(JI)%NJSUP < 1 )        CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_JSUP too small (<1)' )
+    IF ( TOUT_BOXES(JI)%NJSUP > NJMAX_ll ) CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_JSUP too large (>NJMAX)')
+    IF ( TOUT_BOXES(JI)%NJSUP < TOUT_BOXES(JI)%NJINF ) &
+                                           CALL Print_msg( NVERB_ERROR, 'GEN', 'INI_LES_n', 'NOUT_BOX_JSUP < NOUT_BOX_JINF' )
+
+    IF ( TOUT_BOXES(JI)%NKINF < 1 )     CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_KINF too small (<1)' )
+    IF ( TOUT_BOXES(JI)%NKINF > NKMAX ) CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_KINF too large (>NKMAX)' )
+    IF ( TOUT_BOXES(JI)%NKSUP < 1 )     CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_KSUP too small (<1)' )
+    IF ( TOUT_BOXES(JI)%NKSUP > NKMAX ) CALL Print_msg( NVERB_ERROR, 'GEN', 'IO_BOX_PREPARE', 'NOUT_BOX_KSUP too large (>NKMAX)' )
+    IF ( TOUT_BOXES(JI)%NKSUP < TOUT_BOXES(JI)%NKINF ) &
+                                           CALL Print_msg( NVERB_ERROR, 'GEN', 'INI_LES_n', 'NOUT_BOX_KSUP < NOUT_BOX_KINF' )
+  END DO
+END SUBROUTINE IO_BOX_PREPARE
+
 END SUBROUTINE IO_Bakout_struct_prepare
 
 
diff --git a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
index 0efd6aaa2..001b20ad6 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
@@ -21,6 +21,7 @@
 !  P. Wautelet 27/05/2021: improve IO_Mnhname_clean to autocorrect names to be CF compliant
 !  P. Wautelet 08/10/2021: add 2 new dimensions: LW_bands and SW_bands
 !  P. Wautelet 15/02/2024: add time dimension for Lagrangian trajectories
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-----------------------------------------------------------------
 #ifdef MNH_IOCDF4
 module mode_io_tools_nc4
@@ -270,12 +271,16 @@ use modd_field,         only: NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_NI_U, NMNHDIM_NJ_U
                               NMNHDIM_SERIES_LEVEL, NMNHDIM_SERIES_LEVEL_W,                                     &
                               NMNHDIM_SERIES_TIME, NMNHDIM_PROFILER_TIME, NMNHDIM_STATION_TIME,                 &
                               NMNHDIM_PAIR,                                                                     &
+                              NMNHDIM_BOX_FIRST_ENTRY, NMNHDIM_BOX_LAST_ENTRY,                                  &
+                              NMNHDIM_BOX_NI, NMNHDIM_BOX_NJ, NMNHDIM_BOX_NI_U, NMNHDIM_BOX_NJ_U,               &
+                              NMNHDIM_BOX_NI_V, NMNHDIM_BOX_NJ_V, NMNHDIM_BOX_LEVEL, NMNHDIM_BOX_LEVEL_W,       &
                               NMNHDIM_ARAKAWA,                                                                  &
                               NMNHDIM_LASTDIM_NODIACHRO, NMNHDIM_LASTDIM_DIACHRO
 
 use modd_les,           only: lles_pdf, nles_k, npdf, nspectra_k, xles_temp_mean_start, xles_temp_mean_step, xles_temp_mean_end
 use modd_les_n,         only: nles_times, nspectra_ni, nspectra_nj
 use modd_nsv,           only: nsv
+use modd_out_n,         only: nout_nboxes, tout_boxes
 USE MODD_PARAMETERS_ll, ONLY: JPHEXT, JPVEXT
 use modd_param_n,       only: crad
 use modd_profiler_n,    only: lprofiler, tprofilers_time
@@ -292,6 +297,7 @@ CHARACTER(LEN=:),ALLOCATABLE :: YPROGRAM
 integer                      :: iavg, iprof, istation
 integer                      :: ispectra_ni, ispectra_nj
 INTEGER                      :: IIU_ll, IJU_ll, IKU
+integer                      :: jbox
 
 CALL PRINT_MSG(NVERB_DEBUG,'IO','IO_Knowndims_set_nc4','called for '//TRIM(TPFILE%CNAME))
 
@@ -352,6 +358,41 @@ if ( tpfile%ctype == 'MNHDIAG' .and. Trim( yprogram ) == 'DIAG' ) then
   if ( ltraj .and. ntrajstlg > 0 ) call IO_Add_dim_nc4( tpfile, NMNHDIM_TRAJ_TIME, 'time_traj', ntrajstlg )
 end if
 
+!Write dimensions used in boxes (subdomains) for MNHOUTPUT files
+if ( tpfile%ctype == 'MNHOUTPUT' ) then
+  !Allocate arrays
+  Allocate( tpfile%nboxncid  (nout_nboxes) )
+  Allocate( tpfile%tboxncdims(nout_nboxes) )
+
+  !Loop on the boxes
+  do jbox = 1, nout_nboxes
+    !Create the HDF group for each box
+    tpfile%nboxncid(jbox) = IO_Box_group_create_nc4( tpfile, tout_boxes(jbox) )
+
+    !Allocate dimension list
+    tpfile%tboxncdims(jbox)%nmaxdims = NMNHDIM_BOX_LAST_ENTRY - NMNHDIM_BOX_FIRST_ENTRY + 1
+    Allocate( tpfile%tboxncdims(jbox)%tdims(NMNHDIM_BOX_FIRST_ENTRY:NMNHDIM_BOX_LAST_ENTRY) )
+
+    !Write the box dimensions
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_NI,      'box_ni',      tout_boxes(jbox)%nisup-tout_boxes(jbox)%niinf+1 )
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_NI_U,    'box_ni_u',    tout_boxes(jbox)%nisup-tout_boxes(jbox)%niinf+1 )
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_NI_V,    'box_ni_v',    tout_boxes(jbox)%nisup-tout_boxes(jbox)%niinf+1 )
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_NJ,      'box_nj',      tout_boxes(jbox)%njsup-tout_boxes(jbox)%njinf+1 )
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_NJ_U,    'box_nj_u',    tout_boxes(jbox)%njsup-tout_boxes(jbox)%njinf+1 )
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_NJ_V,    'box_nj_v',    tout_boxes(jbox)%njsup-tout_boxes(jbox)%njinf+1 )
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_LEVEL,   'box_level',   tout_boxes(jbox)%nksup-tout_boxes(jbox)%nkinf+1 )
+    call IO_Add_dim_box_nc4( tpfile, jbox, NMNHDIM_BOX_LEVEL_W, 'box_level_w', tout_boxes(jbox)%nksup-tout_boxes(jbox)%nkinf+1 )
+
+    !Write the box attributes
+    call IO_box_attribute_write_nc4( tpfile, jbox, 'min_I_index_in_physical_domain', tout_boxes(jbox)%niinf )
+    call IO_box_attribute_write_nc4( tpfile, jbox, 'max_I_index_in_physical_domain', tout_boxes(jbox)%nisup )
+    call IO_box_attribute_write_nc4( tpfile, jbox, 'min_J_index_in_physical_domain', tout_boxes(jbox)%njinf )
+    call IO_box_attribute_write_nc4( tpfile, jbox, 'max_J_index_in_physical_domain', tout_boxes(jbox)%njsup )
+    call IO_box_attribute_write_nc4( tpfile, jbox, 'min_K_index_in_physical_domain', tout_boxes(jbox)%nkinf )
+    call IO_box_attribute_write_nc4( tpfile, jbox, 'max_K_index_in_physical_domain', tout_boxes(jbox)%nksup )
+  end do
+end if
+
 !Write dimensions used in diachronic files
 if ( tpfile%ctype == 'MNHDIACHRONIC' ) then
   !Dimension of size 2 used for NMNHDIM_COMPLEX
@@ -485,6 +526,96 @@ if ( istatus /= NF90_NOERR ) &
 end subroutine IO_Add_dim_nc4
 
 
+function IO_Box_group_create_nc4( tpfile, tpbox ) result( kgroupid )
+
+  use modd_parameters, only: NMNHNAMELGTMAX
+  use modd_precision,  only: CDFINT
+  use modd_out_n,      only: toutboxmetadata
+
+  use NETCDF,          only: NF90_DEF_GRP, NF90_INQ_NCID, NF90_NOERR
+
+  type(tfiledata),       intent(in) :: tpfile
+  type(toutboxmetadata), intent(in) :: tpbox
+  integer(kind=CDFINT)              :: kgroupid
+
+  character(len=nmnhnamelgtmax) :: ygroupname
+  integer(kind=CDFINT)          :: istatus
+
+  CALL IO_Mnhname_clean( tpbox%cname, ygroupname )
+  istatus = NF90_INQ_NCID( tpfile%nncid, ygroupname, kgroupid )
+  if ( istatus == NF90_NOERR ) then
+    call Print_msg( NVERB_ERROR, 'IO', 'IO_Box_group_create_nc4', Trim(tpfile%cname) // &
+                    ': group ' // Trim(ygroupname) // ' already exists', olocal=.true. )
+  else
+    istatus = NF90_DEF_GRP( tpfile%nncid, ygroupname, kgroupid )
+    if ( istatus /= NF90_NOERR ) &
+      call IO_Err_handle_nc4( istatus, 'IO_Box_group_create_nc4', 'NF90_DEF_GRP', 'for ' // trim( ygroupname ) )
+  end if
+
+end function IO_Box_group_create_nc4
+
+
+subroutine IO_Add_dim_box_nc4( tpfile, kbox, kidx, hdimname, klen )
+  use modd_field, only: NMNHDIM_BOX_FIRST_ENTRY, NMNHDIM_BOX_LAST_ENTRY
+
+  use NETCDF, only: NF90_DEF_DIM, NF90_NOERR
+
+  type(tfiledata),  intent(inout) :: tpfile
+  integer,          intent(in)    :: kbox     !Index of the box
+  integer,          intent(in)    :: kidx     !Position of the dimension in the list
+  character(len=*), intent(in)    :: hdimname !Name of the dimension
+  integer,          intent(in)    :: klen     !Length of the dimension
+
+  character(len=Len(hdimname))  :: ydimname_clean
+  integer(kind=CDFINT)          :: istatus
+
+  call IO_Mnhname_clean( hdimname, ydimname_clean )
+
+  if ( kidx < NMNHDIM_BOX_FIRST_ENTRY .or. kidx > NMNHDIM_BOX_LAST_ENTRY)                                                 &
+    call Print_msg( NVERB_FATAL, 'IO', 'IO_Add_dim_box_nc4', 'index out of range for dimension ' // Trim( ydimname_clean ) // &
+                    ' of file ' //Trim( tpfile%cname ) )
+
+  if ( tpfile%tboxncdims(kbox)%tdims(kidx)%nlen /= -1 .or. tpfile%tboxncdims(kbox)%tdims(kidx)%nid /= -1 ) &
+    call Print_msg( NVERB_WARNING, 'IO', 'IO_Add_dim_box_nc4', 'dimension ' // Trim( ydimname_clean ) //   &
+                    ' already defined for file ' //Trim( tpfile%cname ) )
+
+  tpfile%tboxncdims(kbox)%tdims(kidx)%cname = ydimname_clean
+  tpfile%tboxncdims(kbox)%tdims(kidx)%nlen  = Int( klen, kind = CDFINT )
+
+  istatus = NF90_DEF_DIM( tpfile%nboxncid(kbox), Trim( ydimname_clean ), Int( klen, kind = CDFINT ), &
+                          tpfile%tboxncdims(kbox)%tdims(kidx)%nid )
+  if ( istatus /= NF90_NOERR ) &
+    call IO_Err_handle_nc4( istatus, 'IO_Add_dim_nc4', 'NF90_DEF_DIM', Trim( ydimname_clean ) )
+
+end subroutine IO_Add_dim_box_nc4
+
+
+subroutine IO_box_attribute_write_nc4( tpfile, kbox, hattname, kdata )
+  use NETCDF,            only: NF90_INQUIRE_ATTRIBUTE, NF90_PUT_ATT, NF90_GLOBAL, NF90_NOERR
+
+  use modd_precision,    only: CDFINT
+
+  type(tfiledata),  intent(in) :: tpfile
+  integer,          intent(in) :: kbox     !Index of the box
+  character(len=*), intent(in) :: hattname
+  integer,          intent(in) :: kdata
+
+  character(len=Len(hattname)) :: yattname
+  integer(kind=CDFINT)         :: istatus
+
+  call IO_Mnhname_clean( hattname, yattname )
+
+  istatus = NF90_INQUIRE_ATTRIBUTE( tpfile%nboxncid(kbox), NF90_GLOBAL, yattname )
+  if (istatus == NF90_NOERR ) &
+    call Print_msg( NVERB_ERROR, 'IO', 'IO_box_attribute_write_nc4', 'attribute ' // yattname // ' already exists' )
+
+  istatus = NF90_PUT_ATT( tpfile%nboxncid(kbox), NF90_GLOBAL, yattname, kdata )
+  if (istatus /= NF90_NOERR ) &
+   call IO_Err_handle_nc4( istatus, 'IO_box_attribute_write_nc4', 'NF90_PUT_ATT', Trim( yattname ) )
+
+end subroutine IO_box_attribute_write_nc4
+
+
 SUBROUTINE IO_Iocdf_dealloc_nc4(tpdimsnc)
 TYPE(tdimsnc),  POINTER :: tpdimsnc
 
@@ -501,12 +632,12 @@ END SUBROUTINE IO_Iocdf_dealloc_nc4
 
 SUBROUTINE IO_Vdims_fill_nc4(TPFILE, TPFIELD, KSHAPE, KVDIMS)
 
-use NETCDF, only: NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION
-
 use modd_conf,   only: l2d, lpack
 use modd_field,  only: NMNHDIM_UNKNOWN, NMNHDIM_ONE, NMNHDIM_COMPLEX,                                   &
                        NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_NI_U, NMNHDIM_NJ_U, NMNHDIM_NI_V, NMNHDIM_NJ_V,  &
                        NMNHDIM_LEVEL, NMNHDIM_LEVEL_W, NMNHDIM_TIME,                                    &
+                       NMNHDIM_BOX_NI, NMNHDIM_BOX_NJ, NMNHDIM_BOX_NI_U, NMNHDIM_BOX_NJ_U,              &
+                       NMNHDIM_BOX_NI_V, NMNHDIM_BOX_NJ_V, NMNHDIM_BOX_LEVEL, NMNHDIM_BOX_LEVEL_W,      &
                        NMNHDIM_BUDGET_CART_NI,      NMNHDIM_BUDGET_CART_NJ,  NMNHDIM_BUDGET_CART_NI_U,  &
                        NMNHDIM_BUDGET_CART_NJ_U,    NMNHDIM_BUDGET_CART_NI_V, NMNHDIM_BUDGET_CART_NJ_V, &
                        NMNHDIM_BUDGET_CART_LEVEL,   NMNHDIM_BUDGET_CART_LEVEL_W,                        &
@@ -535,8 +666,6 @@ INTEGER                       :: IGRID
 integer                       :: iidx
 integer                       :: iresp
 INTEGER                       :: JI
-integer(kind=CDFINT)          :: ilen
-integer(kind=CDFINT)          :: istatus
 !
 CALL PRINT_MSG(NVERB_DEBUG,'IO','IO_Vdims_fill_nc4','called for '//TRIM(TPFIELD%CMNHNAME))
 !
@@ -589,18 +718,31 @@ if ( Any( tpfield%ndimlist(:) /= NMNHDIM_UNKNOWN ) ) then
     end if
 
     if ( tpfield%ndimlist(ji) == NMNHDIM_FLYER_TIME ) then
-      istatus = NF90_INQ_DIMID( tpfile%nncid, 'time_flyer', kvdims(ji) )
-      if ( istatus /= NF90_NOERR ) &
-        call IO_Err_handle_nc4( istatus, 'IO_Vdims_fill_nc4','NF90_INQ_DIMID', Trim( tpfield%cmnhname ) )
-
-      istatus = NF90_INQUIRE_DIMENSION( tpfile%nncid, kvdims(ji), len = ilen)
-      if ( istatus /= NF90_NOERR ) &
-        call IO_Err_handle_nc4( istatus, 'IO_Vdims_fill_nc4','NF90_INQUIRE_DIMENSION', Trim( tpfield%cmnhname ) )
+      call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'time_flyer', tpfield%cmnhname, kshape(ji), kvdims(ji) )
+      cycle
+    end if
 
-      if ( kshape(ji) /= ilen ) then
-        call Print_msg( NVERB_FATAL, 'IO', 'IO_Vdims_fill_nc4', &
-                                     'wrong size for dimension '// 'time_flyer' // ' of field ' // Trim( tpfield%cmnhname ) )
-      end if
+    if ( tpfield%ndimlist(ji) >= NMNHDIM_BOX_NI  .AND. tpfield%ndimlist(ji) <= NMNHDIM_BOX_LEVEL_W) then
+      select case ( tpfield%ndimlist(ji) )
+        case ( NMNHDIM_BOX_NI )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_ni',      tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case ( NMNHDIM_BOX_NJ )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_nj',      tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case ( NMNHDIM_BOX_NI_U )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_ni_u',    tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case ( NMNHDIM_BOX_NJ_U )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_nj_u',    tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case ( NMNHDIM_BOX_NI_V )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_ni_v',    tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case ( NMNHDIM_BOX_NJ_V )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_nj_v',    tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case ( NMNHDIM_BOX_LEVEL )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_level',   tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case ( NMNHDIM_BOX_LEVEL_W )
+          call IO_Dim_localgroup_check( tpfile%nncid, tpfield%ndimlist(ji), 'box_level_w', tpfield%cmnhname, kshape(ji), kvdims(ji))
+        case default
+          call Print_msg( NVERB_FATAL, 'IO', 'IO_Vdims_fill_nc4', 'invalid box dimension' )
+      end select
       cycle
     end if
 
@@ -758,6 +900,63 @@ end if
 end subroutine IO_Dim_find_create_nc4
 
 
+subroutine IO_Dim_localgroup_check( kgroupid, kmnhdim, hdimname, hvarname, kdimlen, kdimid )
+!Subroutine to check dimension inside a local HDF group (dimension is not global)
+use NETCDF, only: NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_NOERR
+
+integer(kind=CDFINT), intent(in)  :: kgroupid
+integer,              intent(in)  :: kmnhdim
+character(len=*),     intent(in)  :: hdimname
+character(len=*),     intent(in)  :: hvarname
+integer,              intent(in)  :: kdimlen
+integer(kind=CDFINT), intent(out) :: kdimid
+
+character(len=:),    allocatable :: ygroupname
+integer(kind=CDFINT)             :: ilen
+integer(kind=CDFINT)             :: istatus
+
+istatus = NF90_INQ_DIMID( kgroupid, hdimname, kdimid )
+if ( istatus /= NF90_NOERR ) then
+  ygroupname = Get_groupname()
+  call IO_Err_handle_nc4( istatus, 'IO_Dim_localgroup_check','NF90_INQ_DIMID', &
+                          Trim(ygroupname) // ':' // Trim(hvarname) // ':' // Trim(hdimname) )
+end if
+
+istatus = NF90_INQUIRE_DIMENSION( kgroupid, kdimid, len = ilen)
+if ( istatus /= NF90_NOERR ) then
+  ygroupname = Get_groupname()
+  call IO_Err_handle_nc4( istatus, 'IO_Dim_localgroup_check', 'NF90_INQUIRE_DIMENSION', &
+                          Trim(ygroupname) // ':' // Trim(hvarname) // ':' // Trim(hdimname) )
+end if
+
+if ( kdimlen /= ilen ) then
+  ygroupname = Get_groupname()
+  call Print_msg( NVERB_FATAL, 'IO', 'IO_Dim_localgroup_check', &
+                  'wrong size for dimension '// Trim(hdimname) // ' of field ' // Trim(ygroupname) // ':' // Trim(hvarname) )
+end if
+
+contains
+function Get_groupname()
+  use NETCDF, only: NF90_INQ_GRPNAME_LEN, NF90_INQ_GRPNAME_FULL, NF90_NOERR
+
+  character(len=:), allocatable :: Get_groupname
+
+  integer(kind=CDFINT) :: igrplen
+
+  istatus = NF90_INQ_GRPNAME_LEN( kgroupid, igrplen )
+  if ( istatus /= NF90_NOERR ) then
+    Get_groupname = 'unknown_group'
+    return
+  end if
+
+  allocate( character(len=igrplen) :: Get_groupname )
+  istatus = NF90_INQ_GRPNAME_FULL( kgroupid, igrplen, Get_groupname )
+end function Get_groupname
+
+
+end subroutine IO_Dim_localgroup_check
+
+
 FUNCTION IO_Strdimid_get_nc4(TPFILE,KLEN)
 use modd_netcdf, only: tdimnc
 
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index ee99fe61a..fbf37f88f 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -34,6 +34,7 @@
 !  P. Wautelet 13/01/2023: IO_Coordvar_write_nc4: add optional dummy argument TPDTMODELN to force written model time
 !  P. Wautelet 14/12/2023: add lossy compression for output files
 !  P. Wautelet 15/02/2024: IO_Coordvar_write_nc4: add time dimension for Lagrangian trajectories
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-----------------------------------------------------------------
 #ifdef MNH_IOCDF4
 module mode_io_write_nc4
@@ -1465,6 +1466,7 @@ use modd_les,        only: cles_level_type, cspectra_level_type, nlesn_iinf, nle
 use modd_les_n,      only: nles_dtcount, nles_mean_end, nles_mean_start, nles_mean_step, nles_mean_times, &
                            nles_times, nspectra_ni, nspectra_nj, tles_dates, xles_times
 use modd_netcdf,     only: tdimnc
+use modd_out_n,      only: nout_nboxes, tout_boxes
 use modd_parameters, only: jphext, JPVEXT
 use modd_profiler_n, only: lprofiler, tprofilers_time
 use modd_series,     only: lseries
@@ -1699,6 +1701,15 @@ if ( tzfile%lmaster ) then
   end if
 end if
 
+! Write coordinates for boxes (subdomains) in OUTPUT files
+if ( tzfile%lmaster ) then
+  if ( tzfile%ctype == 'MNHOUTPUT' .and. nout_nboxes > 0 ) then
+    do ji = 1, nout_nboxes
+      call Write_box_coords( ji, tout_boxes(ji) )
+    end do
+  end if
+end if
+
 if ( tzfile%lmaster ) then
   !Write coordinates used in diachronic files
   if ( tzfile%ctype == 'MNHDIACHRONIC' ) then
@@ -2369,6 +2380,61 @@ subroutine Write_flyer_time_coord( tpflyer )
 
 end subroutine Write_flyer_time_coord
 
+
+subroutine Write_box_coords( kbox, tpbox )
+  use modd_field, only: NMNHDIM_BOX_NI,   NMNHDIM_BOX_NJ,   NMNHDIM_BOX_NI_U,  NMNHDIM_BOX_NJ_U,  &
+                        NMNHDIM_BOX_NI_V, NMNHDIM_BOX_NJ_V, NMNHDIM_BOX_LEVEL, NMNHDIM_BOX_LEVEL_W
+
+  use modd_out_n, only: toutboxmetadata
+
+  integer,               intent(in) :: kbox
+  type(toutboxmetadata), intent(in) :: tpbox
+
+  integer(kind=CDFINT) :: incid_root
+
+  ! Go to the box HDF group
+  incid_root = tzfile%nncid
+  tzfile%nncid = tzfile%nboxncid(kbox)
+
+  incid = tzfile%nboxncid(kbox)
+
+  call Write_hor_coord1d( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_NI),   'x-dimension of the box', &
+                          trim(ystdnameprefix)//'_x_coordinate',               'X', 0.,   0, 0,      &
+                          zxhatm_glob(tpbox%niinf + jphext : tpbox%nisup + jphext)                   )
+
+  call Write_hor_coord1d( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_NI_U), 'x-dimension of the box at u location', &
+                          trim(ystdnameprefix)//'_x_coordinate_at_u_location', 'X', -0.5, 0, 0,                    &
+                          zxhat_glob (tpbox%niinf + jphext : tpbox%nisup + jphext)                                 )
+
+  call Write_hor_coord1d( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_NI_V), 'x-dimension of the box at v location', &
+                          trim(ystdnameprefix)//'_x_coordinate_at_v_location', 'X', 0.,   0, 0,                    &
+                          zxhatm_glob(tpbox%niinf + jphext : tpbox%nisup + jphext)                                 )
+
+  call Write_hor_coord1d( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_NJ),   'y-dimension of the box', &
+                          trim(ystdnameprefix)//'_y_coordinate',               'Y', 0.,   0, 0,      &
+                          zyhatm_glob(tpbox%njinf + jphext : tpbox%njsup + jphext)                   )
+
+  call Write_hor_coord1d( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_NJ_U), 'y-dimension of the box at u location', &
+                          trim(ystdnameprefix)//'_y_coordinate_at_u_location', 'Y', 0.,   0, 0,                    &
+                          zyhatm_glob(tpbox%njinf + jphext : tpbox%njsup + jphext)                                 )
+
+  call Write_hor_coord1d( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_NJ_V), 'y-dimension of the box at v location', &
+                          trim(ystdnameprefix)//'_y_coordinate_at_v_location', 'Y', -0.5, 0, 0,                    &
+                          zyhat_glob (tpbox%njinf + jphext : tpbox%njsup + jphext)                                )
+
+  call Write_ver_coord( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_LEVEL),                                             &
+                        'position z in the transformed space of the box',                                             &
+                        '', 'altitude',               0.,   0, 0, zzhatm(tpbox%nkinf + JPVEXT : tpbox%nksup + JPVEXT) )
+
+  call Write_ver_coord( tzfile%tboxncdims(kbox)%tdims(NMNHDIM_BOX_LEVEL_W),                                           &
+                        'position z in the transformed space at w location of the box',                               &
+                        '', 'altitude_at_w_location', -0.5, 0, 0, zzhat (tpbox%nkinf + JPVEXT : tpbox%nksup + JPVEXT) )
+
+  ! Restore the root HDF group
+  tzfile%nncid = incid_root
+
+end subroutine Write_box_coords
+
 END SUBROUTINE IO_Coordvar_write_nc4
 
 
diff --git a/src/MNH/modd_bakout.f90 b/src/MNH/modd_bakout.f90
index 84e1eab13..98032cbc2 100644
--- a/src/MNH/modd_bakout.f90
+++ b/src/MNH/modd_bakout.f90
@@ -37,6 +37,7 @@
 !  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 14/12/2023: add lossy compression for output files
 !  P. Wautelet 07/02/2024: add compression for backup files
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -49,6 +50,7 @@ IMPLICIT NONE
 SAVE
 !
 INTEGER, PARAMETER :: NCOMPRNAMELGTMAX  = 10 ! Maximum length of compression algorithm name
+INTEGER, PARAMETER :: NOUT_BOXMAX       = 20 ! Maximum number of boxes (subdomains) for each model output
 !
 LOGICAL :: LBAK_BEG = .FALSE. ! Force a backup/output at the first timestep
 LOGICAL :: LOUT_BEG = .FALSE. ! of the segment for all models
@@ -89,5 +91,16 @@ CHARACTER(LEN=NMNHNAMELGTMAX), ALLOCATABLE, DIMENSION(:,:) :: COUT_VAR ! Name of
 !
 !Directory names for backups/outputs
 CHARACTER(LEN=NDIRNAMELGTMAX) :: CBAK_DIR='', COUT_DIR=''
-!
+
+! Boxes (subdomains) for outputs
+INTEGER, DIMENSION(JPMODELMAX) :: NOUT_BOXES = 0 ! Number of sub-boxes inside each modelgrid
+CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:,:), ALLOCATABLE :: COUT_BOX_NAME ! Names of the boxes
+
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_IINF ! Box coordinates in physical domain (for each model and for each box)
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_ISUP
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_JINF
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_JSUP
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_KINF
+INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_KSUP
+
 END MODULE MODD_BAKOUT
diff --git a/src/MNH/modd_outn.f90 b/src/MNH/modd_outn.f90
index b48907563..2e3057515 100644
--- a/src/MNH/modd_outn.f90
+++ b/src/MNH/modd_outn.f90
@@ -33,21 +33,33 @@
 !!      Original    20/10/94                      
 !  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 02/02/2024: restructure backups/outputs lists
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
 !             ------------
 !
 !
-USE MODD_PARAMETERS, ONLY: JPMODELMAX, NNEGUNDEF
+USE MODD_PARAMETERS, ONLY: JPMODELMAX, NMNHNAMELGTMAX, NNEGUNDEF
 
 IMPLICIT NONE
 
 SAVE
 
+TYPE TOUTBOXMETADATA
+  CHARACTER(LEN=NMNHNAMELGTMAX):: CNAME = '' ! Name of the box
+  INTEGER :: NIINF = NNEGUNDEF ! Box coordinates in physical domain
+  INTEGER :: NISUP = NNEGUNDEF
+  INTEGER :: NJINF = NNEGUNDEF
+  INTEGER :: NJSUP = NNEGUNDEF
+  INTEGER :: NKINF = NNEGUNDEF
+  INTEGER :: NKSUP = NNEGUNDEF
+END TYPE TOUTBOXMETADATA
+
 TYPE OUT_t
   INTEGER                            :: NBAK_NUMB = 0 ! Number of backups performed by model n
   INTEGER                            :: NOUT_NUMB = 0 ! Number of outputs performed by model n
+  INTEGER                            :: NOUT_NBOXES = 0 ! Number of boxes (sub-domains) to write in output files
   INTEGER                            :: NBAK_PREV_STEP = 1 ! Timestep number of the previous backup
   INTEGER                            :: NBAK_STEPFREQ      = NNEGUNDEF ! Time steps between 2 backups
   INTEGER                            :: NOUT_STEPFREQ      = NNEGUNDEF ! Time steps between 2 outputs
@@ -56,12 +68,14 @@ TYPE OUT_t
   INTEGER, DIMENSION(:), ALLOCATABLE :: NBAK_STEPLIST                  ! List of time steps to do backups (except regular series)
   INTEGER, DIMENSION(:), ALLOCATABLE :: NOUT_STEPLIST                  ! List of time steps to do outputs (except regular series)
   INTEGER, DIMENSION(:), ALLOCATABLE :: NOUT_FIELDLIST                 ! List of fields to write in outputs
+  TYPE(TOUTBOXMETADATA), DIMENSION(:), ALLOCATABLE :: TOUT_BOXES
 END TYPE OUT_t
 
 TYPE(OUT_t), DIMENSION(JPMODELMAX), TARGET :: OUT_MODEL
 
 INTEGER,               POINTER :: NBAK_NUMB          => NULL()
 INTEGER,               POINTER :: NOUT_NUMB          => NULL()
+INTEGER,               POINTER :: NOUT_NBOXES        => NULL()
 INTEGER,               POINTER :: NBAK_PREV_STEP     => NULL()
 INTEGER,               POINTER :: NBAK_STEPFREQ      => NULL()
 INTEGER,               POINTER :: NOUT_STEPFREQ      => NULL()
@@ -70,6 +84,7 @@ INTEGER,               POINTER :: NOUT_STEPFREQFIRST => NULL()
 INTEGER, DIMENSION(:), POINTER :: NBAK_STEPLIST      => NULL()
 INTEGER, DIMENSION(:), POINTER :: NOUT_STEPLIST      => NULL()
 INTEGER, DIMENSION(:), POINTER :: NOUT_FIELDLIST     => NULL()
+TYPE(TOUTBOXMETADATA), DIMENSION(:), POINTER :: TOUT_BOXES  => NULL()
 
 CONTAINS
 
@@ -79,6 +94,7 @@ INTEGER, INTENT(IN) :: KFROM, KTO
 ! Current model is set to model KTO
 NBAK_NUMB          => OUT_MODEL(KTO)%NBAK_NUMB
 NOUT_NUMB          => OUT_MODEL(KTO)%NOUT_NUMB
+NOUT_NBOXES        => OUT_MODEL(KTO)%NOUT_NBOXES
 NBAK_PREV_STEP     => OUT_MODEL(KTO)%NBAK_PREV_STEP
 NBAK_STEPFREQ      => OUT_MODEL(KTO)%NBAK_STEPFREQ
 NOUT_STEPFREQ      => OUT_MODEL(KTO)%NOUT_STEPFREQ
@@ -87,6 +103,7 @@ NOUT_STEPFREQFIRST => OUT_MODEL(KTO)%NOUT_STEPFREQFIRST
 NBAK_STEPLIST      => OUT_MODEL(KTO)%NBAK_STEPLIST
 NOUT_STEPLIST      => OUT_MODEL(KTO)%NOUT_STEPLIST
 NOUT_FIELDLIST     => OUT_MODEL(KTO)%NOUT_FIELDLIST
+TOUT_BOXES         => OUT_MODEL(KTO)%TOUT_BOXES
 
 END SUBROUTINE OUT_GOTO_MODEL
 
diff --git a/src/MNH/modn_output.f90 b/src/MNH/modn_output.f90
index df3fc7f41..5766aafd8 100644
--- a/src/MNH/modn_output.f90
+++ b/src/MNH/modn_output.f90
@@ -35,6 +35,7 @@
 !  P. Wautelet       2016: new structures for outputs/backups
 !  P. Wautelet 02/10/2017: split NAM_OUTPUT in NAM_BACKUP and NAM_OUTPUT
 !  P. Wautelet 14/12/2023: add lossy compression for output files
+!  P. Wautelet 20/03/2024: add boxes for output files
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -52,8 +53,12 @@ NAMELIST/NAM_OUTPUT/LOUT_BEG,LOUT_END,&
                    LOUT_REDUCE_FLOAT_PRECISION, &
                    LOUT_COMPRESS, NOUT_COMPRESS_LEVEL,&
                    LOUT_COMPRESS_LOSSY, COUT_COMPRESS_LOSSY_ALGO, NOUT_COMPRESS_LOSSY_NSD, &
-                   COUT_DIR
-!
+                   COUT_DIR, &
+                   NOUT_BOXES, COUT_BOX_NAME, NOUT_BOX_IINF, NOUT_BOX_ISUP, NOUT_BOX_JINF, NOUT_BOX_JSUP, &
+                   NOUT_BOX_KINF, NOUT_BOX_KSUP
+
+LOGICAL, SAVE, PRIVATE :: LOUTPUT_NML_ALLOCATED = .FALSE.
+
 CONTAINS
 
 SUBROUTINE OUTPUT_NML_ALLOCATE( )
@@ -61,30 +66,69 @@ SUBROUTINE OUTPUT_NML_ALLOCATE( )
   USE MODD_IO,         ONLY: NFILE_NUM_MAX
   USE MODD_PARAMETERS, ONLY: NNEGUNDEF, XNEGUNDEF
 
-  IF ( .NOT.ALLOCATED(XOUT_TIME) ) THEN
-    ALLOCATE( XOUT_TIME(NMODEL, NFILE_NUM_MAX) )
-    XOUT_TIME(:,:) = XNEGUNDEF
-  END IF
+  USE MODE_MSG
 
-  IF ( .NOT.ALLOCATED(NOUT_STEP) ) THEN
-    ALLOCATE( NOUT_STEP(NMODEL, NFILE_NUM_MAX) )
-    NOUT_STEP(:,:) = NNEGUNDEF
-  END IF
+  INTEGER :: JI
 
-  IF ( .NOT.ALLOCATED(COUT_VAR) ) THEN
-    ALLOCATE( COUT_VAR(NMODEL, JPOUTVARMAX) )
-    COUT_VAR(:,:) = ''
-  END IF
+  IF ( LOUTPUT_NML_ALLOCATED ) THEN
+      CALL PRINT_MSG( NVERB_DEBUG, 'GEN', 'OUTPUT_NML_ALLOCATE', 'was altready allocated' )
+      RETURN
+    END IF
+
+  ALLOCATE( XOUT_TIME(NMODEL, NFILE_NUM_MAX) )
+  ALLOCATE( NOUT_STEP(NMODEL, NFILE_NUM_MAX) )
+  ALLOCATE( COUT_VAR(NMODEL, JPOUTVARMAX) )
+
+  ALLOCATE( COUT_BOX_NAME(NMODEL, NOUT_BOXMAX) )
+
+  ALLOCATE( NOUT_BOX_IINF(NMODEL, NOUT_BOXMAX) )
+  ALLOCATE( NOUT_BOX_ISUP(NMODEL, NOUT_BOXMAX) )
+  ALLOCATE( NOUT_BOX_JINF(NMODEL, NOUT_BOXMAX) )
+  ALLOCATE( NOUT_BOX_JSUP(NMODEL, NOUT_BOXMAX) )
+  ALLOCATE( NOUT_BOX_KINF(NMODEL, NOUT_BOXMAX) )
+  ALLOCATE( NOUT_BOX_KSUP(NMODEL, NOUT_BOXMAX) )
+
+  XOUT_TIME(:,:) = XNEGUNDEF
+  NOUT_STEP(:,:) = NNEGUNDEF
+  COUT_VAR(:,:)  = ''
+
+  COUT_BOX_NAME(:,:) = ''
+
+  NOUT_BOX_IINF(:,:) = NNEGUNDEF
+  NOUT_BOX_ISUP(:,:) = NNEGUNDEF
+  NOUT_BOX_JINF(:,:) = NNEGUNDEF
+  NOUT_BOX_JSUP(:,:) = NNEGUNDEF
+  NOUT_BOX_KINF(:,:) = NNEGUNDEF
+  NOUT_BOX_KSUP(:,:) = NNEGUNDEF
+
+  LOUTPUT_NML_ALLOCATED = .TRUE.
 
 END SUBROUTINE OUTPUT_NML_ALLOCATE
 
 
 SUBROUTINE OUTPUT_NML_DEALLOCATE( )
+  USE MODE_MSG
+
+  IF ( .NOT. LOUTPUT_NML_ALLOCATED ) THEN
+    CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'OUTPUT_NML_DEALLOCATE', 'was not cleanly allocated' )
+    RETURN
+  END IF
 
   DEALLOCATE( XOUT_TIME )
   DEALLOCATE( NOUT_STEP )
   DEALLOCATE( COUT_VAR  )
 
+  DEALLOCATE( COUT_BOX_NAME )
+
+  DEALLOCATE( NOUT_BOX_IINF )
+  DEALLOCATE( NOUT_BOX_ISUP )
+  DEALLOCATE( NOUT_BOX_JINF )
+  DEALLOCATE( NOUT_BOX_JSUP )
+  DEALLOCATE( NOUT_BOX_KINF )
+  DEALLOCATE( NOUT_BOX_KSUP )
+
+  LOUTPUT_NML_ALLOCATED = .FALSE.
+
 END SUBROUTINE OUTPUT_NML_DEALLOCATE
 
 END MODULE MODN_OUTPUT
-- 
GitLab