diff --git a/src/LIB/SURCOUCHE/src/modd_field.f90 b/src/LIB/SURCOUCHE/src/modd_field.f90
index 4133b8d7b91099d5e2c5ac72c59e0faa130fffe5..607ea1a695ad5588de2928299be1a2565a872f18 100644
--- a/src/LIB/SURCOUCHE/src/modd_field.f90
+++ b/src/LIB/SURCOUCHE/src/modd_field.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2016-2023 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2016-2024 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.
@@ -16,6 +16,7 @@
 !  P. Wautelet 08/10/2021: add 2 new dimensions: LW_bands (NMNHDIM_NLWB) and SW_bands (NMNHDIM_NSWB)
 !  P. Wautelet 14/10/2021: dynamically allocate tfieldlist (+ reallocate if necessary)
 !  P. Wautelet 04/11/2021: add TFIELDMETADATA type
+!  P. Wautelet 15/02/2024: add time dimension for Lagrangian trajectories
 !-----------------------------------------------------------------
 module modd_field
 
@@ -51,49 +52,51 @@ integer, parameter :: NMNHDIM_ONE                 = 10
 integer, parameter :: NMNHDIM_NSWB                = 11
 integer, parameter :: NMNHDIM_NLWB                = 12
 
-integer, parameter :: NMNHDIM_LASTDIM_NODIACHRO   = 12 ! Index of the last defined dimension for non-diachronic files
+integer, PARAMETER :: NMNHDIM_TRAJ_TIME           = 13
 
-integer, parameter :: NMNHDIM_COMPLEX             = 13
+integer, parameter :: NMNHDIM_LASTDIM_NODIACHRO   = 13 ! Index of the last defined dimension for non-diachronic files
 
-integer, parameter :: NMNHDIM_BUDGET_CART_NI      = 14
-integer, parameter :: NMNHDIM_BUDGET_CART_NJ      = 15
-integer, parameter :: NMNHDIM_BUDGET_CART_NI_U    = 16
-integer, parameter :: NMNHDIM_BUDGET_CART_NJ_U    = 17
-integer, parameter :: NMNHDIM_BUDGET_CART_NI_V    = 18
-integer, parameter :: NMNHDIM_BUDGET_CART_NJ_V    = 19
-integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL   = 20
-integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL_W = 21
+integer, parameter :: NMNHDIM_COMPLEX             = 14
 
-integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL   = 22
-integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL_W = 23
-integer, parameter :: NMNHDIM_BUDGET_MASK_NBUMASK = 24
+integer, parameter :: NMNHDIM_BUDGET_CART_NI      = 15
+integer, parameter :: NMNHDIM_BUDGET_CART_NJ      = 16
+integer, parameter :: NMNHDIM_BUDGET_CART_NI_U    = 17
+integer, parameter :: NMNHDIM_BUDGET_CART_NJ_U    = 18
+integer, parameter :: NMNHDIM_BUDGET_CART_NI_V    = 19
+integer, parameter :: NMNHDIM_BUDGET_CART_NJ_V    = 20
+integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL   = 21
+integer, parameter :: NMNHDIM_BUDGET_CART_LEVEL_W = 22
 
-integer, parameter :: NMNHDIM_BUDGET_TIME         = 25
+integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL   = 23
+integer, parameter :: NMNHDIM_BUDGET_MASK_LEVEL_W = 24
+integer, parameter :: NMNHDIM_BUDGET_MASK_NBUMASK = 25
 
-integer, parameter :: NMNHDIM_BUDGET_LES_TIME     = 26
-integer, parameter :: NMNHDIM_BUDGET_LES_AVG_TIME = 27
-integer, parameter :: NMNHDIM_BUDGET_LES_LEVEL    = 28
-integer, parameter :: NMNHDIM_BUDGET_LES_SV       = 29
-integer, parameter :: NMNHDIM_BUDGET_LES_PDF      = 30
+integer, parameter :: NMNHDIM_BUDGET_TIME         = 26
+
+integer, parameter :: NMNHDIM_BUDGET_LES_TIME     = 27
+integer, parameter :: NMNHDIM_BUDGET_LES_AVG_TIME = 28
+integer, parameter :: NMNHDIM_BUDGET_LES_LEVEL    = 29
+integer, parameter :: NMNHDIM_BUDGET_LES_SV       = 30
+integer, parameter :: NMNHDIM_BUDGET_LES_PDF      = 31
 integer, parameter :: NMNHDIM_BUDGET_LES_MASK     = 100 ! This is not a true dimension
 
-integer, parameter :: NMNHDIM_SPECTRA_2PTS_NI     = 31
-integer, parameter :: NMNHDIM_SPECTRA_2PTS_NJ     = 32
-integer, parameter :: NMNHDIM_SPECTRA_SPEC_NI     = 33
-integer, parameter :: NMNHDIM_SPECTRA_SPEC_NJ     = 34
-integer, parameter :: NMNHDIM_SPECTRA_LEVEL       = 35
+integer, parameter :: NMNHDIM_SPECTRA_2PTS_NI     = 32
+integer, parameter :: NMNHDIM_SPECTRA_2PTS_NJ     = 33
+integer, parameter :: NMNHDIM_SPECTRA_SPEC_NI     = 34
+integer, parameter :: NMNHDIM_SPECTRA_SPEC_NJ     = 35
+integer, parameter :: NMNHDIM_SPECTRA_LEVEL       = 36
 
