From 04c791a5f9fbacba88f84bc76ba163c4226515ef Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Fri, 29 Mar 2024 16:01:21 +0100
Subject: [PATCH] Philippe 29/03/2024: outputs: remove non physical horizontal
 levels by default (possible to write them if needed)

---
 src/LIB/SURCOUCHE/src/mode_io_field_write.f90 |  24 +-
 .../SURCOUCHE/src/mode_io_manage_struct.f90   |  20 +-
 src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90   |  23 +-
 src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90   | 214 ++++++++++++++----
 src/MNH/modd_bakout.f90                       |   1 +
 src/MNH/modd_outn.f90                         |   3 +
 src/MNH/modn_output.f90                       |   2 +-
 7 files changed, 225 insertions(+), 62 deletions(-)

diff --git a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90 b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
index 7c7b86323..b0a8dda19 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
@@ -3961,7 +3961,7 @@ END SUBROUTINE IO_Fieldlist_write
 
 SUBROUTINE IO_Fieldlist_1field_write( TPOUTPUT, KMI, TPFIELD, TPBOX )
 
-USE MODD_DIM_n,      ONLY: NKMAX
+USE MODD_DIM_n,      ONLY: NIMAX_ll, NJMAX_ll, NKMAX
 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
@@ -3973,7 +3973,7 @@ INTEGER,                         INTENT(IN) :: KMI
 TYPE(TFIELDDATA),                INTENT(IN) :: TPFIELD
 TYPE(TOUTBOXMETADATA), OPTIONAL, INTENT(IN) :: TPBOX
 
-INTEGER              :: IKSUP
+INTEGER              :: IISUP, IJSUP, IKSUP
 TYPE(TFIELDMETADATA) :: TZFIELDMD
 
 IF ( PRESENT(TPBOX) ) THEN
@@ -4286,16 +4286,28 @@ NDIMS: SELECT CASE (TPFIELD%NDIMS)
                   TZFIELDMD%NDIMLIST(3) = NMNHDIM_BOX_LEVEL_W
               END SELECT
             END IF
+            IF ( TPBOX%CNAME == CMAINDOMAINNAME .AND. TZFIELDMD%NGRID == 2 ) THEN
+              ! There is one more PHYSICAL horizontal layer for u points (but the same number for physical + unphysical borders)
+              IISUP = MIN( JPHEXT + TPBOX%NISUP + 1, NIMAX_ll + 2*JPHEXT )
+            ELSE
+              IISUP = JPHEXT + TPBOX%NISUP
+            END IF
+            IF ( TPBOX%CNAME == CMAINDOMAINNAME .AND. TZFIELDMD%NGRID == 3 ) THEN
+              ! There is one more PHYSICAL horizontal layer for v points (but the same number for physical + unphysical borders)
+              IJSUP = MIN( JPHEXT + TPBOX%NJSUP + 1, NJMAX_ll + 2*JPHEXT )
+            ELSE
+              IJSUP = JPHEXT + TPBOX%NJSUP
+            END IF
             IF ( TPBOX%CNAME == CMAINDOMAINNAME .AND. TZFIELDMD%NGRID == 4 ) THEN
               ! There is one more PHYSICAL vertical layer for w points (but the same number for physical + unphysical borders)
               IKSUP = MIN( JPVEXT + TPBOX%NKSUP + 1, NKMAX + 2*JPVEXT )
             ELSE
               IKSUP = JPVEXT + TPBOX%NKSUP
             END IF
