diff --git a/src/LIB/SURCOUCHE/src/modd_io.f90 b/src/LIB/SURCOUCHE/src/modd_io.f90
index c8457cc2d32624a020d122c41cb318faa86a1d4c..d51c4dd2b221574ae6b2f4175add47076e170c43 100644
--- a/src/LIB/SURCOUCHE/src/modd_io.f90
+++ b/src/LIB/SURCOUCHE/src/modd_io.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2023 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -12,12 +12,17 @@
 !  P. Wautelet 17/01/2020: add 'BUD' category for Print_msg + corresponding namelist variables
 !  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
 !-----------------------------------------------------------------
 
 #define MNH_REDUCE_DIMENSIONS_IN_FILES 1
 
 MODULE MODD_IO
 !
+#ifdef MNH_IOCDF4
+USE NETCDF, ONLY: NF90_QUANTIZE_GRANULARBR
+#endif
+!
 use modd_netcdf,     only: tdimsnc
 USE MODD_PARAMETERS, ONLY: NDIRNAMELGTMAX, NFILENAMELGTMAX
 use modd_precision,  only: CDFINT, LFIINT
@@ -123,6 +128,9 @@ TYPE TFILEDATA
                                                                 ! instead of double precision
   LOGICAL                :: LNCCOMPRESS = .FALSE. ! Do compression on fields
   INTEGER(KIND=CDFINT)   :: NNCCOMPRESS_LEVEL = 0 ! Compression level
+  LOGICAL                :: LNCCOMPRESS_LOSSY      = .FALSE.                  ! Do lossy compression on float fields
+  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
 #endif
   !
diff --git a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90 b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
index 8b7c8c9089385dec3964b49b76939382e54cac69..54d694804d2b40d7d94df08d376c6227b274c813 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_manage_struct.f90
@@ -22,6 +22,7 @@
 !  P. Wautelet 18/03/2022: minor bugfix in ISTEP_MAX computation + adapt diagnostics messages
 !                          (change verbosity level and remove some unnecessary warnings)
 !  P. Wautelet 13/01/2023: set NMODEL field for backup and output files
+!  P. Wautelet 14/12/2023: add lossy compression for output files
 !-----------------------------------------------------------------
 MODULE MODE_IO_MANAGE_STRUCT
 !
@@ -523,7 +524,13 @@ END SUBROUTINE SORT_ENTRIES
 SUBROUTINE POPULATE_STRUCT(TPFILE_FIRST,TPFILE_LAST,KSTEPS,HFILETYPE,TPBAKOUTN,KMI)
 !#########################################################################
   !
-  USE MODD_CONFZ, ONLY: NB_PROCIO_W
+#ifdef MNH_IOCDF4
+  USE NETCDF, ONLY: NF90_QUANTIZE_BITGROOM, NF90_QUANTIZE_BITROUND, NF90_QUANTIZE_GRANULARBR
+#endif
+  !
+  USE MODD_CONFZ,      ONLY: NB_PROCIO_W
+  !
+  USE MODE_TOOLS,      ONLY: UPCASE
   !
   TYPE(TFILEDATA),           POINTER,INTENT(INOUT) :: TPFILE_FIRST,TPFILE_LAST
   INTEGER,DIMENSION(:),              INTENT(IN)    :: KSTEPS
