diff --git a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90 b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
index ea5ff4f74f7d7c0afb42ae47599d31f862ba8fac..7cda6cf62c7a6cf6e61b7ade694a93aaa338c94d 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
@@ -4931,40 +4931,116 @@ END SUBROUTINE IO_Fieldlist_1field_write
 SUBROUTINE IO_Box_coords_enable( TPFIELD, TPBOX )
   ! This subroutine set to true a coordinate as soon as it is used
   ! Once it is true, it must not be set again to false (it is to know if the coordinate has been used and has to be written)
-  USE MODD_FIELD, ONLY: TFIELDMETADATA
+  USE MODD_FIELD, ONLY: NMNHDIM_BOX_LEVEL, NMNHDIM_BOX_LEVEL_W, NMNHDIM_LEVEL, NMNHDIM_LEVEL_W, TFIELDMETADATA
   USE MODD_OUT_n, ONLY: TOUTBOXMETADATA
 
   CLASS(TFIELDMETADATA),  INTENT(IN)    :: TPFIELD
   CLASS(TOUTBOXMETADATA), INTENT(INOUT) :: TPBOX
 
-  IF ( TPFIELD%NDIMS < 2 ) RETURN
-
-  SELECT CASE( TPFIELD%NGRID )
-    CASE(1) !Mass point
-      TPBOX%LCOORD_LATM = .TRUE.
-      TPBOX%LCOORD_LONM = .TRUE.
-    CASE(2) !u point
-      TPBOX%LCOORD_LATU = .TRUE.
-      TPBOX%LCOORD_LONU = .TRUE.
-    CASE(3) !v point
-      TPBOX%LCOORD_LATV = .TRUE.
-      TPBOX%LCOORD_LONV = .TRUE.
-    CASE(4) !w point
-      TPBOX%LCOORD_LATM = .TRUE.
-      TPBOX%LCOORD_LONM = .TRUE.
-    CASE(5) !xi vorticity point (=f point =uv point)
-      TPBOX%LCOORD_LATF = .TRUE.
-      TPBOX%LCOORD_LONF = .TRUE.
-    CASE(6) !eta vorticity point (=uw point)
-      TPBOX%LCOORD_LATU = .TRUE.
-      TPBOX%LCOORD_LONU = .TRUE.
-    CASE(7) !zeta vorticity point (=vw point)
-      TPBOX%LCOORD_LATV = .TRUE.
-      TPBOX%LCOORD_LONV = .TRUE.
-    CASE(8) !fw point (=uvw point)
-      TPBOX%LCOORD_LATF = .TRUE.
-      TPBOX%LCOORD_LONF = .TRUE.
-  END SELECT
+  IF ( TPFIELD%NDIMS == 1 ) THEN
+    SELECT CASE ( TPFIELD%CDIR )
+      CASE('XX')
+        SELECT CASE( TPFIELD%NGRID )
+          CASE(1) !Mass point
+            TPBOX%LCOORD_NIM = .TRUE.
+          CASE(2) !u point
+            TPBOX%LCOORD_NIU = .TRUE.
+          CASE(3) !v point
+            TPBOX%LCOORD_NIV = .TRUE.
+          CASE(4) !w point
+            TPBOX%LCOORD_NIM = .TRUE.
+          CASE(5) !xi vorticity point (=f point =uv point)
+            TPBOX%LCOORD_NIU = .TRUE.
+          CASE(6) !eta vorticity point (=uw point)
+            TPBOX%LCOORD_NIU = .TRUE.
+          CASE(7) !zeta vorticity point (=vw point)
+            TPBOX%LCOORD_NIV = .TRUE.
+          CASE(8) !fw point (=uvw point)
+            TPBOX%LCOORD_NIU = .TRUE.
+        END SELECT
+      CASE('YY')
+        SELECT CASE( TPFIELD%NGRID )
+          CASE(1) !Mass point
+            TPBOX%LCOORD_NJM = .TRUE.
+          CASE(2) !u point
+            TPBOX%LCOORD_NJU = .TRUE.
+          CASE(3) !v point
+            TPBOX%LCOORD_NJV = .TRUE.
+          CASE(4) !w point
+            TPBOX%LCOORD_NJM = .TRUE.
+          CASE(5) !xi vorticity point (=f point =uv point)
+            TPBOX%LCOORD_NJV = .TRUE.
+          CASE(6) !eta vorticity point (=uw point)
+            TPBOX%LCOORD_NJU = .TRUE.
+          CASE(7) !zeta vorticity point (=vw point)
+            TPBOX%LCOORD_NJV = .TRUE.
+          CASE(8) !fw point (=uvw point)
+            TPBOX%LCOORD_NJV = .TRUE.
+        END SELECT
+      CASE('ZZ')
+        SELECT CASE( TPFIELD%NGRID )
+          CASE(1, 2, 3, 5)
+            TPBOX%LCOORD_LEVELM = .TRUE.
+          CASE(4, 6, 7, 8)
+            TPBOX%LCOORD_LEVELW = .TRUE.
+        END SELECT
+    END SELECT
+  ELSE IF ( TPFIELD%NDIMS >= 2 ) THEN
+    SELECT CASE( TPFIELD%NGRID )
+      CASE(1) !Mass point
+        TPBOX%LCOORD_NIM  = .TRUE.
+        TPBOX%LCOORD_NJM  = .TRUE.
+        TPBOX%LCOORD_LATM = .TRUE.
+        TPBOX%LCOORD_LONM = .TRUE.
+      CASE(2) !u point
+        TPBOX%LCOORD_NIU  = .TRUE.
+        TPBOX%LCOORD_NJU  = .TRUE.
+        TPBOX%LCOORD_LATU = .TRUE.
+        TPBOX%LCOORD_LONU = .TRUE.
+      CASE(3) !v point
+        TPBOX%LCOORD_NIV  = .TRUE.
+        TPBOX%LCOORD_NJV  = .TRUE.
+        TPBOX%LCOORD_LATV = .TRUE.
+        TPBOX%LCOORD_LONV = .TRUE.
+      CASE(4) !w point
+        TPBOX%LCOORD_NIM  = .TRUE.
+        TPBOX%LCOORD_NJM  = .TRUE.
+        TPBOX%LCOORD_LATM = .TRUE.
+        TPBOX%LCOORD_LONM = .TRUE.
+      CASE(5) !xi vorticity point (=f point =uv point)
+        TPBOX%LCOORD_NIU  = .TRUE.
+        TPBOX%LCOORD_NJV  = .TRUE.
+        TPBOX%LCOORD_LATF = .TRUE.
+        TPBOX%LCOORD_LONF = .TRUE.
+      CASE(6) !eta vorticity point (=uw point)
+        TPBOX%LCOORD_NIU  = .TRUE.
+        TPBOX%LCOORD_NJU  = .TRUE.
+        TPBOX%LCOORD_LATU = .TRUE.
+        TPBOX%LCOORD_LONU = .TRUE.
+      CASE(7) !zeta vorticity point (=vw point)
+        TPBOX%LCOORD_NIV  = .TRUE.
+        TPBOX%LCOORD_NJV  = .TRUE.
+        TPBOX%LCOORD_LATV = .TRUE.
+        TPBOX%LCOORD_LONV = .TRUE.
+      CASE(8) !fw point (=uvw point)
+        TPBOX%LCOORD_NIU  = .TRUE.
+        TPBOX%LCOORD_NJV  = .TRUE.
+        TPBOX%LCOORD_LATF = .TRUE.
+        TPBOX%LCOORD_LONF = .TRUE.
+    END SELECT
+
+    ! Treat vertical dimension (special case 1D ZZ already treated)
+    ! It is assumed that the vertical dimension is the 3rd one
+    IF ( TPFIELD%NDIMS >= 3 ) THEN
+      IF ( TPFIELD%NDIMLIST(3) == NMNHDIM_LEVEL .OR. TPFIELD%NDIMLIST(3) == NMNHDIM_BOX_LEVEL ) THEN
+        TPBOX%LCOORD_LEVELM = .TRUE.
+      ELSE IF ( TPFIELD%NDIMLIST(3) == NMNHDIM_LEVEL_W .OR. TPFIELD%NDIMLIST(3) == NMNHDIM_BOX_LEVEL_W ) THEN
+        TPBOX%LCOORD_LEVELW = .TRUE.
+      ELSE
+        call Print_msg( NVERB_WARNING, 'IO', 'IO_Box_coords_enable', trim(tpfield%cmnhname) // ': incorrect vertical dimension' )
+      END IF
+    END IF
+  END IF
 END SUBROUTINE IO_Box_coords_enable
 
 
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index fa5aabd08d4d5e56250dd521cb8b5b285288dbcf..2e9a95b4a925d1c190c6bf021865b6aeee9e8af7 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -1817,7 +1817,7 @@ else
 end if
 
 
-if ( tzfile%ctype /= 'MNHOUTPUT' .or. lout_bigbox_write ) then !Do not write coordinates of main domain if it is not written
+if ( tzfile%ctype /= 'MNHOUTPUT' ) then !Do not write coordinates of main domain if it is not written
   call Write_hor_coord1d( tzdim_ni,   'x-dimension of the grid', &
                           trim(ystdnameprefix)//'_x_coordinate',               'X', 0.,   iim_bdmin, iim_bdmax, &
                           zxhatm_glob(iim_min:iim_max) )
@@ -1836,6 +1836,31 @@ if ( tzfile%ctype /= 'MNHOUTPUT' .or. lout_bigbox_write ) then !Do not write coo
   call Write_hor_coord1d( tzdim_nj_v, 'y-dimension of the grid at v location', &
                           trim(ystdnameprefix)//'_y_coordinate_at_v_location', 'Y', -0.5, ijv_bdmin, ijv_bdmax, &
                           zyhat_glob(ijv_min:ijv_max)  )
+else if ( tzfile%ctype == 'MNHOUTPUT' .and. lout_bigbox_write ) then !Do not write coordinates of main domain if it is not written
+  if ( tout_boxes(0)%lcoord_nim ) &
+    call Write_hor_coord1d( tzdim_ni,   'x-dimension of the grid', &
+                            trim(ystdnameprefix)//'_x_coordinate',               'X', 0.,   iim_bdmin, iim_bdmax, &
+                          zxhatm_glob(iim_min:iim_max) )
+  if ( tout_boxes(0)%lcoord_njm ) &
+    call Write_hor_coord1d( tzdim_nj,   'y-dimension of the grid', &
+                            trim(ystdnameprefix)//'_y_coordinate',               'Y', 0.,   ijm_bdmin, ijm_bdmax, &
+                            zyhatm_glob(ijm_min:ijm_max) )
+  if ( tout_boxes(0)%lcoord_niu ) &
+    call Write_hor_coord1d( tzdim_ni_u, 'x-dimension of the grid at u location', &
+                            trim(ystdnameprefix)//'_x_coordinate_at_u_location', 'X', -0.5, iiu_bdmin, iiu_bdmax, &
+                            zxhat_glob(iiu_min:iiu_max)  )
+  if ( tout_boxes(0)%lcoord_nju ) &
+    call Write_hor_coord1d( tzdim_nj_u, 'y-dimension of the grid at u location', &
+                            trim(ystdnameprefix)//'_y_coordinate_at_u_location', 'Y', 0.,   iju_bdmin, iju_bdmax, &
+                            zyhatm_glob(iju_min:iju_max) )
+  if ( tout_boxes(0)%lcoord_niv ) &
+    call Write_hor_coord1d( tzdim_ni_v, 'x-dimension of the grid at v location', &
+                            trim(ystdnameprefix)//'_x_coordinate_at_v_location', 'X', 0.,   iiv_bdmin, iiv_bdmax, &
+                          zxhatm_glob(iiv_min:iiv_max) )
+  if ( tout_boxes(0)%lcoord_njv ) &
+    call Write_hor_coord1d( tzdim_nj_v, 'y-dimension of the grid at v location', &
+                            trim(ystdnameprefix)//'_y_coordinate_at_v_location', 'Y', -0.5, ijv_bdmin, ijv_bdmax, &
+                            zyhat_glob(ijv_min:ijv_max)  )
 end if
 
 !Warning: the following block has to be reenabled if IO_Coordvar_write_nc4 is again called for Z-split files
@@ -1894,7 +1919,7 @@ if ( .not. lcartesian ) then
 end if
 
 if ( tzfile%lmaster ) then !vertical coordinates in the transformed space are the same on all processes
-  if ( tzfile%ctype /= 'MNHOUTPUT' .or. lout_bigbox_write ) then !Do not write coordinates of main domain if it is not written
+  if ( tzfile%ctype /= 'MNHOUTPUT' ) then !Do not write coordinates of main domain if it is not written
     if ( Trim( yprogram ) /= 'PGD' .and. Trim( yprogram ) /= 'NESPGD' .and. Trim( yprogram ) /= 'ZOOMPG' &
         .and. .not. ( Trim( yprogram ) == 'REAL' .and. cstorage_type == 'SU') ) then !condition to detect prep_surfex
       call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_LEVEL),  'position z in the transformed space',              '', &