-            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 =                IKSUP  )
+            CALL IO_Field_write_BOX( TPOUTPUT, TZFIELDMD, 'OTHER', TPFIELD%TFIELD_X3D(KMI)%DATA, &
+                                     KXOBOX = JPHEXT + TPBOX%NIINF, KXEBOX = IISUP,              &
+                                     KYOBOX = JPHEXT + TPBOX%NJINF, KYEBOX = IJSUP,              &
+                                     KZOBOX = JPVEXT + TPBOX%NKINF, KZEBOX = IKSUP               )
           END IF
         ELSE
           call Print_msg( NVERB_ERROR, 'IO', 'IO_Fieldlist_1field_write', trim(tpfield%cmnhname)// &
diff --git a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90 b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
index fc854cb42..c1a0c3a94 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
@@ -815,6 +815,9 @@ SUBROUTINE IO_BOX_PREPARE( KMI )
     OUT_MODEL(IMI)%LOUT_TAL_REMOVE = .FALSE.
   END IF
 
+  ! Remove the horizontal unphysical points?
+  OUT_MODEL(IMI)%LOUT_HOR_BORDER_REMOVE = LOUT_UNPHYSICAL_HOR_CELLS_REMOVE(IMI)
+
   ! Remove the vertical unphysical points?
   OUT_MODEL(IMI)%LOUT_VER_BORDER_REMOVE = LOUT_UNPHYSICAL_VER_CELLS_REMOVE(IMI)
 
@@ -832,11 +835,18 @@ SUBROUTINE IO_BOX_PREPARE( KMI )
 
   ALLOCATE( TOUT_BOXES(0)%NFIELDLIST_SUPP(0) )
 
-  ! Set boundaries in physical domain coordinates, but must cover all the domain with non-physical values
-  TOUT_BOXES(0)%NIINF = 1 - JPHEXT
-  TOUT_BOXES(0)%NISUP = NIMAX_ll + JPHEXT
-  TOUT_BOXES(0)%NJINF = 1 - JPHEXT
-  TOUT_BOXES(0)%NJSUP = NJMAX_ll + JPHEXT
+  IF ( OUT_MODEL(IMI)%LOUT_HOR_BORDER_REMOVE ) THEN
+    TOUT_BOXES(0)%NIINF = 1
+    TOUT_BOXES(0)%NISUP = NIMAX_ll
+    TOUT_BOXES(0)%NJINF = 1
+    TOUT_BOXES(0)%NJSUP = NJMAX_ll
+  ELSE
+    ! Set boundaries in physical domain coordinates, but must cover all the domain with non-physical values
+    TOUT_BOXES(0)%NIINF = 1 - JPHEXT
+    TOUT_BOXES(0)%NISUP = NIMAX_ll + JPHEXT
+    TOUT_BOXES(0)%NJINF = 1 - JPHEXT
+    TOUT_BOXES(0)%NJSUP = NJMAX_ll + JPHEXT
+  END IF
   IF ( OUT_MODEL(IMI)%LOUT_VER_BORDER_REMOVE ) THEN
     TOUT_BOXES(0)%NKINF = 1
     TOUT_BOXES(0)%NKSUP = NKMAX
diff --git a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
index 1a0eceddb..e22aab89d 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
@@ -296,7 +296,8 @@ CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: HPROGRAM_ORIG !To emulate a file coming
 CHARACTER(LEN=:),ALLOCATABLE :: YPROGRAM
 integer                      :: iavg, iprof, istation
 integer                      :: ispectra_ni, ispectra_nj
-INTEGER                      :: IIU_ll, IJU_ll, IKU, IKU_MAX
+INTEGER                      :: IIU_ll, IJU_ll, IKU
+INTEGER                      :: IIU_MAX, IJU_MAX, IKU_MAX
 integer                      :: jbox
 
 CALL PRINT_MSG(NVERB_DEBUG,'IO','IO_Knowndims_set_nc4','called for '//TRIM(TPFILE%CNAME))
@@ -307,14 +308,18 @@ ELSE
   YPROGRAM = CPROGRAM
 ENDIF
 
-IIU_ll = NIMAX_ll + 2*JPHEXT
-IJU_ll = NJMAX_ll + 2*JPHEXT
 IF ( TPFILE%CTYPE == 'MNHOUTPUT' ) THEN
-  IKU = TOUT_BOXES(0)%NKSUP - TOUT_BOXES(0)%NKINF + 1
+  IIU_ll = TOUT_BOXES(0)%NISUP - TOUT_BOXES(0)%NIINF + 1
+  IJU_ll = TOUT_BOXES(0)%NJSUP - TOUT_BOXES(0)%NJINF + 1
+  IKU    = TOUT_BOXES(0)%NKSUP - TOUT_BOXES(0)%NKINF + 1
 ELSE
-  IKU = NKMAX + 2*JPVEXT
+  IIU_ll = NIMAX_ll + 2*JPHEXT
+  IJU_ll = NJMAX_ll + 2*JPHEXT
+  IKU    = NKMAX    + 2*JPVEXT
 END IF
-IKU_MAX = NKMAX + 2*JPVEXT
+IIU_MAX = NIMAX_ll + 2*JPHEXT
+IJU_MAX = NJMAX_ll + 2*JPHEXT
+IKU_MAX = NKMAX    + 2*JPVEXT
 
 if ( .not.Associated( tpfile%tncdims ) ) then
   call Print_msg( NVERB_ERROR, 'IO', 'IO_Knowndims_set_nc4', 'tncdims not associated for ' // Trim( tpfile%cname ) )
@@ -336,10 +341,12 @@ end if
 
 call IO_Add_dim_nc4( tpfile, NMNHDIM_NI,   'ni',   IIU_ll )
 call IO_Add_dim_nc4( tpfile, NMNHDIM_NJ,   'nj',   IJU_ll )
-call IO_Add_dim_nc4( tpfile, NMNHDIM_NI_U, 'ni_u', IIU_ll )
+! There is one more PHYSICAL horizontal layer for u points (but the same number for physical + unphysical borders)
+call IO_Add_dim_nc4( tpfile, NMNHDIM_NI_U, 'ni_u', MIN( IIU_ll+1, IIU_MAX ) )
 call IO_Add_dim_nc4( tpfile, NMNHDIM_NJ_U, 'nj_u', IJU_ll )
 call IO_Add_dim_nc4( tpfile, NMNHDIM_NI_V, 'ni_v', IIU_ll )
-call IO_Add_dim_nc4( tpfile, NMNHDIM_NJ_V, 'nj_v', IJU_ll )
+! There is one more PHYSICAL horizontal layer for v points (but the same number for physical + unphysical borders)
+call IO_Add_dim_nc4( tpfile, NMNHDIM_NJ_V, 'nj_v', MIN( IJU_ll+1, IJU_MAX ) )
 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 IO_Add_dim_nc4( tpfile, NMNHDIM_LEVEL,   'level',   IKU )
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index 6d7256d09..b43314177 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -1444,7 +1444,7 @@ use modd_budget,     only: cbutype, lbu_icp, lbu_jcp, lbu_kcp, nbuih, nbuil, nbu
 use modd_conf,       only: cprogram, l2d, lcartesian
 use modd_conf_n,     only: cstorage_type
 use modd_diag_flag,  only: ltraj
-use modd_dim_n,      only: nkmax
+use modd_dim_n,      only: nimax_ll, njmax_ll, nkmax
 use modd_dyn_n,      only: nalbot, xtstep
 use modd_field,      only: NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_NI_U, NMNHDIM_NJ_U, NMNHDIM_NI_V, NMNHDIM_NJ_V, &
                            NMNHDIM_LEVEL, NMNHDIM_LEVEL_W, NMNHDIM_TIME, NMNHDIM_TRAJ_TIME,                &
@@ -1489,13 +1489,21 @@ character(len=:),                         allocatable :: ystdnameprefix
 character(len=:),                         allocatable :: yprogram
 integer                                               :: iiu, iju
 integer                                               :: id, iid, iresp
-integer                                               :: ilow, imax
-integer                                               :: iboundlow, iboundmax
 integer                                               :: imi
 integer                                               :: ji
 integer                                               :: jt
 integer                                               :: jtb, jte
 integer(kind=cdfint)                                  :: incid
+integer                                               :: iim_min, iim_max, ijm_min, ijm_max
+integer                                               :: iiu_min, iiu_max, iju_min, iju_max
+integer                                               :: iiv_min, iiv_max, ijv_min, ijv_max
+integer                                               :: ikm_min, ikm_max
+integer                                               :: ikw_min, ikw_max
+integer                                               :: iim_bdmin, iim_bdmax, ijm_bdmin, ijm_bdmax
+integer                                               :: iiu_bdmin, iiu_bdmax, iju_bdmin, iju_bdmax
+integer                                               :: iiv_bdmin, iiv_bdmax, ijv_bdmin, ijv_bdmax
+integer                                               :: ikm_bdmin, ikm_bdmax
+integer                                               :: ikw_bdmin, ikw_bdmax
 logical                                               :: gchangemodel
 logical                                               :: gdealloc
 logical,                         pointer              :: gsleve
@@ -1548,6 +1556,110 @@ incid = tzfile%nncid
 
 call Get_model_number_ll( imi )
 
+! Get domain boundaries (modified in some cases if not all domain is written in file such as for MNHOUTPUT files)
+if ( tzfile%ctype /= 'MNHOUTPUT' ) then
+  ! Coordinates boundaries
+  iim_min = 1
+  iim_max = nimax_ll + 2 * jphext
+  ijm_min = 1
+  ijm_max = njmax_ll + 2 * jphext
+
+  iiu_min = iim_min
+  iiu_max = iim_max
+  iju_min = ijm_min
+  iju_max = ijm_max
+
+  iiv_min = iim_min
+  iiv_max = iim_max
+  ijv_min = ijm_min
+  ijv_max = ijm_max
+
+  ikm_min = 1
+  ikm_max = nkmax + 2 * JPVEXT
+
+  ikw_min = ikm_min
+  ikw_max = ikm_max
+
+  ! Coordinates: number of unphysical values on boundaries
+  iim_bdmin = jphext
+  iim_bdmax = jphext
+  ijm_bdmin = jphext
+  ijm_bdmax = jphext
+
+  iiu_bdmin = jphext
+  iiu_bdmax = jphext - 1
+  iju_bdmin = jphext
+  iju_bdmax = jphext
+
+  iiv_bdmin = jphext
+  iiv_bdmax = jphext
+  ijv_bdmin = jphext
+  ijv_bdmax = jphext - 1
+
+  ikm_bdmin = JPVEXT
+  ikm_bdmax = JPVEXT
+
+  ikw_bdmin = JPVEXT
+  ikw_bdmax = JPVEXT - 1
+else
+  ! Take the selected dimensions (ie if the unphysical points are removed)
+  ! Coordinates boundaries
+  iim_min = tout_boxes(0)%niinf + jphext
+  iim_max = tout_boxes(0)%nisup + jphext
+  ijm_min = tout_boxes(0)%njinf + jphext
+  ijm_max = tout_boxes(0)%njsup + jphext
+
+  iiu_min = iim_min
+  iiu_max = min( iim_max + 1, nimax_ll + 2 * jphext )
+  iju_min = ijm_min
+  iju_max = ijm_max
+
+  iiv_min = iim_min
+  iiv_max = iim_max
+  ijv_min = ijm_min
+  ijv_max = min( ijm_max + 1, njmax_ll + 2 * jphext )
+
+  ikm_min = tout_boxes(0)%nkinf + JPVEXT
+  ikm_max = tout_boxes(0)%nksup + JPVEXT
+
+  ikw_min = ikm_min
+  ikw_max = min( ikm_max + 1, nkmax + 2 * JPVEXT )
+
+  ! Coordinates: number of unphysical values on boundaries
+  iim_bdmin = max( 0, jphext - iim_min + 1 )
+  iim_bdmax = max( 0, iim_max - (nimax_ll + 2*jphext) + 1 )
+  ijm_bdmin = max( 0, jphext - ijm_min + 1 )
+  ijm_bdmax = max( 0, ijm_max - (njmax_ll + 2*jphext) + 1 )
+
+  iiu_bdmin = max( 0, jphext - iiu_min + 1 )
+  iiu_bdmax = max( 0, iim_max - (nimax_ll + 2*jphext) )
+  iju_bdmin = max( 0, jphext - iju_min + 1 )
+  iju_bdmax = max( 0, iju_max - (njmax_ll + 2*jphext) + 1 )
+
+  iiv_bdmin = max( 0, jphext - iiv_min + 1 )
+  iiv_bdmax = max( 0, iiv_max - (nimax_ll + 2*jphext) + 1 )
+  ijv_bdmin = max( 0, jphext - ijv_min + 1 )
+  ijv_bdmax = max( 0, iiv_max - (nimax_ll + 2*jphext) )
+
+  ikm_bdmin = max( 0, JPVEXT - ikm_min )
+  ikm_bdmax = max( 0, ikm_max - (nkmax + 2*JPVEXT) + 1 )
+
+  ikw_bdmin = max( 0, JPVEXT - ikw_min )
+  ikw_bdmax = max( 0, ikm_max - (nkmax + 2*JPVEXT) )
+
+  ! Adapt the dimensions to the selected ones
+  if ( tzfile%lmaster ) then
+    tzfile%tncdims%tdims(NMNHDIM_NI)%nlen      = iim_max - iim_min + 1
+    tzfile%tncdims%tdims(NMNHDIM_NJ)%nlen      = ijm_max - ijm_min + 1
+    tzfile%tncdims%tdims(NMNHDIM_NI_U)%nlen    = iiu_max - iiu_min + 1
+    tzfile%tncdims%tdims(NMNHDIM_NJ_U)%nlen    = iju_max - iju_min + 1
+    tzfile%tncdims%tdims(NMNHDIM_NI_V)%nlen    = iiv_max - iiv_min + 1
+    tzfile%tncdims%tdims(NMNHDIM_NJ_V)%nlen    = ijv_max - ijv_min + 1
+    tzfile%tncdims%tdims(NMNHDIM_LEVEL)%nlen   = ikm_max - ikm_min + 1
+    tzfile%tncdims%tdims(NMNHDIM_LEVEL_W)%nlen = ikw_max - ikw_min + 1
+  end if
+end if
+
 if ( tzfile%nmodel > 0 ) then
   call Find_field_id_from_mnhname( 'XHAT', iid, iresp )
   zxhat => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
@@ -1617,17 +1729,23 @@ else
 end if
 
 call Write_hor_coord1d( tzdim_ni,   'x-dimension of the grid', &
-                        trim(ystdnameprefix)//'_x_coordinate',               'X', 0.,   jphext, jphext, zxhatm_glob )
+                        trim(ystdnameprefix)//'_x_coordinate',               'X', 0.,   iim_bdmin, iim_bdmax, &
+                        zxhatm_glob(iim_min:iim_max) )
 call Write_hor_coord1d( tzdim_nj,   'y-dimension of the grid', &
-                        trim(ystdnameprefix)//'_y_coordinate',               'Y', 0.,   jphext, jphext, zyhatm_glob )
+                        trim(ystdnameprefix)//'_y_coordinate',               'Y', 0.,   ijm_bdmin, ijm_bdmax, &
+                        zyhatm_glob(ijm_min:ijm_max) )
 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, jphext, 0,      zxhat_glob  )
+                        trim(ystdnameprefix)//'_x_coordinate_at_u_location', 'X', -0.5, iiu_bdmin, iiu_bdmax, &
+                        zxhat_glob(iiu_min:iiu_max)  )
 call Write_hor_coord1d( tzdim_nj_u, 'y-dimension of the grid at u location', &
-                        trim(ystdnameprefix)//'_y_coordinate_at_u_location', 'Y', 0.,   jphext, jphext, zyhatm_glob )
+                        trim(ystdnameprefix)//'_y_coordinate_at_u_location', 'Y', 0.,   iju_bdmin, iju_bdmax, &
+                        zyhatm_glob(iju_min:iju_max) )
 call Write_hor_coord1d( tzdim_ni_v, 'x-dimension of the grid at v location', &
-                        trim(ystdnameprefix)//'_x_coordinate_at_v_location', 'X', 0.,   jphext, jphext, zxhatm_glob )
+                        trim(ystdnameprefix)//'_x_coordinate_at_v_location', 'X', 0.,   iiv_bdmin, iiv_bdmax, &
+                        zxhatm_glob(iiv_min:iiv_max) )
 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, jphext, 0,      zyhat_glob  )