@@ -559,16 +566,109 @@ SUBROUTINE POPULATE_STRUCT(TPFILE_FIRST,TPFILE_LAST,KSTEPS,HFILETYPE,TPBAKOUTN,K
         IF (TRIM(HFILETYPE)=='MNHOUTPUT') THEN
           ! Add a "OUT" suffix for output files
           TPBAKOUTN(IPOS)%TFILE%CNAME=ADJUSTL(ADJUSTR(IO_SURF_MNH_MODEL(IMI)%COUTFILE)//'.OUT.'//YNUMBER)
+
+#ifdef MNH_IOCDF4
           !Reduce the float precision if asked
           TPBAKOUTN(IPOS)%TFILE%LNCREDUCE_FLOAT_PRECISION = LOUT_REDUCE_FLOAT_PRECISION(IMI)
+
           !Set compression if asked
           TPBAKOUTN(IPOS)%TFILE%LNCCOMPRESS = LOUT_COMPRESS(IMI)
           IF ( NOUT_COMPRESS_LEVEL(IMI)<0 .OR. NOUT_COMPRESS_LEVEL(IMI)>9 ) THEN
-            CALL PRINT_MSG(NVERB_WARNING,'IO','POPULATE_STRUCT',&
-                           'NOUT_COMPRESS_LEVEL must be in the [0..9] range. Value forced to 4')
+            CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT', &
+                            'NOUT_COMPRESS_LEVEL must be in the [0..9] range' )
             NOUT_COMPRESS_LEVEL(IMI) = 4
           END IF
           TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LEVEL = NOUT_COMPRESS_LEVEL(IMI)
+
+          !Set lossy compression
+          TPBAKOUTN(IPOS)%TFILE%LNCCOMPRESS_LOSSY = LOUT_COMPRESS_LOSSY(IMI)
+          IF ( LOUT_COMPRESS_LOSSY(IMI) ) THEN
+            !Force compression if lossy compression is enabled
+            TPBAKOUTN(IPOS)%TFILE%LNCCOMPRESS = .TRUE.
+
+            !Set lossy compression algorithm
+            SELECT CASE ( UPCASE( COUT_COMPRESS_LOSSY_ALGO(IMI) ) )
+              CASE ( 'BITGROOM' )
+                TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_ALGO = NF90_QUANTIZE_BITGROOM
+              CASE ( 'GRANULARBR' )
+                TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_ALGO = NF90_QUANTIZE_GRANULARBR
+              CASE ( 'BITROUND' )
+                TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_ALGO = NF90_QUANTIZE_BITROUND
+              CASE DEFAULT
+                CMNHMSG(1) = 'invalid COUT_COMPRESS_LOSSY_ALGO'
+                CMNHMSG(2) = 'Accepted algorithms: BITGROOM, GRANULARBR (default choice), BITROUND'
+                CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT' )
+                TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_ALGO = NF90_QUANTIZE_GRANULARBR
+            END SELECT
+
+            !Set number of significant digits/bits for lossy compression algorithm
+#if (MNH_REAL == 4)
+            SELECT CASE ( TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_ALGO )
+              CASE ( NF90_QUANTIZE_BITROUND )
+                ! For 32 bit reals, number of significant bits must be in the 1 to 23 range
+                IF ( NOUT_COMPRESS_LOSSY_NSD(IMI) < 1 .OR. NOUT_COMPRESS_LOSSY_NSD(IMI) > 23 ) THEN
+                  CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT', &
+                                  'NOUT_COMPRESS_LOSSY_NSD must be in the [1..23] range' )
+                  NOUT_COMPRESS_LOSSY_NSD(IMI) = 7
+                END IF
+              CASE ( NF90_QUANTIZE_BITGROOM, NF90_QUANTIZE_GRANULARBR )
+                ! For 32 bit reals, number of significant digits must be in the 1 to 7 range
+                IF ( NOUT_COMPRESS_LOSSY_NSD(IMI) < 1 .OR. NOUT_COMPRESS_LOSSY_NSD(IMI) > 7 ) THEN
+                  CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT', &
+                                  'NOUT_COMPRESS_LOSSY_NSD must be in the [1..7] range' )
+                  NOUT_COMPRESS_LOSSY_NSD(IMI) = 3
+                END IF
+              CASE DEFAULT
+                CALL PRINT_MSG( NVERB_FATAL, 'IO', 'POPULATE_STRUCT', 'invalid NNCCOMPRESS_LOSSY_ALGO (internal fatal error)' )
+            END SELECT
+#elif (MNH_REAL == 8)
+            IF ( TPBAKOUTN(IPOS)%TFILE%LNCREDUCE_FLOAT_PRECISION ) THEN
+              SELECT CASE ( TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_ALGO )
+                CASE ( NF90_QUANTIZE_BITROUND )
+                  ! For 32 bit reals, number of significant bits must be in the 1 to 23 range
+                  IF ( NOUT_COMPRESS_LOSSY_NSD(IMI) < 1 .OR. NOUT_COMPRESS_LOSSY_NSD(IMI) > 23 ) THEN
+                    CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT', &
+                                    'NOUT_COMPRESS_LOSSY_NSD must be in the [1..23] range' )
+                    NOUT_COMPRESS_LOSSY_NSD(IMI) = 7
+                  END IF
+                CASE ( NF90_QUANTIZE_BITGROOM, NF90_QUANTIZE_GRANULARBR )
+                  ! For 32 bit reals, number of significant digits must be in the 1 to 7 range
+                  IF ( NOUT_COMPRESS_LOSSY_NSD(IMI) < 1 .OR. NOUT_COMPRESS_LOSSY_NSD(IMI) > 7 ) THEN
+                    CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT', &
+                                    'NOUT_COMPRESS_LOSSY_NSD must be in the [1..7] range' )
+                    NOUT_COMPRESS_LOSSY_NSD(IMI) = 3
+                  END IF
+                CASE DEFAULT
+                  CALL PRINT_MSG( NVERB_FATAL, 'IO', 'POPULATE_STRUCT', 'invalid NNCCOMPRESS_LOSSY_ALGO (internal fatal error)' )
+              END SELECT
+            ELSE
+              SELECT CASE ( TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_ALGO )
+                CASE ( NF90_QUANTIZE_BITROUND )
+                  ! For 64 bit reals, number of significant bits must be in the 1 to 52 range
+                  IF ( NOUT_COMPRESS_LOSSY_NSD(IMI) < 1 .OR. NOUT_COMPRESS_LOSSY_NSD(IMI) > 52 ) THEN
+                    CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT', &
+                                    'NOUT_COMPRESS_LOSSY_NSD must be in the [1..52] range' )
+                    NOUT_COMPRESS_LOSSY_NSD(IMI) = 7
+                  END IF
+                CASE ( NF90_QUANTIZE_BITGROOM, NF90_QUANTIZE_GRANULARBR )
+                  ! For 64 bit reals, number of significant digits must be in the 1 to 15 range
+                  IF ( NOUT_COMPRESS_LOSSY_NSD(IMI) < 1 .OR. NOUT_COMPRESS_LOSSY_NSD(IMI) > 15 ) THEN
+                    CALL PRINT_MSG( NVERB_ERROR, 'IO', 'POPULATE_STRUCT', &
+                                    'NOUT_COMPRESS_LOSSY_NSD must be in the [1..15] range')
+                    NOUT_COMPRESS_LOSSY_NSD(IMI) = 3
+                  END IF
+                CASE DEFAULT
+                  CALL PRINT_MSG( NVERB_FATAL, 'IO', 'POPULATE_STRUCT', 'invalid NNCCOMPRESS_LOSSY_ALGO (internal fatal error)' )
+              END SELECT
+            END IF
+#else
+#error "Invalid MNH_REAL"
+#endif
+            TPBAKOUTN(IPOS)%TFILE%NNCCOMPRESS_LOSSY_NSD = NOUT_COMPRESS_LOSSY_NSD(IMI)
+          END IF
+#endif
+
+          !Set output directory
           IF (LEN_TRIM(COUT_DIR)>0) THEN
             TPBAKOUTN(IPOS)%TFILE%CDIRNAME=TRIM(COUT_DIR)
           ELSE IF (LEN_TRIM(CIO_DIR)>0) THEN
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index 9559f6b17ca0dcaa3816b5dcb764b2db48b8a678..cdc9e0c43a2ea38b320bb7c7f45295ce43ad81c5 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -32,6 +32,7 @@
 !  P. Wautelet    06/2022: reorganize flyers
 !  P. Wautelet 21/06/2022: bugfix: time_budget was not computed correctly (tdtexp -> tdtseg)
 !  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
 !-----------------------------------------------------------------
 #ifdef MNH_IOCDF4
 module mode_io_write_nc4
@@ -44,9 +45,10 @@ use modd_precision,    only: CDFINT, MNHINT_NF90, MNHREAL32, MNHREAL_MPI, MNHREA
 use mode_io_tools_nc4, only: IO_Mnhname_clean, IO_Vdims_fill_nc4, IO_Dim_find_create_nc4, IO_Strdimid_get_nc4, IO_Err_handle_nc4
 use mode_msg
 
-use NETCDF,            only: NF90_CHAR, NF90_FLOAT, NF90_INT1,                                    &
-                             NF90_GLOBAL, NF90_NOERR,                                             &
-                             NF90_DEF_VAR, NF90_DEF_VAR_DEFLATE, NF90_GET_ATT, NF90_INQ_VARID,    &
+use NETCDF,            only: NF90_CHAR, NF90_FLOAT, NF90_INT1,                          &
+                             NF90_GLOBAL, NF90_NOERR,                                   &
+                             NF90_DEF_VAR, NF90_DEF_VAR_DEFLATE, NF90_DEF_VAR_QUANTIZE, &
+                             NF90_GET_ATT, NF90_INQ_VARID,                              &
                              NF90_INQUIRE_ATTRIBUTE, NF90_PUT_ATT, NF90_PUT_VAR
 
 implicit none
@@ -90,8 +92,6 @@ use modd_parameters, only: jphext
 
 use mode_tools_ll,   only: Get_globaldims_ll
 
-use NETCDF,          only: NF90_FLOAT
-
 type(tfiledata),       intent(in) :: tpfile
 class(tfieldmetadata), intent(in) :: tpfield
 integer,               intent(in) :: knblocks
@@ -377,7 +377,6 @@ END SUBROUTINE IO_Field_attr_write_nc4
 
 
 subroutine IO_Field_create_nc4( tpfile, tpfield, kshape, hcalendar, oiscoord, kvarid, oisempty )
-use NETCDF, only: NF90_CHAR, NF90_FLOAT, NF90_INT1
 
 use modd_field,     only: NMNHDIM_TIME, TYPECHAR, TYPEDATE, TYPEINT, TYPELOG, TYPEREAL, TYPEUNDEF
 use modd_precision, only: MNHINT_NF90, MNHREAL_NF90
@@ -517,6 +516,12 @@ if ( istatus /= NF90_NOERR ) then
   if ( tpfile%lnccompress .and. tpfield%ntype == TYPEREAL .and. tpfield%ndims >= 1 ) then
     istatus = NF90_DEF_VAR_DEFLATE( tpfile%nncid, ivarid, SHUFFLE, DEFLATE, tpfile%nnccompress_level )
     if ( istatus /= NF90_NOERR ) call IO_Err_handle_nc4( istatus, 'IO_Field_create_nc4', 'NF90_DEF_VAR_DEFLATE', Trim( yvarname ) )