-integer, parameter :: NMNHDIM_SERIES_LEVEL        = 36
-integer, parameter :: NMNHDIM_SERIES_LEVEL_W      = 37
-integer, parameter :: NMNHDIM_SERIES_TIME         = 38  ! Time dimension for time series
+integer, parameter :: NMNHDIM_SERIES_LEVEL        = 37
+integer, parameter :: NMNHDIM_SERIES_LEVEL_W      = 38
+integer, parameter :: NMNHDIM_SERIES_TIME         = 39  ! Time dimension for time series
 
-integer, parameter :: NMNHDIM_FLYER_TIME          = 39  ! Time dimension for aircraft/balloon (dimension local to each flyer)
-integer, parameter :: NMNHDIM_PROFILER_TIME       = 40  ! Time dimension for profilers
-integer, parameter :: NMNHDIM_STATION_TIME        = 41  ! Time dimension for stations
+integer, parameter :: NMNHDIM_FLYER_TIME          = 40  ! Time dimension for aircraft/balloon (dimension local to each flyer)
+integer, parameter :: NMNHDIM_PROFILER_TIME       = 41  ! Time dimension for profilers
+integer, parameter :: NMNHDIM_STATION_TIME        = 42  ! Time dimension for stations
 
-integer, parameter :: NMNHDIM_PAIR                = 42  ! For values coming by pair (ie boundaries)
+integer, parameter :: NMNHDIM_PAIR                = 43  ! For values coming by pair (ie boundaries)
 
-integer, parameter :: NMNHDIM_LASTDIM_DIACHRO     = 42  ! Index of the last defined dimension for diachronic files
+integer, parameter :: NMNHDIM_LASTDIM_DIACHRO     = 43  ! Index of the last defined dimension for diachronic files
 
 integer, parameter :: NMNHDIM_BUDGET_NGROUPS      = 101 ! This is not a true dimension
 integer, parameter :: NMNHDIM_FLYER_PROC          = 102 ! This is not a true dimension
diff --git a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
index 7d98869b43caecd3e2c356e92a4d7dbe73ae42a2..2a1c3f8d54d994dec378dd1accf523939f943623 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2023 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2024 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.
@@ -20,6 +20,7 @@
 !  P. Wautelet 04/05/2021: improve IO_Vdims_fill_nc4 if l2d and lpack
 !  P. Wautelet 27/05/2021: improve IO_Mnhname_clean to autocorrect names to be CF compliant
 !  P. Wautelet 08/10/2021: add 2 new dimensions: LW_bands and SW_bands
+!  P. Wautelet 15/02/2024: add time dimension for Lagrangian trajectories
 !-----------------------------------------------------------------
 #ifdef MNH_IOCDF4
 module mode_io_tools_nc4
@@ -250,12 +251,13 @@ use modd_budget,        only: cbutype, lbu_icp, lbu_jcp, lbu_kcp, nbuimax_ll, nb
 use modd_lbc_n,         only: clbcx, clbcy
 USE MODD_CONF,          ONLY: CPROGRAM, l2d, lpack
 USE MODD_CONF_n,        ONLY: CSTORAGE_TYPE
+USE MODD_DIAG_FLAG,     ONLY: LTRAJ
 USE MODD_DIM_n,         ONLY: NIMAX_ll, NJMAX_ll, NKMAX
 use modd_dyn,           only: xseglen
 use modd_dyn_n,         only: dyn_model
 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_ONE,  NMNHDIM_NSWB, NMNHDIM_NLWB, NMNHDIM_COMPLEX,                        &
+                              NMNHDIM_ONE,  NMNHDIM_NSWB, NMNHDIM_NLWB, NMNHDIM_TRAJ_TIME, NMNHDIM_COMPLEX,     &
                               NMNHDIM_BUDGET_CART_NI, NMNHDIM_BUDGET_CART_NJ, NMNHDIM_BUDGET_CART_NI_U,         &
                               NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NI_V, NMNHDIM_BUDGET_CART_NJ_V,     &
                               NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W,                           &
@@ -281,6 +283,7 @@ use modd_radiations_n,  only: nlwb_mnh, nswb_mnh
 use modd_series,        only: lseries
 use modd_series_n,      only: nsnbstept
 use modd_station_n,     only: lstation, tstations_time
+use modd_sto_file,      only: ntrajstlg
 
 TYPE(TFILEDATA),INTENT(INOUT)        :: TPFILE
 CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: HPROGRAM_ORIG !To emulate a file coming from this program
@@ -344,6 +347,11 @@ if ( tpfile%ctype /= 'MNHDIACHRONIC' .and. crad /= 'NONE' ) then
   call IO_Add_dim_nc4( tpfile, NMNHDIM_NLWB, 'LW_bands', nlwb_mnh ) !number of LW bands practically used
 end if
 
+if ( tpfile%ctype == 'MNHDIAG' .and. Trim( yprogram ) == 'DIAG' ) then
+   !Number of times for Lagrangian trajectories
+  if ( ltraj .and. ntrajstlg > 0 ) call IO_Add_dim_nc4( tpfile, NMNHDIM_TRAJ_TIME, 'time_lagrangian', ntrajstlg )
+end if
+
 !Write dimensions used in diachronic files
 if ( tpfile%ctype == 'MNHDIACHRONIC' ) then
   !Dimension of size 2 used for NMNHDIM_COMPLEX
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index 01d333dba15f5db4e8a066d287c0975de597f94f..d13611fe0d687f2ea39732043c93591c8d99a3d3 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -33,6 +33,7 @@
 !  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