+                        trim(ystdnameprefix)//'_y_coordinate_at_v_location', 'Y', -0.5, ijv_bdmin, ijv_bdmax, &
+                        zyhat_glob(ijv_min:ijv_max)  )
 
 !Warning: the following block has to be reenabled if IO_Coordvar_write_nc4 is again called for Z-split files
 #if 0
@@ -1655,13 +1773,13 @@ if ( .not. lcartesian ) then
   end if
 
   ! Mass point
-  call Write_hor_coord2d( zlatm_glob, zlonm_glob, 'latitude',  'longitude')
+  call Write_hor_coord2d( zlatm_glob, zlonm_glob,  'latitude',  'longitude',  iim_min, iim_max, ijm_min, ijm_max )
   ! u point
-  call Write_hor_coord2d( zlatu_glob, zlonu_glob, 'latitude_u','longitude_u')
+  call Write_hor_coord2d( zlatu_glob, zlonu_glob, 'latitude_u','longitude_u', iiu_min, iiu_max, iju_min, iju_max )
   ! v point
-  call Write_hor_coord2d( zlatv_glob, zlonv_glob, 'latitude_v','longitude_v')
+  call Write_hor_coord2d( zlatv_glob, zlonv_glob, 'latitude_v','longitude_v', iiv_min, iiv_max, ijv_min, ijv_max )
   ! xi vorticity point (=f point =uv point)