+
+    if ( tpfile%lnccompress_lossy ) then
+      istatus = NF90_DEF_VAR_QUANTIZE( tpfile%nncid, ivarid, tpfile%nnccompress_lossy_algo, tpfile%nnccompress_lossy_nsd )
+      if ( istatus /= NF90_NOERR ) call IO_Err_handle_nc4( istatus, 'IO_Field_create_nc4', &
+                                                           'NF90_DEF_VAR_QUANTIZE', Trim( yvarname ) )
+    end if
   end if
 else
   gexisted = .true.
@@ -1503,10 +1508,15 @@ real, dimension(:,:), pointer, save :: zlatm_glob  => null(), zlonm_glob  => nul
 real, dimension(:,:), pointer, save :: zlatu_glob  => null(), zlonu_glob  => null()
 real, dimension(:,:), pointer, save :: zlatv_glob  => null(), zlonv_glob  => null()
 real, dimension(:,:), pointer, save :: zlatf_glob  => null(), zlonf_glob  => null()
+type(tfiledata)                                       :: tzfile
 
 
 call Print_msg( NVERB_DEBUG, 'IO', 'IO_Coordvar_write_nc4', 'called for ' // Trim( tpfile%cname ) )
 
+tzfile = tpfile
+! Disable lossy compression to not loose precision for coordinates variables
+tzfile%lnccompress_lossy = .false.
+
 zxhat  => null()
 zyhat  => null()
 zzhat  => null()
@@ -1527,37 +1537,37 @@ else
 endif
 
 ! Get the Netcdf file ID
-incid = tpfile%nncid
+incid = tzfile%nncid
 
 call Get_model_number_ll( imi )
 
-if ( tpfile%nmodel > 0 ) then
+if ( tzfile%nmodel > 0 ) then
   call Find_field_id_from_mnhname( 'XHAT', iid, iresp )
-  zxhat => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zxhat => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'YHAT', iid, iresp )
-  zyhat => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zyhat => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'XHATM', iid, iresp )
-  zxhatm => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zxhatm => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'YHATM', iid, iresp )
-  zyhatm => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zyhatm => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'ZHAT', iid, iresp )
-  zzhat => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zzhat => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'ZHATM', iid, iresp )
-  zzhatm => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zzhatm => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'XHAT_ll', iid, iresp )
-  zxhat_glob => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zxhat_glob => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'YHAT_ll', iid, iresp )
-  zyhat_glob => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zyhat_glob => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'XHATM_ll', iid, iresp )
-  zxhatm_glob => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zxhatm_glob => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'YHATM_ll', iid, iresp )
-  zyhatm_glob => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  zyhatm_glob => tfieldlist(iid)%tfield_x1d(tzfile%nmodel)%data
   call Find_field_id_from_mnhname( 'SLEVE', iid, iresp )
-  gsleve => tfieldlist(iid)%tfield_l0d(tpfile%nmodel)%data
+  gsleve => tfieldlist(iid)%tfield_l0d(tzfile%nmodel)%data
 
-  if ( imi /= tpfile%nmodel ) then
+  if ( imi /= tzfile%nmodel ) then
     !This is necessary to have correct domain sizes (used by Gather_xxfield)
-    call Go_tomodel_ll( tpfile%nmodel, iresp )
+    call Go_tomodel_ll( tzfile%nmodel, iresp )
     gchangemodel = .true.
   end if
 else
@@ -1583,13 +1593,13 @@ else
   ystdnameprefix = 'projection'
 endif
 
-if ( Associated( tpfile%tncdims ) ) then
-  tzdim_ni   => tpfile%tncdims%tdims(NMNHDIM_NI)
-  tzdim_nj   => tpfile%tncdims%tdims(NMNHDIM_NJ)
-  tzdim_ni_u => tpfile%tncdims%tdims(NMNHDIM_NI_U)
-  tzdim_nj_u => tpfile%tncdims%tdims(NMNHDIM_NJ_U)
-  tzdim_ni_v => tpfile%tncdims%tdims(NMNHDIM_NI_V)
-  tzdim_nj_v => tpfile%tncdims%tdims(NMNHDIM_NJ_V)
+if ( Associated( tzfile%tncdims ) ) then
+  tzdim_ni   => tzfile%tncdims%tdims(NMNHDIM_NI)
+  tzdim_nj   => tzfile%tncdims%tdims(NMNHDIM_NJ)
+  tzdim_ni_u => tzfile%tncdims%tdims(NMNHDIM_NI_U)
+  tzdim_nj_u => tzfile%tncdims%tdims(NMNHDIM_NJ_U)
+  tzdim_ni_v => tzfile%tncdims%tdims(NMNHDIM_NI_V)
+  tzdim_nj_v => tzfile%tncdims%tdims(NMNHDIM_NJ_V)
 else
   tzdim_ni   => Null()
   tzdim_nj   => Null()
@@ -1616,9 +1626,9 @@ call Write_hor_coord1d( tzdim_nj_v, 'y-dimension of the grid at v location', &
 #if 0
 !Deallocate only if it is a non Z-split file or the last Z-split subfile
 gdealloc = .false.
-if ( Associated( tpfile%tmainfile ) ) then
-  if ( tpfile%cname == tpfile%tmainfile%tfiles_ioz(tpfile%tmainfile%nsubfiles_ioz)%tfile%cname ) gdealloc = .true.
-else if ( tpfile%nsubfiles_ioz == 0 .and. .not. Associated( tpfile%tmainfile ) ) then
+if ( Associated( tzfile%tmainfile ) ) then
+  if ( tzfile%cname == tzfile%tmainfile%tfiles_ioz(tzfile%tmainfile%nsubfiles_ioz)%tfile%cname ) gdealloc = .true.
+else if ( tzfile%nsubfiles_ioz == 0 .and. .not. Associated( tzfile%tmainfile ) ) then
   gdealloc = .true.
 end if
 #else
@@ -1630,7 +1640,7 @@ if ( .not. lcartesian ) then
   Allocate( zlat(iiu, iju), zlon(iiu, iju) )
 
   !If the file is a Z-split subfile, coordinates are already collected
-  if ( .not. associated( tpfile%tmainfile ) ) then
+  if ( .not. associated( tzfile%tmainfile ) ) then
     call Gather_hor_coord2d( zxhatm, zyhatm, zlatm_glob, zlonm_glob )
     call Gather_hor_coord2d( zxhat,  zyhatm, zlatu_glob, zlonu_glob )
     call Gather_hor_coord2d( zxhatm, zyhat,  zlatv_glob, zlonv_glob )
@@ -1653,78 +1663,78 @@ if ( .not. lcartesian ) then
   if ( gdealloc ) Deallocate( zlatm_glob, zlonm_glob, zlatu_glob, zlonu_glob, zlatv_glob, zlonv_glob, zlatf_glob, zlonf_glob )
 end if
 
-if ( tpfile%lmaster ) then !vertical coordinates in the transformed space are the same on all processes
+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
 
-    call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_LEVEL),  'position z in the transformed space',              '', &
+    call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_LEVEL),  'position z in the transformed space',              '', &
                           'altitude',                0.,  JPVEXT, JPVEXT, ZZHATM )
-    call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_LEVEL_W),'position z in the transformed space at w location','', &
+    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 )
   END IF
 END IF
 
 !Write time scale