+!  P. Wautelet 15/02/2024: IO_Coordvar_write_nc4: add time dimension for Lagrangian trajectories
 !-----------------------------------------------------------------
 #ifdef MNH_IOCDF4
 module mode_io_write_nc4
@@ -1441,10 +1442,11 @@ use modd_budget,     only: cbutype, lbu_icp, lbu_jcp, lbu_kcp, nbuih, nbuil, nbu
                            nbustep, nbutotwrite
 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_dyn_n,      only: 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_LEVEL, NMNHDIM_LEVEL_W, NMNHDIM_TIME, NMNHDIM_TRAJ_TIME,                &
                            NMNHDIM_BUDGET_CART_NI,   NMNHDIM_BUDGET_CART_NJ,   NMNHDIM_BUDGET_CART_NI_U,   &
                            NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NI_V, NMNHDIM_BUDGET_CART_NJ_V,   &
                            NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W,                         &
@@ -1688,6 +1690,14 @@ if ( tzfile%lmaster ) then !time scale is the same on all processes
   end if
 end if
 
+! Write time scale for lagrangian variables
+if ( tzfile%lmaster ) then
+  if ( Trim( yprogram ) == 'DIAG' ) then
+    if ( ltraj .and. ntrajstlg > 0 ) &
+      call Write_time_coord( tzfile%tncdims%tdims(NMNHDIM_TRAJ_TIME), 'time axis for Lagrangian trajectories', tlagr_dates )
+  end if
+end if
+
 if ( tzfile%lmaster ) then
   !Write coordinates used in diachronic files
   if ( tzfile%ctype == 'MNHDIACHRONIC' ) then
diff --git a/src/MNH/compute_r00.f90 b/src/MNH/compute_r00.f90
index 727b230ffd52a9b7afa1b20d5748bc45e34876b4..7c5246cccb502b4a09d854dd216108f52b002238 100644
--- a/src/MNH/compute_r00.f90
+++ b/src/MNH/compute_r00.f90
@@ -1,26 +1,98 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2024 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.
 !-----------------------------------------------------------------
-!     ###############################
-      MODULE MODI_COMPUTE_R00
-!     ###############################
-!
-INTERFACE
-SUBROUTINE COMPUTE_R00(TPFILE)
-!
-USE MODD_IO, ONLY: TFILEDATA
+!###############################
+MODULE MODE_COMPUTE_R00
+! ###############################
+
+  IMPLICIT NONE
+
+  PRIVATE
+
+  PUBLIC :: INI_COMPUTE_R00
+  PUBLIC :: COMPUTE_R00
+
+  REAL, PARAMETER :: NSPVAL = -1.E+11
+
+  INTEGER                 :: NFILES
+  INTEGER                 :: NNBR_START
+  INTEGER, DIMENSION(100) :: NBRFILES
+
+!PW:TODO: DOCTOR
+!PW: si JF=1: reprendre comportement orig DIAG
+
+
+CONTAINS
+
+!#################################
+SUBROUTINE INI_COMPUTE_R00()
+!#################################
+
+  USE MODD_LUNIT_n,    ONLY: CINIFILE
+  USE MODD_PARAMETERS, ONLY: NUNDEF
+  USE MODD_STO_FILE,   ONLY: CFILES, NSTART_SUPP, NTRAJSTLG
+
+  USE MODE_MSG
+
+  IMPLICIT NONE
+
+  INTEGER :: IFILECUR, JFILECUR
+  INTEGER :: JLOOP
+  INTEGER :: JF
+
+  JF=1
+  DO WHILE (LEN_TRIM(CFILES(JF))/=0)
+    JF=JF+1
+  END DO
+  !
+  IF (JF/=1) THEN
+    IF (CINIFILE==CFILES(JF-1)) THEN
+      CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'INI_COMPUTE_R00', 'initial file not treated' )
+    END IF
+  END IF
 !
-TYPE(TFILEDATA),   INTENT(IN) :: TPFILE ! Output file
 !
-END SUBROUTINE COMPUTE_R00
-END INTERFACE
-END MODULE MODI_COMPUTE_R00
+!*       2.0    FIND THE FILE TO BE TREATED AND THE INIT-SV FILES
+!               -------------------------------------------------
 !
-!     ###############################
-      SUBROUTINE COMPUTE_R00(TPFILE)
-!     ###############################
+! Search the number of the file to be treated
+  IFILECUR=0
+  DO JFILECUR=1,100
+    IF (CINIFILE==CFILES(JFILECUR)) THEN
+      IFILECUR=JFILECUR
+      EXIT
+    END IF
+  END DO
+  !
+  IF ( IFILECUR == 0) CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'COMPUTE_R00', 'problem with CINIFILE:' // TRIM(CINIFILE) )
+  !
+  ! Search the number of the files(NFILES), where the Lagrangian tracers
+  !have been reinitialized
+  NFILES=0
+  DO JFILECUR=IFILECUR+1,100
+    IF (LEN_TRIM(CFILES(JFILECUR)) /= 0) THEN
+      NFILES= NFILES +1
+      NBRFILES(NFILES)=JFILECUR       ! contains the number of the files where
+                                      ! the Lag. tracers have been restarted
+    ENDIF
+  END DO
+  !
+  ! compute the number of supplementary cumulative starts
+  NNBR_START=1
+  DO JLOOP=1,NFILES-1
+    IF (NSTART_SUPP(JLOOP)/=NUNDEF .AND. NSTART_SUPP(JLOOP)> IFILECUR ) THEN
+      NNBR_START=NNBR_START+1
+    END IF
+  END DO
+
+  NTRAJSTLG = NNBR_START
+END SUBROUTINE INI_COMPUTE_R00
+
+!#############################
+SUBROUTINE COMPUTE_R00(TPFILE)
+!#############################
 !
 !!**** 
 !!