-  call Write_hor_coord2d( zlatf_glob, zlonf_glob, 'latitude_f','longitude_f')
+  call Write_hor_coord2d( zlatf_glob, zlonf_glob, 'latitude_f','longitude_f', iiu_min, iiu_max, ijv_min, ijv_max )
 
   Deallocate( zlat, zlon )
 
@@ -1673,28 +1791,10 @@ end if
 if ( tzfile%lmaster ) then !vertical coordinates in the transformed space are the same on all processes
   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
-    if ( tzfile%ctype == 'MNHOUTPUT' ) then
-      ! Write only selected levels
-      ilow = tout_boxes(0)%nkinf + JPVEXT
-      imax = tout_boxes(0)%nksup + JPVEXT
-      ! Determine the dynamic range (physical limits)
-      ! If we are fully in the physical domain => iboundlow/iboundhigh = 0,
-      ! else remove of the dynamic range the non-physical element
-      iboundlow = MAX( 0, JPVEXT - ilow + 1 )
-      iboundmax = MAX( 0, imax - (nkmax + 2*JPVEXT) + 1 )
-      call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_LEVEL),  'position z in the transformed space',              '', &
-                            'altitude',                0.,  iboundlow, iboundmax, ZZHATM(ilow:imax) )
-      ! There is one more PHYSICAL vertical layer for w points (but the same number for physical + unphysical borders)
-      imax = min( imax+1, nkmax + 2*JPVEXT )
-      iboundmax = max( iboundmax-1, 0 )
-      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, iboundlow, iboundmax, ZZHAT (ilow:imax) )
-    else
-      call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_LEVEL),  'position z in the transformed space',              '', &
-                            'altitude',                0.,  JPVEXT, JPVEXT,   ZZHATM )
-      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, JPVEXT, JPVEXT-1, ZZHAT  )
-    end if
+    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) )
+    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
 