-if ( tpfile%lmaster ) then !time scale is the same on all processes
+if ( tzfile%lmaster ) then !time scale is 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 ( tpfile%ctype /= 'MNHDIACHRONIC' ) then
+    if ( tzfile%ctype /= 'MNHDIACHRONIC' ) then
       if ( Present( tpdtmodeln ) ) then
-        call Write_time_coord( tpfile%tncdims%tdims(nmnhdim_time), 'time axis', [ tpdtmodeln ] )
+        call Write_time_coord( tzfile%tncdims%tdims(nmnhdim_time), 'time axis', [ tpdtmodeln ] )
       else if ( Associated( tdtcur ) ) then
-        call Write_time_coord( tpfile%tncdims%tdims(nmnhdim_time), 'time axis', [ tdtcur ] )
+        call Write_time_coord( tzfile%tncdims%tdims(nmnhdim_time), 'time axis', [ tdtcur ] )
       end if
     end if
   end if
 end if
 
-if ( tpfile%lmaster ) then
+if ( tzfile%lmaster ) then
   !Write coordinates used in diachronic files
-  if ( tpfile%ctype == 'MNHDIACHRONIC' ) then
+  if ( tzfile%ctype == 'MNHDIACHRONIC' ) then
     if ( cbutype == 'CART' .or. cbutype == 'SKIP' ) then
       !Coordinates for the budgets in cartesian boxes
       if ( .not. lbu_icp )                                                                                                  &
-        call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NI), 'x-dimension of the budget cartesian box',    &
+        call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NI), 'x-dimension of the budget cartesian box',    &
                          trim(ystdnameprefix)//'_x_coordinate', 'X', 0., 0, 0, zxhatm_glob(nbuil + jphext : nbuih + jphext) )
       if ( .not. lbu_jcp )                                                                                                  &
-        call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NJ), 'y-dimension of the budget cartesian box',    &
+        call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NJ), 'y-dimension of the budget cartesian box',    &
                          trim(ystdnameprefix)//'_y_coordinate', 'Y', 0., 0, 0, zyhatm_glob(nbujl + jphext : nbujh + jphext) )
       if ( .not. lbu_icp )                                                                    &
-        call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NI_U),               &
+        call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NI_U),               &
                                 'x-dimension of the budget cartesian box at u location',      &
                                 trim(ystdnameprefix)//'_x_coordinate_at_u_location',          &
                                 'X', -0.5, 0, 0, zxhat_glob (nbuil + jphext : nbuih + jphext) )
       if ( .not. lbu_jcp )                                                                    &
-        call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NJ_U),               &
+        call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NJ_U),               &
                                 'y-dimension of the budget cartesian box at u location',      &
                                 trim(ystdnameprefix)//'_y_coordinate_at_u_location',          &
                                 'Y', 0.,   0, 0, zyhatm_glob(nbujl + jphext : nbujh + jphext) )
       if ( .not. lbu_icp )                                                                    &
-        call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NI_V),               &
+        call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NI_V),               &
                                 'x-dimension of the budget cartesian box at v location',      &
                                 trim(ystdnameprefix)//'_x_coordinate_at_v_location',          &
                                 'X', 0.,   0, 0, zxhatm_glob(nbuil + jphext : nbuih + jphext) )
       if ( .not. lbu_jcp )                                                                    &
-        call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NJ_V),               &
+        call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_NJ_V),               &
                                 'y-dimension of the budget cartesian box at v location',      &
                                 trim(ystdnameprefix)//'_y_coordinate_at_v_location',          &
                                 'Y', -0.5, 0, 0, zyhat_glob (nbujl + jphext : nbujh + jphext) )
       if ( .not. lbu_kcp )                                                                                      &
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_LEVEL),                                  &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_LEVEL),                                  &
                               'position z in the transformed space of the budget cartesian box',                &
                               '', 'altitude',               0.,   0, 0, zzhatm(nbukl + JPVEXT : nbukh + JPVEXT) )
       if ( .not. lbu_kcp )                                                                                     &
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_CART_LEVEL_W),                               &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_CART_LEVEL_W),                               &
                               'position z in the transformed space at w location of the budget cartesian box', &
                               '', 'altitude_at_w_location', -0.5, 0, 0, zzhat (nbukl + JPVEXT : nbukh + JPVEXT) )
     else if ( cbutype == 'MASK' ) then
       !Coordinates for the budgets masks
       if ( nbukmax > 0 )                                                                                        &
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_MASK_LEVEL),                                  &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_MASK_LEVEL),                                  &
                               'position z in the transformed space of the budget mask',                         &
                               '', 'altitude',               0.,   0, 0, zzhatm(nbukl + JPVEXT : nbukh + JPVEXT) )
       if ( nbukmax > 0 )                                                                                        &
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_MASK_LEVEL_W),                                &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_MASK_LEVEL_W),                                &
                               'position z in the transformed space at w location of the budget mask',           &
                               '', 'altitude',               -0.5, 0, 0, zzhat (nbukl + JPVEXT : nbukh + JPVEXT) )
 
@@ -1753,7 +1763,7 @@ if ( tpfile%lmaster ) then
         tzdates_bound(2, jt)%xtime  = tdtseg%xtime + nbustep * jt * xtstep
       end do
 
-      call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_TIME), 'time axis for budgets', tzdates, tzdates_bound )
+      call Write_time_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_TIME), 'time axis for budgets', tzdates, tzdates_bound )
 
       Deallocate( tzdates_bound )
       Deallocate( tzdates )
@@ -1761,7 +1771,7 @@ if ( tpfile%lmaster ) then
 
     !Coordinates for the number of LES budget time samplings
     if ( nles_times > 0 ) &
-      call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_LES_TIME), 'time axis for LES budgets', tles_dates )
+      call Write_time_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_LES_TIME), 'time axis for LES budgets', tles_dates )
 
     !Coordinates for the number of LES budget time averages
     !Condition also on nles_times to not create this coordinate when not used (no time average if nles_times=0)
@@ -1793,7 +1803,7 @@ if ( tpfile%lmaster ) then
         tzdates_bound(2, jt)%nday   = tdtseg%nday
         tzdates_bound(2, jt)%xtime  = tdtseg%xtime + xles_times(jte)
       end do
-      call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_LES_AVG_TIME), 'time axis for LES budget time averages', &
+      call Write_time_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_LES_AVG_TIME), 'time axis for LES budget time averages', &
                              tzdates, tzdates_bound )
 
       Deallocate( tzdates_bound )
@@ -1807,12 +1817,12 @@ if ( tpfile%lmaster ) then
         do ji = 1, nles_k
           zles_levels(ji) = zzhatm(nles_levels(ji) + JPVEXT)
         end do
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_LES_LEVEL),           &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_LES_LEVEL),           &
                               'position z in the transformed space of the LES budgets', &
                               '', 'altitude', 0., 0, 0, zles_levels(:)                  )
         Deallocate( zles_levels )
       else if ( cles_level_type == 'Z' ) then
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_BUDGET_LES_LEVEL), 'altitude levels for the LES budgets', &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_BUDGET_LES_LEVEL), 'altitude levels for the LES budgets', &
                               'altitude', '', 0., 0, 0, xles_altitudes(1:nles_k) )
       else
         call Print_msg( NVERB_ERROR, 'IO', 'IO_Coordvar_write_nc4','invalid cles_level_type' )