@@ -58,6 +130,7 @@ END MODULE MODI_COMPUTE_R00
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 07/02/2019: force TYPE to a known value for IO_File_add2list
 !  P. Wautelet 11/04/2019: bugfix: nullify TZTRACFILE when appropriate
+!  P. Wautelet 15/02/2024: add time dimension for Lagrangian trajectories
 !-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -65,7 +138,7 @@ END MODULE MODI_COMPUTE_R00
 !
 USE MODD_CONF
 USE MODD_GRID_n
-use modd_field,            only: tfieldmetadata, TYPEREAL
+use modd_field,            only: NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_LEVEL, NMNHDIM_TRAJ_TIME, tfieldmetadata, TYPEREAL
 USE MODD_FIELD_n
 USE MODD_IO,               ONLY: TFILEDATA
 USE MODD_LUNIT_n
@@ -76,7 +149,7 @@ USE MODD_TYPE_DATE
 USE MODD_VAR_ll
 !
 USE MODE_IO_FIELD_READ,    only: IO_Field_read
-USE MODE_IO_FIELD_WRITE,   only: IO_Field_write
+USE MODE_IO_FIELD_WRITE,   only: IO_Field_create, IO_Field_write
 USE MODE_IO_FILE,          only: IO_File_close, IO_File_open
 USE MODE_IO_MANAGE_STRUCT, ONLY: IO_File_add2list
 USE MODE_ll
@@ -92,10 +165,10 @@ TYPE(TFILEDATA),   INTENT(IN) :: TPFILE ! Output file
 !
 !*       0.2   declarations of local variables
 !
-INTEGER                            :: IFILECUR,JFILECUR,NIU,NJU,NKU
-INTEGER                            :: NFILES,JLOOP
+INTEGER                            :: NIU,NJU,NKU
+INTEGER                            :: JFILECUR
+INTEGER                            :: JLOOP
 REAL                               :: ZXOR,ZYOR,ZDX,ZDY
-REAL                               :: ZSPVAL
 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZX0, ZY0, ZZ0        ! origin of the 
        ! particules colocated with the mesh-grid points read in the file
 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZX00, ZY00, ZZ00, ZZL ! cumulative
@@ -106,15 +179,11 @@ REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZRV0          ! same fields
        ! for Rv as for the coordinates of the origin
 REAL, ALLOCATABLE, DIMENSION(:,:,:):: ZWORK1,ZWORK2,ZWORK3       
 TYPE(DATE_TIME)                    :: TDTCUR_START
-CHARACTER(LEN=NMNHNAMELGTMAX)      :: YMNHNAME
-CHARACTER(LEN=24)                  :: YDATE 
-INTEGER                            :: IHOUR, IMINUTE
-REAL                               :: ZSECOND, ZREMAIN
 LOGICAL                            :: GSTART
-INTEGER                            :: INBR_START
 REAL                               :: ZXMAX,ZYMAX,ZZMAX  ! domain extrema
-INTEGER, DIMENSION(100)            :: NBRFILES 
 TYPE(TFIELDMETADATA)               :: TZFIELD
+TYPE(TFIELDMETADATA)               :: TZFIELD_X0, TZFIELD_Y0, TZFIELD_Z0
+TYPE(TFIELDMETADATA)               :: TZFIELD_TH0, TZFIELD_RV0
 TYPE(TFILEDATA),POINTER            :: TZTRACFILE
 !
 !-------------------------------------------------------------------------------
@@ -123,47 +192,6 @@ TYPE(TFILEDATA),POINTER            :: TZTRACFILE
 !               --------------
 !
 TZTRACFILE => NULL()
-ZSPVAL=-1.E+11
-!
-!-------------------------------------------------------------------------------
-!
-!*       2.0    FIND THE FILE TO BE TREATED AND THE INIT-SV FILES
-!               -------------------------------------------------
-!
-! Search the number of the file to be treated
-IFILECUR=0
-DO JFILECUR=1,100
-  IF (CINIFILE==CFILES(JFILECUR)) THEN
-    IFILECUR=JFILECUR
-    EXIT
-  END IF
-END DO
-!
-IF (IFILECUR==0) THEN
-  PRINT*,'PROBLEM WITH THE FOLLOWING FILE: ',CINIFILE
-  PRINT*,CFILES
-!callabortstop
-  CALL PRINT_MSG(NVERB_FATAL,'GEN','COMPUTE_R00','')
-ENDIF
-!
-! Search the number of the files(NFILES), where the Lagrangian tracers 
-!have been reinitialized 
-NFILES=0
-DO JFILECUR=IFILECUR+1,100
-  IF (LEN_TRIM(CFILES(JFILECUR)) /= 0) THEN
-    NFILES= NFILES +1 
-    NBRFILES(NFILES)=JFILECUR       ! contains the number of the files where
-                                    ! the Lag. tracers have been restarted
-  ENDIF
-END DO
-!
-! compute the number of supplementary cumulative starts
-INBR_START=1
-DO JLOOP=1,NFILES-1
-  IF (NSTART_SUPP(JLOOP)/=NUNDEF .AND. NSTART_SUPP(JLOOP)> IFILECUR ) THEN
-    INBR_START=INBR_START+1
-  END IF
-END DO
 !
 !-------------------------------------------------------------------------------
 !
