diff --git a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90 b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
index 4e3a02ea91c4dbaa22f7b0692f50111c012f474c..7c7b8632372d3024acf6f0befb6db3053b6afdde 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_field_write.f90
@@ -3961,6 +3961,7 @@ END SUBROUTINE IO_Fieldlist_write
 
 SUBROUTINE IO_Fieldlist_1field_write( TPOUTPUT, KMI, TPFIELD, TPBOX )
 
+USE MODD_DIM_n,      ONLY: 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
@@ -3972,6 +3973,7 @@ INTEGER,                         INTENT(IN) :: KMI
 TYPE(TFIELDDATA),                INTENT(IN) :: TPFIELD
 TYPE(TOUTBOXMETADATA), OPTIONAL, INTENT(IN) :: TPBOX
 
+INTEGER              :: IKSUP
 TYPE(TFIELDMETADATA) :: TZFIELDMD
 
 IF ( PRESENT(TPBOX) ) THEN
@@ -4284,10 +4286,16 @@ NDIMS: SELECT CASE (TPFIELD%NDIMS)
                   TZFIELDMD%NDIMLIST(3) = NMNHDIM_BOX_LEVEL_W
               END SELECT
             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 = JPVEXT + TPBOX%NKSUP  )
+                                     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 5aabfbcce72378ad159c0a46aeb3e66dcf16eb76..fc854cb42193a46a715da7b20933d749b792b808 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
@@ -815,36 +815,39 @@ SUBROUTINE IO_BOX_PREPARE( KMI )
     OUT_MODEL(IMI)%LOUT_TAL_REMOVE = .FALSE.
   END IF
 
+  ! Remove the vertical unphysical points?
+  OUT_MODEL(IMI)%LOUT_VER_BORDER_REMOVE = LOUT_UNPHYSICAL_VER_CELLS_REMOVE(IMI)
+
   OUT_MODEL(IMI)%NOUT_NBOXES = NOUT_BOXES(IMI)
 
   ! Allocate boxes
-  IF ( OUT_MODEL(IMI)%LOUT_BIGBOX_WRITE ) THEN
     ! Allocate also a special box for the main domain (box number 0)
     ! This is useful to store its boundaries (ie used if we remove unphysical boundaries or top absorbing layer)
-    ALLOCATE( OUT_MODEL(IMI)%TOUT_BOXES(0:NOUT_BOXES(IMI)) )
-  ELSE
-    ALLOCATE( OUT_MODEL(IMI)%TOUT_BOXES(NOUT_BOXES(IMI)) )
-  END IF
+  ALLOCATE( OUT_MODEL(IMI)%TOUT_BOXES(0:NOUT_BOXES(IMI)) )
 
   TOUT_BOXES => OUT_MODEL(IMI)%TOUT_BOXES
 
   ! Treat special box for main domain
-  IF ( OUT_MODEL(IMI)%LOUT_BIGBOX_WRITE ) THEN
-    TOUT_BOXES(0)%CNAME = CMAINDOMAINNAME
-
-    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
+  TOUT_BOXES(0)%CNAME = CMAINDOMAINNAME
+
+  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_VER_BORDER_REMOVE ) THEN
+    TOUT_BOXES(0)%NKINF = 1
+    TOUT_BOXES(0)%NKSUP = NKMAX
+  ELSE
     TOUT_BOXES(0)%NKINF = 1 - JPVEXT
-    IF ( OUT_MODEL(IMI)%LOUT_TAL_REMOVE ) THEN
-      TOUT_BOXES(0)%NKSUP = NALBOT - JPVEXT
-    ELSE
-      TOUT_BOXES(0)%NKSUP = NKMAX + JPVEXT
-    END IF
+    TOUT_BOXES(0)%NKSUP = NKMAX + JPVEXT
+  END IF
+
+  IF ( OUT_MODEL(IMI)%LOUT_TAL_REMOVE ) THEN
+    ! Min to manage case when the TAL is in the unphysical domain (should not happen)
+    TOUT_BOXES(0)%NKSUP = MIN( TOUT_BOXES(0)%NKSUP, NALBOT - JPVEXT )
   END IF
 
   ! Treat boxes
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index 3c6fa953c83d140dd69559725fac1d1b22180127..6d7256d0976322ecd82350bf5e6ea6e0df679686 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -1489,6 +1489,8 @@ 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
@@ -1671,17 +1673,27 @@ 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' .and. lout_tal_remove ) then
-      ! Remove levels above the Top Absorbing Layer
+    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.,  JPVEXT, 0,      ZZHATM(1:NALBOT) )
+                            '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, JPVEXT, 0,      ZZHAT (1:NALBOT) )
+                            '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 )
+                            '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, 0,      ZZHAT )
+                            'altitude_at_w_location', -0.5, JPVEXT, JPVEXT-1, ZZHAT  )
     end if
   end if
 end if