@@ -1823,11 +1833,11 @@ if ( tpfile%lmaster ) then
 
     !Coordinates for the number of horizontal wavelengths for non-local LES budgets (2 points correlations)
     if ( nspectra_ni > 0 ) &
-      call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_SPECTRA_2PTS_NI), 'x-dimension of the LES budget cartesian box',    &
+      call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_SPECTRA_2PTS_NI), 'x-dimension of the LES budget cartesian box',    &
                             trim(ystdnameprefix)//'_x_coordinate', 'X', 0., 0, 0,                                              &
                             zxhatm_glob(nlesn_iinf(imi) + jphext : nlesn_isup(imi) + jphext) )
     if ( nspectra_nj > 0 .and. .not. l2d ) &
-      call Write_hor_coord1d( tpfile%tncdims%tdims(NMNHDIM_SPECTRA_2PTS_NJ), 'y-dimension of the LES budget cartesian box',    &
+      call Write_hor_coord1d( tzfile%tncdims%tdims(NMNHDIM_SPECTRA_2PTS_NJ), 'y-dimension of the LES budget cartesian box',    &
                               trim(ystdnameprefix)//'_y_coordinate', 'Y', 0., 0, 0,                                            &
                               zyhatm_glob(nlesn_jinf(imi) + jphext : nlesn_jsup(imi) + jphext) )
 
@@ -1841,12 +1851,12 @@ if ( tpfile%lmaster ) then
         do ji = 1, nspectra_k
           zspectra_levels(ji) = zzhatm(nspectra_levels(ji) + JPVEXT)
         end do
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_SPECTRA_LEVEL),                        &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_SPECTRA_LEVEL),                        &
                               'position z in the transformed space of the non-local LES budgets', &
                               '', 'altitude', 0., 0, 0, zspectra_levels(:)                        )
         Deallocate( zspectra_levels )
       else if ( cspectra_level_type == 'Z' ) then
-        call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_SPECTRA_LEVEL), 'altitude levels for the non-local LES budgets', &
+        call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_SPECTRA_LEVEL), 'altitude levels for the non-local LES budgets', &
                               'altitude', '', 0., 0, 0, xspectra_altitudes(1 : nspectra_k) )
       else
         call Print_msg( NVERB_ERROR, 'IO', 'IO_Coordvar_write_nc4','invalid cspectra_level_type' )
@@ -1855,21 +1865,21 @@ if ( tpfile%lmaster ) then
 
     !Coordinates for the number of profiler times
     if ( lprofiler ) &
-      call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_PROFILER_TIME), 'time axis for profilers', tprofilers_time%tpdates )
+      call Write_time_coord( tzfile%tncdims%tdims(NMNHDIM_PROFILER_TIME), 'time axis for profilers', tprofilers_time%tpdates )
 
     !Coordinates for the number of station times
     if ( lstation ) &
-      call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_STATION_TIME), 'time axis for stations', tstations_time%tpdates )
+      call Write_time_coord( tzfile%tncdims%tdims(NMNHDIM_STATION_TIME), 'time axis for stations', tstations_time%tpdates )
 
     !Dimension for the number of series times
     if ( lseries .and. nsnbstept > 0 ) then
-      call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_SERIES_LEVEL),                   &
+      call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_SERIES_LEVEL),                   &
                             'position z in the transformed space of the temporal series', &
                             '', 'altitude', 0., 0, 0, zzhatm(1 + JPVEXT : nkmax + JPVEXT) )
-      call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_SERIES_LEVEL_W),                                 &
+      call Write_ver_coord( tzfile%tncdims%tdims(NMNHDIM_SERIES_LEVEL_W),                                 &
                             'position z in the transformed space at w location of the temporal series',   &
                             '', 'altitude_at_w_location', -0.5, 0, 0, zzhat (1 + JPVEXT : nkmax + JPVEXT) )
-      call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_SERIES_TIME), 'time axis for temporal series', tpsdates )
+      call Write_time_coord( tzfile%tncdims%tdims(NMNHDIM_SERIES_TIME), 'time axis for temporal series', tpsdates )
     end if
 
     if ( lflyer ) then
@@ -1927,10 +1937,10 @@ subroutine Gather_hor_coord1d( haxis, pcoords_loc, pcoords_glob )
   ! Allocate pcoords_glob
   if ( gsmonoproc ) then ! sequential execution
     allocate( pcoords_glob( size( pcoords_loc) ) )
-  else if ( tpfile%nsubfiles_ioz > 0 ) then
+  else if ( tzfile%nsubfiles_ioz > 0 ) then
     !If there are Z-split subfiles, all subfile writers need the coordinates
     call Allocbuffer_ll( pcoords_glob, pcoords_loc, ydir, galloc )
-  else if ( .not. tpfile%lmaster ) then
+  else if ( .not. tzfile%lmaster ) then
     allocate( pcoords_glob(0 ) ) !to prevent false positive with valgrind
   else !Master process
     call Allocbuffer_ll( pcoords_glob, pcoords_loc, ydir, galloc )
@@ -1940,12 +1950,12 @@ subroutine Gather_hor_coord1d( haxis, pcoords_loc, pcoords_glob )
   if ( gsmonoproc ) then ! sequential execution
       pcoords_glob(: ) = pcoords_loc(: )
   else ! multiprocesses execution
-      call Gather_xxfield( ydir, pcoords_loc, pcoords_glob, tpfile%nmaster_rank, tpfile%nmpicomm )
+      call Gather_xxfield( ydir, pcoords_loc, pcoords_glob, tzfile%nmaster_rank, tzfile%nmpicomm )
   endif
 
   !If the file has Z-split subfiles, broadcast the coordinates to all processes
-  if ( tpfile%nsubfiles_ioz > 0 ) &
-    call MPI_BCAST( pcoords_glob, size( pcoords_glob ), MNHREAL_MPI, tpfile%nmaster_rank - 1,  tpfile%nmpicomm, ierr )
+  if ( tzfile%nsubfiles_ioz > 0 ) &
+    call MPI_BCAST( pcoords_glob, size( pcoords_glob ), MNHREAL_MPI, tzfile%nmaster_rank - 1,  tzfile%nmpicomm, ierr )
 
 end subroutine Gather_hor_coord1d
 #endif
@@ -1972,11 +1982,11 @@ subroutine Gather_hor_coord2d( px, py, plat_glob, plon_glob )
   if ( gsmonoproc ) then ! sequential execution
     allocate( plat_glob( size( zlat, 1 ), size( zlat, 2 ) ) )
     allocate( plon_glob( size( zlon, 1 ), size( zlon, 2 ) ) )
-  else if ( tpfile%nsubfiles_ioz > 0 ) then
+  else if ( tzfile%nsubfiles_ioz > 0 ) then
     !If there are Z-split subfiles, all subfile writers need the coordinates
     call Allocbuffer_ll( plat_glob, zlat, 'XY', galloc1 )
     call Allocbuffer_ll( plon_glob, zlon, 'XY', galloc2 )
-  else if ( .not. tpfile%lmaster ) then
+  else if ( .not. tzfile%lmaster ) then
     allocate( plat_glob( 0, 0 ), plon_glob( 0, 0 ) ) !to prevent false positive with valgrind
   else !Master process
     call Allocbuffer_ll( plat_glob, zlat, 'XY', galloc1 )
@@ -1988,14 +1998,14 @@ subroutine Gather_hor_coord2d( px, py, plat_glob, plon_glob )
       plat_glob(:, : ) = zlat(:, : )
       plon_glob(:, : ) = zlon(:, : )
   else ! multiprocesses execution