@@ -187,6 +215,9 @@ ALLOCATE(ZZ00(NIU,NJU,NKU))
 ALLOCATE(ZZL(NIU,NJU,NKU))
 ALLOCATE(ZTH0(NIU,NJU,NKU))
 ALLOCATE(ZRV0(NIU,NJU,NKU))
+
+ALLOCATE( TLAGR_DATES(NTRAJSTLG) )
+
 ! initial values
 ZXOR=0.5 * (XXHAT(2)+XXHAT(3)) 
 ZYOR=0.5 * (XYHAT(2)+XYHAT(3))
@@ -214,19 +245,85 @@ ZZ00(:,:,:)=XSVT(:,:,:,NSV_LGEND)*1.E-3   ! ZZ0 in km
 IF (L2D) THEN
   WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
           ZZ00>ZZMAX)
-    ZX00=ZSPVAL
-    ZZ00=ZSPVAL
+    ZX00=NSPVAL
+    ZZ00=NSPVAL
   END WHERE
 ELSE
   WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
           ZY00<ZYOR .OR. ZY00>ZYMAX .OR. &
 	      ZZ00>ZZMAX)
-    ZX00=ZSPVAL
-    ZY00=ZSPVAL
-    ZZ00=ZSPVAL
+    ZX00=NSPVAL
+    ZY00=NSPVAL
+    ZZ00=NSPVAL
   END WHERE
 END IF
 !
+! Create the metadata of the fields (has to be done only once)
+TZFIELD_X0 = TFIELDMETADATA( &
+  CMNHNAME   = 'X0',         &
+  CSTDNAME   = '',           &
+  CLONGNAME  = 'X0',         &
+  CUNITS     = 'km',         &
+  CDIR       = 'XY',         &
+  CCOMMENT   = 'X_Y_Z_X0',   &
+  NGRID      = 1,            &
+  NTYPE      = TYPEREAL,     &
+  NDIMS      = 4,            &
+  NDIMLIST   = [ NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_LEVEL, NMNHDIM_TRAJ_TIME ] )
+CALL IO_FIELD_CREATE( TPFILE, TZFIELD_X0 )
+
+TZFIELD_Y0 = TFIELDMETADATA( &
+  CMNHNAME   = 'Y0',         &
+  CSTDNAME   = '',           &
+  CLONGNAME  = 'Y0',         &
+  CUNITS     = 'km',         &
+  CDIR       = 'XY',         &
+  CCOMMENT   = 'X_Y_Z_Y0',   &
+  NGRID      = 1,            &
+  NTYPE      = TYPEREAL,     &
+  NDIMS      = 4,            &
+  NDIMLIST   = [ NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_LEVEL, NMNHDIM_TRAJ_TIME ] )
+CALL IO_FIELD_CREATE( TPFILE, TZFIELD_Y0 )
+
+TZFIELD_Z0 = TFIELDMETADATA( &
+  CMNHNAME   = 'Z0',         &
+  CSTDNAME   = '',           &
+  CLONGNAME  = 'Z0',         &
+  CUNITS     = 'km',         &
+  CDIR       = 'XY',         &
+  CCOMMENT   = 'X_Y_Z_Z0',   &
+  NGRID      = 1,            &
+  NTYPE      = TYPEREAL,     &
+  NDIMS      = 4,            &
+  NDIMLIST   = [ NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_LEVEL, NMNHDIM_TRAJ_TIME ] )
+CALL IO_FIELD_CREATE( TPFILE, TZFIELD_Z0 )
+
+TZFIELD_TH0 = TFIELDMETADATA(   &
+  CMNHNAME   = 'TH0',           &
+  CSTDNAME   = '',              &
+  CLONGNAME  = 'TH0',           &
+  CUNITS     = 'K',             &
+  CDIR       = 'XY',            &
+  CCOMMENT   = 'X_Y_Z_'//'TH0', &
+  NGRID      = 1,               &
+  NTYPE      = TYPEREAL,        &
+  NDIMS      = 4,               &
+  NDIMLIST   = [ NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_LEVEL, NMNHDIM_TRAJ_TIME ] )
+CALL IO_FIELD_CREATE( TPFILE, TZFIELD_TH0 )
+
+TZFIELD_RV0 = TFIELDMETADATA(   &
+  CMNHNAME   = 'RV0',           &
+  CSTDNAME   = '',              &
+  CLONGNAME  = 'RV0',           &
+  CUNITS     = 'g kg-1',        &
+  CDIR       = 'XY',            &
+  CCOMMENT   = 'X_Y_Z_'//'RV0', &
+  NGRID      = 1,               &
+  NTYPE      = TYPEREAL,        &
+  NDIMS      = 4,               &
+  NDIMLIST   = [ NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_LEVEL, NMNHDIM_TRAJ_TIME ] )
+CALL IO_FIELD_CREATE( TPFILE, TZFIELD_RV0 )
+
 !-------------------------------------------------------------------------------
 !
 !*       4.0    COMPUTE THE ORIGIN STEP BY STEP
@@ -246,134 +343,49 @@ DO JFILECUR=1,NFILES
   GSTART=.FALSE.
   DO JLOOP=1,NFILES
     IF (NBRFILES(JFILECUR)==NSTART_SUPP(JLOOP) .OR. JFILECUR==NFILES) THEN
-      INBR_START=INBR_START-1
+      NNBR_START=NNBR_START-1
       GSTART=.TRUE.
       EXIT
     END IF
   ENDDO
+
+  IF (GSTART) THEN
 !
 !*       4.2 read the potential temp or the water vapor at the start instant      
 !