@@ -2117,17 +2217,47 @@ subroutine Write_hor_coord1d(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUNDLOW,KBOU
 end subroutine Write_hor_coord1d
 
 
-subroutine Write_hor_coord2d( plat, plon, hlat, hlon )
-  real,dimension(:,:), intent(in) :: plat
-  real,dimension(:,:), intent(in) :: plon
-  character(len=*),    intent(in) :: hlat
-  character(len=*),    intent(in) :: hlon
+subroutine Write_hor_coord2d( plat, plon, hlat, hlon, kimin, kimax, kjmin, kjmax )
+  real,dimension(:,:),           intent(in) :: plat
+  real,dimension(:,:),           intent(in) :: plon
+  character(len=*),              intent(in) :: hlat
+  character(len=*),              intent(in) :: hlon
+  integer,             optional, intent(in) :: kimin
+  integer,             optional, intent(in) :: kimax
+  integer,             optional, intent(in) :: kjmin
+  integer,             optional, intent(in) :: kjmax
+
+  integer :: iimin, iimax, ijmin, ijmax
 
   if ( tzfile%lmaster ) then
+    if ( present(kimin) ) then
+      iimin = kimin
+    else
+      iimin = LBOUND( plat, dim=1 )
+    end if
+
+    if ( present(kimax) ) then
+      iimax = kimax
+    else
+      iimax = UBOUND( plat, dim=1 )
+    end if
+
+    if ( present(kjmin) ) then
+      ijmin = kjmin
+    else
+      ijmin = LBOUND( plat, dim=2 )
+    end if
+
+    if ( present(kjmax) ) then
+      ijmax = kjmax
+    else
+      ijmax = UBOUND( plat, dim=2 )
+    end if
+
     call Find_field_id_from_mnhname( hlat, id, iresp )
-    call IO_Field_write_nc4_x2( tzfile, tfieldlist(id ), plat, iresp, oiscoord = .true. )
+    call IO_Field_write_nc4_x2( tzfile, tfieldlist(id ), plat(iimin:iimax,ijmin:ijmax), iresp, oiscoord = .true. )
     call Find_field_id_from_mnhname( hlon, id, iresp )
-    call IO_Field_write_nc4_x2( tzfile, tfieldlist(id ), plon, iresp, oiscoord = .true. )
+    call IO_Field_write_nc4_x2( tzfile, tfieldlist(id ), plon(iimin:iimax,ijmin:ijmax), iresp, oiscoord = .true. )
   end if
 end subroutine Write_hor_coord2d
 
diff --git a/src/MNH/modd_bakout.f90 b/src/MNH/modd_bakout.f90
index 8a987a1d8..c5c133683 100644
--- a/src/MNH/modd_bakout.f90
+++ b/src/MNH/modd_bakout.f90
@@ -100,6 +100,7 @@ CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:,:,:), ALLOCATABLE :: COUT_BOX_VAR_SUP
 LOGICAL, DIMENSION(JPMODELMAX) :: LOUT_MAINDOMAIN_WRITE = .FALSE. ! True to write the main domain
                                                                   ! (automatically forced to .TRUE. if NOUT_BOXES=0)
 LOGICAL, DIMENSION(JPMODELMAX) :: LOUT_TOP_ABSORBING_LAYER_REMOVE  = .TRUE. ! Remove the top absorbing layer