-      call Gather_xyfield( zlat, plat_glob, tpfile%nmaster_rank, tpfile%nmpicomm )
-      call Gather_xyfield( zlon, plon_glob, tpfile%nmaster_rank, tpfile%nmpicomm )
+      call Gather_xyfield( zlat, plat_glob, tzfile%nmaster_rank, tzfile%nmpicomm )
+      call Gather_xyfield( zlon, plon_glob, tzfile%nmaster_rank, tzfile%nmpicomm )
   endif
 
   !If the file has Z-split subfiles, broadcast the coordinates to all processes
-  if ( tpfile%nsubfiles_ioz > 0 ) then
-    call MPI_BCAST( plat_glob, size( plat_glob ), MNHREAL_MPI, tpfile%nmaster_rank - 1,  tpfile%nmpicomm, ierr )
-    call MPI_BCAST( plon_glob, size( plon_glob ), MNHREAL_MPI, tpfile%nmaster_rank - 1,  tpfile%nmpicomm, ierr )
+  if ( tzfile%nsubfiles_ioz > 0 ) then
+    call MPI_BCAST( plat_glob, size( plat_glob ), MNHREAL_MPI, tzfile%nmaster_rank - 1,  tzfile%nmpicomm, ierr )
+    call MPI_BCAST( plon_glob, size( plon_glob ), MNHREAL_MPI, tzfile%nmaster_rank - 1,  tzfile%nmpicomm, ierr )
   end if
 end subroutine Gather_hor_coord2d
 