-  IF (GSTART) THEN
-    !
     CALL IO_Field_read(TZTRACFILE,'DTCUR',TDTCUR_START)
-    IHOUR   = INT(TDTCUR_START%xtime/3600.)
-    ZREMAIN = MOD(TDTCUR_START%xtime,3600.)
-    IMINUTE = INT(ZREMAIN/60.)
-    ZSECOND = MOD(ZREMAIN,60.)
-    WRITE(YDATE,FMT='(1X,I4.4,I2.2,I2.2,2X,I2.2,"H",I2.2,"M", &
-         & F5.2,"S")') TDTCUR_START%nyear, TDTCUR_START%nmonth, TDTCUR_START%nday, IHOUR,IMINUTE,ZSECOND
+    TLAGR_DATES(NNBR_START+1) = TDTCUR_START
     !
     CALL IO_Field_read(TZTRACFILE,'THT',ZTH0(:,:,:))
     !
     CALL IO_Field_read(TZTRACFILE,'RVT',ZRV0(:,:,:))
     ZRV0(:,:,:)=ZRV0(:,:,:)*1.E+3  ! ZRV0 in g/kg
     !
-  END IF
 !
 !*       4.3  store the X0,Y0,Z0 field for the current start before 
 !             computing the new origin
 !
-  IF (GSTART) THEN
-    PRINT *,'INBR_START',INBR_START,' NBRFILES(JFILECUR)',NBRFILES(JFILECUR)
-    WRITE(YMNHNAME,'(A2,I2.2)')'X0',INBR_START
-    TZFIELD = TFIELDMETADATA(                       &
-      CMNHNAME   = TRIM( YMNHNAME ),                &
-      CSTDNAME   = '',                              &
-      CLONGNAME  = TRIM( YMNHNAME ),                &
-      CUNITS     = 'km',                            &
-      CDIR       = 'XY',                            &
-      CCOMMENT   = 'X_Y_Z_'//TRIM(YMNHNAME)//YDATE, &
-      NGRID      = 1,                               &
-      NTYPE      = TYPEREAL,                        &
-      NDIMS      = 3,                               &
-      LTIMEDEP   = .TRUE.                           )
-    PRINT *,'YCOMMENT = ',TRIM(TZFIELD%CCOMMENT)
-    CALL IO_Field_write(TPFILE,TZFIELD,ZX00(:,:,:))
-    !
-    WRITE(YMNHNAME,'(A2,I2.2)')'Y0',INBR_START
-    TZFIELD = TFIELDMETADATA(                       &
-      CMNHNAME   = TRIM( YMNHNAME ),                &
-      CSTDNAME   = '',                              &
-      CLONGNAME  = TRIM(TZFIELD%CMNHNAME),          &
-      CUNITS     = 'km',                            &
-      CDIR       = 'XY',                            &
-      CCOMMENT   = 'X_Y_Z_'//TRIM(YMNHNAME)//YDATE, &
-      NGRID      = 1,                               &
-      NTYPE      = TYPEREAL,                        &
-      NDIMS      = 3,                               &
-      LTIMEDEP   = .TRUE.                           )
-    PRINT *,'YCOMMENT = ',TRIM(TZFIELD%CCOMMENT)
-    CALL IO_Field_write(TPFILE,TZFIELD,ZY00(:,:,:))
-    !
-    WRITE(YMNHNAME,'(A2,I2.2)')'Z0',INBR_START
-    TZFIELD = TFIELDMETADATA(                       &
-      CMNHNAME   = TRIM( YMNHNAME ),                &
-      CSTDNAME   = '',                              &
-      CLONGNAME  = TRIM( YMNHNAME ),                &
-      CUNITS     = 'km',                            &
-      CDIR       = 'XY',                            &
-      CCOMMENT   = 'X_Y_Z_'//TRIM(YMNHNAME)//YDATE, &
-      NGRID      = 1,                               &
-      NTYPE      = TYPEREAL,                        &
-      NDIMS      = 3,                               &
-      LTIMEDEP   = .TRUE.                           )
-    PRINT *,'YCOMMENT = ',TRIM(TZFIELD%CCOMMENT)
-    CALL IO_Field_write(TPFILE,TZFIELD,ZZ00(:,:,:))
-  END IF
-!
+    PRINT *,'INBR_START',NNBR_START,' NBRFILES(JFILECUR)',NBRFILES(JFILECUR)
+    CALL IO_Field_write( TPFILE, TZFIELD_X0, RESHAPE( ZX00(:,:,:), [ SHAPE(ZX00), 1 ] ), KOFFSET = [0, 0, 0, NNBR_START ] )
+    CALL IO_Field_write( TPFILE, TZFIELD_Y0, RESHAPE( ZY00(:,:,:), [ SHAPE(ZY00), 1 ] ), KOFFSET = [0, 0, 0, NNBR_START ] )
+    CALL IO_Field_write( TPFILE, TZFIELD_Z0, RESHAPE( ZZ00(:,:,:), [ SHAPE(ZZ00), 1 ] ), KOFFSET = [0, 0, 0, NNBR_START ] )
 !
 !*       4.6   compute and store potential temp and water vapor at the origin
 !