@@ -1902,6 +1927,13 @@ if ( tzfile%lmaster ) then !vertical coordinates in the transformed space are th
       call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_LEVEL_W),'position z in the transformed space at w location','', &
                             'altitude_at_w_location', -0.5, ikw_bdmin, ikw_bdmax, ZZHAT (ikw_min:ikw_max) )
     end if
+  else if ( tzfile%ctype == 'MNHOUTPUT' .and. lout_bigbox_write ) then !Do not write coordinates of main domain if it is not written
+    if ( tout_boxes(0)%lcoord_levelm ) &
+      call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_LEVEL),  'position z in the transformed space',              '', &
+                            'altitude',                0.,  ikm_bdmin, ikm_bdmax, ZZHATM(ikm_min:ikm_max) )
+    if ( tout_boxes(0)%lcoord_levelw ) &
+      call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_LEVEL_W),'position z in the transformed space at w location','', &
+                            'altitude_at_w_location', -0.5, ikw_bdmin, ikw_bdmax, ZZHAT (ikw_min:ikw_max) )
   end if
 end if
 
@@ -2659,37 +2691,45 @@ subroutine Write_box_coords( kbox, tpbox )
 
   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) )
+  if ( tpbox%lcoord_nim ) &
+    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)                   )
+
+  if ( tpbox%lcoord_njm ) &
+    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)                   )
+
+  if ( tpbox%lcoord_niu ) &
+    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)                                 )
+
+  if ( tpbox%lcoord_nju ) &
+    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)                                 )
+
+  if ( tpbox%lcoord_niv ) &
+    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)                                 )
+
+  if ( tpbox%lcoord_njv ) &
+    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)                                )
+
+  if ( tpbox%lcoord_levelm ) &
+    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) )
+
+  if ( tpbox%lcoord_levelw ) &
+    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) )
 
   if ( .not. lcartesian ) then
     ! Mass point