@@ -2021,7 +2031,7 @@ subroutine Write_hor_coord1d(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUNDLOW,KBOU
   INTEGER(KIND=CDFINT)          :: IVDIM
   INTEGER(KIND=CDFINT)          :: ISTATUS
 
-  IF (TPFILE%LMASTER) THEN
+  IF (tzfile%LMASTER) THEN
     isize    = tdim%nlen
     yvarname = Trim( tdim%cname )
     ivdim    = tdim%nid
@@ -2072,11 +2082,11 @@ subroutine Write_hor_coord2d( plat, plon, hlat, hlon )
   character(len=*),    intent(in) :: hlat
   character(len=*),    intent(in) :: hlon
 
-  if ( tpfile%lmaster ) then
+  if ( tzfile%lmaster ) then
     call Find_field_id_from_mnhname( hlat, id, iresp )
-    call IO_Field_write_nc4_x2( tpfile, tfieldlist(id ), plat, iresp, oiscoord = .true. )
+    call IO_Field_write_nc4_x2( tzfile, tfieldlist(id ), plat, iresp, oiscoord = .true. )
     call Find_field_id_from_mnhname( hlon, id, iresp )
-    call IO_Field_write_nc4_x2( tpfile, tfieldlist(id ), plon, iresp, oiscoord = .true. )
+    call IO_Field_write_nc4_x2( tzfile, tfieldlist(id ), plon, iresp, oiscoord = .true. )
   end if
 end subroutine Write_hor_coord2d
 
@@ -2257,7 +2267,7 @@ subroutine Write_time_coord( tdim, hlongname, tpdates, tpdates_bound )
     istatus = NF90_INQ_VARID( incid, yvarname, ivarid )
     if ( istatus /= NF90_NOERR ) then
       ! Define the coordinate variable
-      ivdims(1) = tpfile%tncdims%tdims(NMNHDIM_PAIR)%nid
+      ivdims(1) = tzfile%tncdims%tdims(NMNHDIM_PAIR)%nid
       ivdims(2) = tdim%nid
       istatus = NF90_DEF_VAR( incid, yvarname, MNHREAL_NF90, ivdims, ivarid )
       if ( istatus /= NF90_NOERR ) call IO_Err_handle_nc4( istatus, 'Write_time_coord', 'NF90_DEF_VAR', Trim( yvarname ) )
@@ -2306,10 +2316,10 @@ subroutine Write_flyer_time_coord( tpflyer )
   if ( Count( tpflyer%xx /= XUNDEF) > 1 ) then
     Allocate( tzdim )
 
-    istatus = NF90_INQ_NCID( tpfile%nncid, 'Flyers', icatid )
+    istatus = NF90_INQ_NCID( tzfile%nncid, 'Flyers', icatid )
     if ( istatus /= NF90_NOERR ) then
       call Print_msg( NVERB_ERROR, 'IO', 'Write_flyer_time_coord', &
-                      Trim( tpfile%cname ) // ': group Flyers not found' )
+                      Trim( tzfile%cname ) // ': group Flyers not found' )
     end if
 
     call Aircraft_balloon_longtype_get( tpflyer, ytype )
@@ -2317,19 +2327,19 @@ subroutine Write_flyer_time_coord( tpflyer )
     istatus = NF90_INQ_NCID( icatid, Trim( ytype_clean ), isubcatid )
     if ( istatus /= NF90_NOERR ) then
       call Print_msg( NVERB_ERROR, 'IO', 'Write_flyer_time_coord', &
-                      Trim( tpfile%cname ) // ': group ' // Trim( ytype_clean ) // ' not found' )
+                      Trim( tzfile%cname ) // ': group ' // Trim( ytype_clean ) // ' not found' )
     end if
 
     istatus = NF90_INQ_NCID( isubcatid, Trim( tpflyer%cname ), incid )
     if ( istatus /= NF90_NOERR ) then
       call Print_msg( NVERB_ERROR, 'IO', 'Write_flyer_time_coord', &
-                      Trim( tpfile%cname ) // ': group '// Trim( tpflyer%cname ) // ' not found' )
+                      Trim( tzfile%cname ) // ': group '// Trim( tpflyer%cname ) // ' not found' )
     end if
 
     istatus = NF90_INQ_DIMID( incid, 'time_flyer', idimid )
     if ( istatus /= NF90_NOERR ) then
       call Print_msg( NVERB_ERROR, 'IO', 'Write_flyer_time_coord', &
-                      Trim( tpfile%cname ) // ': group ' // Trim( tpflyer%cname ) // ' time_flyer dimension not found' )
+                      Trim( tzfile%cname ) // ': group ' // Trim( tpflyer%cname ) // ' time_flyer dimension not found' )
     end if
 
     tzdim%cname = 'time_flyer'
@@ -2342,7 +2352,7 @@ subroutine Write_flyer_time_coord( tpflyer )
     Deallocate( tzdim )
 
     !Restore file identifier to root group
-    incid = tpfile%nncid
+    incid = tzfile%nncid
   end if
   end if
 
diff --git a/src/LIB/hdf5-1.14.0.tar.gz b/src/LIB/hdf5-1.14.0.tar.gz
deleted file mode 100644
index 14467a66342048527655995b2cd28292bcbc3272..0000000000000000000000000000000000000000
--- a/src/LIB/hdf5-1.14.0.tar.gz
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:a571cc83efda62e1a51a0a912dd916d01895801c5025af91669484a1575a6ef4
-size 19285771
diff --git a/src/LIB/hdf5-1.14.2.tar.gz b/src/LIB/hdf5-1.14.2.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..2a264607b301231525dd537b0087a0000f8cedf9
--- /dev/null
+++ b/src/LIB/hdf5-1.14.2.tar.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1c342e634008284a8c2794c8e7608e2eaf26d01d445fb3dfd7f33cb2fb51ac53
+size 19815645
diff --git a/src/LIB/libaec-1.1.2.tar.gz b/src/LIB/libaec-1.1.2.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..03a033ad9eb4d6ca7632218b31becc11638d70ad
--- /dev/null
+++ b/src/LIB/libaec-1.1.2.tar.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a5d03a242468014b92b1b196de397f1e9a7927782db1a0423ea01db78421e373
+size 3138136
diff --git a/src/LIB/libaec-v1.0.6.tar.gz b/src/LIB/libaec-v1.0.6.tar.gz
deleted file mode 100644
index 79f28267131696d2d501ece5120aae3278233ea4..0000000000000000000000000000000000000000
--- a/src/LIB/libaec-v1.0.6.tar.gz
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:c5409f7903dd84e75fb48d8d5e603b2a06e82737a23ac86297bedd6218c20d43
-size 3629633
diff --git a/src/LIB/netcdf-c-4.9.0.tar.gz b/src/LIB/netcdf-c-4.9.0.tar.gz
deleted file mode 100644
index 6ee47e081789963c760d947169817a78c299d41b..0000000000000000000000000000000000000000
--- a/src/LIB/netcdf-c-4.9.0.tar.gz
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:4c956022b79c08e5e14eee8df51b13c28e6121c2b7e7faadc21b375949400b49
-size 7103958
diff --git a/src/LIB/netcdf-c-4.9.2.tar.gz b/src/LIB/netcdf-c-4.9.2.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..7d35658a6e6b13bb086bc46364b7b6adfc3fd771
--- /dev/null
+++ b/src/LIB/netcdf-c-4.9.2.tar.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cf11babbbdb9963f09f55079e0b019f6d0371f52f8e1264a5ba8e9fdab1a6c48
+size 7142536
diff --git a/src/LIB/netcdf-fortran-4.6.0.tar.gz b/src/LIB/netcdf-fortran-4.6.0.tar.gz
deleted file mode 100644
index acc4addf308e9f42bfcaeac1876fa79edb2342c9..0000000000000000000000000000000000000000
--- a/src/LIB/netcdf-fortran-4.6.0.tar.gz
+++ /dev/null
@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:198bff6534cc85a121adc9e12f1c4bc53406c403bda331775a1291509e7b2f23
-size 1110214
diff --git a/src/LIB/netcdf-fortran-4.6.1.tar.gz b/src/LIB/netcdf-fortran-4.6.1.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..022946e1c589e3bb3f08623bc0bff4b5bd1ff982
--- /dev/null
+++ b/src/LIB/netcdf-fortran-4.6.1.tar.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b50b0c72b8b16b140201a020936aa8aeda5c79cf265c55160986cd637807a37a
+size 1112270
diff --git a/src/MNH/modd_bakout.f90 b/src/MNH/modd_bakout.f90
index 487cd44265ae2ff81768fae25de1bc3549e8a078..92bee25c9396dc96d5d2d81097bc8152d6ae1e0a 100644
--- a/src/MNH/modd_bakout.f90
+++ b/src/MNH/modd_bakout.f90
@@ -1,6 +1,6 @@
-!MNH_LIC Copyright 1996-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1996-2023 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 !     ##################
@@ -16,7 +16,7 @@
 !       (compression and precision reduction) of the outputs and backups realized by
 !       all models.
 !       Introduced to facilitate the output FM-file managment in case of nesting
-!     
+!
 !!**  IMPLICIT ARGUMENTS
 !!    ------------------
 !!
@@ -26,15 +26,16 @@
 !!    REFERENCE
 !!    ---------
 !!      Book2 of Meso-NH documentation (module MODD_BAKOUT)
-!!          
+!!
 !!    AUTHOR
 !!    ------
 !!	J.P. Lafore      *Meteo France*
 !!
 !!    MODIFICATIONS
 !!    -------------
-!!      Original    26/07/96                     
-!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!!      Original    26/07/96
+!  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
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -44,31 +45,43 @@ USE MODD_PARAMETERS
 !
 IMPLICIT NONE
 !
-LOGICAL,SAVE :: LBAK_BEG = .FALSE. ! Force a backup/output at the first timestep
-LOGICAL,SAVE :: LOUT_BEG = .FALSE. ! of the segment for all models
-LOGICAL,SAVE :: LBAK_END = .FALSE. ! Force a backup/output at the last timestep
-LOGICAL,SAVE :: LOUT_END = .FALSE. ! of the segment for all models
-LOGICAL,SAVE,DIMENSION(JPMODELMAX) :: LOUT_REDUCE_FLOAT_PRECISION = .FALSE.
+SAVE
+!
+INTEGER, PARAMETER :: NCOMPRNAMELGTMAX  = 10 ! Maximum length of compression algorithm name
+!
+LOGICAL :: LBAK_BEG = .FALSE. ! Force a backup/output at the first timestep
+LOGICAL :: LOUT_BEG = .FALSE. ! of the segment for all models
+LOGICAL :: LBAK_END = .FALSE. ! Force a backup/output at the last timestep
+LOGICAL :: LOUT_END = .FALSE. ! of the segment for all models
+
+! Compression
+LOGICAL, DIMENSION(JPMODELMAX) :: LOUT_REDUCE_FLOAT_PRECISION = .FALSE.
 ! Reduce the precision of floats to single precision instead of double precision (for netCDF)
-LOGICAL,SAVE,DIMENSION(JPMODELMAX) :: LOUT_COMPRESS = .FALSE. ! Compress (float) arrays (for netCDF)
-INTEGER,SAVE,DIMENSION(JPMODELMAX) :: NOUT_COMPRESS_LEVEL = 4 ! Compression level (for netCDF)
-REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)  ::   XBAK_TIME, XOUT_TIME
-! XBAK_TIME(m,i) / XOUT_TIME(m,i) array of 
+LOGICAL, DIMENSION(JPMODELMAX) :: LOUT_COMPRESS = .FALSE. ! Compress (float) arrays (for netCDF)
+INTEGER, DIMENSION(JPMODELMAX) :: NOUT_COMPRESS_LEVEL = 4 ! Compression level (for netCDF)
+
+! Lossy compression
+LOGICAL, DIMENSION(JPMODELMAX) :: LOUT_COMPRESS_LOSSY = .FALSE. ! Lossy compression of (float) arrays (for netCDF)
+CHARACTER(LEN=NCOMPRNAMELGTMAX), DIMENSION(JPMODELMAX) :: COUT_COMPRESS_LOSSY_ALGO = 'GRANULARBR' ! Quantization algorithm
+INTEGER, DIMENSION(JPMODELMAX) :: NOUT_COMPRESS_LOSSY_NSD  = 3  ! Number of Significant Digits (or Bits)
+
+REAL, ALLOCATABLE, DIMENSION(:,:)  ::   XBAK_TIME, XOUT_TIME
+! XBAK_TIME(m,i) / XOUT_TIME(m,i) array of
 ! the increments in seconds from the beginning of the segment to the
 ! instant where the i-th fields output on FM-files is realized by model "m"
-INTEGER,SAVE,ALLOCATABLE,DIMENSION(:,:)  ::   NBAK_STEP, NOUT_STEP
+INTEGER, ALLOCATABLE, DIMENSION(:,:)  ::   NBAK_STEP, NOUT_STEP
 ! NBAK_STEP(m,i) / NOUT_STEP(m,i) array of
 ! the increments in steps from the beginning of the segment to the
 ! step where the i-th fields output on FM-files is realized by model "m"
-INTEGER,SAVE,DIMENSION(JPMODELMAX) :: NBAK_STEP_FREQ = NNEGUNDEF, NOUT_STEP_FREQ = NNEGUNDEF
+INTEGER, DIMENSION(JPMODELMAX) :: NBAK_STEP_FREQ = NNEGUNDEF, NOUT_STEP_FREQ = NNEGUNDEF
 ! Number of timesteps between 2 backups/outputs for each model
-INTEGER,SAVE,DIMENSION(JPMODELMAX) :: NBAK_STEP_FREQ_FIRST = 1, NOUT_STEP_FREQ_FIRST = 1
+INTEGER, DIMENSION(JPMODELMAX) :: NBAK_STEP_FREQ_FIRST = 1, NOUT_STEP_FREQ_FIRST = 1
 ! First timestep numbers between 2 backups/outputs for each model (if NBAK/OUT_STEP_FREQ is set)
-REAL,SAVE,DIMENSION(JPMODELMAX) :: XBAK_TIME_FREQ = XNEGUNDEF, XOUT_TIME_FREQ = XNEGUNDEF
+REAL, DIMENSION(JPMODELMAX) :: XBAK_TIME_FREQ = XNEGUNDEF, XOUT_TIME_FREQ = XNEGUNDEF
 ! Time between 2 backups/outputs for each model
-REAL,SAVE,DIMENSION(JPMODELMAX) :: XBAK_TIME_FREQ_FIRST = 0., XOUT_TIME_FREQ_FIRST = 0.
+REAL, DIMENSION(JPMODELMAX) :: XBAK_TIME_FREQ_FIRST = 0., XOUT_TIME_FREQ_FIRST = 0.
 ! Time for first backup/output for each model (if XBAK/OUT_TIME_FREQ is set)
-CHARACTER(LEN=NMNHNAMELGTMAX),SAVE,ALLOCATABLE,DIMENSION(:,:) :: COUT_VAR ! Name of the fields to output
+CHARACTER(LEN=NMNHNAMELGTMAX), ALLOCATABLE, DIMENSION(:,:) :: COUT_VAR ! Name of the fields to output
 !
 !Directory names for backups/outputs
 CHARACTER(LEN=NDIRNAMELGTMAX) :: CBAK_DIR='', COUT_DIR=''
diff --git a/src/MNH/modn_output.f90 b/src/MNH/modn_output.f90
index 0ac56ebddb71d76fd4159373bbf140ad6ac062ca..64429041d9ca7aa91ec9767e7ef55c72663f731a 100644
--- a/src/MNH/modn_output.f90
+++ b/src/MNH/modn_output.f90
@@ -1,6 +1,6 @@
-!MNH_LIC Copyright 1996-2017 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1996-2023 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 !     ##################
@@ -24,16 +24,17 @@
 !!    REFERENCE
 !!    ---------
 !!      Book2 of Meso-NH documentation (module MODD_BAKOUT)
-!!          
+!!
 !!    AUTHOR
 !!    ------
 !!	J.P. Lafore      *Meteo France*
 !!
 !!    MODIFICATIONS
 !!    -------------
-!!      Original    26/07/96                      
-!!      Ph. Wautelet 2016       new structures for outputs/backups
-!!      Ph. Wautelet 02/10/2017 split NAM_OUTPUT in NAM_BACKUP and NAM_OUTPUT
+!!      Original    26/07/96
+!  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
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -50,6 +51,7 @@ NAMELIST/NAM_OUTPUT/LOUT_BEG,LOUT_END,&
                    COUT_VAR, &
                    LOUT_REDUCE_FLOAT_PRECISION, &
                    LOUT_COMPRESS, NOUT_COMPRESS_LEVEL,&
+                   LOUT_COMPRESS_LOSSY, COUT_COMPRESS_LOSSY_ALGO, NOUT_COMPRESS_LOSSY_NSD, &
                    COUT_DIR
 !
 END MODULE MODN_OUTPUT
diff --git a/src/Makefile b/src/Makefile
index b13e1975e741ee65bd1cf087ad21d8916d83d71b..d903b69fdec9669451f95961dc7fafd4030c6d96 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -352,7 +352,7 @@ $(CDF_MOD) :
 	cd ${DIR_HDF} && ./configure --enable-fortran --disable-shared --prefix=${CDF_PATH} --libdir=${CDF_PATH}/lib64 --with-szlib=${CDF_PATH}/include,${CDF_PATH}/lib64 \
 	CC="$(CC)" CFLAGS="$(HDF_OPT)" ${HDF_CONF} FC="$(FC)" FCFLAGS="$(NETCDF_OPT)" LDFLAGS="-L${CDF_PATH}/lib64" LIBS="-lsz -laec -lz" && \
 	$(MAKE) && $(MAKE) install && $(MAKE) clean
-	cd ${DIR_CDFC} && ./configure --disable-shared --prefix=${CDF_PATH} --libdir=${CDF_PATH}/lib64 --disable-dap \
+	cd ${DIR_CDFC} && ./configure --disable-shared --prefix=${CDF_PATH} --libdir=${CDF_PATH}/lib64 --disable-dap --disable-byterange \
 	CC="$(CC)" CFLAGS="$(NETCDF_OPT)" CPPFLAGS="${INC_NETCDF}" ${CDF_CONF} LDFLAGS="-L${CDF_PATH}/lib64" LIBS="-lhdf5_hl -lhdf5 -lsz -laec -lz -ldl " && \
 	$(MAKE) && $(MAKE) install && $(MAKE) clean
 ifdef MNH_FOREFIRE
diff --git a/src/configure b/src/configure
index 1c765c1827274ca722e16863ece4c15ba4b2ed10..1ae2a41366d3ea557cf211a90ec137a42e11e1b4 100755
--- a/src/configure
+++ b/src/configure
@@ -12,11 +12,11 @@ export VERSION_MASTER=${VERSION_MASTER:-MNH-V5-6}
 export VERSION_BUG=${VERSION_BUG:-2}
 export VERSION_XYZ=${VERSION_XYZ:-${VERSION_MASTER}-${VERSION_BUG}${VER_OASIS:+-${VER_OASIS}}}
 export VERSION_DATE=${VERSION_DATE:-"28/11/2023"}
-export VERSION_LIBAEC=${VERSION_LIBAEC:-"v1.0.6"}
-export VERSION_HDF=${VERSION_HDF:-"1.14.0"}
-export VERSION_CDFC=${VERSION_CDFC:-"4.9.0"}
+export VERSION_LIBAEC=${VERSION_LIBAEC:-"1.1.2"}
+export VERSION_HDF=${VERSION_HDF:-"1.14.2"}
+export VERSION_CDFC=${VERSION_CDFC:-"4.9.2"}
 export VERSION_CDFCXX=${VERSION_CDFCXX:-"4.3.1"}
-export VERSION_CDFF=${VERSION_CDFF:-"4.6.0"}
+export VERSION_CDFF=${VERSION_CDFF:-"4.6.1"}
 export VERSION_GRIBAPI=${VERSION_GRIBAPI:-"1.26.0-Source"}
 export VERSION_ECCODES=${VERSION_ECCODES:-"2.18.0"}
 export ECCODES_DEFINITION_PATH=${ECCODES_DEFINITION_PATH:-${SRC_MESONH}/src/LIB/eccodes-${VERSION_ECCODES}/definitions/}