-  IF (GSTART) THEN
-    !
-    CALL INTERPXYZ(ZX00,ZY00,ZZ00,     &
-                   ZTH0,ZWORK1         )
-    !
-    CALL INTERPXYZ(ZX00,ZY00,ZZ00,     &
-                   ZRV0,ZWORK2         )
-    !
+    CALL INTERPXYZ( ZX00, ZY00, ZZ00, ZTH0, ZWORK1 )
+    CALL INTERPXYZ( ZX00, ZY00, ZZ00, ZRV0, ZWORK2 )
+
+    CALL IO_Field_write( TPFILE, TZFIELD_TH0, RESHAPE( ZWORK1(:,:,:), [ SHAPE(ZWORK1), 1 ] ), KOFFSET = [0, 0, 0, NNBR_START ] )
+    CALL IO_Field_write( TPFILE, TZFIELD_RV0, RESHAPE( ZWORK2(:,:,:), [ SHAPE(ZWORK2), 1 ] ), KOFFSET = [0, 0, 0, NNBR_START ] )
   END IF
 !
-  IF (GSTART) THEN
-    !
-    WRITE(YMNHNAME,'(A3,I2.2)')'TH0',INBR_START
-    TZFIELD = TFIELDMETADATA(                       &
-      CMNHNAME   = TRIM( YMNHNAME ),                &
-      CSTDNAME   = '',                              &
-      CLONGNAME  = TRIM( YMNHNAME ),                &
-      CUNITS     = 'K',                             &
-      CDIR       = 'XY',                            &
-      CCOMMENT   = 'X_Y_Z_'//TRIM(YMNHNAME)//YDATE, &
-      NGRID      = 1,                               &
-      NTYPE      = TYPEREAL,                        &
-      NDIMS      = 3,                               &
-      LTIMEDEP   = .TRUE.                           )
-    PRINT *,'YCOMMENT = ',TRIM(TZFIELD%CCOMMENT)
-    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK1(:,:,:))
-    !
-    WRITE(YMNHNAME,'(A3,I2.2)')'RV0',INBR_START
-    TZFIELD = TFIELDMETADATA(                       &
-      CMNHNAME   = TRIM( YMNHNAME ),                &
-      CSTDNAME   = '',                              &
-      CLONGNAME  = TRIM( YMNHNAME ),                &
-      CUNITS     = 'g kg-1',                        &
-      CDIR       = 'XY',                            &
-      CCOMMENT   = 'X_Y_Z_'//TRIM(YMNHNAME)//YDATE, &
-      NGRID      = 1,                               &
-      NTYPE      = TYPEREAL,                        &
-      NDIMS      = 3,                               &
-      LTIMEDEP   = .TRUE.                           )
-    PRINT *,'YCOMMENT = ',TRIM(TZFIELD%CCOMMENT)
-    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK2(:,:,:))
-  ENDIF
 !*       4.4   compute the origin of the particules using one more segment
 !
   IF (JFILECUR /= NFILES) THEN
     TZFIELD = TFIELDMETADATA(&
-      CMNHNAME   = 'LGXT',   &
+      CMNHNAME   = 'LGX',    &
       CSTDNAME   = '',       &
-      CLONGNAME  = 'LGXT',   &
+      CLONGNAME  = 'LGX',    &
       CUNITS     = 'm',      &
       CDIR       = 'XY',     &
       CCOMMENT   = '',       &
@@ -384,13 +396,13 @@ DO JFILECUR=1,NFILES
     CALL IO_Field_read(TZTRACFILE,TZFIELD,ZX0)
     ZX0(:,:,:)=ZX0(:,:,:)*1.E-3   ! ZX0 in km
     !
-    TZFIELD%CMNHNAME   = 'LGYT'
-    TZFIELD%CLONGNAME  = 'LGYT'
+    TZFIELD%CMNHNAME   = 'LGY'
+    TZFIELD%CLONGNAME  = 'LGY'
     CALL IO_Field_read(TZTRACFILE,TZFIELD,ZY0)
     ZY0(:,:,:)=ZY0(:,:,:)*1.E-3   ! ZY0 in km
     !
-    TZFIELD%CMNHNAME   = 'LGZT'
-    TZFIELD%CLONGNAME  = 'LGZT'
+    TZFIELD%CMNHNAME   = 'LGZ'
+    TZFIELD%CLONGNAME  = 'LGZ'
     CALL IO_Field_read(TZTRACFILE,TZFIELD,ZZ0)
     ZZ0(:,:,:)=ZZ0(:,:,:)*1.E-3   ! ZZ0 in km
     !
@@ -410,16 +422,16 @@ DO JFILECUR=1,NFILES
     IF (L2D) THEN
       WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
               ZZ00>ZZMAX)
-        ZX00=ZSPVAL
-        ZZ00=ZSPVAL
+        ZX00=NSPVAL
+        ZZ00=NSPVAL
       END WHERE
     ELSE
       WHERE ( ZX00<ZXOR .OR. ZX00>ZXMAX .OR. &
               ZY00<ZYOR .OR. ZY00>ZYMAX .OR. &
               ZZ00>ZZMAX)
-        ZX00=ZSPVAL
-        ZY00=ZSPVAL
-        ZZ00=ZSPVAL
+        ZX00=NSPVAL
+        ZY00=NSPVAL
+        ZZ00=NSPVAL
       END WHERE
     END IF
     !
@@ -479,17 +491,17 @@ DO JK=1,NKU
       !
       ! remove external points
       IF (L2D) THEN
-        GEXT=(ZX==ZSPVAL).OR.(ZZ==ZSPVAL)
+        GEXT=(ZX==NSPVAL).OR.(ZZ==NSPVAL)
       ELSE
-        GEXT=(ZX==ZSPVAL).OR.(ZY==ZSPVAL).OR.(ZZ==ZSPVAL)
+        GEXT=(ZX==NSPVAL).OR.(ZY==NSPVAL).OR.(ZZ==NSPVAL)
       END IF
       IF (GEXT) THEN