diff --git a/src/MNH/modd_outn.f90 b/src/MNH/modd_outn.f90
index 423b60b7cbc6fee847f41434fa5cc7a156bde7c6..0125394fe5fd54b82b963ad7e430da0b631b8b90 100644
--- a/src/MNH/modd_outn.f90
+++ b/src/MNH/modd_outn.f90
@@ -60,6 +60,16 @@ TYPE TOUTBOXMETADATA
 
   ! Variables to decide if some coordinates are used and are to be written
   ! Set to true if the coordinate is necessary
+  ! 1D coordinates
+  LOGICAL :: LCOORD_NIM    = .FALSE. ! ni
+  LOGICAL :: LCOORD_NJM    = .FALSE. ! ni
+  LOGICAL :: LCOORD_NIU    = .FALSE. ! ni_u
+  LOGICAL :: LCOORD_NJU    = .FALSE. ! ni_u
+  LOGICAL :: LCOORD_NIV    = .FALSE. ! ni_v
+  LOGICAL :: LCOORD_NJV    = .FALSE. ! ni_v
+  LOGICAL :: LCOORD_LEVELM = .FALSE. ! level
+  LOGICAL :: LCOORD_LEVELW = .FALSE. ! level_w
+  ! 2D coordinates
   LOGICAL :: LCOORD_LATM = .FALSE. ! Latitude  at mass points
   LOGICAL :: LCOORD_LONM = .FALSE. ! Longitude at mass points
   LOGICAL :: LCOORD_LATU = .FALSE. ! Latitude  at u points