+LOGICAL, DIMENSION(JPMODELMAX) :: LOUT_UNPHYSICAL_HOR_CELLS_REMOVE = .TRUE. ! Remove the JPHEXT horizontal external points
 LOGICAL, DIMENSION(JPMODELMAX) :: LOUT_UNPHYSICAL_VER_CELLS_REMOVE = .TRUE. ! Remove the JPVEXT vertical external points
 
 INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_IINF ! Box coordinates in physical domain (for each model and for each box)
diff --git a/src/MNH/modd_outn.f90 b/src/MNH/modd_outn.f90
index e6dc62428..31f28032c 100644
--- a/src/MNH/modd_outn.f90
+++ b/src/MNH/modd_outn.f90
@@ -73,6 +73,7 @@ TYPE OUT_t
   INTEGER, DIMENSION(:), ALLOCATABLE :: NOUT_FIELDLIST                 ! List of fields to write in outputs
   LOGICAL                            :: LOUT_BIGBOX_WRITE      = .TRUE. ! Write the main box/domain (if there are boxes)
   LOGICAL                            :: LOUT_TAL_REMOVE        = .TRUE. ! Remove the top absorbing layer
+  LOGICAL                            :: LOUT_HOR_BORDER_REMOVE = .TRUE. ! Remove the JPHEXT horizontal external points
   LOGICAL                            :: LOUT_VER_BORDER_REMOVE = .TRUE. ! Remove the JPVEXT vertical external points
 
   TYPE(TOUTBOXMETADATA), DIMENSION(:), ALLOCATABLE :: TOUT_BOXES