diff --git a/src/MNH/modd_bakout.f90 b/src/MNH/modd_bakout.f90
index e9283580b3d71329b71da02a85bb0db121a40392..8a987a1d8cab2dc0e6c7cfdabc5dac864bae86c9 100644
--- a/src/MNH/modd_bakout.f90
+++ b/src/MNH/modd_bakout.f90
@@ -99,7 +99,8 @@ CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:,:,:), ALLOCATABLE :: COUT_BOX_VAR_SUP
                                                                                   ! in the different boxes (added to the COUT_VAR)
 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_TOP_ABSORBING_LAYER_REMOVE  = .TRUE. ! Remove the top absorbing layer
+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)
 INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOUT_BOX_ISUP
diff --git a/src/MNH/modd_outn.f90 b/src/MNH/modd_outn.f90
index 862faa7dcbfed8dc1e0f44f0a03dbab3b1842de0..e6dc62428f39876fb84e960fd5b583af69a78609 100644
--- a/src/MNH/modd_outn.f90
+++ b/src/MNH/modd_outn.f90
@@ -71,27 +71,29 @@ 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
-  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_BIGBOX_WRITE      = .TRUE. ! Write the main box/domain (if there are boxes)
+  LOGICAL                            :: LOUT_TAL_REMOVE        = .TRUE. ! Remove the top absorbing layer
+  LOGICAL                            :: LOUT_VER_BORDER_REMOVE = .TRUE. ! Remove the JPVEXT vertical external points
 
   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()
-INTEGER,               POINTER :: NBAK_STEPFREQFIRST => NULL()
-INTEGER,               POINTER :: NOUT_STEPFREQFIRST => NULL()
-INTEGER, DIMENSION(:), POINTER :: NBAK_STEPLIST      => NULL()
-INTEGER, DIMENSION(:), POINTER :: NOUT_STEPLIST      => NULL()
-INTEGER, DIMENSION(:), POINTER :: NOUT_FIELDLIST     => NULL()
-LOGICAL,               POINTER :: LOUT_BIGBOX_WRITE  => NULL()
-LOGICAL,               POINTER :: LOUT_TAL_REMOVE    => NULL()
+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()
+INTEGER,               POINTER :: NBAK_STEPFREQFIRST        => NULL()
+INTEGER,               POINTER :: NOUT_STEPFREQFIRST        => NULL()
+INTEGER, DIMENSION(:), POINTER :: NBAK_STEPLIST             => NULL()
+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_VER_BORDER_REMOVE    => NULL()
 TYPE(TOUTBOXMETADATA), DIMENSION(:), POINTER :: TOUT_BOXES  => NULL()
 
 CONTAINS
@@ -100,20 +102,21 @@ SUBROUTINE OUT_GOTO_MODEL(KFROM, KTO)
 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
-NBAK_STEPFREQFIRST => OUT_MODEL(KTO)%NBAK_STEPFREQFIRST
-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
-LOUT_BIGBOX_WRITE  => OUT_MODEL(KTO)%LOUT_BIGBOX_WRITE
-LOUT_TAL_REMOVE    => OUT_MODEL(KTO)%LOUT_TAL_REMOVE
-TOUT_BOXES         => OUT_MODEL(KTO)%TOUT_BOXES
+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
+NBAK_STEPFREQFIRST     => OUT_MODEL(KTO)%NBAK_STEPFREQFIRST
+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
+LOUT_BIGBOX_WRITE      => OUT_MODEL(KTO)%LOUT_BIGBOX_WRITE
+LOUT_TAL_REMOVE        => OUT_MODEL(KTO)%LOUT_TAL_REMOVE
+LOUT_VER_BORDER_REMOVE => OUT_MODEL(KTO)%LOUT_VER_BORDER_REMOVE
+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 56db2dbe85424307d64c26d71f8b96fef279bcef..842c16d9d2e32c5f543869a298e80dcbdbc5803e 100644
--- a/src/MNH/modn_output.f90
+++ b/src/MNH/modn_output.f90
@@ -55,7 +55,8 @@ NAMELIST/NAM_OUTPUT/LOUT_BEG,LOUT_END,&
                    LOUT_COMPRESS_LOSSY, COUT_COMPRESS_LOSSY_ALGO, NOUT_COMPRESS_LOSSY_NSD, &
                    COUT_DIR, &
                    NOUT_BOXES, COUT_BOX_NAME, COUT_BOX_VAR_SUPP, &
-                   LOUT_MAINDOMAIN_WRITE, LOUT_TOP_ABSORBING_LAYER_REMOVE, &
+                   LOUT_MAINDOMAIN_WRITE, &
+                   LOUT_TOP_ABSORBING_LAYER_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.