-        POUT1(JI,JJ,JK) = ZSPVAL
+        POUT1(JI,JJ,JK) = NSPVAL
         IF (PRESENT(PIN2)) THEN
-          POUT2(JI,JJ,JK) = ZSPVAL
+          POUT2(JI,JJ,JK) = NSPVAL
         END IF
         IF (PRESENT(PIN3)) THEN
-          POUT3(JI,JJ,JK) = ZSPVAL
+          POUT3(JI,JJ,JK) = NSPVAL
         ENDIF
         !
         CYCLE
@@ -574,3 +586,5 @@ END SUBROUTINE INTERPXYZ
 !-------------------------------------------------------------------------------
 !
 END SUBROUTINE COMPUTE_R00 
+
+END MODULE MODE_COMPUTE_R00
diff --git a/src/MNH/diag.f90 b/src/MNH/diag.f90
index fd25ea531d56d9849962cdda587f46d2bebd84c6..143696420466c1ab512d17a54920d53649bf2beb 100644
--- a/src/MNH/diag.f90
+++ b/src/MNH/diag.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1999-2023 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1999-2024 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.
@@ -138,6 +138,7 @@ USE MODD_TIME_n
 USE MODD_TURB_n
 USE MODD_VAR_ll
 !
+USE MODE_COMPUTE_R00
 USE MODE_DATETIME
 USE MODE_FINALIZE_MNH,     only: FINALIZE_MNH
 USE MODE_IO_FILE,          only: IO_File_close, IO_File_open
@@ -154,7 +155,6 @@ USE MODE_POS
 USE MODE_TIME
 !
 USE MODI_CH_MONITOR_n
-USE MODI_COMPUTE_R00
 USE MODI_DIAG_SURF_ATM_N
 USE MODI_INIT_MNH
 USE MODI_PHYS_PARAM_n
@@ -187,7 +187,6 @@ REAL(kind=MNHTIME), DIMENSION(2) :: ZSTART, ZINIT, ZWRIT, ZPHYS, ZSURF, ZWRITS,
 INTEGER(KIND=LFIINT) :: INPRAR ! number of articles predicted  in the LFIFM file
 INTEGER :: ILUNAM      ! Logical unit numbers for the namelist file
                        ! and for output_listing file
-INTEGER        :: JF =0   !  loop index
 LOGICAL :: GFOUND         ! Return code when searching namelist
 LOGICAL, DIMENSION(:,:),ALLOCATABLE     :: GMASKkids ! kids domains mask
 LOGICAL:: GCLOUD_ONLY          ! conditionnal radiation computations for
@@ -406,21 +405,7 @@ CALL IO_File_close(TZNMLFILE)
 CINIFILE = YINIFILE(1)
 CINIFILEPGD = YINIFILEPGD(1)
 !!
-IF (LTRAJ) THEN
-  JF=1
-  DO WHILE (LEN_TRIM(CFILES(JF))/=0)
-    JF=JF+1
-  END DO
-  !
-  IF (JF/=1) THEN
-    IF (CINIFILE==CFILES(JF-1)) THEN
-!callabortstop
-      CALL PRINT_MSG(NVERB_FATAL,'GEN','DIAG','initial file not treated')
-    END IF
-  END IF
-!
-END IF
-!
+IF ( LTRAJ ) CALL INI_COMPUTE_R00()
 !
 !-------------------------------------------------------------------------------
 !
@@ -655,7 +640,7 @@ ZTIME1=ZTIME2
 !
 !*       8.0    Initial positions computation (back into simulation segments)
 !
-IF (LTRAJ .AND. JF/=1) CALL COMPUTE_R00(TOUTDATAFILE)
+IF ( LTRAJ ) CALL COMPUTE_R00( TOUTDATAFILE )
 !
 CALL SECOND_MNH2(ZTIME2)
 ZTRAJ =ZTIME2-ZTIME1
diff --git a/src/MNH/modd_sto_file.f90 b/src/MNH/modd_sto_file.f90
index 0278b8e0ee9a272935331d629b05ae77c53ec1a6..34c7587017518c3ab388ebce58a72b1c7c02e5ab 100644
--- a/src/MNH/modd_sto_file.f90
+++ b/src/MNH/modd_sto_file.f90
@@ -1,20 +1,25 @@
-!MNH_LIC Copyright 1994-2023 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2024 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.
 !-----------------------------------------------------------------
 ! Modifications:
-!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
+!  P. Wautelet 15/02/2024: add time dimension for Lagrangian trajectories
 !-----------------------------------------------------------------
-!     ######spl
+!###################
 MODULE MODD_STO_FILE
-!################
+!###################
 
 USE MODD_PARAMETERS, ONLY: NFILENAMELGTMAX
+USE MODD_TYPE_DATE,  ONLY: DATE_TIME
 
 SAVE
 
-CHARACTER(LEN=NFILENAMELGTMAX) :: CFILES(100)       ! names of the files to be treated
+CHARACTER(LEN=NFILENAMELGTMAX) :: CFILES(100)      ! names of the files to be treated
 INTEGER                        :: NSTART_SUPP(100) ! supplementary starts for the lagrangian trajectories
-!
+
+INTEGER                                    :: NTRAJSTLG = 0 ! Number of time starts for Lagrangian trajectories
+TYPE(DATE_TIME), DIMENSION(:), ALLOCATABLE :: TLAGR_DATES   ! Times for Lagrangian trajectories
+
 END MODULE MODD_STO_FILE