@@ -93,6 +94,7 @@ INTEGER, DIMENSION(:), POINTER :: NOUT_STEPLIST             => NULL()
 INTEGER, DIMENSION(:), POINTER :: NOUT_FIELDLIST            => NULL()
 LOGICAL,               POINTER :: LOUT_BIGBOX_WRITE         => NULL()
 LOGICAL,               POINTER :: LOUT_TAL_REMOVE           => NULL()
+LOGICAL,               POINTER :: LOUT_HOR_BORDER_REMOVE    => NULL()
 LOGICAL,               POINTER :: LOUT_VER_BORDER_REMOVE    => NULL()
 TYPE(TOUTBOXMETADATA), DIMENSION(:), POINTER :: TOUT_BOXES  => NULL()
 
@@ -115,6 +117,7 @@ NOUT_STEPLIST          => OUT_MODEL(KTO)%NOUT_STEPLIST
 NOUT_FIELDLIST         => OUT_MODEL(KTO)%NOUT_FIELDLIST
 LOUT_BIGBOX_WRITE      => OUT_MODEL(KTO)%LOUT_BIGBOX_WRITE
 LOUT_TAL_REMOVE        => OUT_MODEL(KTO)%LOUT_TAL_REMOVE
+LOUT_HOR_BORDER_REMOVE => OUT_MODEL(KTO)%LOUT_HOR_BORDER_REMOVE
 LOUT_VER_BORDER_REMOVE => OUT_MODEL(KTO)%LOUT_VER_BORDER_REMOVE
 TOUT_BOXES             => OUT_MODEL(KTO)%TOUT_BOXES
 
diff --git a/src/MNH/modn_output.f90 b/src/MNH/modn_output.f90
index 842c16d9d..94b7fc108 100644
--- a/src/MNH/modn_output.f90
+++ b/src/MNH/modn_output.f90
@@ -56,7 +56,7 @@ NAMELIST/NAM_OUTPUT/LOUT_BEG,LOUT_END,&
                    COUT_DIR, &
                    NOUT_BOXES, COUT_BOX_NAME, COUT_BOX_VAR_SUPP, &
                    LOUT_MAINDOMAIN_WRITE, &
-                   LOUT_TOP_ABSORBING_LAYER_REMOVE, LOUT_UNPHYSICAL_VER_CELLS_REMOVE, &
+                   LOUT_TOP_ABSORBING_LAYER_REMOVE, LOUT_UNPHYSICAL_HOR_CELLS_REMOVE, LOUT_UNPHYSICAL_VER_CELLS_REMOVE, &
                    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.
-- 
GitLab