From 33d772457b3a15f00c931c7b9ec20f373a584f1e Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Fri, 29 Apr 2022 15:15:00 +0200
Subject: [PATCH] Philippe 29/04/2022: restructure profilers for better
 performance, reduce memory usage and correct some problems/bugs + share code
 with stations

---
 src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90 |    4 +-
 src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 |    4 +-
 src/MNH/goto_model_wrapper.f90              |    2 -
 src/MNH/ini_modeln.f90                      |    4 +-
 src/MNH/ini_posprofilern.f90                |  329 ++---
 src/MNH/ini_profilern.f90                   |  176 ---
 src/MNH/ini_surfstationn.f90                |   26 +-
 src/MNH/modd_profilern.f90                  |   45 +-
 src/MNH/modd_stationn.f90                   |    5 +-
 src/MNH/modd_sub_profilern.f90              |  109 --
 src/MNH/modd_type_profiler.f90              |  115 --
 src/MNH/modd_type_statprof.f90              |  159 ++-
 src/MNH/modeln.f90                          |   13 +-
 src/MNH/modn_stationn.f90                   |    1 -
 src/MNH/profilern.f90                       | 1242 +++++++------------
 src/MNH/stationn.f90                        |   84 +-
 src/MNH/statprof_reader.f90                 |   99 +-
 src/MNH/statprof_tools.f90                  |  542 ++++++--
 src/MNH/write_profilern.f90                 |  525 +++++---
 src/MNH/write_stationn.f90                  |   56 +-
 20 files changed, 1659 insertions(+), 1881 deletions(-)
 delete mode 100644 src/MNH/ini_profilern.f90
 delete mode 100644 src/MNH/modd_sub_profilern.f90
 delete mode 100644 src/MNH/modd_type_profiler.f90

diff --git a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
index e462d6f19..09d1c225a 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_tools_nc4.f90
@@ -278,7 +278,7 @@ use modd_les_n,         only: nles_times, nspectra_ni, nspectra_nj
 use modd_nsv,           only: nsv
 USE MODD_PARAMETERS_ll, ONLY: JPHEXT, JPVEXT
 use modd_param_n,       only: crad
-use modd_profiler_n,    only: lprofiler, tprofiler
+use modd_profiler_n,    only: lprofiler, tprofilers_time
 use modd_radiations_n,  only: nlwb_mnh, nswb_mnh
 use modd_series,        only: lseries
 use modd_series_n,      only: nsnbstept
@@ -423,7 +423,7 @@ if ( tpfile%ctype == 'MNHDIACHRONIC' ) then
 
   !Dimension for the number of profiler times
   if ( lprofiler ) then
-    iprof = Nint ( ( xseglen - dyn_model(1)%xtstep ) / tprofiler%step ) + 1
+    iprof = Nint ( ( xseglen - dyn_model(1)%xtstep ) / tprofilers_time%xtstep ) + 1
     call IO_Add_dim_nc4( tpfile, NMNHDIM_PROFILER_TIME, 'time_profiler', iprof )
   end if
 
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index 223246cca..c554d7da1 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -1456,7 +1456,7 @@ use modd_les_n,      only: nles_dtcount, nles_mean_end, nles_mean_start, nles_me
                            nles_times, nspectra_ni, nspectra_nj, tles_dates, xles_times
 use modd_netcdf,     only: tdimnc
 use modd_parameters, only: jphext, JPVEXT
-use modd_profiler_n, only: lprofiler, tprofiler
+use modd_profiler_n, only: lprofiler, tprofilers_time
 use modd_series,     only: lseries
 use modd_series_n,   only: nsnbstept, tpsdates
 use modd_station_n,  only: lstation, tstations_time
@@ -1837,7 +1837,7 @@ 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', tprofiler%tpdates )
+      call Write_time_coord( tpfile%tncdims%tdims(NMNHDIM_PROFILER_TIME), 'time axis for profilers', tprofilers_time%tpdates )
 
     !Coordinates for the number of station times
     if ( lstation ) &
diff --git a/src/MNH/goto_model_wrapper.f90 b/src/MNH/goto_model_wrapper.f90
index c0e7e24dd..33dbefdf2 100644
--- a/src/MNH/goto_model_wrapper.f90
+++ b/src/MNH/goto_model_wrapper.f90
@@ -109,7 +109,6 @@ USE MODD_SUB_ELEC_n
 USE MODD_SUB_MODEL_n
 USE MODD_SUB_PASPOL_n
 USE MODD_SUB_PHYS_PARAM_n
-USE MODD_SUB_PROFILER_n
 USE MODD_TIMEZ
 USE MODD_TURB_n
 !
@@ -205,7 +204,6 @@ CALL SUB_CH_FIELD_VALUE_GOTO_MODEL(KFROM, KTO)
 CALL SUB_CH_MONITOR_GOTO_MODEL(KFROM, KTO)
 CALL SUB_MODEL_GOTO_MODEL(KFROM, KTO)
 CALL SUB_PHYS_PARAM_GOTO_MODEL(KFROM, KTO)
-CALL SUB_PROFILER_GOTO_MODEL(KFROM, KTO)
 CALL SUB_PASPOL_GOTO_MODEL(KFROM, KTO)
 CALL SUB_ELEC_GOTO_MODEL(KFROM, KTO)
 !CALL TIME_GOTO_MODEL(KFROM, KTO)
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 6ee23d80a..27d610835 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -2572,9 +2572,7 @@ CALL INI_SURFSTATION_n( )
 !*      25.     PROFILER initializations
 !              ------------------------
 !
-CALL INI_POSPROFILER_n(XTSTEP, XSEGLEN, NRR, NSV,  &
-                       CTURB=="TKEL",              &
-                       XLATORI, XLONORI            )
+CALL INI_POSPROFILER_n( )
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/ini_posprofilern.f90 b/src/MNH/ini_posprofilern.f90
index 00279c8a8..ad6340211 100644
--- a/src/MNH/ini_posprofilern.f90
+++ b/src/MNH/ini_posprofilern.f90
@@ -9,19 +9,7 @@ MODULE MODI_INI_POSPROFILER_n
 !
 INTERFACE
 !
-      SUBROUTINE INI_POSPROFILER_n(PTSTEP, PSEGLEN,          &
-                                   KRR, KSV, OUSETKE,        &
-                                   PLATOR, PLONOR            )
-!
-REAL,               INTENT(IN) :: PTSTEP  ! time step
-REAL,               INTENT(IN) :: PSEGLEN ! segment length
-INTEGER,            INTENT(IN) :: KRR     ! number of moist variables
-INTEGER,            INTENT(IN) :: KSV     ! number of scalar variables
-LOGICAL,            INTENT(IN) :: OUSETKE ! flag to use tke
-REAL,               INTENT(IN) :: PLATOR  ! latitude of origine point
-REAL,               INTENT(IN) :: PLONOR  ! longitude of origine point
-!
-!-------------------------------------------------------------------------------
+      SUBROUTINE INI_POSPROFILER_n( )
 !
 END SUBROUTINE INI_POSPROFILER_n
 !
@@ -29,11 +17,9 @@ END INTERFACE
 !
 END MODULE MODI_INI_POSPROFILER_n
 !
-!     ########################################################
-      SUBROUTINE INI_POSPROFILER_n(PTSTEP, PSEGLEN,          &
-                                   KRR, KSV, OUSETKE,        &
-                                   PLATOR, PLONOR            )
-!     ########################################################
+!     ###############################
+      SUBROUTINE INI_POSPROFILER_n( )
+!     ###############################
 !
 !
 !!****  *INI_POSPROFILER_n* - 
@@ -66,219 +52,144 @@ END MODULE MODI_INI_POSPROFILER_n
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
 !  M. Taufour  05/07/2021: modify RARE for hydrometeors containing ice and add bright band calculation for RARE
+!  P. Wautelet    04/2022: restructure profilers for better performance, reduce memory usage and correct some problems/bugs
 !! --------------------------------------------------------------------------
-!       
+!
 !*      0. DECLARATIONS
 !          ------------
 !
-USE MODD_CONF
-USE MODD_DYN_n
-USE MODD_GRID
-USE MODD_GRID_n
-USE MODD_LUNIT_n,      ONLY: TLUOUT
-USE MODD_PARAMETERS
-USE MODD_PROFILER_n
-USE MODD_RADIATIONS_n, ONLY: NAER
-USE MODD_TYPE_PROFILER
-!
-USE MODE_GRIDPROJ
-USE MODE_ll
+USE MODD_ALLPROFILER_n
+USE MODD_CONF,           ONLY: LCARTESIAN
+USE MODD_DYN,            ONLY: XSEGLEN
+USE MODD_DYN_n,          ONLY: DYN_MODEL, XTSTEP
+USE MODD_GRID_n,         ONLY: XXHAT, XYHAT
+USE MODD_PARAMETERS,     ONLY: JPHEXT, JPVEXT
+USE MODD_PROFILER_n,     ONLY: LPROFILER, NUMBPROFILER_LOC, TPROFILERS, TPROFILERS_TIME
+USE MODD_TYPE_STATPROF,  ONLY: TPROFILERDATA
+
+USE MODE_ALLOCBUFFER_ll,  ONLY: ALLOCBUFFER_ll
+USE MODE_GATHER_ll,       ONLY: GATHERALL_FIELD_ll
 USE MODE_MSG
-!
-USE MODI_INI_PROFILER_N
-!
+USE MODE_STATPROF_READER, ONLY: STATPROF_CSV_READ
+USE MODE_STATPROF_TOOLS,  ONLY: PROFILER_ADD, PROFILER_ALLOCATE, STATPROF_INI_INTERP, STATPROF_POSITION
+
 IMPLICIT NONE
 !
 !
 !*      0.1  declarations of arguments
 !
-!
-REAL,               INTENT(IN) :: PTSTEP  ! time step
-REAL,               INTENT(IN) :: PSEGLEN ! segment length
-INTEGER,            INTENT(IN) :: KRR     ! number of moist variables
-INTEGER,            INTENT(IN) :: KSV     ! number of scalar variables
-LOGICAL,            INTENT(IN) :: OUSETKE ! flag to use tke
-REAL,               INTENT(IN) :: PLATOR  ! latitude of origine point
-REAL,               INTENT(IN) :: PLONOR  ! longitude of origine point
+! NONE
 !
 !-------------------------------------------------------------------------------
 !
 !       0.2  declaration of local variables
 !
-INTEGER :: ISTORE   ! number of storage instants
-INTEGER :: ILUOUT   ! logical unit
-INTEGER :: IKU      !
+INTEGER :: IERR
+INTEGER :: IIU
+INTEGER :: IJU
+INTEGER :: INUMBPROF                        ! Total number of profilers (inside physical domain of model)
+INTEGER :: ISTORE                           ! number of storage instants
+INTEGER :: JI
+LOGICAL :: GALLOCX, GALLOCY
+LOGICAL :: GINSIDE                          ! True if profiler is inside physical domain of model
+LOGICAL :: GPRESENT                         ! True if profiler is present on the current process
+REAL    :: ZXHATM_PHYS_MIN, ZYHATM_PHYS_MIN ! Minimum X coordinate of mass points in the physical domain
+REAL    :: ZXHATM_PHYS_MAX, ZYHATM_PHYS_MAX ! Minimum X coordinate of mass points in the physical domain
+REAL, DIMENSION(SIZE(XXHAT)) :: ZXHATM      ! mass point coordinates
+REAL, DIMENSION(SIZE(XYHAT)) :: ZYHATM      ! mass point coordinates
+REAL, DIMENSION(:), POINTER  :: ZXHAT_GLOB
+REAL, DIMENSION(:), POINTER  :: ZYHAT_GLOB
+TYPE(TPROFILERDATA)          :: TZPROFILER
 !
 !----------------------------------------------------------------------------
-ILUOUT = TLUOUT%NLU
-!----------------------------------------------------------------------------
-!
-!*      1.   Default values
-!            --------------
-IKU =  SIZE(XZZ,3)     ! nombre de niveaux verticaux
-!
-CALL DEFAULT_PROFILER_n(TPROFILER)
-!
-!
-!*      3.   Stations initialization
-!            -----------------------
-!
-CALL INI_PROFILER_n
-LPROFILER = (NUMBPROFILER>0)
-!
-!----------------------------------------------------------------------------
-!
-!*      4.   Allocations of storage arrays
-!            -----------------------------
-!
-IF(NUMBPROFILER>0) THEN
-  CALL ALLOCATE_PROFILER_n(TPROFILER)
-  CALL INI_INTERP_PROFILER_n(TPROFILER)
+
+TPROFILERS_TIME%XTSTEP = XSTEP_PROF
+
+if ( tprofilers_time%xtstep < xtstep ) then
+  call Print_msg( NVERB_WARNING, 'GEN', 'INI_POSPROFILER_n', 'Timestep for profilers was smaller than model timestep' )
+  tprofilers_time%xtstep = xtstep
+end if
+
+ISTORE = NINT ( ( XSEGLEN - DYN_MODEL(1)%XTSTEP ) / TPROFILERS_TIME%XTSTEP ) + 1
+
+allocate( tprofilers_time%tpdates(istore) )
+!
+! Prepare positioning data
+!
+IF ( CFILE_PROF /= "NO_INPUT_CSV" .OR. NNUMB_PROF > 0 ) THEN
+  IIU = SIZE( XXHAT )
+  IJU = SIZE( XYHAT )
+
+  ! Get global XHAT and YHAT (needed by PROFILER_POSITION)
+  CALL ALLOCBUFFER_ll( ZXHAT_GLOB, XXHAT, 'XX', GALLOCX )
+  CALL ALLOCBUFFER_ll( ZYHAT_GLOB, XYHAT, 'YY', GALLOCY )
+  CALL GATHERALL_FIELD_ll( 'XX', XXHAT, ZXHAT_GLOB, IERR )
+  CALL GATHERALL_FIELD_ll( 'YY', XYHAT, ZYHAT_GLOB, IERR )
+
+  ! Interpolations of model variables to mass points
+  ZXHATM(1:IIU-1) = 0.5 * XXHAT(1:IIU-1) + 0.5 * XXHAT(2:IIU  )
+  ZXHATM(  IIU  ) = 1.5 * XXHAT(  IIU  ) - 0.5 * XXHAT(  IIU-1)
+
+  ZYHATM(1:IJU-1) = 0.5 * XYHAT(1:IJU-1) + 0.5 * XYHAT(2:IJU  )
+  ZYHATM(  IJU  ) = 1.5 * XYHAT(  IJU  ) - 0.5 * XYHAT(  IJU-1)
+
+  ZXHATM_PHYS_MIN = 0.5 * ( ZXHAT_GLOB(1+JPHEXT) + ZXHAT_GLOB(2+JPHEXT) )
+  ZXHATM_PHYS_MAX = 0.5 * ( ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT) + ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT+1) )
+  ZYHATM_PHYS_MIN = 0.5 * ( ZYHAT_GLOB(1+JPHEXT) + ZYHAT_GLOB(2+JPHEXT) )
+  ZYHATM_PHYS_MAX = 0.5 * ( ZYHAT_GLOB(UBOUND(ZYHAT_GLOB,1)-JPHEXT) + ZYHAT_GLOB(UBOUND(ZYHAT_GLOB,1)-JPHEXT+1) )
 END IF
-!----------------------------------------------------------------------------
 !
-CONTAINS
-!
-!----------------------------------------------------------------------------
-SUBROUTINE DEFAULT_PROFILER_n(TPROFILER)
-!
-TYPE(PROFILER), INTENT(INOUT) :: TPROFILER
-!
-NUMBPROFILER     = 0
-TPROFILER%T_CUR  = XUNDEF
-TPROFILER%N_CUR  = 0
-TPROFILER%STEP   = XTSTEP 
-!
-END SUBROUTINE DEFAULT_PROFILER_n
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-SUBROUTINE ALLOCATE_PROFILER_n(TPROFILER)
-!
-TYPE(PROFILER), INTENT(INOUT) :: TPROFILER
-!
-ISTORE = NINT( ( PSEGLEN - DYN_MODEL(1)%XTSTEP ) / TPROFILER%STEP ) + 1
-!
-allocate( tprofiler%tpdates( istore ) )
-ALLOCATE(TPROFILER%ERROR (NUMBPROFILER))
-ALLOCATE(TPROFILER%X     (NUMBPROFILER))
-ALLOCATE(TPROFILER%Y     (NUMBPROFILER))
-ALLOCATE(TPROFILER%ZON   (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%MER   (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%FF    (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%DD    (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%W     (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%P     (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%ZZ    (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%TH    (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%THV   (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%RHOD  (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%VISI  (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%VISIKUN(ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%CRARE  (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%CRARE_ATT(ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%LWCZ  (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%IWCZ  (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%CIZ  (ISTORE,IKU,NUMBPROFILER))
-ALLOCATE(TPROFILER%R     (ISTORE,IKU,NUMBPROFILER,KRR))
-ALLOCATE(TPROFILER%SV    (ISTORE,IKU,NUMBPROFILER,KSV))
-ALLOCATE(TPROFILER%AER   (ISTORE,IKU,NUMBPROFILER,NAER))
-IF (OUSETKE) THEN
-  ALLOCATE(TPROFILER%TKE (ISTORE,IKU,NUMBPROFILER))
+! Profilers initialization
+!
+NUMBPROFILER_LOC = 0
+
+IF (CFILE_PROF=="NO_INPUT_CSV") THEN
+  ! Treat namelist
+  INUMBPROF = 0
+  IF ( NNUMB_PROF > 0 ) THEN
+    DO JI = 1, NNUMB_PROF
+      IF ( LCARTESIAN ) THEN
+        TZPROFILER%XX = XX_PROF(JI)
+        TZPROFILER%XY = XY_PROF(JI)
+      ELSE
+        TZPROFILER%XLAT = XLAT_PROF(JI)
+        TZPROFILER%XLON = XLON_PROF(JI)
+        CALL STATPROF_INI_INTERP( TZPROFILER )
+      END IF
+      TZPROFILER%XZ    = XZ_PROF(JI)
+      TZPROFILER%CNAME = CNAME_PROF(JI)
+
+      CALL STATPROF_POSITION( TZPROFILER, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,                 &
+                              ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX, ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
+                              GINSIDE, GPRESENT                                                   )
+
+      IF ( GINSIDE ) THEN
+        INUMBPROF = INUMBPROF + 1
+        TZPROFILER%NID = INUMBPROF
+      END IF
+
+      IF ( GPRESENT ) CALL PROFILER_ADD( TZPROFILER )
+    END DO
+  END IF
 ELSE
-  ALLOCATE(TPROFILER%TKE (0,IKU,0))
+  !Treat CSV datafile
+  CALL STATPROF_CSV_READ( TZPROFILER, CFILE_PROF, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,    &
+                          ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX,ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
+                          INUMBPROF                                                          )
 END IF
-ALLOCATE(TPROFILER%T2M     (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%Q2M     (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%HU2M    (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%ZON10M  (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%MER10M  (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%RN      (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%H       (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%LE      (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%LEI     (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%GFLUX   (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%LWD     (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%LWU     (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%SWD     (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%SWU     (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%IWV   (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%ZTD   (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%ZWD   (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%ZHD   (ISTORE,NUMBPROFILER))
-ALLOCATE(TPROFILER%TKE_DISS(ISTORE,IKU,NUMBPROFILER))
-!
-!
-TPROFILER%ERROR= .FALSE.
-TPROFILER%ZON  = XUNDEF
-TPROFILER%MER  = XUNDEF
-TPROFILER%FF   = XUNDEF
-TPROFILER%DD   = XUNDEF
-TPROFILER%W    = XUNDEF
-TPROFILER%P    = XUNDEF
-TPROFILER%ZZ   = XUNDEF
-TPROFILER%TH   = XUNDEF
-TPROFILER%THV  = XUNDEF
-TPROFILER%RHOD = XUNDEF
-TPROFILER%VISI = XUNDEF
-TPROFILER%VISIKUN = XUNDEF
-TPROFILER%CRARE = XUNDEF
-TPROFILER%CRARE_ATT = XUNDEF
-TPROFILER%LWCZ = XUNDEF
-TPROFILER%IWCZ = XUNDEF
-TPROFILER%CIZ = XUNDEF
-TPROFILER%IWV  = XUNDEF
-TPROFILER%ZTD  = XUNDEF
-TPROFILER%ZWD  = XUNDEF
-TPROFILER%ZHD  = XUNDEF
-TPROFILER%R    = XUNDEF
-TPROFILER%SV   = XUNDEF
-TPROFILER%AER  = XUNDEF
-TPROFILER%TKE  = XUNDEF
-TPROFILER%T2M      = XUNDEF
-TPROFILER%Q2M      = XUNDEF
-TPROFILER%HU2M     = XUNDEF
-TPROFILER%ZON10M   = XUNDEF
-TPROFILER%MER10M   = XUNDEF
-TPROFILER%RN       = XUNDEF
-TPROFILER%H        = XUNDEF
-TPROFILER%LE       = XUNDEF
-TPROFILER%GFLUX    = XUNDEF
-TPROFILER%LEI      = XUNDEF
-TPROFILER%LWD      = XUNDEF
-TPROFILER%LWU      = XUNDEF
-TPROFILER%SWD      = XUNDEF
-TPROFILER%SWU      = XUNDEF
-TPROFILER%TKE_DISS = XUNDEF
-!
-END SUBROUTINE ALLOCATE_PROFILER_n
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-SUBROUTINE INI_INTERP_PROFILER_n(TPROFILER)
-!
-TYPE(PROFILER), INTENT(INOUT) :: TPROFILER
-INTEGER :: III
-INTEGER :: IIU, IJU
-!
-DO III=1,NUMBPROFILER
-  CALL GET_DIM_EXT_ll ('B',IIU,IJU)
-  CALL SM_XYHAT(PLATOR,PLONOR,                          &
-                TPROFILER%LAT(III), TPROFILER%LON(III), &
-                TPROFILER%X(III),   TPROFILER%Y(III)    )
-ENDDO
-!
-IF ( ANY(TPROFILER%LAT(:)==XUNDEF) .OR. ANY(TPROFILER%LON(:)==XUNDEF) ) THEN
-  WRITE(ILUOUT,*) 'Error in station position '
-  WRITE(ILUOUT,*) 'either LATitude or LONgitude segment'
-  WRITE(ILUOUT,*) 'definiton is not complete.'
-!callabortstop
-  CALL PRINT_MSG(NVERB_FATAL,'GEN','INI_INTERP_PROFILER_n','')
+
+LPROFILER = ( INUMBPROF > 0 )
+
+DO JI = 1, NUMBPROFILER_LOC
+  CALL PROFILER_ALLOCATE( TPROFILERS(JI), ISTORE )
+END DO
+!
+! Clean positioning data
+!
+IF ( CFILE_PROF /= "NO_INPUT_CSV" .OR. NNUMB_PROF > 0 ) THEN
+  IF ( GALLOCX ) DEALLOCATE( ZXHAT_GLOB )
+  IF ( GALLOCY ) DEALLOCATE( ZYHAT_GLOB )
 END IF
-!
-TPROFILER%STEP  = MAX ( PTSTEP, TPROFILER%STEP )
-!
-!
-END SUBROUTINE INI_INTERP_PROFILER_n
-!----------------------------------------------------------------------------
 !----------------------------------------------------------------------------
 !
 END SUBROUTINE INI_POSPROFILER_n
diff --git a/src/MNH/ini_profilern.f90 b/src/MNH/ini_profilern.f90
deleted file mode 100644
index 598980099..000000000
--- a/src/MNH/ini_profilern.f90
+++ /dev/null
@@ -1,176 +0,0 @@
-!MNH_LIC Copyright 2002-2020 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.
-!-----------------------------------------------------------------
-!     #######################
-      SUBROUTINE INI_PROFILER_n
-!     #######################
-!
-!
-!!****  *INI_PROFILER_n* - user initializes the station location
-!!
-!!    PURPOSE
-!!    -------
-!
-!
-!!**  METHOD
-!!    ------
-!!    
-!!   Must be defined (for each aircraft):
-!!   ---------------
-!!
-!!  No default exist for these variables.
-!!  ************************************
-!!
-!!  1) Number of stations
-!!  2) the model in which these stations are
-!!     if NOT initialized, the stations are NOT used.
-!!
-!!  3) the (LAT, LON, ALT) latitude,longitude and altitude of the station location.
-!!  4) the station name
-!!
-!!
-!!
-!!   Can be defined  (for each aircraft):
-!!   --------------
-!!
-!!
-!!  9) the time step for data storage.
-!!    default is 60s
-!!
-!! 10) the name or title describing the balloon (8 characters)
-!!     default is the balloon type (6 characters) + the balloon numbers (2 characters)
-!!
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!      Pierre Tulet             * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!     Original 15/01/2002
-!!              July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
-!!                                      aircraft, ballon and profiler
-!!
-!! --------------------------------------------------------------------------
-!       
-!*      0. DECLARATIONS
-!          ------------
-!
-use modd_conf,       only: lcartesian
-USE MODD_PARAMETERS
-USE MODD_PROFILER_n
-
-use mode_msg
-
-IMPLICIT NONE
-!
-!
-!*      0.1  declarations of arguments
-!
-!
-!-------------------------------------------------------------------------------
-!
-!       0.2  declaration of local variables
-!
-!
-!----------------------------------------------------------------------------
-!
-!*      1.   Nameliste 
-!            ---------
-NUMBPROFILER             = 0
-
-if ( numbprofiler > 0 .and. lcartesian ) &
-  call Print_msg( NVERB_FATAL, 'GEN', 'INI_PROFILER_n', 'profilers are not available if LCARTESIAN=T' )
-
-IF (NUMBPROFILER > 0) THEN
-ALLOCATE(TPROFILER%LAT (NUMBPROFILER))
-ALLOCATE(TPROFILER%LON (NUMBPROFILER))
-ALLOCATE(TPROFILER%ALT (NUMBPROFILER))
-ALLOCATE(TPROFILER%NAME(NUMBPROFILER))
-ALLOCATE(TPROFILER%TYPE(NUMBPROFILER))
-!
-TPROFILER%LON = XUNDEF
-TPROFILER%LAT = XUNDEF
-TPROFILER%ALT = XUNDEF
-TPROFILER%NAME = "    "
-TPROFILER%TYPE = "         "
-!
-TPROFILER%STEP = 900.
-!
-!* location (latitude, longitude, altitude)
-!
-!
-TPROFILER%LAT                = (/ 43.3000, 43.3300, 43.6200, 43.6550, 43.3400, &
-                                  43.3000, 43.3500, 43.8128, 44.1711, 44.1689, &
-                                  44.0833, 43.6200, 43.7164, 43.5333, 43.4833, &
-                                  43.4800, 43.4856,                            &
-                                  43.5423, 43.3000, 43.5000, 43.2568, 43.5000, &
-                                  43.3333, 43.9070, 43.5430, 43.5300, 43.5000, &
-                                  43.4300,                                     &
-                                  43.2780, 43.3277, 43.3396, 43.3230, 43.3559, &
-                                  43.3060, 43.3724, 43.2568, 43.3834, 43.3417, &
-                                  43.3165, 43.3873, 43.2833, 43.2900, 43.2842, &
-                                  43.3597, 43.2788, 43.4999 /)
-!
-TPROFILER%LON                = (/ 5.3790, 4.8200, 5.2000, 6.0986, 5.4100, &
-                                  5.3833, 5.4000, 5.1506, 5.2875, 5.0692, &
-                                  5.0500, 5.4000, 5.7653, 5.0667, 5.2833, &
-                                  5.3200, 5.3405,                         &
-                                  4.9010, 5.4000, 4.9300, 5.4054, 5.3700, &
-                                  5.4167, 4.8988, 5.0740, 5.0700, 5.3700, &
-                                  5.2300,                                 &
-                                  5.5762, 5.4782, 5.5873, 5.3668, 5.4550, &
-                                  5.3951, 5.4842, 5.4054, 5.4253, 5.4113, &
-                                  5.4369, 5.3893, 5.5114, 5.3788, 5.4611, &
-                                  5.3994, 5.3538, 5.3674 /)
-TPROFILER%ALT                = (/ 2.,2.,2.,2.,2.,&
-                                  2.,2.,2.,2.,2.,&
-                                  2.,2.,2.,2.,2.,&
-                                  2.,2.,         &
-                                  2.,2.,2.,2.,2.,&
-                                  2.,2.,2.,2.,2.,&
-                                  2.,            &
-                                  2.,2.,2.,2.,2.,&
-                                  2.,2.,2.,2.,2.,&
-                                  2.,2.,2.,2.,2.,&
-                                  2.,2.,2.       /)                          
-!
-TPROFILER%NAME               = (/ 'CAAM    ', 'CRAU    ', 'BARD    ', 'MONT    ', 'IUTF    ', &
-                                  'OBSF    ', 'VALF    ', 'LUBE    ', 'VENT    ', 'DENT    ', &
-                                  'CARP    ', 'DUPA    ', 'VINO    ', 'CHAM    ', 'REA1    ', &
-                                  'REA2    ', 'REA3    ',                                     &
-                                  'SOD_4M  ', 'UHF_4M  ', 'VHF_STM ', 'SOD_CO  ', 'UHF_DEG ', &
-                                  'SOD_EC  ', 'SOD_E1  ', 'SOD_E2  ', 'UHF_ED  ', 'VHF_LS  ', &
-                                  'UHF_DEP ',                                                 &
-                                  'AGRI    ', 'ALLA    ', 'BOYE    ', 'CANE    ', 'CHAT    ', &
-                                  'CINQ    ', 'COTE    ', 'DR12    ', 'ETOI    ', 'IUT1    ', &
-                                  'JULI    ', 'ONYX    ', 'PENN    ', 'TRIB    ', 'VALB    ', &
-                                  'VALL    ', 'MARS    ', 'AIRB    ' /)                                               
-!
-TPROFILER%TYPE               = (/ 'flux    ', 'flux    ', 'flux    ', 'flux    ', 'flux    ', &
-                                  'flux    ', 'flux    ', 'flux    ', 'flux    ', 'flux    ', &
-                                  'flux    ', 'flux    ', 'flux    ', 'flux    ', 'flux    ', &
-                                  'flux    ', 'flux    ',                                     &
-                                  'sodar   ', 'uhf     ', 'vhf     ', 'sodar   ', 'uhf     ', &
-                                  'sodar   ', 'sodar   ', 'sodar   ', 'uhf     ', 'vhf     ', &
-                                  'uhf     ',                                                 &
-                                  'gps     ', 'gps     ', 'gps     ', 'gps     ', 'gps     ', &
-                                  'gps     ', 'gps     ', 'gps     ', 'gps     ', 'gps     ', &
-                                  'gps     ', 'gps     ', 'gps     ', 'gps     ', 'gps     ', &
-                                  'gps     ', 'gps     ', 'gps     ' /)                     
-!
-!----------------------------------------------------------------------------
-ENDIF
-!
-END SUBROUTINE INI_PROFILER_n
diff --git a/src/MNH/ini_surfstationn.f90 b/src/MNH/ini_surfstationn.f90
index 8811fa566..5452a5f0e 100644
--- a/src/MNH/ini_surfstationn.f90
+++ b/src/MNH/ini_surfstationn.f90
@@ -68,11 +68,11 @@ USE MODD_PARAMETERS,     ONLY: JPHEXT, JPVEXT
 USE MODD_STATION_n
 USE MODD_TYPE_STATPROF
 !
-USE MODE_MSG
-USE MODE_STATPROF_READER
-USE MODE_STATPROF_TOOLS,  ONLY: STATION_ADD, STATION_ALLOCATE, STATION_INI_INTERP, STATION_POSITION
 USE MODE_ALLOCBUFFER_ll,  ONLY: ALLOCBUFFER_ll
 USE MODE_GATHER_ll,       ONLY: GATHERALL_FIELD_ll
+USE MODE_MSG
+USE MODE_STATPROF_READER, ONLY: STATPROF_CSV_READ
+USE MODE_STATPROF_TOOLS,  ONLY: STATION_ADD, STATION_ALLOCATE, STATPROF_INI_INTERP, STATPROF_POSITION
 !
 IMPLICIT NONE
 !
@@ -88,6 +88,7 @@ IMPLICIT NONE
 INTEGER :: IERR
 INTEGER :: IIU
 INTEGER :: IJU
+INTEGER :: INUMBSTAT                        ! Total number of stations (inside physical domain of model)
 INTEGER :: ISTORE                           ! number of storage instants
 INTEGER :: JI
 LOGICAL :: GALLOCX, GALLOCY
@@ -120,7 +121,7 @@ IF ( CFILE_STAT /= "NO_INPUT_CSV" .OR. NNUMB_STAT > 0 ) THEN
   IIU = SIZE( XXHAT )
   IJU = SIZE( XYHAT )
 
-  ! Get global XHAT and YHAT (needed by STATION_POSITION)
+  ! Get global XHAT and YHAT (needed by STATPROF_POSITION)
   CALL ALLOCBUFFER_ll( ZXHAT_GLOB, XXHAT, 'XX', GALLOCX )
   CALL ALLOCBUFFER_ll( ZYHAT_GLOB, XYHAT, 'YY', GALLOCY )
   CALL GATHERALL_FIELD_ll( 'XX', XXHAT, ZXHAT_GLOB, IERR )
@@ -145,7 +146,7 @@ NUMBSTAT_LOC = 0
 
 IF (CFILE_STAT=="NO_INPUT_CSV") THEN
   ! Treat namelist
-  NUMBSTAT = 0
+  INUMBSTAT = 0
   IF ( NNUMB_STAT > 0 ) THEN
     DO JI = 1, NNUMB_STAT
       IF ( LCARTESIAN ) THEN
@@ -154,18 +155,18 @@ IF (CFILE_STAT=="NO_INPUT_CSV") THEN
       ELSE
         TZSTATION%XLAT = XLAT_STAT(JI)
         TZSTATION%XLON = XLON_STAT(JI)
-        CALL STATION_INI_INTERP( TZSTATION )
+        CALL STATPROF_INI_INTERP( TZSTATION )
       END IF
       TZSTATION%XZ    = XZ_STAT(JI)
       TZSTATION%CNAME = CNAME_STAT(JI)
 
-      CALL STATION_POSITION( TZSTATION, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,                  &
+      CALL STATPROF_POSITION( TZSTATION, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,                 &
                              ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX, ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
                              GINSIDE, GPRESENT                                                   )
 
       IF ( GINSIDE ) THEN
-        NUMBSTAT = NUMBSTAT + 1
-        TZSTATION%NID = NUMBSTAT
+        INUMBSTAT = INUMBSTAT + 1
+        TZSTATION%NID = INUMBSTAT
       END IF
 
       IF ( GPRESENT ) CALL STATION_ADD( TZSTATION )
@@ -173,11 +174,12 @@ IF (CFILE_STAT=="NO_INPUT_CSV") THEN
   END IF
 ELSE
   !Treat CSV datafile
-  CALL READ_CSV_STATION( CFILE_STAT, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,               &
-                         ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX,ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX )
+  CALL STATPROF_CSV_READ( TZSTATION, CFILE_STAT, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,     &
+                          ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX,ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
+                          INUMBSTAT                                                          )
 END IF
 
-LSTATION = ( NUMBSTAT > 0 )
+LSTATION = ( INUMBSTAT > 0 )
 
 DO JI = 1, NUMBSTAT_LOC
   CALL STATION_ALLOCATE( TSTATIONS(JI), ISTORE )
diff --git a/src/MNH/modd_profilern.f90 b/src/MNH/modd_profilern.f90
index d8fb2469f..3abfe6611 100644
--- a/src/MNH/modd_profilern.f90
+++ b/src/MNH/modd_profilern.f90
@@ -1,23 +1,18 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-2022 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.
 !-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 modd 2006/05/18 13:07:25
-!-----------------------------------------------------------------
 !     ############################
       MODULE MODD_PROFILER_n
 !     ############################
 !
-!!****  *MODD_PROFILER* - declaration of stations
+!!****  *MODD_PROFILER* - declaration of profilers
 !!
 !!    PURPOSE
 !!    -------
 !       The purpose of this declarative module is to define
-!      the different stations types.
+!      the different profilers types.
 !
 !!
 !!**  IMPLICIT ARGUMENTS
@@ -34,32 +29,42 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original    15/01/02
+!  P. Wautelet    04/2022: restructure profilers for better performance, reduce memory usage and correct some problems/bugs
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
 !             ------------
 !
 !
-USE MODD_TYPE_PROFILER
-USE MODD_PARAMETERS, ONLY: JPMODELMAX
+USE MODD_PARAMETERS,    ONLY: JPMODELMAX
+USE MODD_TYPE_STATPROF, ONLY: TPROFILERDATA, TSTATPROFTIME
+
 IMPLICIT NONE
 
+PRIVATE
+
+PUBLIC :: LPROFILER, NUMBPROFILER_LOC, TPROFILERS_TIME, TPROFILERS
+
+PUBLIC :: PROFILER_GOTO_MODEL
+
 TYPE PROFILER_t
 !
 !-------------------------------------------------------------------------------------------
 !
-  LOGICAL                          :: LPROFILER    ! flag to use stations
-  INTEGER                          :: NUMBPROFILER    ! number of stations
+  LOGICAL                          :: LPROFILER    ! flag to use profilers
+  INTEGER                          :: NUMBPROFILER_LOC = 0 ! number of profilers on this process
 !
-  TYPE(PROFILER) :: TPROFILER ! characteristics and records of an aircraft
+  TYPE(TSTATPROFTIME) :: TPROFILERS_TIME
+  TYPE(TPROFILERDATA), DIMENSION(:), POINTER :: TPROFILERS ! characteristics and records of the profilers
 !
 END TYPE PROFILER_t
 
 TYPE(PROFILER_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: PROFILER_MODEL
 
 LOGICAL, POINTER :: LPROFILER=>NULL()
-INTEGER, POINTER :: NUMBPROFILER=>NULL()
-TYPE(PROFILER), POINTER :: TPROFILER=>NULL()
+INTEGER, POINTER :: NUMBPROFILER_LOC=>NULL()
+TYPE(TSTATPROFTIME),               POINTER :: TPROFILERS_TIME => NULL()
+TYPE(TPROFILERDATA), DIMENSION(:), POINTER :: TPROFILERS      => NULL()
 
 CONTAINS
 
@@ -67,11 +72,13 @@ SUBROUTINE PROFILER_GOTO_MODEL(KFROM, KTO)
 INTEGER, INTENT(IN) :: KFROM, KTO
 !
 ! Save current state for allocated arrays
+PROFILER_MODEL(KFROM)%TPROFILERS => TPROFILERS
 !
 ! Current model is set to model KTO
-LPROFILER=>PROFILER_MODEL(KTO)%LPROFILER
-NUMBPROFILER=>PROFILER_MODEL(KTO)%NUMBPROFILER
-TPROFILER=>PROFILER_MODEL(KTO)%TPROFILER
+LPROFILER        => PROFILER_MODEL(KTO)%LPROFILER
+NUMBPROFILER_LOC => PROFILER_MODEL(KTO)%NUMBPROFILER_LOC
+TPROFILERS_TIME  => PROFILER_MODEL(KTO)%TPROFILERS_TIME
+TPROFILERS       => PROFILER_MODEL(KTO)%TPROFILERS
 
 END SUBROUTINE PROFILER_GOTO_MODEL
 
diff --git a/src/MNH/modd_stationn.f90 b/src/MNH/modd_stationn.f90
index bd86a1595..406d14909 100644
--- a/src/MNH/modd_stationn.f90
+++ b/src/MNH/modd_stationn.f90
@@ -43,7 +43,7 @@ IMPLICIT NONE
 
 PRIVATE
 
-PUBLIC :: LSTATION, NUMBSTAT, NUMBSTAT_LOC, TSTATIONS_TIME, TSTATIONS
+PUBLIC :: LSTATION, NUMBSTAT_LOC, TSTATIONS_TIME, TSTATIONS
 
 PUBLIC :: STATION_GOTO_MODEL
 
@@ -52,7 +52,6 @@ TYPE STATION_t
 !-------------------------------------------------------------------------------------------
 !
   LOGICAL                          :: LSTATION    ! flag to use stations
-  INTEGER                          :: NUMBSTAT    ! number of stations
   INTEGER                          :: NUMBSTAT_LOC = 0 ! number of stations on this process
 !
   TYPE(TSTATPROFTIME) :: TSTATIONS_TIME
@@ -63,7 +62,6 @@ END TYPE STATION_t
 TYPE(STATION_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: STATION_MODEL
 
 LOGICAL, POINTER :: LSTATION=>NULL()
-INTEGER, POINTER :: NUMBSTAT=>NULL()
 INTEGER, POINTER :: NUMBSTAT_LOC=>NULL()
 TYPE(TSTATPROFTIME),              POINTER :: TSTATIONS_TIME => NULL()
 TYPE(TSTATIONDATA), DIMENSION(:), POINTER :: TSTATIONS      => NULL()
@@ -78,7 +76,6 @@ STATION_MODEL(KFROM)%TSTATIONS => TSTATIONS
 !
 ! Current model is set to model KTO
 LSTATION       => STATION_MODEL(KTO)%LSTATION
-NUMBSTAT       => STATION_MODEL(KTO)%NUMBSTAT
 NUMBSTAT_LOC   => STATION_MODEL(KTO)%NUMBSTAT_LOC
 TSTATIONS_TIME => STATION_MODEL(KTO)%TSTATIONS_TIME
 TSTATIONS      => STATION_MODEL(KTO)%TSTATIONS
diff --git a/src/MNH/modd_sub_profilern.f90 b/src/MNH/modd_sub_profilern.f90
deleted file mode 100644
index 14582b74b..000000000
--- a/src/MNH/modd_sub_profilern.f90
+++ /dev/null
@@ -1,109 +0,0 @@
-!MNH_LIC Copyright 1994-2014 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.
-!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 modd 2006/06/27 12:30:56
-!-----------------------------------------------------------------
-!     ############################
-      MODULE MODD_SUB_PROFILER_n
-!     ############################
-!
-!!****  *MODD_PROFILER* - declaration of stations
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this declarative module is to define
-!      the different stations types.
-!
-!!
-!!**  IMPLICIT ARGUMENTS
-!!    ------------------
-!!      NONE 
-!!
-!!    REFERENCE
-!!    --------- 
-!!       
-!!    AUTHOR
-!!    ------
-!!	P. Tulet   *Meteo France*
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    15/01/02
-!-------------------------------------------------------------------------------
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-!
-USE MODD_PARAMETERS, ONLY: JPMODELMAX
-IMPLICIT NONE
-
-TYPE SUB_PROFILER_t
-!
-!-------------------------------------------------------------------------------------------
-!
-  LOGICAL :: GPROFILERFIRSTCALL  = .TRUE.
-!
-  INTEGER,DIMENSION(:), POINTER :: II=>NULL()       ! mass lidar position (x index)
-  INTEGER,DIMENSION(:), POINTER :: IJ=>NULL()       ! mass lidar position (y index)
-  INTEGER,DIMENSION(:), POINTER :: IU=>NULL()       ! U flux point lidar position (x index)
-  INTEGER,DIMENSION(:), POINTER :: IV=>NULL()       ! V flux point lidar position (y index)
-!
-  REAL, DIMENSION(:), POINTER :: ZTHIS_PROCS=>NULL()
-!
-  REAL,DIMENSION(:), POINTER :: ZXCOEF=>NULL()   ! X direction interpolation coefficient
-  REAL,DIMENSION(:), POINTER :: ZUCOEF=>NULL()   ! X direction interpolation coefficient (for U)
-  REAL,DIMENSION(:), POINTER :: ZYCOEF=>NULL()   ! Y direction interpolation coefficient
-  REAL,DIMENSION(:), POINTER :: ZVCOEF=>NULL()   ! Y direction interpolation coefficient (for V)
-
-END TYPE SUB_PROFILER_t
-
-TYPE(SUB_PROFILER_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: SUB_PROFILER_MODEL
-
-LOGICAL, POINTER :: GPROFILERFIRSTCALL=>NULL()
-INTEGER,DIMENSION(:), POINTER :: II=>NULL()
-INTEGER,DIMENSION(:), POINTER :: IJ=>NULL()
-INTEGER,DIMENSION(:), POINTER :: IU=>NULL()
-INTEGER,DIMENSION(:), POINTER :: IV=>NULL()
-REAL, DIMENSION(:), POINTER :: ZTHIS_PROCS=>NULL()
-REAL,DIMENSION(:), POINTER :: ZXCOEF=>NULL()
-REAL,DIMENSION(:), POINTER :: ZUCOEF=>NULL()
-REAL,DIMENSION(:), POINTER :: ZYCOEF=>NULL()
-REAL,DIMENSION(:), POINTER :: ZVCOEF=>NULL()
-
-CONTAINS
-
-SUBROUTINE SUB_PROFILER_GOTO_MODEL(KFROM, KTO)
-INTEGER, INTENT(IN) :: KFROM, KTO
-!
-! Save current state for allocated arrays
-SUB_PROFILER_MODEL(KFROM)%II=>II
-SUB_PROFILER_MODEL(KFROM)%IJ=>IJ
-SUB_PROFILER_MODEL(KFROM)%IU=>IU
-SUB_PROFILER_MODEL(KFROM)%IV=>IV
-SUB_PROFILER_MODEL(KFROM)%ZTHIS_PROCS=>ZTHIS_PROCS
-SUB_PROFILER_MODEL(KFROM)%ZXCOEF=>ZXCOEF
-SUB_PROFILER_MODEL(KFROM)%ZUCOEF=>ZUCOEF
-SUB_PROFILER_MODEL(KFROM)%ZYCOEF=>ZYCOEF
-SUB_PROFILER_MODEL(KFROM)%ZVCOEF=>ZVCOEF
-!
-! Current model is set to model KTO
-GPROFILERFIRSTCALL=>SUB_PROFILER_MODEL(KTO)%GPROFILERFIRSTCALL
-II=>SUB_PROFILER_MODEL(KTO)%II
-IJ=>SUB_PROFILER_MODEL(KTO)%IJ
-IU=>SUB_PROFILER_MODEL(KTO)%IU
-IV=>SUB_PROFILER_MODEL(KTO)%IV
-ZTHIS_PROCS=>SUB_PROFILER_MODEL(KTO)%ZTHIS_PROCS
-ZXCOEF=>SUB_PROFILER_MODEL(KTO)%ZXCOEF
-ZUCOEF=>SUB_PROFILER_MODEL(KTO)%ZUCOEF
-ZYCOEF=>SUB_PROFILER_MODEL(KTO)%ZYCOEF
-ZVCOEF=>SUB_PROFILER_MODEL(KTO)%ZVCOEF
-
-END SUBROUTINE SUB_PROFILER_GOTO_MODEL
-
-END MODULE MODD_SUB_PROFILER_n
diff --git a/src/MNH/modd_type_profiler.f90 b/src/MNH/modd_type_profiler.f90
deleted file mode 100644
index b5fedbf60..000000000
--- a/src/MNH/modd_type_profiler.f90
+++ /dev/null
@@ -1,115 +0,0 @@
-!MNH_LIC Copyright 2002-2021 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 MODD_TYPE_PROFILER
-!     ############################
-!
-!!****  *MODD_PROFILER* - declaration of stations
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this declarative module is to define
-!      the different stations types.
-!
-!!
-!!**  IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!
-!!    REFERENCE
-!!    --------- 
-!!       
-!!    AUTHOR
-!!    ------
-!!	P. Tulet   *Meteo France*
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original    15/01/02
-!!     C.Lac 10/2016  Add visibility diagnostic
-!  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
-!  M. Taufour  05/07/2021: modify RARE for hydrometeors containing ice and add bright band calculation for RARE
-!-------------------------------------------------------------------------------
-!
-!*       0.   DECLARATIONS
-!             ------------
-!
-use modd_type_date, only: date_time
-
-implicit none
-
-TYPE PROFILER
-!
-!
-!* general information
-!
-!
-!* storage monitoring
-!
-REAL                          :: T_CUR  ! current time since last storage
-INTEGER                       :: N_CUR  ! current step of storage
-REAL                          :: STEP   ! storage time step
-!
-!* data records
-!
-CHARACTER(LEN=8),DIMENSION(:), POINTER   :: NAME=>NULL()   ! station name
-CHARACTER(LEN=8),DIMENSION(:), POINTER   :: TYPE=>NULL()   ! station type
-!
-type(date_time), dimension(:), pointer :: tpdates => NULL() ! dates(n) (n: recording instants)
-LOGICAL, DIMENSION(:),    POINTER :: ERROR=>NULL()
-REAL, DIMENSION(:),       POINTER :: X=>NULL()        ! X(n)
-REAL, DIMENSION(:),       POINTER :: Y=>NULL()        ! Y(n)
-REAL, DIMENSION(:),       POINTER :: LON=>NULL()      ! longitude(n)
-REAL, DIMENSION(:),       POINTER :: LAT=>NULL()      ! latitude (n)
-REAL, DIMENSION(:),       POINTER :: ALT=>NULL()      ! altitude (n)
-REAL, DIMENSION(:,:,:),   POINTER :: ZON=>NULL()      ! zonal wind(n)
-REAL, DIMENSION(:,:,:),   POINTER :: MER=>NULL()      ! meridian wind(n)
-REAL, DIMENSION(:,:,:),   POINTER :: FF=>NULL()       ! wind intensity
-REAL, DIMENSION(:,:,:),   POINTER :: DD=>NULL()       ! wind direction
-REAL, DIMENSION(:,:,:),   POINTER :: W=>NULL()        ! w(n)  (air vertical speed)
-REAL, DIMENSION(:,:,:),   POINTER :: P=>NULL()        ! p(n)
-REAL, DIMENSION(:,:,:),   POINTER :: ZZ=>NULL()       ! altitude(n)
-REAL, DIMENSION(:,:,:),   POINTER :: TKE=>NULL()      ! tke(n)
-REAL, DIMENSION(:,:,:),   POINTER :: TH=>NULL()       ! th(n)
-REAL, DIMENSION(:,:,:),   POINTER :: THV=>NULL()      ! thv(n)
-REAL, DIMENSION(:,:,:),   POINTER :: VISI=>NULL()     ! VISI(n)
-REAL, DIMENSION(:,:,:),   POINTER :: VISIKUN=>NULL()  ! VISI KUNKEL(n)
-REAL, DIMENSION(:,:,:),   POINTER :: CRARE=>NULL()     ! radar reflectivity (n)
-REAL, DIMENSION(:,:,:),   POINTER :: CRARE_ATT=>NULL() ! radar attenuated reflectivity (n)
-REAL, DIMENSION(:,:,:),   POINTER :: CIZ=>NULL()       ! Ice number concentration ICE3 (n)
-REAL, DIMENSION(:,:,:),   POINTER :: LWCZ=>NULL()      ! liquid water content (n)
-REAL, DIMENSION(:,:,:),   POINTER :: IWCZ=>NULL()      ! ice water content (n)
-REAL, DIMENSION(:,:,:),   POINTER :: RHOD=>NULL()     ! density of dry air/moist air
-REAL, DIMENSION(:,:,:,:), POINTER :: R=>NULL()        ! r*(n)
-REAL, DIMENSION(:,:,:,:), POINTER :: SV=>NULL()       ! Sv*(n)
-REAL, DIMENSION(:,:,:,:), POINTER :: AER=>NULL()      ! AER*(n) aerosol extinction
-!
-REAL, DIMENSION(:,:),     POINTER :: T2M=>NULL()      ! 2 m air temperature (°C)
-REAL, DIMENSION(:,:),     POINTER :: Q2M=>NULL()      ! 2 m humidity (kg/kg)
-REAL, DIMENSION(:,:),     POINTER :: HU2M=>NULL()     ! 2 m relative humidity (%)
-REAL, DIMENSION(:,:),     POINTER :: ZON10M=>NULL()   ! 10 m zonal wind (m/s)  
-REAL, DIMENSION(:,:),     POINTER :: MER10M=>NULL()   ! 10 m merid. wind (m/s)
-REAL, DIMENSION(:,:),     POINTER :: RN=>NULL()       ! net radiation (W m2)
-REAL, DIMENSION(:,:),     POINTER :: H=>NULL()        ! sensible heat flux (W m2)
-REAL, DIMENSION(:,:),     POINTER :: LE=>NULL()       ! Total latent heat flux (W m2)
-REAL, DIMENSION(:,:),     POINTER :: LEI=>NULL()      ! Solid latent heat flux (W m2)
-REAL, DIMENSION(:,:),     POINTER :: GFLUX=>NULL()    ! storage heat flux (W m2)
-REAL, DIMENSION(:,:),     POINTER :: LWD=>NULL()       ! IR downward radiation (W m2)
-REAL, DIMENSION(:,:),     POINTER :: LWU=>NULL()       ! IR upward radiation (W m2)
-REAL, DIMENSION(:,:),     POINTER :: SWD=>NULL()       ! solar downward radiation (W m2)
-REAL, DIMENSION(:,:),     POINTER :: SWU=>NULL()       ! solar upward radiation (W m2)
-REAL, DIMENSION(:,:),   POINTER :: IWV=>NULL()      ! integrated water vpour(n)
-REAL, DIMENSION(:,:),   POINTER :: ZTD=>NULL()      ! GPS zenith tropo delay(n)
-REAL, DIMENSION(:,:),   POINTER :: ZWD=>NULL()      ! GPS zenith wet delay(n)
-REAL, DIMENSION(:,:),   POINTER :: ZHD=>NULL()      ! GPS zenith hydro delay(n)
-!
-REAL, DIMENSION(:,:,:),   POINTER :: TKE_DISS=>NULL() ! TKE dissipation rate
-!
-!
-END TYPE PROFILER
-!
-END MODULE MODD_TYPE_PROFILER
-
diff --git a/src/MNH/modd_type_statprof.f90 b/src/MNH/modd_type_statprof.f90
index 12150549a..c456ae947 100644
--- a/src/MNH/modd_type_statprof.f90
+++ b/src/MNH/modd_type_statprof.f90
@@ -30,7 +30,7 @@
 !!    -------------
 !!      Original    15/01/02
 !  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
-!  P. Wautelet    04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
+!  P. Wautelet    04/2022: restructure stations/profilers for better performance, reduce memory usage and correct some problems/bugs
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -43,7 +43,8 @@ implicit none
 
 private
 
-public :: TSTATPROFTIME, TSTATIONDATA
+public :: TSTATPROFTIME
+public :: TPROFILERDATA, TSTATIONDATA, TSTATPROFDATA
 
 TYPE :: TSTATPROFTIME
   REAL                                       :: XTIME_CUR = XUNDEF  ! current time since last storage
@@ -52,63 +53,105 @@ TYPE :: TSTATPROFTIME
   type(date_time), dimension(:), ALLOCATABLE :: tpdates             ! dates(n) (n: recording instants)
 END TYPE TSTATPROFTIME
 
-TYPE TSTATIONDATA
-! Type to store all the data of 1 station
-CHARACTER(LEN=NSTATPROFNAMELGTMAX) :: CNAME = ''  ! station name
-
-
-INTEGER :: NID = 0 ! Global identification number of the station (from 1 to total number of stations of the model)
-
-REAL :: XX   = XUNDEF  ! X(n)
-REAL :: XY   = XUNDEF  ! Y(n)
-REAL :: XZ   = XUNDEF  ! Z(n)
-REAL :: XLON = XUNDEF  ! longitude(n)
-REAL :: XLAT = XUNDEF  ! latitude (n)
-REAL :: XZS  = XUNDEF  ! zs(n)
-
-! Position in the mesh
-INTEGER :: NI_M = NNEGUNDEF ! X position for mass-point axis (between this one and the next one)
-INTEGER :: NJ_M = NNEGUNDEF ! Y position for mass-point axis (between this one and the next one)
-INTEGER :: NI_U = NNEGUNDEF ! X position for u-point axis (between this one and the next one)
-INTEGER :: NJ_V = NNEGUNDEF ! Y position for v-point axis (between this one and the next one)
-
-! Coefficient to interpolate values (stations are usually not exactly on mesh points)
-REAL :: XXMCOEF = XUNDEF ! Interpolation coefficient for X (mass-point)
-REAL :: XYMCOEF = XUNDEF ! Interpolation coefficient for Y (mass-point)
-REAL :: XXUCOEF = XUNDEF ! Interpolation coefficient for X (U-point)
-REAL :: XYVCOEF = XUNDEF ! Interpolation coefficient for Y (V-point)
-
-INTEGER :: NK = NNEGUNDEF ! Model level for altitude comparisons
-
-REAL, DIMENSION(:),   ALLOCATABLE :: XZON    ! zonal wind(n)
-REAL, DIMENSION(:),   ALLOCATABLE :: XMER    ! meridian wind(n)
-REAL, DIMENSION(:),   ALLOCATABLE :: XW      ! w(n)  (air vertical speed)
-REAL, DIMENSION(:),   ALLOCATABLE :: XP      ! p(n)
-REAL, DIMENSION(:),   ALLOCATABLE :: XTKE    ! tke(n)
-REAL, DIMENSION(:),   ALLOCATABLE :: XTH     ! th(n)
-REAL, DIMENSION(:,:), ALLOCATABLE :: XR      ! r*(n)
-REAL, DIMENSION(:,:), ALLOCATABLE :: XSV     ! Sv*(n)
-REAL, DIMENSION(:),   ALLOCATABLE :: XTSRAD  ! Ts(n)
-
-REAL, DIMENSION(:),   ALLOCATABLE :: XT2M
-REAL, DIMENSION(:),   ALLOCATABLE :: XQ2M
-REAL, DIMENSION(:),   ALLOCATABLE :: XHU2M
-REAL, DIMENSION(:),   ALLOCATABLE :: XZON10M
-REAL, DIMENSION(:),   ALLOCATABLE :: XMER10M
-REAL, DIMENSION(:),   ALLOCATABLE :: XRN
-REAL, DIMENSION(:),   ALLOCATABLE :: XH
-REAL, DIMENSION(:),   ALLOCATABLE :: XLE
-REAL, DIMENSION(:),   ALLOCATABLE :: XLEI
-REAL, DIMENSION(:),   ALLOCATABLE :: XGFLUX
-REAL, DIMENSION(:),   ALLOCATABLE :: XSWD
-REAL, DIMENSION(:),   ALLOCATABLE :: XSWU
-REAL, DIMENSION(:),   ALLOCATABLE :: XLWD
-REAL, DIMENSION(:),   ALLOCATABLE :: XLWU
-REAL, DIMENSION(:),   ALLOCATABLE :: XSWDIR
-REAL, DIMENSION(:),   ALLOCATABLE :: XSWDIFF
-REAL, DIMENSION(:),   ALLOCATABLE :: XDSTAOD ! Dust Aerosol Optical Depth
-REAL, DIMENSION(:),   ALLOCATABLE :: XSFCO2  ! CO2 surface flux
+TYPE :: TSTATPROFDATA
+  ! Type to store data common to stations and profilers
+  ! It is used as a basis for the TSTATIONDATA and TPROFILERDATA
+  ! and for common procedures for these 2 types
+  CHARACTER(LEN=NSTATPROFNAMELGTMAX) :: CNAME = ''  ! Station/profiler name
 
+  INTEGER :: NID = 0 ! Global identification number of the station/profiler (from 1 to total number)
+
+  REAL :: XX   = XUNDEF  ! X(n)
+  REAL :: XY   = XUNDEF  ! Y(n)
+  REAL :: XZ   = XUNDEF  ! Z(n)
+  REAL :: XLON = XUNDEF  ! longitude(n)
+  REAL :: XLAT = XUNDEF  ! latitude (n)
+
+  ! Position in the mesh
+  INTEGER :: NI_M = NNEGUNDEF ! X position for mass-point axis (between this one and the next one)
+  INTEGER :: NJ_M = NNEGUNDEF ! Y position for mass-point axis (between this one and the next one)
+  INTEGER :: NI_U = NNEGUNDEF ! X position for u-point axis (between this one and the next one)
+  INTEGER :: NJ_V = NNEGUNDEF ! Y position for v-point axis (between this one and the next one)
+
+  ! Coefficient to interpolate values (stations are usually not exactly on mesh points)
+  REAL :: XXMCOEF = XUNDEF ! Interpolation coefficient for X (mass-point)
+  REAL :: XYMCOEF = XUNDEF ! Interpolation coefficient for Y (mass-point)
+  REAL :: XXUCOEF = XUNDEF ! Interpolation coefficient for X (U-point)
+  REAL :: XYVCOEF = XUNDEF ! Interpolation coefficient for Y (V-point)
+
+  ! Dimension corresponds to recording instants
+  REAL, DIMENSION(:),   ALLOCATABLE :: XT2M    ! 2 m air temperature (C)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XQ2M    ! 2 m humidity (kg/kg)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XHU2M   ! 2 m relative humidity (%)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XZON10M ! 10 m zonal wind (m/s)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XMER10M ! 10 m merid. wind (m/s)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XRN     ! net radiation (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XH      ! sensible heat flux (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XLE     ! Total latent heat flux (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XLEI    ! Solid latent heat flux (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XGFLUX  ! storage heat flux (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XSWD    ! IR downward radiation (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XSWU    ! IR upward radiation (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XLWD    ! solar downward radiation (W m2)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XLWU    ! solar upward radiation (W m2)
+END TYPE
+
+TYPE, EXTENDS( TSTATPROFDATA ) ::  TSTATIONDATA
+  ! Type to store all the data of 1 station
+  INTEGER :: NK = NNEGUNDEF ! Model level for altitude comparisons
+
+  REAL :: XZS  = XUNDEF  ! zs(n)
+
+  ! (n: recording instants)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XZON    ! zonal wind(n)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XMER    ! meridian wind(n)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XW      ! w(n)  (air vertical speed)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XP      ! p(n)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XTKE    ! tke(n)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XTH     ! th(n)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: XR      ! r*(n)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: XSV     ! Sv*(n)
+  REAL, DIMENSION(:),   ALLOCATABLE :: XTSRAD  ! Ts(n)
+
+  REAL, DIMENSION(:),   ALLOCATABLE :: XSWDIR
+  REAL, DIMENSION(:),   ALLOCATABLE :: XSWDIFF
+  REAL, DIMENSION(:),   ALLOCATABLE :: XDSTAOD ! Dust Aerosol Optical Depth
+  REAL, DIMENSION(:),   ALLOCATABLE :: XSFCO2  ! CO2 surface flux
 END TYPE TSTATIONDATA
 
+TYPE, EXTENDS( TSTATPROFDATA ) ::  TPROFILERDATA
+  ! Type to store all the data of 1 profiler
+  CHARACTER(LEN=NSTATPROFNAMELGTMAX) :: CTYPE = ''  ! Profiler type
+
+  ! (n: recording instants)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XZON       ! zonal wind(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XMER       ! meridian wind(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XFF        ! wind intensity
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XDD        ! wind direction
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XW         ! w(n)  (air vertical speed)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XP         ! p(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XZZ        ! altitude(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XTKE       ! tke(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XTH        ! th(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XTHV       ! thv(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XVISI      ! VISI(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XVISIKUN   ! VISI KUNKEL(n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XCRARE     ! radar reflectivity (n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XCRARE_ATT ! radar attenuated reflectivity (n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XCIZ       ! Ice number concentration ICE3 (n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XLWCZ      ! liquid water content (n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XIWCZ      ! ice water content (n)
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: XRHOD      ! density of dry air/moist air
+  REAL, DIMENSION(:,:,:), ALLOCATABLE :: XR         ! r*(n)
+  REAL, DIMENSION(:,:,:), ALLOCATABLE :: XSV        ! Sv*(n)
+  REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAER       ! AER*(n) aerosol extinction
+
+  REAL, DIMENSION(:), ALLOCATABLE :: XIWV ! integrated water vpour(n)
+  REAL, DIMENSION(:), ALLOCATABLE :: XZTD ! GPS zenith tropo delay(n)
+  REAL, DIMENSION(:), ALLOCATABLE :: XZWD ! GPS zenith wet delay(n)
+  REAL, DIMENSION(:), ALLOCATABLE :: XZHD ! GPS zenith hydro delay(n)
+
+  REAL, DIMENSION(:,:), ALLOCATABLE :: XTKE_DISS ! TKE dissipation rate
+END TYPE
+
 END MODULE MODD_TYPE_STATPROF
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index 61c50c483..c778b2e1e 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -2116,21 +2116,20 @@ IF (LFLYER)                                                                   &
 !*       24.2    STATION (observation diagnostic)
 !               --------------------------------
 !
-IF (LSTATION)                                                            &
-  CALL STATION_n(XTSTEP,                                                 &
-                 XXHAT, XYHAT, XZZ,                                      &
-                 XUT, XVT, XWT, XTHT, XRT, XSVT, XTKET, XTSRAD, XPABST   )
+IF ( LSTATION )                                                        &
+  CALL STATION_n(XTSTEP, XZZ,                                          &
+                 XUT, XVT, XWT, XTHT, XRT, XSVT, XTKET, XTSRAD, XPABST )
 !
 !---------------------------------------------------------
 !
 !*       24.3    PROFILER (observation diagnostic)
 !               ---------------------------------
 !
-IF (LPROFILER)                                                           &
+IF ( LPROFILER )                                                         &
   CALL PROFILER_n(XTSTEP,                                                &
-                  XXHAT, XYHAT, XZZ,XRHODREF,                            &
+                  XZZ, XRHODREF,                                         &
                   XUT, XVT, XWT, XTHT, XRT, XSVT, XTKET, XTSRAD, XPABST, &
-                  XAER, XCLDFR, XCIT,PSEA=ZSEA(:,:))
+                  XAER, XCIT, PSEA=ZSEA(:,:))
 !
 !
 CALL SECOND_MNH2(ZTIME2)
diff --git a/src/MNH/modn_stationn.f90 b/src/MNH/modn_stationn.f90
index 9b8737c03..094b3dbb1 100644
--- a/src/MNH/modn_stationn.f90
+++ b/src/MNH/modn_stationn.f90
@@ -37,7 +37,6 @@ USE MODD_ALLSTATION_n, ONLY:&
         CFILE_STAT_n    =>CFILE_STAT    ,&
         LDIAG_SURFRAD_n =>LDIAG_SURFRAD
 USE MODD_PARAMETERS, ONLY: NFILENAMELGTMAX, NSTATPROFNAMELGTMAX
-USE MODD_STATION_n
 !
 !-----------------------------------------------------------------------------
 !
diff --git a/src/MNH/profilern.f90 b/src/MNH/profilern.f90
index 4424b816f..ccd26a860 100644
--- a/src/MNH/profilern.f90
+++ b/src/MNH/profilern.f90
@@ -9,27 +9,24 @@ MODULE MODI_PROFILER_n
 !
 INTERFACE
 !
-      SUBROUTINE PROFILER_n(PTSTEP,                               &
-                            PXHAT, PYHAT, PZ,PRHODREF,            &
-                            PU, PV, PW, PTH, PR, PSV, PTKE,       &
-                            PTS, PP, PAER, PCLDFR, PCIT, PSEA)
+      SUBROUTINE PROFILER_n( PTSTEP,                         &
+                             PZ, PRHODREF,                   &
+                             PU, PV, PW, PTH, PR, PSV, PTKE, &
+                             PTS, PP, PAER, PCIT, PSEA )
 !
 REAL,                     INTENT(IN)     :: PTSTEP ! time step
-REAL, DIMENSION(:),       INTENT(IN)     :: PXHAT  ! x coordinate
-REAL, DIMENSION(:),       INTENT(IN)     :: PYHAT  ! y coordinate
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PZ     ! z array
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PU     ! horizontal wind X component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PV     ! horizontal wind Y component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PW     ! vertical wind
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PTH    ! potential temperature
-REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PRHODREF                            
+REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PRHODREF
 REAL, DIMENSION(:,:,:,:), INTENT(IN)     :: PR     ! water mixing ratios
 REAL, DIMENSION(:,:,:,:), INTENT(IN)     :: PSV    ! Scalar variables
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PTKE   ! turbulent kinetic energy
 REAL, DIMENSION(:,:),     INTENT(IN)     :: PTS    ! surface temperature
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PP     ! pressure
 REAL, DIMENSION(:,:,:,:), INTENT(IN)     :: PAER   ! aerosol extinction
-REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PCLDFR ! cloud fraction
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PCIT   ! ice concentration
 REAL, DIMENSION(:,:),     INTENT(IN)     :: PSEA   ! for radar 
 !
@@ -41,11 +38,11 @@ END INTERFACE
 !
 END MODULE MODI_PROFILER_n
 !
-!  ########################################################
-      SUBROUTINE PROFILER_n(PTSTEP,                               &
-                            PXHAT, PYHAT, PZ,PRHODREF,            &
-                            PU, PV, PW, PTH, PR, PSV, PTKE,       &
-                            PTS, PP, PAER, PCLDFR, PCIT, PSEA)
+!     ########################################################
+      SUBROUTINE PROFILER_n( PTSTEP,                         &
+                             PZ, PRHODREF,                   &
+                             PU, PV, PW, PTH, PR, PSV, PTKE, &
+                             PTS, PP, PAER, PCIT, PSEA )
 !     ########################################################
 !
 !
@@ -88,54 +85,50 @@ END MODULE MODI_PROFILER_n
 !  P. Wautelet 09/02/2022: add message when some variables not computed
 !                          + bugfix: put values in variables in this case
 !                          + move some operations outside a do loop
+!  P. Wautelet    04/2022: restructure profilers for better performance, reduce memory usage and correct some problems/bugs
 ! --------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
 !          ------------
 !
-USE MODD_CONF
-USE MODD_CST
+USE MODD_CST,              ONLY: XCPD, XG, XLIGHTSPEED, XP00, XPI, XRD, XRHOLW, XRV, XTT
 USE MODD_DIAG_IN_RUN
-USE MODD_GRID
-USE MODD_SUB_PROFILER_n
-USE MODD_NSV
-USE MODD_PARAMETERS
-USE MODD_PARAM_n,        ONLY: CCLOUD, CRAD
+USE MODD_GRID,             ONLY: XBETA, XLON0, XRPK
+USE MODD_NSV,              ONLY: NSV_C2R2BEG, NSV_C2R2END, NSV_LIMA_NC, NSV_LIMA_NI, NSV_LIMA_NR
+USE MODD_PARAMETERS,       ONLY: JPVEXT, XUNDEF
+USE MODD_PARAM_LIMA,       ONLY: XALPHAR_L => XALPHAR, XNUR_L => XNUR, XALPHAS_L => XALPHAS, XNUS_L => XNUS, &
+                                 XALPHAG_L => XALPHAG, XNUG_L => XNUG, XALPHAI_L => XALPHAI, XNUI_L => XNUI, &
+                                 XRTMIN_L => XRTMIN, XALPHAC_L => XALPHAC, XNUC_L => XNUC
+USE MODD_PARAM_LIMA_COLD,  ONLY: XDI_L => XDI, XLBEXI_L => XLBEXI, XLBI_L => XLBI, XAI_L => XAI, XBI_L => XBI, XC_I_L => XC_I, &
+                                 XLBEXS_L => XLBEXS, XLBS_L => XLBS, XCCS_L => XCCS,                                           &
+                                 XAS_L => XAS, XBS_L => XBS, XCXS_L => XCXS
+USE MODD_PARAM_LIMA_MIXED, ONLY: XDG_L => XDG, XLBEXG_L => XLBEXG, XLBG_L => XLBG, XCCG_L => XCCG, &
+                                 XAG_L => XAG, XBG_L => XBG, XCXG_L => XCXG, XCG_L => XCG
+USE MODD_PARAM_LIMA_WARM,  ONLY: XLBEXR_L => XLBEXR, XLBR_L => XLBR, XBR_L => XBR, XAR_L => XAR, &
+                                 XBC_L => XBC, XAC_L => XAC
+USE MODD_PARAM_n,          ONLY: CCLOUD, CRAD
 USE MODD_PROFILER_n
-USE MODD_TIME,           only: tdtexp
-USE MODD_TIME_n,         only: tdtcur
-!
-USE MODE_ll
+USE MODD_RAIN_ICE_DESCR,   ONLY: XALPHAR_I => XALPHAR, XNUR_I => XNUR, XLBEXR_I => XLBEXR,                                 &
+                                 XLBR_I => XLBR, XCCR_I => XCCR, XBR_I => XBR, XAR_I => XAR,                               &
+                                 XALPHAC_I => XALPHAC, XNUC_I => XNUC,                                                     &
+                                 XLBC_I => XLBC, XBC_I => XBC, XAC_I => XAC,                                               &
+                                 XALPHAC2_I => XALPHAC2, XNUC2_I => XNUC2,                                                 &
+                                 XALPHAS_I => XALPHAS, XNUS_I => XNUS, XLBEXS_I => XLBEXS,                                 &
+                                 XLBS_I => XLBS, XCCS_I => XCCS, XAS_I => XAS, XBS_I => XBS, XCXS_I => XCXS,               &
+                                 XALPHAG_I => XALPHAG, XNUG_I => XNUG, XDG_I => XDG, XLBEXG_I => XLBEXG,                   &
+                                 XLBG_I => XLBG, XCCG_I => XCCG, XAG_I => XAG, XBG_I => XBG, XCXG_I => XCXG, XCG_I => XCG, &
+                                 XALPHAI_I => XALPHAI, XNUI_I => XNUI, XDI_I => XDI, XLBEXI_I => XLBEXI,                   &
+                                 XLBI_I => XLBI, XAI_I => XAI, XBI_I => XBI, XC_I_I => XC_I,                               &
+                                 XRTMIN_I => XRTMIN, XCONC_LAND, XCONC_SEA
+USE MODD_TIME_n,           only: tdtcur
+!
+USE MODE_FGAU,             ONLY: GAULAG
+USE MODE_FSCATTER,         ONLY: BHMIE, QEPSI, QEPSW, MG, MOMG
 USE MODE_MSG
+USE MODE_STATPROF_TOOLS,   ONLY: STATPROF_INTERP_2D, STATPROF_INTERP_3D, STATPROF_INTERP_3D_U, STATPROF_INTERP_3D_V
 !
 USE MODI_GPS_ZENITH_GRID
-USE MODI_LIDAR
-USE MODI_RADAR_RAIN_ICE
 USE MODI_WATER_SUM
-USE MODE_FGAU,             ONLY : GAULAG
-USE MODE_FSCATTER,         ONLY: QEPSW,QEPSI,BHMIE,MOMG,MG
-USE MODD_PARAM_LIMA,       ONLY: XALPHAR_L=>XALPHAR,XNUR_L=>XNUR,XALPHAS_L=>XALPHAS,XNUS_L=>XNUS,&
-                                 XALPHAG_L=>XALPHAG,XNUG_L=>XNUG, XALPHAI_L=>XALPHAI,XNUI_L=>XNUI,&
-                                 XRTMIN_L=>XRTMIN,XALPHAC_L=>XALPHAC,XNUC_L=>XNUC
-USE MODD_PARAM_LIMA_COLD,  ONLY: XDI_L=>XDI,XLBEXI_L=>XLBEXI,XLBI_L=>XLBI,XAI_L=>XAI,XBI_L=>XBI,XC_I_L=>XC_I,&
-                                 XLBEXS_L=>XLBEXS,XLBS_L=>XLBS,XCCS_L=>XCCS,&
-                                 XAS_L=>XAS,XBS_L=>XBS,XCXS_L=>XCXS
-USE MODD_PARAM_LIMA_MIXED, ONLY: XDG_L=>XDG,XLBEXG_L=>XLBEXG,XLBG_L=>XLBG,XCCG_L=>XCCG,&
-                                 XAG_L=>XAG,XBG_L=>XBG,XCXG_L=>XCXG,XCG_L=>XCG
-USE MODD_PARAM_LIMA_WARM,  ONLY: XLBEXR_L=>XLBEXR,XLBR_L=>XLBR,XBR_L=>XBR,XAR_L=>XAR,&
-                                 XBC_L=>XBC,XAC_L=>XAC
-USE MODD_RAIN_ICE_DESCR,   ONLY: XALPHAR_I=>XALPHAR,XNUR_I=>XNUR,XLBEXR_I=>XLBEXR,&
-                                 XLBR_I=>XLBR,XCCR_I=>XCCR,XBR_I=>XBR,XAR_I=>XAR,&
-                                 XALPHAC_I=>XALPHAC,XNUC_I=>XNUC,&
-                                 XLBC_I=>XLBC,XBC_I=>XBC,XAC_I=>XAC,&
-                                 XALPHAC2_I=>XALPHAC2,XNUC2_I=>XNUC2,&
-                                 XALPHAS_I=>XALPHAS,XNUS_I=>XNUS,XLBEXS_I=>XLBEXS,&
-                                 XLBS_I=>XLBS,XCCS_I=>XCCS,XAS_I=>XAS,XBS_I=>XBS,XCXS_I=>XCXS,&
-                                 XALPHAG_I=>XALPHAG,XNUG_I=>XNUG,XDG_I=>XDG,XLBEXG_I=>XLBEXG,&
-                                 XLBG_I=>XLBG,XCCG_I=>XCCG,XAG_I=>XAG,XBG_I=>XBG,XCXG_I=>XCXG,XCG_I=>XCG,&
-                                 XALPHAI_I=>XALPHAI,XNUI_I=>XNUI,XDI_I=>XDI,XLBEXI_I=>XLBEXI,&
-                                 XLBI_I=>XLBI,XAI_I=>XAI,XBI_I=>XBI,XC_I_I=>XC_I,&
-                                 XRTMIN_I=>XRTMIN,XCONC_LAND,XCONC_SEA
 !
 !
 IMPLICIT NONE
@@ -145,51 +138,40 @@ IMPLICIT NONE
 !
 !
 REAL,                     INTENT(IN)     :: PTSTEP ! time step
-REAL, DIMENSION(:),       INTENT(IN)     :: PXHAT  ! x coordinate
-REAL, DIMENSION(:),       INTENT(IN)     :: PYHAT  ! y coordinate
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PZ     ! z array
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PU     ! horizontal wind X component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PV     ! horizontal wind Y component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PW     ! vertical wind
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PTH    ! potential temperature
-REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PRHODREF                            
+REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PRHODREF
 REAL, DIMENSION(:,:,:,:), INTENT(IN)     :: PR     ! water mixing ratios
 REAL, DIMENSION(:,:,:,:), INTENT(IN)     :: PSV    ! Scalar variables
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PTKE   ! turbulent kinetic energy
 REAL, DIMENSION(:,:),     INTENT(IN)     :: PTS    ! surface temperature
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PP     ! pressure
 REAL, DIMENSION(:,:,:,:), INTENT(IN)     :: PAER   ! aerosol extinction
-REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PCLDFR ! cloud fraction
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PCIT   ! ice concentration
-REAL, DIMENSION(:,:),     INTENT(IN)     :: PSEA   ! for radar 
+REAL, DIMENSION(:,:),     INTENT(IN)     :: PSEA   ! for radar
 !
 !-------------------------------------------------------------------------------
 !
 !       0.2  declaration of local variables
 !
 !
-INTEGER :: IIB        ! current processor domain sizes
-INTEGER :: IJB
+INTEGER, PARAMETER :: JPTS_GAULAG = 9 ! number of points for Gauss-Laguerre quadrature
+!
 INTEGER :: IKB
-INTEGER :: IIE
-INTEGER :: IJE
 INTEGER :: IKE
-INTEGER :: IIU
-INTEGER :: IJU
 INTEGER :: IKU
 !
 !
-REAL, DIMENSION(SIZE(PXHAT))        :: ZXHATM ! mass point coordinates
-REAL, DIMENSION(SIZE(PYHAT))        :: ZYHATM ! mass point coordinates
-!
 REAL, DIMENSION(SIZE(PSV,1),SIZE(PSV,2),SIZE(PSV,3),SIZE(PSV,4))  :: ZWORK  
 REAL, DIMENSION(SIZE(PSV,1),SIZE(PSV,2),SIZE(PSV,3),SIZE(PAER,4))  :: ZWORK2  
 !
-LOGICAL :: GSTORE   ! storage occurs at this time step
-!
 INTEGER :: IN     ! time index
 INTEGER :: JSV    ! loop counter
 INTEGER :: JK     ! loop
+INTEGER :: JP     ! loop for profilers
 INTEGER :: IKRAD
 !
 REAL,DIMENSION(SIZE(PZ,3)) :: ZU_PROFILER ! horizontal wind speed profile at station location (along x)
@@ -219,34 +201,25 @@ REAL                       :: ZZWD_PROFILER ! ZWD at station location
 REAL                       :: ZZHDR         ! ZHD correction at station location
 REAL                       :: ZZWDR         ! ZWD correction at station location
 !
-INTEGER                    :: IINFO_ll    ! return code
-INTEGER                    :: ILUOUT      ! logical unit
-INTEGER                    :: IRESP       ! return code
-INTEGER                    :: I           ! loop for stations
-!
 REAL,DIMENSION(SIZE(PTH,1),SIZE(PTH,2))              :: ZZTD,ZZHD,ZZWD
-REAL,DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3))  :: ZTEMP,ZRARE,ZTHV,ZTEMPV
-REAL,DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3))  :: ZWORK32,ZWORK33,ZWORK34
+REAL,DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3))  :: ZTEMP,ZTHV,ZTEMPV
 REAL,DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3))  :: ZVISI,ZVISIKUN
 REAL ::  ZK1,ZK2,ZK3            ! k1, k2 and K3 atmospheric refractivity constants
 REAL  :: ZRDSRV                 ! XRD/XRV
 !
 ! specific to cloud radar
-INTEGER                        :: JLOOP,JLOOP2    ! loop counter
+INTEGER                        :: JLOOP    ! loop counter
 REAL, DIMENSION(SIZE(PR,3))    :: ZTEMPZ! vertical profile of temperature
 REAL, DIMENSION(SIZE(PR,3))    :: ZRHODREFZ ! vertical profile of dry air density of the reference state
 REAL, DIMENSION(SIZE(PR,3))    :: ZCIT     ! pristine ice concentration
 REAL, DIMENSION(SIZE(PR,3))    :: ZCCI,ZCCR,ZCCC     ! ICE,RAIN CLOUD concentration (LIMA)
-REAL, DIMENSION(SIZE(PR,1),SIZE(PR,2),SIZE(PR,3))    :: ZR   
 REAL, DIMENSION(SIZE(PR,3),SIZE(PR,4)+1) :: ZRZ  ! vertical profile of hydrometeor mixing ratios
 REAL                           :: ZA,ZB,ZCC,ZCX,ZALPHA,ZNU,ZLB,ZLBEX,ZRHOHYD,XLAM_CRAD   ! generic microphysical parameters
 INTEGER                        :: JJ    ! loop counter for quadrature
 COMPLEX                        :: QMW,QMI,QM,QB,QEPSIW,QEPSWI   ! dielectric parameter
 REAL                           :: ZAETOT,ZAETMP,ZREFLOC,ZQSCA,ZQBACK,ZQEXT ! temporary scattering parameters
 REAL,DIMENSION(:),ALLOCATABLE  :: ZAELOC,ZZMZ ! temporary arrays
-INTEGER                        :: JPTS_GAULAG=9 ! number of points for Gauss-Laguerre quadrature
 REAL                           :: ZLBDA   ! slope distribution parameter
-REAL                           :: ZFRAC_ICE  ! ice water fraction
 REAL                           :: ZDELTA_EQUIV ! mass-equivalent Gauss-Laguerre point
 REAL                           :: ZFW ! liquid fraction
 REAL                           :: ZFPW ! weight for mixed-phase reflectivity
@@ -270,150 +243,31 @@ XLAM_CRAD        = 3.154E-3 ! (in m) <=> 95.04 GHz = Rasta cloud radar frequency
 !*      2.1  Indices
 !            -------
 !
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
 IKU =   SIZE(PZ,3)     ! nombre de niveaux sur la verticale
 IKB = JPVEXT+1
 IKE = IKU-JPVEXT
 !
-!
-!*      2.2  Interpolations of model variables to mass points
-!            ------------------------------------------------
-!
-IIU=SIZE(PXHAT)
-IJU=SIZE(PYHAT)
-!
-ZXHATM(1:IIU-1)=0.5*PXHAT(1:IIU-1)+0.5*PXHAT(2:IIU  )
-ZXHATM(  IIU  )=1.5*PXHAT(  IIU  )-0.5*PXHAT(  IIU-1)
-!
-ZYHATM(1:IJU-1)=0.5*PYHAT(1:IJU-1)+0.5*PYHAT(2:IJU  )
-ZYHATM(  IJU  )=1.5*PYHAT(  IJU  )-0.5*PYHAT(  IJU-1)
-!
 !----------------------------------------------------------------------------
 !
 !
 !*      3.4  instant of storage
 !            ------------------
 !
-IF ( TPROFILER%T_CUR == XUNDEF ) TPROFILER%T_CUR = TPROFILER%STEP - PTSTEP
+
+IF ( TPROFILERS_TIME%XTIME_CUR == XUNDEF ) TPROFILERS_TIME%XTIME_CUR = TPROFILERS_TIME%XTSTEP - PTSTEP
 !
-TPROFILER%T_CUR = TPROFILER%T_CUR + PTSTEP
+TPROFILERS_TIME%XTIME_CUR = TPROFILERS_TIME%XTIME_CUR + PTSTEP
 !
-IF ( TPROFILER%T_CUR >= TPROFILER%STEP - 1.E-10 ) THEN
-  GSTORE = .TRUE.
-  TPROFILER%T_CUR = TPROFILER%T_CUR - TPROFILER%STEP
-  TPROFILER%N_CUR = TPROFILER%N_CUR + 1
-  IN = TPROFILER%N_CUR
+IF ( TPROFILERS_TIME%XTIME_CUR >= TPROFILERS_TIME%XTSTEP - 1.E-10 ) THEN
+  TPROFILERS_TIME%XTIME_CUR = TPROFILERS_TIME%XTIME_CUR - TPROFILERS_TIME%XTSTEP
+  TPROFILERS_TIME%N_CUR = TPROFILERS_TIME%N_CUR + 1
+  IN = TPROFILERS_TIME%N_CUR
+  tprofilers_time%tpdates(in) = tdtcur
 ELSE
-  GSTORE = .FALSE.
-END IF
-!
-IF (GSTORE) THEN
-#if 0
-  tprofiler%tpdates(in)%date%year  = tdtexp%date%year
-  tprofiler%tpdates(in)%date%month = tdtexp%date%month
-  tprofiler%tpdates(in)%date%day   = tdtexp%date%day
-  tprofiler%tpdates(in)%xtime      = tdtexp%xtime + ( in - 1 ) * tprofiler%step
-#else
-  tprofiler%tpdates(in) = tdtcur
-#endif
+  !No profiler storage at this time step
+  RETURN
 END IF
 !
-!
-!----------------------------------------------------------------------------
-!
-!*      4.   PROFILER POSITION
-!            --------------
-!
-!*      4.0  initialization of processor test
-!            --------------------------------
-IF (GPROFILERFIRSTCALL) THEN
-GPROFILERFIRSTCALL=.FALSE.
-!
- IF (.NOT.(ASSOCIATED(ZTHIS_PROCS))) ALLOCATE(ZTHIS_PROCS(NUMBPROFILER))
-!
-IF (.NOT.(ASSOCIATED(II))) ALLOCATE(II(NUMBPROFILER))
-IF (.NOT.(ASSOCIATED(IJ))) ALLOCATE(IJ(NUMBPROFILER))
-IF (.NOT.(ASSOCIATED(IV))) ALLOCATE(IV(NUMBPROFILER))
-IF (.NOT.(ASSOCIATED(IU))) ALLOCATE(IU(NUMBPROFILER))
-IF (.NOT.(ASSOCIATED(ZXCOEF))) ALLOCATE(ZXCOEF(NUMBPROFILER))
-IF (.NOT.(ASSOCIATED(ZUCOEF))) ALLOCATE(ZUCOEF(NUMBPROFILER))
-IF (.NOT.(ASSOCIATED(ZYCOEF))) ALLOCATE(ZYCOEF(NUMBPROFILER))
-IF (.NOT.(ASSOCIATED(ZVCOEF))) ALLOCATE(ZVCOEF(NUMBPROFILER))
-!
-ZXCOEF(:)=XUNDEF
-ZUCOEF(:)=XUNDEF
-ZYCOEF(:)=XUNDEF
-ZVCOEF(:)=XUNDEF
-!
-DO I=1,NUMBPROFILER
-
-ZTHIS_PROCS(I)=0.
-!
-!*      4.1  X position
-!            ----------
-!
-IU(I)=COUNT( PXHAT (:)<=TPROFILER%X(I) )
-II(I)=COUNT( ZXHATM(:)<=TPROFILER%X(I) )
-!
-IF (II(I)<=IIB-1   .AND. LWEST_ll() .AND. .NOT. L1D) TPROFILER%ERROR(I)=.TRUE.
-IF (II(I)>=IIE     .AND. LEAST_ll() .AND. .NOT. L1D) TPROFILER%ERROR(I)=.TRUE.
-!
-!
-!*      4.2  Y position
-!            ----------
-!
-IV(I)=COUNT( PYHAT (:)<=TPROFILER%Y(I) )
-IJ(I)=COUNT( ZYHATM(:)<=TPROFILER%Y(I) )
-!
-IF (IJ(I)<=IJB-1   .AND. LSOUTH_ll() .AND. .NOT. L1D) TPROFILER%ERROR(I)=.TRUE.
-IF (IJ(I)>=IJE     .AND. LNORTH_ll() .AND. .NOT. L1D) TPROFILER%ERROR(I)=.TRUE.
-!
-!
-!*      4.3  Position of station according to processors
-!            -------------------------------------------
-!
-IF (IU(I)>=IIB .AND. IU(I)<=IIE .AND. IV(I)>=IJB .AND. IV(I)<=IJE) ZTHIS_PROCS(I)=1.
-IF (L1D) ZTHIS_PROCS(I)=1.
-!
-!
-!*      4.4  Computations only on correct processor
-!            --------------------------------------
-ZXCOEF(I) = 0.
-ZYCOEF(I) = 0.
-ZUCOEF(I) = 0.         
-ZVCOEF(I) = 0.
-IF (ZTHIS_PROCS(I) >0. .AND. .NOT. L1D) THEN
-!
-!*      6.1  Interpolation coefficient for X
-!            -------------------------------
-!
-  ZXCOEF(I) = (TPROFILER%X(I) - ZXHATM(II(I))) / (ZXHATM(II(I)+1) - ZXHATM(II(I)))
-!
-!
-!
-!*      6.2  Interpolation coefficient for y
-!            -------------------------------
-!
-  ZYCOEF(I) = (TPROFILER%Y(I) - ZYHATM(IJ(I))) / (ZYHATM(IJ(I)+1) - ZYHATM(IJ(I)))
-!
-!----------------------------------------------------------------------------
-!
-!*      7.   INITIALIZATIONS FOR INTERPOLATIONS OF U AND V
-!            ---------------------------------------------
-!
-!*      7.1  Interpolation coefficient for X (for U)
-!            -------------------------------
-!
-  ZUCOEF(I) = (TPROFILER%X(I) - PXHAT(IU(I))) / (PXHAT(IU(I)+1) - PXHAT(IU(I)))
-!
-!*      7.2  Interpolation coefficient for y (for V)
-!            -------------------------------
-!
-  ZVCOEF(I) = (TPROFILER%Y(I) - PYHAT(IV(I))) / (PYHAT(IV(I)+1) - PYHAT(IV(I)))
-!
-END IF
-ENDDO
-END IF
 !----------------------------------------------------------------------------
 !
 !*      8.   DATA RECORDING
@@ -439,610 +293,406 @@ IF ((SIZE(PR,4) >= 2) .AND. NSV_C2R2END /= 0 ) THEN
   END WHERE
 END IF
 !
-IF (GSTORE) THEN
-  DO I=1,NUMBPROFILER
-    IF ((ZTHIS_PROCS(I)==1.).AND.(.NOT. TPROFILER%ERROR(I))) THEN
-      !
-      ZZ(:)                  = PROFILER_INTERP(PZ)
-      ZRHOD(:)               = PROFILER_INTERP(PRHODREF)
-      ZPRES(:)               = PROFILER_INTERP(PP)
-      ZU_PROFILER(:)         = PROFILER_INTERP_U(PU)
-      ZV_PROFILER(:)         = PROFILER_INTERP_V(PV)
-      ZGAM                   = (XRPK * (TPROFILER%LON(I) - XLON0) - XBETA)*(XPI/180.)
-      ZFF(:)                 = SQRT(ZU_PROFILER(:)**2 + ZV_PROFILER(:)**2)
-      DO JK=1,IKU
-       IF (ZU_PROFILER(JK) >=0. .AND. ZV_PROFILER(JK) > 0.) &
-         ZDD(JK) = ATAN(ABS(ZU_PROFILER(JK)/ZV_PROFILER(JK))) * 180./XPI + 180.
-       IF (ZU_PROFILER(JK) >0. .AND. ZV_PROFILER(JK) <= 0.) &
-         ZDD(JK) = ATAN(ABS(ZV_PROFILER(JK)/ZU_PROFILER(JK))) * 180./XPI + 270.
-       IF (ZU_PROFILER(JK) <=0. .AND. ZV_PROFILER(JK) < 0.) &
-         ZDD(JK) = ATAN(ABS(ZU_PROFILER(JK)/ZV_PROFILER(JK))) * 180./XPI 
-       IF (ZU_PROFILER(JK) <0. .AND. ZV_PROFILER(JK) >= 0.) &
-         ZDD(JK) = ATAN(ABS(ZV_PROFILER(JK)/ZU_PROFILER(JK))) * 180./XPI + 90.
-       IF (ZU_PROFILER(JK) == 0. .AND. ZV_PROFILER(JK) == 0.) &
-         ZDD(JK) = XUNDEF
+PROFILER: DO JP = 1, NUMBPROFILER_LOC
+  ZZ(:)                  = STATPROF_INTERP_3D( TPROFILERS(JP), PZ )
+  ZRHOD(:)               = STATPROF_INTERP_3D( TPROFILERS(JP), PRHODREF )
+  ZPRES(:)               = STATPROF_INTERP_3D( TPROFILERS(JP), PP )
+  ZU_PROFILER(:)         = STATPROF_INTERP_3D_U( TPROFILERS(JP), PU )
+  ZV_PROFILER(:)         = STATPROF_INTERP_3D_V( TPROFILERS(JP), PV )
+  ZGAM                   = (XRPK * (TPROFILERS(JP)%XLON - XLON0) - XBETA)*(XPI/180.)
+  ZFF(:)                 = SQRT(ZU_PROFILER(:)**2 + ZV_PROFILER(:)**2)
+  DO JK=1,IKU
+    IF (ZU_PROFILER(JK) >=0. .AND. ZV_PROFILER(JK) > 0.) &
+      ZDD(JK) = ATAN(ABS(ZU_PROFILER(JK)/ZV_PROFILER(JK))) * 180./XPI + 180.
+    IF (ZU_PROFILER(JK) >0. .AND. ZV_PROFILER(JK) <= 0.) &
+      ZDD(JK) = ATAN(ABS(ZV_PROFILER(JK)/ZU_PROFILER(JK))) * 180./XPI + 270.
+    IF (ZU_PROFILER(JK) <=0. .AND. ZV_PROFILER(JK) < 0.) &
+      ZDD(JK) = ATAN(ABS(ZU_PROFILER(JK)/ZV_PROFILER(JK))) * 180./XPI
+    IF (ZU_PROFILER(JK) <0. .AND. ZV_PROFILER(JK) >= 0.) &
+      ZDD(JK) = ATAN(ABS(ZV_PROFILER(JK)/ZU_PROFILER(JK))) * 180./XPI + 90.
+    IF (ZU_PROFILER(JK) == 0. .AND. ZV_PROFILER(JK) == 0.) &
+      ZDD(JK) = XUNDEF
+  END DO
+  ! GPS IWV and ZTD
+  XZS_GPS=TPROFILERS(JP)%XZ
+  IF ( ABS( ZZ(IKB)-XZS_GPS ) < 150 ) THEN ! distance between real and model orography ok
+    ZRV(:)                 = STATPROF_INTERP_3D( TPROFILERS(JP), PR(:,:,:,1) )
+    ZT(:)                  = STATPROF_INTERP_3D( TPROFILERS(JP), ZTEMP )
+    ZE(:)                  = ZPRES(:)*ZRV(:)/(ZRDSRV+ZRV(:))
+    ZTV(:)                 = STATPROF_INTERP_3D( TPROFILERS(JP), ZTEMPV )
+    ZZTD_PROFILER          = STATPROF_INTERP_2D( TPROFILERS(JP), ZZTD )
+    ZZHD_PROFILER          = STATPROF_INTERP_2D( TPROFILERS(JP), ZZHD )
+    ZZWD_PROFILER          = STATPROF_INTERP_2D( TPROFILERS(JP), ZZWD )
+    ZIWV = 0.
+    DO JK=IKB,IKE
+      ZIWV=ZIWV+ZRHOD(JK)*ZRV(JK)*(ZZ(JK+1)-ZZ(JK))
+    END DO
+    IF (ZZ(IKB) < XZS_GPS) THEN ! station above the model orography
+      DO JK=IKB+1,IKE
+        IF ( ZZ(JK) < XZS_GPS) THEN ! whole layer to remove
+          ZZHDR=( 1.E-6 * ZK1 * ZPRES(JK-1) * ( ZZ(JK) - ZZ(JK-1) ) / ZTV(JK-1))
+          ZZWDR=( 1.E-6 *  ( (ZK2-ZRDSRV*ZK1) + ( ZK3/ZT(JK-1) ) ) * &
+             ZE(JK-1)* ( ZZ(JK) - ZZ(JK-1) ) / ZT(JK-1) )
+          ZZHD_PROFILER=ZZHD_PROFILER-ZZHDR
+          ZZWD_PROFILER=ZZWD_PROFILER-ZZWDR
+          ZZTD_PROFILER=ZZTD_PROFILER-ZZHDR-ZZWDR
+        ELSE                       ! partial layer to remove
+          ZZHDR=( 1.E-6 * ZK1 * ZPRES(JK-1) * ( XZS_GPS - ZZ(JK-1) ) / ZTV(JK-1))
+          ZZWDR=( 1.E-6 *  ( (ZK2-ZRDSRV*ZK1) + ( ZK3/ZT(JK-1) ) ) * &
+             ZE(JK-1)* ( XZS_GPS - ZZ(JK-1) ) / ZT(JK-1) )
+          ZZHD_PROFILER=ZZHD_PROFILER-ZZHDR
+          ZZWD_PROFILER=ZZWD_PROFILER-ZZWDR
+          ZZTD_PROFILER=ZZTD_PROFILER-ZZHDR-ZZWDR
+          EXIT
+        END IF
       END DO
-      ! GPS IWV and ZTD
-      XZS_GPS=TPROFILER%ALT(I)
-      IF ( ABS( ZZ(IKB)-XZS_GPS ) < 150 ) THEN ! distance between real and model orography ok
-        ZRV(:)                 = PROFILER_INTERP(PR(:,:,:,1))
-        ZT(:)                  = PROFILER_INTERP(ZTEMP)
-        ZE(:)                  = ZPRES(:)*ZRV(:)/(ZRDSRV+ZRV(:))
-        ZTV(:)                 = PROFILER_INTERP(ZTEMPV)
-        ZZTD_PROFILER          = PROFILER_INTERP_2D(ZZTD)
-        ZZHD_PROFILER          = PROFILER_INTERP_2D(ZZHD)
-        ZZWD_PROFILER          = PROFILER_INTERP_2D(ZZWD)
-        ZIWV = 0.
-        DO JK=IKB,IKE
-         ZIWV=ZIWV+ZRHOD(JK)*ZRV(JK)*(ZZ(JK+1)-ZZ(JK))
-        END DO
-        IF (ZZ(IKB) < XZS_GPS) THEN ! station above the model orography     
-          DO JK=IKB+1,IKE
-            IF ( ZZ(JK) < XZS_GPS) THEN ! whole layer to remove
-              ZZHDR=( 1.E-6 * ZK1 * ZPRES(JK-1) * ( ZZ(JK) - ZZ(JK-1) ) / ZTV(JK-1))
-              ZZWDR=( 1.E-6 *  ( (ZK2-ZRDSRV*ZK1) + ( ZK3/ZT(JK-1) ) ) * &
-                 ZE(JK-1)* ( ZZ(JK) - ZZ(JK-1) ) / ZT(JK-1) )
-              ZZHD_PROFILER=ZZHD_PROFILER-ZZHDR
-              ZZWD_PROFILER=ZZWD_PROFILER-ZZWDR
-              ZZTD_PROFILER=ZZTD_PROFILER-ZZHDR-ZZWDR
-            ELSE                       ! partial layer to remove
-              ZZHDR=( 1.E-6 * ZK1 * ZPRES(JK-1) * ( XZS_GPS - ZZ(JK-1) ) / ZTV(JK-1)) 
-              ZZWDR=( 1.E-6 *  ( (ZK2-ZRDSRV*ZK1) + ( ZK3/ZT(JK-1) ) ) * &
-                 ZE(JK-1)* ( XZS_GPS - ZZ(JK-1) ) / ZT(JK-1) ) 
-              ZZHD_PROFILER=ZZHD_PROFILER-ZZHDR
-              ZZWD_PROFILER=ZZWD_PROFILER-ZZWDR
-              ZZTD_PROFILER=ZZTD_PROFILER-ZZHDR-ZZWDR
-              EXIT
-            END IF
-          END DO 
-        ELSE ! station below the model orography
+    ELSE ! station below the model orography
 ! Extrapolate variables below the model orography assuming constant T&Tv gradients,
 ! constant rv and hydrostatic law
-          ZZHATM(:)=0.5*(ZZ(1:IKU-1)+ZZ(2:IKU))
-          ZZM_STAT=0.5*(XZS_GPS+ZZ(IKB))
-          ZTM_STAT=ZT(IKB) + ( (ZZM_STAT-ZZHATM(IKB))*&
-             ( ZT(IKB)- ZT(IKB+1) )/(ZZHATM(IKB)-ZZHATM(IKB+1)) )
-          ZTV_STAT=ZTV(IKB) + ( (ZZM_STAT-ZZHATM(IKB))*&
-             ( ZTV(IKB)- ZTV(IKB+1) )/(ZZHATM(IKB)-ZZHATM(IKB+1)) )
-          ZPM_STAT = ZPRES(IKB) * EXP(XG *(ZZM_STAT-ZZHATM(IKB))&
-             /(XRD* 0.5 *(ZTV_STAT+ZTV(IKB))))
-          ZEM_STAT = ZPM_STAT * ZRV(IKB) / ( ZRDSRV + ZRV(IKB) )
-! add contribution below the model orography        
-          ZZHDR=( 1.E-6 * ZK1 * ZPM_STAT * ( ZZ(IKB) - XZS_GPS ) / ZTV_STAT )
-          ZZWDR=( 1.E-6 * ( (ZK2-ZRDSRV*ZK1) + (ZK3/ZTM_STAT) )&
-             * ZEM_STAT* ( ZZ(IKB) - XZS_GPS ) / ZTM_STAT )
-          ZZHD_PROFILER=ZZHD_PROFILER+ZZHDR
-          ZZWD_PROFILER=ZZWD_PROFILER+ZZWDR
-          ZZTD_PROFILER=ZZTD_PROFILER+ZZHDR+ZZWDR
-        END IF
-        TPROFILER%IWV(IN,I)= ZIWV
-        TPROFILER%ZTD(IN,I)= ZZTD_PROFILER
-        TPROFILER%ZWD(IN,I)= ZZWD_PROFILER
-        TPROFILER%ZHD(IN,I)= ZZHD_PROFILER
-      ELSE
-        CMNHMSG(1) = 'altitude of profiler ' // TRIM( TPROFILER%NAME(I) ) // ' is too far from orography'
-        CMNHMSG(2) = 'some variables are therefore not computed (IWV, ZTD, ZWD, ZHD)'
-        CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'PROFILER_n' )
-        TPROFILER%IWV(IN,I)= XUNDEF
-        TPROFILER%ZTD(IN,I)= XUNDEF
-        TPROFILER%ZWD(IN,I)= XUNDEF
-        TPROFILER%ZHD(IN,I)= XUNDEF
-      END IF
-      TPROFILER%ZON (IN,:,I) = ZU_PROFILER(:) * COS(ZGAM) + ZV_PROFILER(:) * SIN(ZGAM)
-      TPROFILER%MER (IN,:,I) = - ZU_PROFILER(:) * SIN(ZGAM) + ZV_PROFILER(:) * COS(ZGAM)
-      TPROFILER%FF  (IN,:,I) = ZFF(:)                                                     
-      TPROFILER%DD  (IN,:,I) = ZDD(:)   
-      TPROFILER%W   (IN,:,I) = PROFILER_INTERP(PW)
-      TPROFILER%TH  (IN,:,I) = PROFILER_INTERP(PTH)
-      TPROFILER%THV (IN,:,I) = PROFILER_INTERP(ZTHV)
-      TPROFILER%VISI(IN,:,I) = PROFILER_INTERP(ZVISI)
-      TPROFILER%VISIKUN(IN,:,I) = PROFILER_INTERP(ZVISIKUN)
-      TPROFILER%ZZ  (IN,:,I) = ZZ(:)
-      TPROFILER%RHOD(IN,:,I) = ZRHOD(:)
-      TPROFILER%CIZ(IN,:,I) = PROFILER_INTERP(PCIT)
+      ZZHATM(:)=0.5*(ZZ(1:IKU-1)+ZZ(2:IKU))
+      ZZM_STAT=0.5*(XZS_GPS+ZZ(IKB))
+      ZTM_STAT=ZT(IKB) + ( (ZZM_STAT-ZZHATM(IKB))*&
+         ( ZT(IKB)- ZT(IKB+1) )/(ZZHATM(IKB)-ZZHATM(IKB+1)) )
+      ZTV_STAT=ZTV(IKB) + ( (ZZM_STAT-ZZHATM(IKB))*&
+         ( ZTV(IKB)- ZTV(IKB+1) )/(ZZHATM(IKB)-ZZHATM(IKB+1)) )
+      ZPM_STAT = ZPRES(IKB) * EXP(XG *(ZZM_STAT-ZZHATM(IKB))&
+         /(XRD* 0.5 *(ZTV_STAT+ZTV(IKB))))
+      ZEM_STAT = ZPM_STAT * ZRV(IKB) / ( ZRDSRV + ZRV(IKB) )
+! add contribution below the model orography
+      ZZHDR=( 1.E-6 * ZK1 * ZPM_STAT * ( ZZ(IKB) - XZS_GPS ) / ZTV_STAT )
+      ZZWDR=( 1.E-6 * ( (ZK2-ZRDSRV*ZK1) + (ZK3/ZTM_STAT) )&
+         * ZEM_STAT* ( ZZ(IKB) - XZS_GPS ) / ZTM_STAT )
+      ZZHD_PROFILER=ZZHD_PROFILER+ZZHDR
+      ZZWD_PROFILER=ZZWD_PROFILER+ZZWDR
+      ZZTD_PROFILER=ZZTD_PROFILER+ZZHDR+ZZWDR
+    END IF
+    TPROFILERS(JP)%XIWV(IN)= ZIWV
+    TPROFILERS(JP)%XZTD(IN)= ZZTD_PROFILER
+    TPROFILERS(JP)%XZWD(IN)= ZZWD_PROFILER
+    TPROFILERS(JP)%XZHD(IN)= ZZHD_PROFILER
+  ELSE
+    CMNHMSG(1) = 'altitude of profiler ' // TRIM( TPROFILERS(JP)%CNAME ) // ' is too far from orography'
+    CMNHMSG(2) = 'some variables are therefore not computed (IWV, ZTD, ZWD, ZHD)'
+    CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'PROFILER_n' )
+    TPROFILERS(JP)%XIWV(IN)= XUNDEF
+    TPROFILERS(JP)%XZTD(IN)= XUNDEF
+    TPROFILERS(JP)%XZWD(IN)= XUNDEF
+    TPROFILERS(JP)%XZHD(IN)= XUNDEF
+  END IF
+  TPROFILERS(JP)%XZON (IN,:) = ZU_PROFILER(:) * COS(ZGAM) + ZV_PROFILER(:) * SIN(ZGAM)
+  TPROFILERS(JP)%XMER (IN,:) = - ZU_PROFILER(:) * SIN(ZGAM) + ZV_PROFILER(:) * COS(ZGAM)
+  TPROFILERS(JP)%XFF  (IN,:) = ZFF(:)
+  TPROFILERS(JP)%XDD  (IN,:) = ZDD(:)
+  TPROFILERS(JP)%XW   (IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), PW )
+  TPROFILERS(JP)%XTH  (IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), PTH )
+  TPROFILERS(JP)%XTHV (IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), ZTHV )
+  TPROFILERS(JP)%XVISI(IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), ZVISI )
+  TPROFILERS(JP)%XVISIKUN(IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), ZVISIKUN )
+  TPROFILERS(JP)%XZZ  (IN,:) = ZZ(:)
+  TPROFILERS(JP)%XRHOD(IN,:) = ZRHOD(:)
+  IF ( CCLOUD == 'ICE3' .OR. CCLOUD == 'ICE4' ) &
+    TPROFILERS(JP)%XCIZ(IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), PCIT )
 ! add RARE
-        ! initialization CRARE and CRARE_ATT + LWC and IWC
-      TPROFILER%CRARE(IN,:,I) = 0.
-      TPROFILER%CRARE_ATT(IN,:,I) = 0.
-      TPROFILER%LWCZ  (IN,:,I) = 0.
-      TPROFILER%IWCZ  (IN,:,I) = 0.
-      IF (CCLOUD=="LIMA" .OR. CCLOUD=="ICE3" ) THEN ! only for ICE3 and LIMA
-       TPROFILER%LWCZ  (IN,:,I) = PROFILER_INTERP((PR(:,:,:,2)+PR(:,:,:,3))*PRHODREF(:,:,:))
-       TPROFILER%IWCZ  (IN,:,I) = PROFILER_INTERP((PR(:,:,:,4)+PR(:,:,:,5)+PR(:,:,:,6))*PRHODREF(:,:,:))
-       ZTEMPZ(:)=PROFILER_INTERP(ZTEMP(:,:,:))
-       ZRHODREFZ(:)=PROFILER_INTERP(PRHODREF(:,:,:))
-       ZCIT(:)=PROFILER_INTERP(PCIT(:,:,:))
-       IF (CCLOUD=="LIMA") THEN
-          ZCCI(:)=PROFILER_INTERP(PSV(:,:,:,NSV_LIMA_NI))
-          ZCCR(:)=PROFILER_INTERP(PSV(:,:,:,NSV_LIMA_NR))
-          ZCCC(:)=PROFILER_INTERP(PSV(:,:,:,NSV_LIMA_NC))
-       ENDIF
-       DO JLOOP=3,6
-          ZRZ(:,JLOOP)=PROFILER_INTERP(PR(:,:,:,JLOOP))
-       END DO
-       DO JK=1,IKU
-          ZRZ(JK,2)=PROFILER_INTERP_2D(PR(:,:,JK,2)*PSEA(:,:))       ! becomes cloud mixing ratio over sea
-          ZRZ(JK,7)=PROFILER_INTERP_2D(PR(:,:,JK,2)*(1.-PSEA(:,:)))  ! becomes cloud mixing ratio over land
-       END DO
-       ALLOCATE(ZAELOC(IKU))
-        !
-       ZAELOC(:)=0.
-       ! initialization of quadrature points and weights
-       ALLOCATE(ZX(JPTS_GAULAG),ZW(JPTS_GAULAG))
-       CALL GAULAG(JPTS_GAULAG,ZX,ZW) ! for integration over diameters
-        ! initialize minimum values
-       ALLOCATE(ZRTMIN(SIZE(PR,4)+1))
-       IF (CCLOUD == 'LIMA') THEN
-          ZRTMIN(2)=XRTMIN_L(2) ! cloud water over sea
-          ZRTMIN(3)=XRTMIN_L(3)
-          ZRTMIN(4)=XRTMIN_L(4)
-          ZRTMIN(5)=1E-10
-          ZRTMIN(6)=XRTMIN_L(6)
-          ZRTMIN(7)=XRTMIN_L(2) ! cloud water over land
-       ELSE
-          ZRTMIN(2)=XRTMIN_I(2) ! cloud water over sea
-          ZRTMIN(3)=XRTMIN_I(3)
-          ZRTMIN(4)=XRTMIN_I(4)
-          ZRTMIN(5)=1E-10
-          ZRTMIN(6)=XRTMIN_I(6)
-          ZRTMIN(7)=XRTMIN_I(2) ! cloud water over land
-       ENDIF
-        ! compute cloud radar reflectivity from vertical profiles of temperature
-        ! and mixing ratios
-       DO JK=1,IKU
-         QMW=SQRT(QEPSW(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
-         QMI=SQRT(QEPSI(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
-         DO JLOOP=2,7
-           IF (CCLOUD == 'LIMA') THEN
-              GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCCI(JK)>0.).AND.&
-                    (JLOOP.NE.3.OR.ZCCR(JK)>0.).AND.((JLOOP.NE.2.AND.JLOOP.NE.7).OR.ZCCC(JK)>0.))
-           ELSE
-              GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCIT(JK)>0.))
-           ENDIF
-           IF(GCALC) THEN
-              SELECT CASE(JLOOP)
-                CASE(2) ! cloud water over sea
-                  IF (CCLOUD == 'LIMA') THEN
-                    ZA=XAC_L
-                    ZB=XBC_L
-                    ZCC=ZCCC(JK)*ZRHODREFZ(JK)
-                    ZCX=0.
-                    ZALPHA=XALPHAC_L
-                    ZNU=XNUC_L
-                    ZLBEX=1.0/(ZCX-ZB)
-                    ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
-                  ELSE
-                    ZA=XAC_I
-                    ZB=XBC_I
-                    ZCC=XCONC_SEA
-                    ZCX=0.
-                    ZALPHA=XALPHAC2_I
-                    ZNU=XNUC2_I
-                    ZLBEX=1.0/(ZCX-ZB)
-                    ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
-                  ENDIF
-                CASE(3) ! rain water
-                  IF (CCLOUD == 'LIMA') THEN
-                    ZA=XAR_L
-                    ZB=XBR_L
-                    ZCC=ZCCR(JK)*ZRHODREFZ(JK)
-                    ZCX=0.
-                    ZALPHA=XALPHAR_L
-                    ZNU=XNUR_L
-                    ZLBEX=1.0/(ZCX-ZB)
-                    ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
-                  ELSE
-                    ZA=XAR_I
-                    ZB=XBR_I
-                    ZCC=XCCR_I
-                    ZCX=-1.
-                    ZALPHA=XALPHAR_I
-                    ZNU=XNUR_I
-                    ZLB=XLBR_I
-                    ZLBEX=XLBEXR_I
-                  ENDIF
-                CASE(4) ! pristine ice
-                  IF (CCLOUD == 'LIMA') THEN
-                    ZA=XAI_L
-                    ZB=XBI_L
-                    ZCC=ZCCI(JK)*ZRHODREFZ(JK)
-                    ZCX=0.
-                    ZALPHA=XALPHAI_L
-                    ZNU=XNUI_L
-                    ZLBEX=1.0/(ZCX-ZB)
-                    ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX) ! because ZCC not included in XLBI
-                    ZFW=0
-                  ELSE
-                    ZA=XAI_I
-                    ZB=XBI_I
-                    ZCC=ZCIT(JK)
-                    ZCX=0.
-                    ZALPHA=XALPHAI_I
-                    ZNU=XNUI_I
-                    ZLBEX=XLBEXI_I
-                    ZLB=XLBI_I*ZCC**(-ZLBEX) ! because ZCC not included in XLBI
-                    ZFW=0
-                  ENDIF
-                CASE(5) ! snow
-                  IF (CCLOUD == 'LIMA') THEN
-                    ZA=XAS_L
-                    ZB=XBS_L
-                    ZCC=XCCS_L
-                    ZCX=XCXS_L
-                    ZALPHA=XALPHAS_L
-                    ZNU=XNUS_L
-                    ZLB=XLBS_L
-                    ZLBEX=XLBEXS_L
-                    ZFW=0
-                  ELSE
-                    ZA=XAS_I
-                    ZB=XBS_I
-                    ZCC=XCCS_I
-                    ZCX=XCXS_I
-                    ZALPHA=XALPHAS_I
-                    ZNU=XNUS_I
-                    ZLB=XLBS_I
-                    ZLBEX=XLBEXS_I
-                    ZFW=0
-                  ENDIF
-                CASE(6) ! graupel
-                  !If temperature between -10 and 10B0C and Mr and Mg over min
-                  !threshold: melting graupel
-                  ! with liquid water fraction Fw=Mr/(Mr+Mg) else dry graupel
-                  ! (Fw=0)    
-                  IF( ZTEMPZ(JK) > XTT-10 .AND. ZTEMPZ(JK) < XTT+10 &
-                    .AND. ZRZ(JK,3) > ZRTMIN(3) ) THEN
-                    ZFW=ZRZ(JK,3)/(ZRZ(JK,3)+ZRZ(JK,JLOOP))
-                  ELSE
-                    ZFW=0
-                  ENDIF
-                  IF (CCLOUD == 'LIMA') THEN
-                    ZA=XAG_L
-                    ZB=XBG_L
-                    ZCC=XCCG_L
-                    ZCX=XCXG_L
-                    ZALPHA=XALPHAG_L
-                    ZNU=XNUG_L
-                    ZLB=XLBG_L
-                    ZLBEX=XLBEXG_L
-                  ELSE
-                    ZA=XAG_I
-                    ZB=XBG_I
-                    ZCC=XCCG_I
-                    ZCX=XCXG_I
-                    ZALPHA=XALPHAG_I
-                    ZNU=XNUG_I
-                    ZLB=XLBG_I
-                    ZLBEX=XLBEXG_I
-                  ENDIF
-                CASE(7) ! cloud water over land
-                  IF (CCLOUD == 'LIMA') THEN
-                    ZA=XAC_L
-                    ZB=XBC_L
-                    ZCC=ZCCC(JK)*ZRHODREFZ(JK)
-                    ZCX=0.
-                    ZALPHA=XALPHAC_L
-                    ZNU=XNUC_L
-                    ZLBEX=1.0/(ZCX-ZB)
-                    ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
-                  ELSE
-                    ZA=XAC_I
-                    ZB=XBC_I
-                    ZCC=XCONC_LAND
-                    ZCX=0.
-                    ZALPHA=XALPHAC_I
-                    ZNU=XNUC_I
-                    ZLBEX=1.0/(ZCX-ZB)
-                    ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
-                  ENDIF
+    ! initialization CRARE and CRARE_ATT + LWC and IWC
+  TPROFILERS(JP)%XCRARE(IN,:) = 0.
+  TPROFILERS(JP)%XCRARE_ATT(IN,:) = 0.
+  TPROFILERS(JP)%XLWCZ  (IN,:) = 0.
+  TPROFILERS(JP)%XIWCZ  (IN,:) = 0.
+  IF (CCLOUD=="LIMA" .OR. CCLOUD=="ICE3" ) THEN ! only for ICE3 and LIMA
+    TPROFILERS(JP)%XLWCZ  (IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), (PR(:,:,:,2)+PR(:,:,:,3))*PRHODREF(:,:,:) )
+    TPROFILERS(JP)%XIWCZ  (IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), (PR(:,:,:,4)+PR(:,:,:,5)+PR(:,:,:,6))*PRHODREF(:,:,:) )
+    ZTEMPZ(:)=STATPROF_INTERP_3D( TPROFILERS(JP), ZTEMP(:,:,:) )
+    ZRHODREFZ(:)=STATPROF_INTERP_3D( TPROFILERS(JP), PRHODREF(:,:,:) )
+    ZCIT(:)=STATPROF_INTERP_3D( TPROFILERS(JP), PCIT(:,:,:) )
+    IF (CCLOUD=="LIMA") THEN
+      ZCCI(:)=STATPROF_INTERP_3D( TPROFILERS(JP), PSV(:,:,:,NSV_LIMA_NI) )
+      ZCCR(:)=STATPROF_INTERP_3D( TPROFILERS(JP), PSV(:,:,:,NSV_LIMA_NR) )
+      ZCCC(:)=STATPROF_INTERP_3D( TPROFILERS(JP), PSV(:,:,:,NSV_LIMA_NC) )
+    END IF
+    DO JLOOP=3,6
+      ZRZ(:,JLOOP)=STATPROF_INTERP_3D( TPROFILERS(JP), PR(:,:,:,JLOOP) )
+    END DO
+    DO JK=1,IKU
+      ZRZ(JK,2)=STATPROF_INTERP_2D( TPROFILERS(JP), PR(:,:,JK,2)*PSEA(:,:) )       ! becomes cloud mixing ratio over sea
+      ZRZ(JK,7)=STATPROF_INTERP_2D( TPROFILERS(JP), PR(:,:,JK,2)*(1.-PSEA(:,:)) )  ! becomes cloud mixing ratio over land
+    END DO
+    ALLOCATE(ZAELOC(IKU))
+    !
+    ZAELOC(:)=0.
+    ! initialization of quadrature points and weights
+    ALLOCATE(ZX(JPTS_GAULAG),ZW(JPTS_GAULAG))
+    CALL GAULAG(JPTS_GAULAG,ZX,ZW) ! for integration over diameters
+    ! initialize minimum values
+    ALLOCATE(ZRTMIN(SIZE(PR,4)+1))
+    IF (CCLOUD == 'LIMA') THEN
+      ZRTMIN(2)=XRTMIN_L(2) ! cloud water over sea
+      ZRTMIN(3)=XRTMIN_L(3)
+      ZRTMIN(4)=XRTMIN_L(4)
+      ZRTMIN(5)=1E-10
+      ZRTMIN(6)=XRTMIN_L(6)
+      ZRTMIN(7)=XRTMIN_L(2) ! cloud water over land
+    ELSE
+      ZRTMIN(2)=XRTMIN_I(2) ! cloud water over sea
+      ZRTMIN(3)=XRTMIN_I(3)
+      ZRTMIN(4)=XRTMIN_I(4)
+      ZRTMIN(5)=1E-10
+      ZRTMIN(6)=XRTMIN_I(6)
+      ZRTMIN(7)=XRTMIN_I(2) ! cloud water over land
+    END IF
+    ! compute cloud radar reflectivity from vertical profiles of temperature
+    ! and mixing ratios
+    DO JK=1,IKU
+      QMW=SQRT(QEPSW(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
+      QMI=SQRT(QEPSI(ZTEMPZ(JK),XLIGHTSPEED/XLAM_CRAD))
+      DO JLOOP=2,7
+        IF (CCLOUD == 'LIMA') THEN
+          GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCCI(JK)>0.).AND.&
+                (JLOOP.NE.3.OR.ZCCR(JK)>0.).AND.((JLOOP.NE.2.AND.JLOOP.NE.7).OR.ZCCC(JK)>0.))
+        ELSE
+          GCALC=(ZRZ(JK,JLOOP)>ZRTMIN(JLOOP).AND.(JLOOP.NE.4.OR.ZCIT(JK)>0.))
+        END IF
+        IF (GCALC) THEN
+          SELECT CASE(JLOOP)
+            CASE(2) ! cloud water over sea
+              IF (CCLOUD == 'LIMA') THEN
+                ZA=XAC_L
+                ZB=XBC_L
+                ZCC=ZCCC(JK)*ZRHODREFZ(JK)
+                ZCX=0.
+                ZALPHA=XALPHAC_L
+                ZNU=XNUC_L
+                ZLBEX=1.0/(ZCX-ZB)
+                ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+              ELSE
+                ZA=XAC_I
+                ZB=XBC_I
+                ZCC=XCONC_SEA
+                ZCX=0.
+                ZALPHA=XALPHAC2_I
+                ZNU=XNUC2_I
+                ZLBEX=1.0/(ZCX-ZB)
+                ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+              END IF
+            CASE(3) ! rain water
+              IF (CCLOUD == 'LIMA') THEN
+                ZA=XAR_L
+                ZB=XBR_L
+                ZCC=ZCCR(JK)*ZRHODREFZ(JK)
+                ZCX=0.
+                ZALPHA=XALPHAR_L
+                ZNU=XNUR_L
+                ZLBEX=1.0/(ZCX-ZB)
+                ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+              ELSE
+                ZA=XAR_I
+                ZB=XBR_I
+                ZCC=XCCR_I
+                ZCX=-1.
+                ZALPHA=XALPHAR_I
+                ZNU=XNUR_I
+                ZLB=XLBR_I
+                ZLBEX=XLBEXR_I
+              END IF
+            CASE(4) ! pristine ice
+              IF (CCLOUD == 'LIMA') THEN
+                ZA=XAI_L
+                ZB=XBI_L
+                ZCC=ZCCI(JK)*ZRHODREFZ(JK)
+                ZCX=0.
+                ZALPHA=XALPHAI_L
+                ZNU=XNUI_L
+                ZLBEX=1.0/(ZCX-ZB)
+                ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX) ! because ZCC not included in XLBI
+                ZFW=0
+              ELSE
+                ZA=XAI_I
+                ZB=XBI_I
+                ZCC=ZCIT(JK)
+                ZCX=0.
+                ZALPHA=XALPHAI_I
+                ZNU=XNUI_I
+                ZLBEX=XLBEXI_I
+                ZLB=XLBI_I*ZCC**(-ZLBEX) ! because ZCC not included in XLBI
+                ZFW=0
+              END IF
+            CASE(5) ! snow
+              IF (CCLOUD == 'LIMA') THEN
+                ZA=XAS_L
+                ZB=XBS_L
+                ZCC=XCCS_L
+                ZCX=XCXS_L
+                ZALPHA=XALPHAS_L
+                ZNU=XNUS_L
+                ZLB=XLBS_L
+                ZLBEX=XLBEXS_L
+                ZFW=0
+              ELSE
+                ZA=XAS_I
+                ZB=XBS_I
+                ZCC=XCCS_I
+                ZCX=XCXS_I
+                ZALPHA=XALPHAS_I
+                ZNU=XNUS_I
+                ZLB=XLBS_I
+                ZLBEX=XLBEXS_I
+                ZFW=0
+              END IF
+            CASE(6) ! graupel
+              !If temperature between -10 and 10B0C and Mr and Mg over min
+              !threshold: melting graupel
+              ! with liquid water fraction Fw=Mr/(Mr+Mg) else dry graupel
+              ! (Fw=0)
+              IF( ZTEMPZ(JK) > XTT-10 .AND. ZTEMPZ(JK) < XTT+10 &
+                .AND. ZRZ(JK,3) > ZRTMIN(3) ) THEN
+                ZFW=ZRZ(JK,3)/(ZRZ(JK,3)+ZRZ(JK,JLOOP))
+              ELSE
+                ZFW=0
+              END IF
+              IF (CCLOUD == 'LIMA') THEN
+                ZA=XAG_L
+                ZB=XBG_L
+                ZCC=XCCG_L
+                ZCX=XCXG_L
+                ZALPHA=XALPHAG_L
+                ZNU=XNUG_L
+                ZLB=XLBG_L
+                ZLBEX=XLBEXG_L
+              ELSE
+                ZA=XAG_I
+                ZB=XBG_I
+                ZCC=XCCG_I
+                ZCX=XCXG_I
+                ZALPHA=XALPHAG_I
+                ZNU=XNUG_I
+                ZLB=XLBG_I
+                ZLBEX=XLBEXG_I
+              END IF
+            CASE(7) ! cloud water over land
+              IF (CCLOUD == 'LIMA') THEN
+                ZA=XAC_L
+                ZB=XBC_L
+                ZCC=ZCCC(JK)*ZRHODREFZ(JK)
+                ZCX=0.
+                ZALPHA=XALPHAC_L
+                ZNU=XNUC_L
+                ZLBEX=1.0/(ZCX-ZB)
+                ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+              ELSE
+                ZA=XAC_I
+                ZB=XBC_I
+                ZCC=XCONC_LAND
+                ZCX=0.
+                ZALPHA=XALPHAC_I
+                ZNU=XNUC_I
+                ZLBEX=1.0/(ZCX-ZB)
+                ZLB=( ZA*ZCC*MOMG(ZALPHA,ZNU,ZB) )**(-ZLBEX)
+              END IF
+          END SELECT
+          ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
+          ZREFLOC=0.
+          ZAETMP=0.
+          DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
+            ZDELTA_EQUIV=ZX(JJ)**(1./ZALPHA)/ZLBDA
+            SELECT CASE(JLOOP)
+              CASE(2,3,7)
+                QM=QMW
+              CASE(4,5,6)
+                ! pristine ice, snow, dry graupel
+                ZRHOHYD=MIN(6.*ZA*ZDELTA_EQUIV**(ZB-3.)/XPI,.92*XRHOLW)
+                QM=sqrt(MG(QMI**2,CMPLX(1,0),ZRHOHYD/.92/XRHOLW))
+                ! water inclusions in ice in air
+                QEPSWI=MG(QMW**2,QM**2,ZFW)
+                ! ice in air inclusions in water
+                QEPSIW=MG(QM**2,QMW**2,1.-ZFW)
+                !MG weighted rule (Matrosov 2008)
+                IF(ZFW .LT. 0.37) THEN
+                  ZFPW=0
+                ELSE IF(ZFW .GT. 0.63) THEN
+                  ZFPW=1
+                ELSE
+                  ZFPW=(ZFW-0.37)/(0.63-0.37)
+                ENDIF
+                QM=sqrt(QEPSWI*(1.-ZFPW)+QEPSIW*ZFPW)
             END SELECT
-            ZLBDA=ZLB*(ZRHODREFZ(JK)*ZRZ(JK,JLOOP))**ZLBEX
-            ZREFLOC=0.
-            ZAETMP=0.
-            DO JJ=1,JPTS_GAULAG ! Gauss-Laguerre quadrature
-               ZDELTA_EQUIV=ZX(JJ)**(1./ZALPHA)/ZLBDA
-               SELECT CASE(JLOOP)
-                  CASE(2,3,7)
-                    QM=QMW
-                  CASE(4,5,6)
-                    ! pristine ice, snow, dry graupel
-                    ZRHOHYD=MIN(6.*ZA*ZDELTA_EQUIV**(ZB-3.)/XPI,.92*XRHOLW)
-                    QM=sqrt(MG(QMI**2,CMPLX(1,0),ZRHOHYD/.92/XRHOLW))
-                    ! water inclusions in ice in air
-                    QEPSWI=MG(QMW**2,QM**2,ZFW)
-                    ! ice in air inclusions in water
-                    QEPSIW=MG(QM**2,QMW**2,1.-ZFW)
-                    !MG weighted rule (Matrosov 2008)
-                    IF(ZFW .LT. 0.37) THEN
-                      ZFPW=0
-                    ELSE IF(ZFW .GT. 0.63) THEN
-                      ZFPW=1
-                    ELSE
-                      ZFPW=(ZFW-0.37)/(0.63-0.37)
-                    ENDIF
-                    QM=sqrt(QEPSWI*(1.-ZFPW)+QEPSIW*ZFPW)
-               END SELECT
-               CALL BHMIE(XPI/XLAM_CRAD*ZDELTA_EQUIV,QM,ZQEXT,ZQSCA,ZQBACK)
-               ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
-               ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
-             END DO
-             ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
-             ZAETMP=ZAETMP  *           XPI    *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
-             TPROFILER%CRARE(IN,JK,I)=TPROFILER%CRARE(IN,JK,I)+ZREFLOC
-             ZAELOC(JK)=ZAELOC(JK)+ZAETMP
-           END IF
-         END DO
-       END DO
-     ! apply attenuation
-       ALLOCATE(ZZMZ(IKU))
-       ZZMZ = ZZ(:) ! PROFILER_INTERP(ZZM(:,:,:))
+            CALL BHMIE(XPI/XLAM_CRAD*ZDELTA_EQUIV,QM,ZQEXT,ZQSCA,ZQBACK)
+            ZREFLOC=ZREFLOC+ZQBACK*ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
+            ZAETMP =ZAETMP +ZQEXT *ZX(JJ)**(ZNU-1)*ZDELTA_EQUIV**2*ZW(JJ)
+          END DO
+          ZREFLOC=ZREFLOC*(XLAM_CRAD/XPI)**4*ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU)*.93)
+          ZAETMP=ZAETMP  *           XPI    *ZCC*ZLBDA**ZCX/(4.*GAMMA(ZNU))
+          TPROFILERS(JP)%XCRARE(IN,JK)=TPROFILERS(JP)%XCRARE(IN,JK)+ZREFLOC
+          ZAELOC(JK)=ZAELOC(JK)+ZAETMP
+        END IF
+      END DO
+    END DO
+ ! apply attenuation
+    ALLOCATE(ZZMZ(IKU))
+    ZZMZ = ZZ(:) ! STATPROF_INTERP_3D( TPROFILERS(JP), ZZM(:,:,:) )
 !        ZZMZ(1)=ZZM_STAT
-        ! zenith
-       ZAETOT=1.
-       DO JK = 2,IKU
-            ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
-        ZAETOT=ZAETOT*EXP(-(ZAELOC(JK-1)+ZAELOC(JK))*(ZZMZ(JK)-ZZMZ(JK-1)))
-        TPROFILER%CRARE_ATT(IN,JK,I)=TPROFILER%CRARE(IN,JK,I)*ZAETOT
-       END DO
-!        TPROFILER%ZZ  (IN,:,I) = ZZMZ(:)
-       DEALLOCATE(ZZMZ,ZAELOC)
-        ! m^3 bmm^6/m^3 bdBZ
-       WHERE(TPROFILER%CRARE(IN,:,I)>0)
-          TPROFILER%CRARE(IN,:,I)=10.*LOG10(1.E18*TPROFILER%CRARE(IN,:,I))
-       ELSEWHERE
-          TPROFILER%CRARE(IN,:,I)=XUNDEF
-       END WHERE
-       WHERE(TPROFILER%CRARE_ATT(IN,:,I)>0)
-          TPROFILER%CRARE_ATT(IN,:,I)=10.*LOG10(1.E18*TPROFILER%CRARE_ATT(IN,:,I))
-       ELSEWHERE
-          TPROFILER%CRARE_ATT(IN,:,I)=XUNDEF
-       END WHERE
-       DEALLOCATE(ZX,ZW,ZRTMIN)
-     END IF ! end LOOP ICE3
+    ! zenith
+    ZAETOT=1.
+    DO JK = 2,IKU
+         ! attenuation from ZAELOC(JK) and ZAELOC(JK-1)
+      ZAETOT=ZAETOT*EXP(-(ZAELOC(JK-1)+ZAELOC(JK))*(ZZMZ(JK)-ZZMZ(JK-1)))
+      TPROFILERS(JP)%XCRARE_ATT(IN,JK)=TPROFILERS(JP)%XCRARE(IN,JK)*ZAETOT
+    END DO
+    DEALLOCATE(ZZMZ,ZAELOC)
+     ! m^3 b mm^6/m^3 b dBZ
+    WHERE(TPROFILERS(JP)%XCRARE(IN,:)>0)
+      TPROFILERS(JP)%XCRARE(IN,:)=10.*LOG10(1.E18*TPROFILERS(JP)%XCRARE(IN,:))
+    ELSEWHERE
+      TPROFILERS(JP)%XCRARE(IN,:)=XUNDEF
+    END WHERE
+    WHERE(TPROFILERS(JP)%XCRARE_ATT(IN,:)>0)
+      TPROFILERS(JP)%XCRARE_ATT(IN,:)=10.*LOG10(1.E18*TPROFILERS(JP)%XCRARE_ATT(IN,:))
+    ELSEWHERE
+      TPROFILERS(JP)%XCRARE_ATT(IN,:)=XUNDEF
+    END WHERE
+    DEALLOCATE(ZX,ZW,ZRTMIN)
+  END IF ! end LOOP ICE3
 ! end add RARE
 !!
-      IF (.NOT. L1D) THEN
-        TPROFILER%P   (IN,:,I) = PROFILER_INTERP(PP(II(I):II(I)+1,IJ(I):IJ(I)+1,:))
-      ELSE
-        TPROFILER%P   (IN,:,I) = PROFILER_INTERP(PP)
-      END IF
-      !
-      DO JSV=1,SIZE(PR,4)
-        TPROFILER%R   (IN,:,I,JSV) = PROFILER_INTERP(PR(:,:,:,JSV))
-      END DO
-        ZWORK(:,:,:,:)=PSV(:,:,:,:)
-        ZWORK(:,:,1,:)=PSV(:,:,2,:)
-      DO JSV=1,SIZE(PSV,4)
-        TPROFILER%SV  (IN,:,I,JSV) = PROFILER_INTERP(ZWORK(:,:,:,JSV))
-      END DO
-      ZWORK2(:,:,:,:) = 0.
-      DO JK=IKB,IKE
-        IKRAD = JK - JPVEXT
-        ZWORK2(:,:,JK,:)=PAER(:,:,IKRAD,:)
-      ENDDO
-      DO JSV=1,SIZE(PAER,4)
-        TPROFILER%AER(IN,:,I,JSV) = PROFILER_INTERP(ZWORK2(:,:,:,JSV))
-      ENDDO
-      IF (SIZE(PTKE)>0) TPROFILER%TKE  (IN,:,I) = PROFILER_INTERP(PTKE)
-      !
-      IF (LDIAG_IN_RUN) THEN
-        TPROFILER%T2M   (IN,I)     = PROFILER_INTERP_2D(XCURRENT_T2M   )
-        TPROFILER%Q2M   (IN,I)     = PROFILER_INTERP_2D(XCURRENT_Q2M   )
-        TPROFILER%HU2M  (IN,I)     = PROFILER_INTERP_2D(XCURRENT_HU2M  )
-        TPROFILER%ZON10M(IN,I)     = PROFILER_INTERP_2D(XCURRENT_ZON10M)
-        TPROFILER%MER10M(IN,I)     = PROFILER_INTERP_2D(XCURRENT_MER10M)
-        TPROFILER%RN    (IN,I)     = PROFILER_INTERP_2D(XCURRENT_RN    )
-        TPROFILER%H     (IN,I)     = PROFILER_INTERP_2D(XCURRENT_H     )
-        TPROFILER%LE    (IN,I)     = PROFILER_INTERP_2D(XCURRENT_LE    )
-        TPROFILER%LEI   (IN,I)     = PROFILER_INTERP_2D(XCURRENT_LEI   )        
-        TPROFILER%GFLUX (IN,I)     = PROFILER_INTERP_2D(XCURRENT_GFLUX )
-       IF (CRAD /= 'NONE') THEN
-        TPROFILER%SWD   (IN,I)     = PROFILER_INTERP_2D(XCURRENT_SWD   )
-        TPROFILER%SWU   (IN,I)     = PROFILER_INTERP_2D(XCURRENT_SWU   )
-        TPROFILER%LWD   (IN,I)     = PROFILER_INTERP_2D(XCURRENT_LWD   )
-        TPROFILER%LWU   (IN,I)     = PROFILER_INTERP_2D(XCURRENT_LWU   )
-       END IF
-        TPROFILER%TKE_DISS(IN,:,I) = PROFILER_INTERP(XCURRENT_TKE_DISS)
-      ENDIF
-    ENDIF
-!
-!----------------------------------------------------------------------------
-!
-!*     11.   EXCHANGE OF INFORMATION BETWEEN PROCESSORS
-!            ------------------------------------------
-!
-!*     11.2  data stored
-!            -----------
-!
-  CALL DISTRIBUTE_PROFILER(TPROFILER%X   (I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%Y   (I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%LON (I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%LAT (I))
-  !
-  IF (LDIAG_IN_RUN) THEN
-    CALL DISTRIBUTE_PROFILER(TPROFILER%T2M   (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%Q2M   (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%HU2M  (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%ZON10M(IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%MER10M(IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%RN    (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%H     (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%LE    (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%LEI   (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%GFLUX (IN,I))
-   IF (CRAD /= 'NONE') THEN
-    CALL DISTRIBUTE_PROFILER(TPROFILER%LWD   (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%LWU   (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%SWD   (IN,I))
-    CALL DISTRIBUTE_PROFILER(TPROFILER%SWU   (IN,I))
-   ENDIF
-  ENDIF
- DO JK=1,IKU
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZON (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%MER (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%FF  (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%DD  (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%W   (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%P   (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZZ  (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%TH  (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%THV (IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%VISI(IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%VISIKUN(IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%RHOD(IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%CRARE(IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%CRARE_ATT(IN,JK,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%CIZ(IN,JK,I))
-
-  !
-  IF (LDIAG_IN_RUN) CALL DISTRIBUTE_PROFILER(TPROFILER%TKE_DISS(IN,JK,I))
+  TPROFILERS(JP)%XP   (IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), PP )
   !
   DO JSV=1,SIZE(PR,4)
-    CALL DISTRIBUTE_PROFILER(TPROFILER%R   (IN,JK,I,JSV))
+    TPROFILERS(JP)%XR   (IN,:,JSV) = STATPROF_INTERP_3D( TPROFILERS(JP), PR(:,:,:,JSV) )
   END DO
+    ZWORK(:,:,:,:)=PSV(:,:,:,:)
+    ZWORK(:,:,1,:)=PSV(:,:,2,:)
   DO JSV=1,SIZE(PSV,4)
-    CALL DISTRIBUTE_PROFILER(TPROFILER%SV  (IN,JK,I,JSV))
+    TPROFILERS(JP)%XSV  (IN,:,JSV) = STATPROF_INTERP_3D( TPROFILERS(JP), ZWORK(:,:,:,JSV) )
   END DO
-  IF (SIZE(PTKE)>0) CALL DISTRIBUTE_PROFILER(TPROFILER%TKE  (IN,JK,I))
- ENDDO
-
-  CALL DISTRIBUTE_PROFILER(TPROFILER%IWV(IN,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZTD(IN,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZHD(IN,I))
-  CALL DISTRIBUTE_PROFILER(TPROFILER%ZWD(IN,I))
-ENDDO
-!
-END IF
-!
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-!
-CONTAINS
-!
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-FUNCTION PROFILER_INTERP_2D(PA) RESULT(PB)
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PA
-REAL                             :: PB
-!
-INTEGER :: JI, JJ
-!
-IF (SIZE(PA,1)==2) THEN
-  JI=1
-  JJ=1
-ELSEIF (L1D) THEN
-     JI=2
-     JJ=2
-ELSE
-  JI=II(I)
-  JJ=IJ(I)
-END IF
-!
-!
-PB = (1.-ZYCOEF(I)) * (1.-ZXCOEF(I)) *  PA(JI,JJ)    + &
-     (1.-ZYCOEF(I)) *    (ZXCOEF(I)) *  PA(JI+1,JJ)  + &
-     (   ZYCOEF(I)) * (1.-ZXCOEF(I)) *  PA(JI,JJ+1)  + &
-     (   ZYCOEF(I)) *    (ZXCOEF(I)) *  PA(JI+1,JJ+1)
-!
-END FUNCTION PROFILER_INTERP_2D
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-FUNCTION PROFILER_INTERP(PA) RESULT(PB)
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PA
-REAL, DIMENSION(SIZE(PA,3))        :: PB
-!
-INTEGER :: JI, JJ,JK
-!
-IF (SIZE(PA,1)==2) THEN
-  JI=1
-  JJ=1
-ELSEIF (L1D) THEN
-     JI=2
-     JJ=2
-ELSE
-  JI=II(I)
-  JJ=IJ(I)
-END IF
-!
-!
-DO JK=1,SIZE(PA,3)
- IF ( (PA(JI,JJ,JK) /= XUNDEF) .AND. (PA(JI+1,JJ,JK) /= XUNDEF) .AND. &
-      (PA(JI,JJ+1,JK) /= XUNDEF) .AND. (PA(JI+1,JJ+1,JK) /= XUNDEF) ) THEN
-    PB(JK) = (1.-ZYCOEF(I)) * (1.-ZXCOEF(I)) *  PA(JI,JJ,JK) + &
-            (1.-ZYCOEF(I)) * (ZXCOEF(I)) *  PA(JI+1,JJ,JK)  + &
-            (ZYCOEF(I)) * (1.-ZXCOEF(I)) *  PA(JI,JJ+1,JK)  + &
-            (ZYCOEF(I)) * (ZXCOEF(I)) *  PA(JI+1,JJ+1,JK) 
- ELSE
-    PB(JK) = XUNDEF 
- END IF
-END DO
-!
-END FUNCTION PROFILER_INTERP
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-FUNCTION PROFILER_INTERP_U(PA) RESULT(PB)
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PA
-REAL, DIMENSION(SIZE(PA,3))        :: PB
-!
-INTEGER :: JI, JJ
-!
-IF (SIZE(PA,1)==2) THEN
-  JI=1
-  JJ=1
-ELSEIF (L1D) THEN
-     JI=2
-     JJ=2
-ELSE
-  JI=IU(I)
-  JJ=IJ(I)
-END IF
-!
-PB(:) = (1.- ZYCOEF(I)) * (1.-ZUCOEF(I)) * PA(JI  ,JJ  ,:) &
-   +    (1.- ZYCOEF(I)) * (   ZUCOEF(I)) * PA(JI+1,JJ  ,:) &
-   + (    ZYCOEF(I)) * (1.-ZUCOEF(I)) *    PA(JI  ,JJ+1,:) &
-   + (    ZYCOEF(I)) * (   ZUCOEF(I)) *    PA(JI+1,JJ+1,:)
-!
-END FUNCTION PROFILER_INTERP_U
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-FUNCTION PROFILER_INTERP_V(PA) RESULT(PB)
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PA
-REAL, DIMENSION(SIZE(PA,3))        :: PB
-!
-INTEGER :: JI, JJ
-!
-IF (SIZE(PA,1)==2) THEN
-  JI=1
-  JJ=1
-ELSEIF (L1D) THEN
-     JI=2
-     JJ=2
-ELSE
-  JI=II(I)
-  JJ=IV(I)
-END IF
-!
-PB(:) = (1.- ZVCOEF(I)) * (1.-ZXCOEF(I)) * PA(JI  ,JJ  ,:) &
-   + (1.- ZVCOEF(I)) * (   ZXCOEF(I)) *  PA(JI+1,JJ  ,:)   &
-   + (    ZVCOEF(I)) * (1.-ZXCOEF(I)) *  PA(JI  ,JJ+1,:)   &
-   + (    ZVCOEF(I)) * (   ZXCOEF(I)) *  PA(JI+1,JJ+1,:) 
-!
-END FUNCTION PROFILER_INTERP_V
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-SUBROUTINE DISTRIBUTE_PROFILER(PAS)
-!
-REAL, INTENT(INOUT) :: PAS
-!
-PAS = PAS * ZTHIS_PROCS(I)
-
-CALL REDUCESUM_ll(PAS,IINFO_ll)
+  ZWORK2(:,:,:,:) = 0.
+  DO JK=IKB,IKE
+    IKRAD = JK - JPVEXT
+    ZWORK2(:,:,JK,:)=PAER(:,:,IKRAD,:)
+  END DO
+  DO JSV=1,SIZE(PAER,4)
+    TPROFILERS(JP)%XAER(IN,:,JSV) = STATPROF_INTERP_3D( TPROFILERS(JP), ZWORK2(:,:,:,JSV) )
+  END DO
+  IF (SIZE(PTKE)>0) TPROFILERS(JP)%XTKE  (IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), PTKE )
+  !
+  IF (LDIAG_IN_RUN) THEN
+    TPROFILERS(JP)%XT2M   (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_T2M    )
+    TPROFILERS(JP)%XQ2M   (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_Q2M    )
+    TPROFILERS(JP)%XHU2M  (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_HU2M   )
+    TPROFILERS(JP)%XZON10M(IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_ZON10M )
+    TPROFILERS(JP)%XMER10M(IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_MER10M )
+    TPROFILERS(JP)%XRN    (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_RN     )
+    TPROFILERS(JP)%XH     (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_H      )
+    TPROFILERS(JP)%XLE    (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_LE     )
+    TPROFILERS(JP)%XLEI   (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_LEI    )
+    TPROFILERS(JP)%XGFLUX (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_GFLUX  )
+    IF (CRAD /= 'NONE') THEN
+      TPROFILERS(JP)%XSWD   (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_SWD    )
+      TPROFILERS(JP)%XSWU   (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_SWU    )
+      TPROFILERS(JP)%XLWD   (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_LWD    )
+      TPROFILERS(JP)%XLWU   (IN)     = STATPROF_INTERP_2D( TPROFILERS(JP), XCURRENT_LWU    )
+    END IF
+    TPROFILERS(JP)%XTKE_DISS(IN,:) = STATPROF_INTERP_3D( TPROFILERS(JP), XCURRENT_TKE_DISS )
+  END IF
+END DO PROFILER
 !
-END SUBROUTINE DISTRIBUTE_PROFILER
 !----------------------------------------------------------------------------
 !
 END SUBROUTINE PROFILER_n
diff --git a/src/MNH/stationn.f90 b/src/MNH/stationn.f90
index 41e95dc53..b10c060a5 100644
--- a/src/MNH/stationn.f90
+++ b/src/MNH/stationn.f90
@@ -9,14 +9,11 @@ MODULE MODI_STATION_n
 !
 INTERFACE
 !
-      SUBROUTINE STATION_n(PTSTEP,                               &
-                           PXHAT, PYHAT, PZ,                     &
-                           PU, PV, PW, PTH, PR, PSV, PTKE,       &
-                           PTS,PP ) 
+      SUBROUTINE STATION_n( PTSTEP, PZ,                     &
+                            PU, PV, PW, PTH, PR, PSV, PTKE, &
+                            PTS, PP )
 !
 REAL,                     INTENT(IN)     :: PTSTEP ! time step
-REAL, DIMENSION(:),       INTENT(IN)     :: PXHAT  ! x coordinate
-REAL, DIMENSION(:),       INTENT(IN)     :: PYHAT  ! y coordinate
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PZ     ! z array
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PU     ! horizontal wind X component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PV     ! horizontal wind Y component
@@ -36,12 +33,11 @@ END INTERFACE
 !
 END MODULE MODI_STATION_n
 !
-!     ########################################################
-      SUBROUTINE STATION_n(PTSTEP,                           &
-                       PXHAT, PYHAT, PZ,                     &
-                       PU, PV, PW, PTH, PR, PSV, PTKE,       &
-                       PTS, PP )
-!     ########################################################
+!     #######################################################
+      SUBROUTINE STATION_n( PTSTEP, PZ,                     &
+                            PU, PV, PW, PTH, PR, PSV, PTKE, &
+                            PTS, PP )
+!     #######################################################
 !
 !
 !!****  *STATION_n* - (advects and) stores 
@@ -97,7 +93,7 @@ USE MODD_STATION_n
 USE MODD_ALLSTATION_n,  ONLY: LDIAG_SURFRAD
 USE MODD_TIME_n,        ONLY: tdtcur
 !
-USE MODE_STATPROF_TOOLS, ONLY: STATION_INTERP_2D, STATION_INTERP_2D_U, STATION_INTERP_2D_V
+USE MODE_STATPROF_TOOLS, ONLY: STATPROF_INTERP_2D, STATPROF_INTERP_2D_U, STATPROF_INTERP_2D_V
 !
 !
 IMPLICIT NONE
@@ -107,8 +103,6 @@ IMPLICIT NONE
 !
 !
 REAL,                     INTENT(IN)     :: PTSTEP ! time step
-REAL, DIMENSION(:),       INTENT(IN)     :: PXHAT  ! x coordinate
-REAL, DIMENSION(:),       INTENT(IN)     :: PYHAT  ! y coordinate
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PZ     ! z array
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PU     ! horizontal wind X component
 REAL, DIMENSION(:,:,:),   INTENT(IN)     :: PV     ! horizontal wind Y component
@@ -163,52 +157,52 @@ STATION: DO JS = 1,NUMBSTAT_LOC
   JK = TSTATIONS(JS)%NK
 
   IF (LCARTESIAN) THEN
-    TSTATIONS(JS)%XZON(IN) =   STATION_INTERP_2D_U( TSTATIONS(JS), PU(:,:,JK) )
-    TSTATIONS(JS)%XMER(IN) =   STATION_INTERP_2D_V( TSTATIONS(JS), PV(:,:,JK) )
+    TSTATIONS(JS)%XZON(IN) =   STATPROF_INTERP_2D_U( TSTATIONS(JS), PU(:,:,JK) )
+    TSTATIONS(JS)%XMER(IN) =   STATPROF_INTERP_2D_V( TSTATIONS(JS), PV(:,:,JK) )
   ELSE
-    ZU_STAT                = STATION_INTERP_2D_U( TSTATIONS(JS), PU(:,:,JK) )
-    ZV_STAT                = STATION_INTERP_2D_V( TSTATIONS(JS), PV(:,:,JK) )
+    ZU_STAT                = STATPROF_INTERP_2D_U( TSTATIONS(JS), PU(:,:,JK) )
+    ZV_STAT                = STATPROF_INTERP_2D_V( TSTATIONS(JS), PV(:,:,JK) )
     ZGAM                   = (XRPK * (TSTATIONS(JS)%XLON - XLON0) - XBETA)*(XPI/180.)
     TSTATIONS(JS)%XZON(IN) =   ZU_STAT * COS(ZGAM) + ZV_STAT * SIN(ZGAM)
     TSTATIONS(JS)%XMER(IN) = - ZU_STAT * SIN(ZGAM) + ZV_STAT * COS(ZGAM)
   END IF
-  TSTATIONS(JS)%XW (IN) = STATION_INTERP_2D( TSTATIONS(JS), PW(:,:,JK) )
-  TSTATIONS(JS)%XTH(IN) = STATION_INTERP_2D( TSTATIONS(JS), PTH(:,:,JK) )
-  TSTATIONS(JS)%XP (IN) = STATION_INTERP_2D( TSTATIONS(JS), PP(:,:,JK) )
+  TSTATIONS(JS)%XW (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), PW(:,:,JK) )
+  TSTATIONS(JS)%XTH(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), PTH(:,:,JK) )
+  TSTATIONS(JS)%XP (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), PP(:,:,JK) )
 
   DO JSV=1,SIZE(PR,4)
-    TSTATIONS(JS)%XR(IN,JSV) = STATION_INTERP_2D( TSTATIONS(JS), PR(:,:,JK,JSV) )
+    TSTATIONS(JS)%XR(IN,JSV) = STATPROF_INTERP_2D( TSTATIONS(JS), PR(:,:,JK,JSV) )
   END DO
 
   DO JSV=1,SIZE(PSV,4)
-    TSTATIONS(JS)%XSV(IN,JSV) = STATION_INTERP_2D( TSTATIONS(JS), PSV(:,:,JK,JSV) )
+    TSTATIONS(JS)%XSV(IN,JSV) = STATPROF_INTERP_2D( TSTATIONS(JS), PSV(:,:,JK,JSV) )
   END DO
 
-  IF (SIZE(PTKE)>0) TSTATIONS(JS)%XTKE(IN) = STATION_INTERP_2D( TSTATIONS(JS), PTKE(:,:,JK) )
-  IF ( CRAD /= 'NONE' ) TSTATIONS(JS)%XTSRAD(IN) = STATION_INTERP_2D( TSTATIONS(JS), PTS )
-  TSTATIONS(JS)%XZS  = STATION_INTERP_2D( TSTATIONS(JS), PZ(:,:,1+JPVEXT))
+  IF (SIZE(PTKE)>0) TSTATIONS(JS)%XTKE(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), PTKE(:,:,JK) )
+  IF ( CRAD /= 'NONE' ) TSTATIONS(JS)%XTSRAD(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), PTS )
+  TSTATIONS(JS)%XZS  = STATPROF_INTERP_2D( TSTATIONS(JS), PZ(:,:,1+JPVEXT))
 
   IF ( LDIAG_SURFRAD ) THEN
-    TSTATIONS(JS)%XZON10M(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_ZON10M )
-    TSTATIONS(JS)%XMER10M(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_MER10M )
-    TSTATIONS(JS)%XT2M   (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_T2M    )
-    TSTATIONS(JS)%XQ2M   (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_Q2M    )
-    TSTATIONS(JS)%XHU2M  (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_HU2M   )
-    TSTATIONS(JS)%XRN    (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_RN     )
-    TSTATIONS(JS)%XH     (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_H      )
-    TSTATIONS(JS)%XLE    (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LE     )
-    TSTATIONS(JS)%XLEI   (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LEI    )
-    TSTATIONS(JS)%XGFLUX (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_GFLUX  )
+    TSTATIONS(JS)%XZON10M(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_ZON10M )
+    TSTATIONS(JS)%XMER10M(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_MER10M )
+    TSTATIONS(JS)%XT2M   (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_T2M    )
+    TSTATIONS(JS)%XQ2M   (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_Q2M    )
+    TSTATIONS(JS)%XHU2M  (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_HU2M   )
+    TSTATIONS(JS)%XRN    (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_RN     )
+    TSTATIONS(JS)%XH     (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_H      )
+    TSTATIONS(JS)%XLE    (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_LE     )
+    TSTATIONS(JS)%XLEI   (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_LEI    )
+    TSTATIONS(JS)%XGFLUX (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_GFLUX  )
     IF ( CRAD /= 'NONE' ) THEN
-      TSTATIONS(JS)%XSWD   (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWD    )
-      TSTATIONS(JS)%XSWU   (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWU    )
-      TSTATIONS(JS)%XLWD   (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LWD    )
-      TSTATIONS(JS)%XLWU   (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_LWU    )
-      TSTATIONS(JS)%XSWDIR (IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWDIR  )
-      TSTATIONS(JS)%XSWDIFF(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SWDIFF )
-      TSTATIONS(JS)%XDSTAOD(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_DSTAOD )
+      TSTATIONS(JS)%XSWD   (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_SWD    )
+      TSTATIONS(JS)%XSWU   (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_SWU    )
+      TSTATIONS(JS)%XLWD   (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_LWD    )
+      TSTATIONS(JS)%XLWU   (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_LWU    )
+      TSTATIONS(JS)%XSWDIR (IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_SWDIR  )
+      TSTATIONS(JS)%XSWDIFF(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_SWDIFF )
+      TSTATIONS(JS)%XDSTAOD(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_DSTAOD )
     END IF
-    TSTATIONS(JS)%XSFCO2(IN) = STATION_INTERP_2D( TSTATIONS(JS), XCURRENT_SFCO2 )
+    TSTATIONS(JS)%XSFCO2(IN) = STATPROF_INTERP_2D( TSTATIONS(JS), XCURRENT_SFCO2 )
   END IF
 END DO STATION
 !
diff --git a/src/MNH/statprof_reader.f90 b/src/MNH/statprof_reader.f90
index ed645d7ab..888cfeac3 100644
--- a/src/MNH/statprof_reader.f90
+++ b/src/MNH/statprof_reader.f90
@@ -11,14 +11,14 @@ IMPLICIT NONE
 
 PRIVATE
 
-PUBLIC :: READ_CSV_STATION
+PUBLIC :: STATPROF_CSV_READ
 
 INTEGER, PARAMETER :: NMAXLINELGT = 400
 
 CONTAINS
 !-------------------------------------------------------------------
 !
-!!****  *READ_CSV_STATION* -
+!!****  *STATPROF_CSV_READ* -
 !!
 !!    PURPOSE
 !!    -------
@@ -31,39 +31,55 @@ CONTAINS
 !!    MODIFICATIONS
 !!    -------------
 !!     03/2020      Original
-!  P. Wautelet    04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
+!  P. Wautelet    04/2022: restructure stations/profilers for better performance, reduce memory usage and correct some problems/bugs
 !---------------------------------------------------------------
 !
-!###############################################################################################
-SUBROUTINE READ_CSV_STATION( HFILE, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,                    &
-                             PXHATM_PHYS_MIN, PXHATM_PHYS_MAX,PYHATM_PHYS_MIN, PYHATM_PHYS_MAX )
-!###############################################################################################
+!#################################################################################################
+SUBROUTINE STATPROF_CSV_READ( TPSTATPROF, HFILE, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,         &
+                              PXHATM_PHYS_MIN, PXHATM_PHYS_MAX,PYHATM_PHYS_MIN, PYHATM_PHYS_MAX, &
+                              KNUMBSTATPROF                                                      )
+!#################################################################################################
 
 USE MODD_CONF,          ONLY: LCARTESIAN
-USE MODD_STATION_n,     ONLY: NUMBSTAT
-USE MODD_TYPE_STATPROF, ONLY: TSTATIONDATA
+USE MODD_TYPE_STATPROF, ONLY: TPROFILERDATA, TSTATIONDATA, TSTATPROFDATA
 
 USE MODE_MSG
-USE MODE_STATPROF_TOOLS, ONLY: STATION_ADD, STATION_INI_INTERP, STATION_POSITION
-
-CHARACTER(LEN=*),   INTENT(IN) :: HFILE ! file to read
-REAL, DIMENSION(:), INTENT(IN) :: PXHAT_GLOB
-REAL, DIMENSION(:), INTENT(IN) :: PYHAT_GLOB
-REAL, DIMENSION(:), INTENT(IN) :: PXHATM ! mass point coordinates
-REAL, DIMENSION(:), INTENT(IN) :: PYHATM ! mass point coordinates
-REAL,               INTENT(IN) :: PXHATM_PHYS_MIN, PYHATM_PHYS_MIN  ! Minimum X coordinate of mass points in the physical domain
-REAL,               INTENT(IN) :: PXHATM_PHYS_MAX, PYHATM_PHYS_MAX  ! Minimum X coordinate of mass points in the physical domain
+USE MODE_STATPROF_TOOLS, ONLY: PROFILER_ADD, STATION_ADD, STATPROF_INI_INTERP, STATPROF_POSITION
+
+CLASS(TSTATPROFDATA), INTENT(IN)  :: TPSTATPROF ! Used only to identify datatype
+CHARACTER(LEN=*),     INTENT(IN)  :: HFILE ! file to read
+REAL, DIMENSION(:),   INTENT(IN)  :: PXHAT_GLOB
+REAL, DIMENSION(:),   INTENT(IN)  :: PYHAT_GLOB
+REAL, DIMENSION(:),   INTENT(IN)  :: PXHATM ! mass point coordinates
+REAL, DIMENSION(:),   INTENT(IN)  :: PYHATM ! mass point coordinates
+REAL,                 INTENT(IN)  :: PXHATM_PHYS_MIN, PYHATM_PHYS_MIN  ! Minimum X coordinate of mass points in the physical domain
+REAL,                 INTENT(IN)  :: PXHATM_PHYS_MAX, PYHATM_PHYS_MAX  ! Minimum X coordinate of mass points in the physical domain
+INTEGER,              INTENT(OUT) :: KNUMBSTATPROF ! Total number of stations/profilers (inside physical domain of model)
 !
 CHARACTER(LEN=NMAXLINELGT) :: YSTRING
-INTEGER            :: ILU      ! logical unit of the file
-INTEGER            :: INBLINE  ! Nb of lines in csv file
-INTEGER            :: JI
-LOGICAL            :: GINSIDE  ! True if station is inside physical domain of model
-LOGICAL            :: GPRESENT ! True if station is present on the current process
-TYPE(TSTATIONDATA) :: TZSTATION
+INTEGER             :: ILU      ! logical unit of the file
+INTEGER             :: INBLINE  ! Nb of lines in csv file
+INTEGER             :: JI
+LOGICAL             :: GINSIDE  ! True if station/profiler is inside physical domain of model
+LOGICAL             :: GPRESENT ! True if station/profiler is present on the current process
+TYPE(TSTATIONDATA),  TARGET :: TZSTATION
+TYPE(TPROFILERDATA), TARGET :: TZPROFILER
 
-INBLINE  = 0 !Number of stations found in the file
-NUMBSTAT = 0 !Number of stations found in the file AND inside the model domain
+CLASS(TSTATPROFDATA), POINTER :: TZSTATPROF
+
+SELECT TYPE( TPSTATPROF )
+  TYPE IS( TPROFILERDATA )
+    TZSTATPROF => TZPROFILER
+
+  TYPE IS( TSTATIONDATA )
+    TZSTATPROF => TZSTATION
+
+  CLASS DEFAULT
+    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_CSV_READ', 'unknown type for TPSTATPROF' )
+END SELECT
+
+INBLINE       = 0 !Number of stations/profilers found in the file
+KNUMBSTATPROF = 0 !Number of stations/profilers found in the file AND inside the model domain
 
 ! Open file
 OPEN( NEWUNIT = ILU, FILE = HFILE, FORM = 'formatted' )
@@ -71,29 +87,40 @@ OPEN( NEWUNIT = ILU, FILE = HFILE, FORM = 'formatted' )
 READ( ILU, END = 101, FMT = '(A)' ) YSTRING ! Reading of header (skip it)
 
 DO
-  ! Read station coordinates
+  ! Read station/profiler coordinates
   READ( ILU, END = 101, FMT = '(A)' ) YSTRING
 
   ! Check if record is written in French convention
   CALL FRENCH_TO_ENGLISH( YSTRING )
 
   IF ( LCARTESIAN ) THEN
-    READ( YSTRING, * ) TZSTATION%CNAME, TZSTATION%XX,   TZSTATION%XY,   TZSTATION%XZ
+    READ( YSTRING, * ) TZSTATPROF%CNAME, TZSTATPROF%XX,   TZSTATPROF%XY,   TZSTATPROF%XZ
   ELSE
-    READ( YSTRING, * ) TZSTATION%CNAME, TZSTATION%XLAT, TZSTATION%XLON, TZSTATION%XZ
+    READ( YSTRING, * ) TZSTATPROF%CNAME, TZSTATPROF%XLAT, TZSTATPROF%XLON, TZSTATPROF%XZ
   END IF
 
-  IF ( .NOT. LCARTESIAN ) CALL STATION_INI_INTERP( TZSTATION )
-  CALL STATION_POSITION( TZSTATION, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,                  &
+  IF ( .NOT. LCARTESIAN ) CALL STATPROF_INI_INTERP( TZSTATPROF )
+  CALL STATPROF_POSITION( TZSTATPROF, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,                 &
                          PXHATM_PHYS_MIN, PXHATM_PHYS_MAX, PYHATM_PHYS_MIN, PYHATM_PHYS_MAX, &
                          GINSIDE, GPRESENT                                                   )
 
   IF ( GINSIDE ) THEN
-    NUMBSTAT = NUMBSTAT + 1
-    TZSTATION%NID = NUMBSTAT
+    KNUMBSTATPROF = KNUMBSTATPROF + 1
+    TZSTATPROF%NID = KNUMBSTATPROF
   END IF
 
-  IF ( GPRESENT ) CALL STATION_ADD( TZSTATION )
+  IF ( GPRESENT ) THEN
+    SELECT TYPE( TZSTATPROF )
+      TYPE IS( TPROFILERDATA )
+        CALL PROFILER_ADD( TZSTATPROF )
+
+      TYPE IS( TSTATIONDATA )
+        CALL STATION_ADD( TZSTATPROF )
+
+      CLASS DEFAULT
+        CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_CSV_READ', 'unknown type for TPSTATPROF' )
+    END SELECT
+  END IF
 
   INBLINE = INBLINE + 1
 END DO
@@ -102,9 +129,9 @@ END DO
 
 CLOSE( ILU )
 
-IF ( INBLINE == 0 ) CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'READ_CSV_STATION', 'Data not found in file ' // TRIM( HFILE ) )
+IF ( INBLINE == 0 ) CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_CSV_READ', 'Data not found in file ' // TRIM( HFILE ) )
 
-END SUBROUTINE READ_CSV_STATION
+END SUBROUTINE STATPROF_CSV_READ
 
 !#########################################################
 SUBROUTINE FRENCH_TO_ENGLISH(HSTRING)
diff --git a/src/MNH/statprof_tools.f90 b/src/MNH/statprof_tools.f90
index 0a5e466c3..ca8ff988c 100644
--- a/src/MNH/statprof_tools.f90
+++ b/src/MNH/statprof_tools.f90
@@ -13,20 +13,146 @@
 MODULE MODE_STATPROF_TOOLS
 !      ###################
 
-USE MODD_TYPE_STATPROF, ONLY: TSTATIONDATA
+USE MODD_TYPE_STATPROF, ONLY: TPROFILERDATA, TSTATIONDATA, TSTATPROFDATA
 
 IMPLICIT NONE
 
 PRIVATE
 
-PUBLIC :: STATION_ALLOCATE
-PUBLIC :: STATION_INI_INTERP
-PUBLIC :: STATION_POSITION
-PUBLIC :: STATION_ADD
-PUBLIC :: STATION_INTERP_2D, STATION_INTERP_2D_U, STATION_INTERP_2D_V
+PUBLIC :: PROFILER_ALLOCATE, STATION_ALLOCATE
+PUBLIC :: STATPROF_INI_INTERP
+PUBLIC :: STATPROF_POSITION
+PUBLIC :: PROFILER_ADD, STATION_ADD
+PUBLIC :: STATPROF_INTERP_2D, STATPROF_INTERP_2D_U, STATPROF_INTERP_2D_V
+PUBLIC :: STATPROF_INTERP_3D, STATPROF_INTERP_3D_U, STATPROF_INTERP_3D_V
 
 CONTAINS
 
+! ################################################
+SUBROUTINE PROFILER_ALLOCATE( TPPROFILER, KSTORE )
+! ################################################
+
+!   USE MODD_ALLSTATION_n, ONLY: LDIAG_SURFRAD
+  USE MODD_CONF_n,       ONLY: NRR
+  USE MODD_DIAG_IN_RUN,  ONLY: LDIAG_IN_RUN
+  USE MODD_DIM_n,        ONLY: NKMAX
+  USE MODD_NSV,          ONLY: NSV
+  USE MODD_PARAMETERS,   ONLY: JPVEXT, XUNDEF
+  USE MODD_PARAM_n,      ONLY: CCLOUD, CRAD, CTURB
+  USE MODD_RADIATIONS_n, ONLY: NAER
+
+  IMPLICIT NONE
+
+  TYPE(TPROFILERDATA), INTENT(INOUT) :: TPPROFILER
+  INTEGER,             INTENT(IN)    :: KSTORE  ! number of moments to store
+
+  INTEGER :: IKU
+
+  IKU = NKMAX + 2 * JPVEXT
+  ALLOCATE( TPPROFILER%XZON      (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XMER      (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XFF       (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XDD       (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XW        (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XP        (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XZZ       (KSTORE, IKU) )
+  IF ( CTURB == 'TKEL' ) THEN
+    ALLOCATE( TPPROFILER%XTKE    (KSTORE, IKU) )
+  ELSE
+    ALLOCATE( TPPROFILER%XTKE    (0, 0) )
+  END IF
+  ALLOCATE( TPPROFILER%XTH       (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XTHV      (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XVISI     (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XVISIKUN  (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XCRARE    (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XCRARE_ATT(KSTORE, IKU) )
+  IF ( CCLOUD == 'ICE3' .OR. CCLOUD == 'ICE4' ) THEN
+    ALLOCATE( TPPROFILER%XCIZ    (KSTORE, IKU) )
+  ELSE
+    ALLOCATE( TPPROFILER%XCIZ    (0, 0) )
+  END IF
+  ALLOCATE( TPPROFILER%XLWCZ     (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XIWCZ     (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XRHOD     (KSTORE, IKU) )
+  ALLOCATE( TPPROFILER%XR        (KSTORE, IKU, NRR) )
+  ALLOCATE( TPPROFILER%XSV       (KSTORE, IKU, NSV) )
+  ALLOCATE( TPPROFILER%XAER      (KSTORE, IKU, NAER) )
+
+  ALLOCATE( TPPROFILER%XIWV(KSTORE) )
+  ALLOCATE( TPPROFILER%XZTD(KSTORE) )
+  ALLOCATE( TPPROFILER%XZWD(KSTORE) )
+  ALLOCATE( TPPROFILER%XZHD(KSTORE) )
+
+!   IF ( LDIAG_IN_RUN ) THEN
+    ALLOCATE( TPPROFILER%XT2M   (KSTORE) )
+    ALLOCATE( TPPROFILER%XQ2M   (KSTORE) )
+    ALLOCATE( TPPROFILER%XHU2M  (KSTORE) )
+    ALLOCATE( TPPROFILER%XZON10M(KSTORE) )
+    ALLOCATE( TPPROFILER%XMER10M(KSTORE) )
+    ALLOCATE( TPPROFILER%XRN    (KSTORE) )
+    ALLOCATE( TPPROFILER%XH     (KSTORE) )
+    ALLOCATE( TPPROFILER%XLE    (KSTORE) )
+    ALLOCATE( TPPROFILER%XLEI   (KSTORE) )
+    ALLOCATE( TPPROFILER%XGFLUX (KSTORE) )
+    IF ( CRAD /= 'NONE' ) THEN
+      ALLOCATE( TPPROFILER%XSWD   (KSTORE) )
+      ALLOCATE( TPPROFILER%XSWU   (KSTORE) )
+      ALLOCATE( TPPROFILER%XLWD   (KSTORE) )
+      ALLOCATE( TPPROFILER%XLWU   (KSTORE) )
+    END IF
+    ALLOCATE( TPPROFILER%XTKE_DISS(KSTORE, IKU) )
+!   END IF
+
+  TPPROFILER%XZON      (:,:) = XUNDEF
+  TPPROFILER%XMER      (:,:) = XUNDEF
+  TPPROFILER%XFF       (:,:) = XUNDEF
+  TPPROFILER%XDD       (:,:) = XUNDEF
+  TPPROFILER%XW        (:,:) = XUNDEF
+  TPPROFILER%XP        (:,:) = XUNDEF
+  TPPROFILER%XZZ       (:,:) = XUNDEF
+  IF ( CTURB == 'TKEL' ) TPPROFILER%XTKE(:,:) = XUNDEF
+  TPPROFILER%XTH       (:,:) = XUNDEF
+  TPPROFILER%XTHV      (:,:) = XUNDEF
+  TPPROFILER%XVISI     (:,:) = XUNDEF
+  TPPROFILER%XVISIKUN  (:,:) = XUNDEF
+  TPPROFILER%XCRARE    (:,:) = XUNDEF
+  TPPROFILER%XCRARE_ATT(:,:) = XUNDEF
+  IF ( CCLOUD == 'ICE3' .OR. CCLOUD == 'ICE4' ) TPPROFILER%XCIZ      (:,:) = XUNDEF
+  TPPROFILER%XLWCZ     (:,:) = XUNDEF
+  TPPROFILER%XIWCZ     (:,:) = XUNDEF
+  TPPROFILER%XRHOD     (:,:) = XUNDEF
+  TPPROFILER%XR        (:,:,:) = XUNDEF
+  TPPROFILER%XSV       (:,:,:) = XUNDEF
+  TPPROFILER%XAER      (:,:,:) = XUNDEF
+
+  TPPROFILER%XIWV(:) = XUNDEF
+  TPPROFILER%XZTD(:) = XUNDEF
+  TPPROFILER%XZWD(:) = XUNDEF
+  TPPROFILER%XZHD(:) = XUNDEF
+
+!   IF ( LDIAG_IN_RUN ) THEN
+    TPPROFILER%XT2M   (:) = XUNDEF
+    TPPROFILER%XQ2M   (:) = XUNDEF
+    TPPROFILER%XHU2M  (:) = XUNDEF
+    TPPROFILER%XZON10M(:) = XUNDEF
+    TPPROFILER%XMER10M(:) = XUNDEF
+    TPPROFILER%XRN    (:) = XUNDEF
+    TPPROFILER%XH     (:) = XUNDEF
+    TPPROFILER%XLE    (:) = XUNDEF
+    TPPROFILER%XLEI   (:) = XUNDEF
+    TPPROFILER%XGFLUX (:) = XUNDEF
+    IF ( CRAD /= 'NONE' ) THEN
+      TPPROFILER%XSWD   (:) = XUNDEF
+      TPPROFILER%XSWU   (:) = XUNDEF
+      TPPROFILER%XLWD   (:) = XUNDEF
+      TPPROFILER%XLWU   (:) = XUNDEF
+    END IF
+    TPPROFILER%XTKE_DISS(:,:) = XUNDEF
+!   END IF
+
+END SUBROUTINE PROFILER_ALLOCATE
+
 ! ##############################################
 SUBROUTINE STATION_ALLOCATE( TPSTATION, KSTORE )
 ! ##############################################
@@ -116,9 +242,9 @@ SUBROUTINE STATION_ALLOCATE( TPSTATION, KSTORE )
 
 END SUBROUTINE STATION_ALLOCATE
 
-! ########################################
-SUBROUTINE STATION_INI_INTERP( TPSTATION )
-! ########################################
+! ##########################################
+SUBROUTINE STATPROF_INI_INTERP( TPSTATPROF )
+! ##########################################
 
   USE MODD_GRID,         ONLY: XLATORI, XLONORI
   USE MODD_PARAMETERS,   ONLY: XUNDEF
@@ -128,28 +254,28 @@ SUBROUTINE STATION_INI_INTERP( TPSTATION )
 
   IMPLICIT NONE
 
-  TYPE(TSTATIONDATA), INTENT(INOUT) :: TPSTATION
+  CLASS(TSTATPROFDATA), INTENT(INOUT) :: TPSTATPROF
 
-  IF ( TPSTATION%XLAT == XUNDEF .OR. TPSTATION%XLON == XUNDEF ) THEN
-    CMNHMSG(1) = 'Error in station position'
+  IF ( TPSTATPROF%XLAT == XUNDEF .OR. TPSTATPROF%XLON == XUNDEF ) THEN
+    CMNHMSG(1) = 'Error in station or profiler position'
     CMNHMSG(2) = 'either LATitude or LONgitude segment'
     CMNHMSG(3) = 'or I and J segment'
     CMNHMSG(4) = 'definition is not complete.'
-    CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'STATION_INI_INTERP' )
+    CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'STATPROF_INI_INTERP' )
   END IF
 
-  CALL SM_XYHAT( XLATORI,        XLONORI,        &
-                 TPSTATION%XLAT, TPSTATION%XLON, &
-                 TPSTATION%XX,   TPSTATION%XY    )
+  CALL SM_XYHAT( XLATORI,         XLONORI,         &
+                 TPSTATPROF%XLAT, TPSTATPROF%XLON, &
+                 TPSTATPROF%XX,   TPSTATPROF%XY    )
 
-END SUBROUTINE STATION_INI_INTERP
+END SUBROUTINE STATPROF_INI_INTERP
 
 ! ###############################################################################################
-SUBROUTINE STATION_POSITION( TPSTATION, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,                 &
+SUBROUTINE STATPROF_POSITION( TPSTATPROF, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,               &
                              PXHATM_PHYS_MIN, PXHATM_PHYS_MAX,PYHATM_PHYS_MIN, PYHATM_PHYS_MAX, &
                              OINSIDE, OPRESENT                                                  )
 ! ###############################################################################################
-! Subroutine to determine the position of a station on the model grid
+! Subroutine to determine the position of a station/profiler on the model grid
 ! and set the useful coefficients for data interpolation
 
   USE MODD_CONF,           ONLY: L1D
@@ -162,15 +288,17 @@ SUBROUTINE STATION_POSITION( TPSTATION, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,
 
   IMPLICIT NONE
 
-  TYPE(TSTATIONDATA), INTENT(INOUT) :: TPSTATION
-  REAL, DIMENSION(:), INTENT(IN)    :: PXHAT_GLOB
-  REAL, DIMENSION(:), INTENT(IN)    :: PYHAT_GLOB
-  REAL, DIMENSION(:), INTENT(IN)    :: PXHATM ! mass point coordinates
-  REAL, DIMENSION(:), INTENT(IN)    :: PYHATM ! mass point coordinates
-  REAL,               INTENT(IN)    :: PXHATM_PHYS_MIN, PYHATM_PHYS_MIN  ! Minimum X coordinate of mass points in the physical domain
-  REAL,               INTENT(IN)    :: PXHATM_PHYS_MAX, PYHATM_PHYS_MAX  ! Minimum X coordinate of mass points in the physical domain
-  LOGICAL,            INTENT(OUT)   :: OINSIDE  ! True if station is inside physical domain of model
-  LOGICAL,            INTENT(OUT)   :: OPRESENT ! True if station is present on the current process
+  CLASS(TSTATPROFDATA), INTENT(INOUT) :: TPSTATPROF
+  REAL, DIMENSION(:),   INTENT(IN)    :: PXHAT_GLOB
+  REAL, DIMENSION(:),   INTENT(IN)    :: PYHAT_GLOB
+  REAL, DIMENSION(:),   INTENT(IN)    :: PXHATM ! mass point coordinates
+  REAL, DIMENSION(:),   INTENT(IN)    :: PYHATM ! mass point coordinates
+  REAL,                 INTENT(IN)    :: PXHATM_PHYS_MIN  ! Minimum X coordinate of mass points in the physical domain
+  REAL,                 INTENT(IN)    :: PYHATM_PHYS_MIN  ! Minimum Y coordinate of mass points in the physical domain
+  REAL,                 INTENT(IN)    :: PXHATM_PHYS_MAX  ! Maximum X coordinate of mass points in the physical domain
+  REAL,                 INTENT(IN)    :: PYHATM_PHYS_MAX  ! Minimum Y coordinate of mass points in the physical domain
+  LOGICAL,              INTENT(OUT)   :: OINSIDE  ! True if station/profiler is inside physical domain of model
+  LOGICAL,              INTENT(OUT)   :: OPRESENT ! True if station/profiler is present on the current process
 
   INTEGER :: IIB ! domain sizes of current process
   INTEGER :: IJB !
@@ -185,39 +313,39 @@ SUBROUTINE STATION_POSITION( TPSTATION, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,
 
   CALL GET_INDICE_ll( IIB, IJB, IIE, IJE )
 
-  IF (       TPSTATION%XX >= PXHAT_GLOB(JPHEXT+1) .AND. TPSTATION%XX <= PXHAT_GLOB(UBOUND(PXHAT_GLOB,1)-JPHEXT+1) &
-       .AND. TPSTATION%XY >= PYHAT_GLOB(JPHEXT+1) .AND. TPSTATION%XY <= PYHAT_GLOB(UBOUND(PYHAT_GLOB,1)-JPHEXT+1) ) THEN
+  IF (       TPSTATPROF%XX >= PXHAT_GLOB(JPHEXT+1) .AND. TPSTATPROF%XX <= PXHAT_GLOB(UBOUND(PXHAT_GLOB,1)-JPHEXT+1) &
+       .AND. TPSTATPROF%XY >= PYHAT_GLOB(JPHEXT+1) .AND. TPSTATPROF%XY <= PYHAT_GLOB(UBOUND(PYHAT_GLOB,1)-JPHEXT+1) ) THEN
     OINSIDE = .TRUE.
   ELSE
     CALL GET_MODEL_NUMBER_ll(IMI)
-    WRITE( CMNHMSG(1), "( 'station ', A, ' is outside of physical domain of model', I3 )" ) TRIM(TPSTATION%CNAME), IMI
-    CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'STATION_POSITION' )
+    WRITE( CMNHMSG(1), "( 'station or profiler ', A, ' is outside of physical domain of model', I3 )" ) TRIM(TPSTATPROF%CNAME), IMI
+    CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'STATPROF_POSITION' )
  END IF
 
   ! X position
-  TPSTATION%NI_U = COUNT( XXHAT (:) <= TPSTATION%XX )
-  TPSTATION%NI_M = COUNT( PXHATM(:) <= TPSTATION%XX )
+  TPSTATPROF%NI_U = COUNT( XXHAT (:) <= TPSTATPROF%XX )
+  TPSTATPROF%NI_M = COUNT( PXHATM(:) <= TPSTATPROF%XX )
 
   ! Y position
-  TPSTATION%NJ_V = COUNT( XYHAT (:) <= TPSTATION%XY )
-  TPSTATION%NJ_M = COUNT( PYHATM(:) <= TPSTATION%XY )
+  TPSTATPROF%NJ_V = COUNT( XYHAT (:) <= TPSTATPROF%XY )
+  TPSTATPROF%NJ_M = COUNT( PYHATM(:) <= TPSTATPROF%XY )
 
-  ! Position of station according to process
-  IF (       TPSTATION%NI_U >= IIB .AND. TPSTATION%NI_U <= IIE &
-       .AND. TPSTATION%NJ_V >= IJB .AND. TPSTATION%NJ_V <= IJE ) OPRESENT = .TRUE.
+  ! Position of station/profiler according to process
+  IF (       TPSTATPROF%NI_U >= IIB .AND. TPSTATPROF%NI_U <= IIE &
+       .AND. TPSTATPROF%NJ_V >= IJB .AND. TPSTATPROF%NJ_V <= IJE ) OPRESENT = .TRUE.
   IF ( L1D ) OPRESENT = .TRUE.
 
-  ! Check if station is too near of physical domain border (outside of physical domain for mass points)
+  ! Check if station/profiler is too near of physical domain border (outside of physical domain for mass points)
   IF ( OINSIDE .AND. .NOT. L1D ) THEN
-    IF (      TPSTATION%XX < PXHATM_PHYS_MIN .OR. TPSTATION%XX > PXHATM_PHYS_MAX &
-         .OR. TPSTATION%XY < PYHATM_PHYS_MIN .OR. TPSTATION%XY > PYHATM_PHYS_MAX ) THEN
+    IF (      TPSTATPROF%XX < PXHATM_PHYS_MIN .OR. TPSTATPROF%XX > PXHATM_PHYS_MAX &
+         .OR. TPSTATPROF%XY < PYHATM_PHYS_MIN .OR. TPSTATPROF%XY > PYHATM_PHYS_MAX ) THEN
       CALL GET_MODEL_NUMBER_ll(IMI)
-      WRITE( CMNHMSG(1), "( 'station ', A, ' is outside of mass-points physical domain of model', I3 )" ) &
-             TRIM(TPSTATION%CNAME), IMI
+      WRITE( CMNHMSG(1), "( 'station or profiler ', A, ' is outside of mass-points physical domain of model', I3 )" ) &
+             TRIM(TPSTATPROF%CNAME), IMI
       CMNHMSG(2) = 'but is inside of flux-points physical domain.'
-      CMNHMSG(3) = 'Meaning: station is too close to the boundaries of physical domain.'
-      CMNHMSG(4) = '=> station disabled (not computed)'
-      CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'STATION_POSITION' )
+      CMNHMSG(3) = 'Meaning: station or profiler is too close to the boundaries of physical domain.'
+      CMNHMSG(4) = '=> station or profiler disabled (not computed)'
+      CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'STATPROF_POSITION' )
       OPRESENT = .FALSE.
       OINSIDE = .FALSE.
     END IF
@@ -226,29 +354,83 @@ SUBROUTINE STATION_POSITION( TPSTATION, PXHAT_GLOB, PYHAT_GLOB, PXHATM, PYHATM,
   ! Computations only on correct process
   IF ( OPRESENT .AND. .NOT. L1D ) THEN
     ! Interpolation coefficient for X (mass-point)
-    TPSTATION%XXMCOEF = ( TPSTATION%XX - PXHATM(TPSTATION%NI_M) ) / ( PXHATM(TPSTATION%NI_M+1) - PXHATM(TPSTATION%NI_M) )
+    TPSTATPROF%XXMCOEF = ( TPSTATPROF%XX - PXHATM(TPSTATPROF%NI_M) ) / ( PXHATM(TPSTATPROF%NI_M+1) - PXHATM(TPSTATPROF%NI_M) )
     ! Interpolation coefficient for Y (mass-point)
-    TPSTATION%XYMCOEF = ( TPSTATION%XY - PYHATM(TPSTATION%NJ_M) ) / ( PYHATM(TPSTATION%NJ_M+1) - PYHATM(TPSTATION%NJ_M) )
+    TPSTATPROF%XYMCOEF = ( TPSTATPROF%XY - PYHATM(TPSTATPROF%NJ_M) ) / ( PYHATM(TPSTATPROF%NJ_M+1) - PYHATM(TPSTATPROF%NJ_M) )
     ! Interpolation coefficient for X (U-point)
-    TPSTATION%XXUCOEF = ( TPSTATION%XX - XXHAT(TPSTATION%NI_U) )  / ( XXHAT(TPSTATION%NI_U+1)  - XXHAT(TPSTATION%NI_U) )
+    TPSTATPROF%XXUCOEF = ( TPSTATPROF%XX - XXHAT(TPSTATPROF%NI_U) )  / ( XXHAT(TPSTATPROF%NI_U+1)  - XXHAT(TPSTATPROF%NI_U) )
     ! Interpolation coefficient for Y (V-point)
-    TPSTATION%XYVCOEF = ( TPSTATION%XY - XYHAT(TPSTATION%NJ_V) )  / ( XYHAT(TPSTATION%NJ_V+1)  - XYHAT(TPSTATION%NJ_V) )
+    TPSTATPROF%XYVCOEF = ( TPSTATPROF%XY - XYHAT(TPSTATPROF%NJ_V) )  / ( XYHAT(TPSTATPROF%NJ_V+1)  - XYHAT(TPSTATPROF%NJ_V) )
   END IF
 
   IF ( OPRESENT ) THEN
-    ! The closest K-level to the station altitude is chosen
-    JK = JPVEXT + 1
-    DO WHILE ( ( STATION_INTERP_2D( TPSTATION, XZZ(:,:,JK) ) - STATION_INTERP_2D( TPSTATION, XZZ(:,:,JPVEXT+1) ) ) < TPSTATION%XZ)
-      JK = JK + 1
-    END DO
-    ZLOW  = STATION_INTERP_2D( TPSTATION, XZZ(:,:,JK-1) ) - STATION_INTERP_2D( TPSTATION, XZZ(:,:,JPVEXT+1) )
-    ZHIGH = STATION_INTERP_2D( TPSTATION, XZZ(:,:,JK  ) ) - STATION_INTERP_2D( TPSTATION, XZZ(:,:,JPVEXT+1) )
-    !If the station is nearer from the lower level, select it
-    IF ( ( ZHIGH - TPSTATION%XZ ) > ( TPSTATION%XZ - ZLOW ) ) JK = JK - 1
-    TPSTATION%NK = JK
+    SELECT TYPE( TPSTATPROF )
+      TYPE IS( TPROFILERDATA )
+        ! Nothing to do
+
+      TYPE IS( TSTATIONDATA )
+        ! The closest K-level to the station altitude is chosen
+        JK = JPVEXT + 1
+        DO WHILE ( ( STATPROF_INTERP_2D( TPSTATPROF, XZZ(:,:,JK) ) - STATPROF_INTERP_2D( TPSTATPROF, XZZ(:,:,JPVEXT+1) ) ) &
+                   < TPSTATPROF%XZ)
+          JK = JK + 1
+        END DO
+        ZLOW  = STATPROF_INTERP_2D( TPSTATPROF, XZZ(:,:,JK-1) ) - STATPROF_INTERP_2D( TPSTATPROF, XZZ(:,:,JPVEXT+1) )
+        ZHIGH = STATPROF_INTERP_2D( TPSTATPROF, XZZ(:,:,JK  ) ) - STATPROF_INTERP_2D( TPSTATPROF, XZZ(:,:,JPVEXT+1) )
+        !If the station/profiler is nearer from the lower level, select it
+        IF ( ( ZHIGH - TPSTATPROF%XZ ) > ( TPSTATPROF%XZ - ZLOW ) ) JK = JK - 1
+        TPSTATPROF%NK = JK
+
+      CLASS DEFAULT
+        CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_POSITION', 'unknown type for TPSTATPROF' )
+    END SELECT
   END IF
 
-END SUBROUTINE STATION_POSITION
+END SUBROUTINE STATPROF_POSITION
+
+! ###################################
+SUBROUTINE PROFILER_ADD( TPPROFILER )
+! ###################################
+! Subroutine to add a station to the local list of profilers
+  USE MODD_PROFILER_n, ONLY: NUMBPROFILER_LOC, TPROFILERS
+
+  IMPLICIT NONE
+
+  CLASS(TSTATPROFDATA), INTENT(IN) :: TPPROFILER
+
+  INTEGER :: JS
+  TYPE(TPROFILERDATA), DIMENSION(:), POINTER :: TZPROFILERS
+
+  NUMBPROFILER_LOC = NUMBPROFILER_LOC + 1
+
+  ALLOCATE( TZPROFILERS( NUMBPROFILER_LOC ) )
+  DO JS = 1, NUMBPROFILER_LOC - 1
+    TZPROFILERS(JS) = TPROFILERS(JS)
+  END DO
+
+  !Copy fields available in TSTATPROFDATA
+  !other fields are not yet set
+  TZPROFILERS(NUMBPROFILER_LOC)%CNAME   = TPPROFILER%CNAME
+  TZPROFILERS(NUMBPROFILER_LOC)%NID     = TPPROFILER%NID
+  TZPROFILERS(NUMBPROFILER_LOC)%XX      = TPPROFILER%XX
+  TZPROFILERS(NUMBPROFILER_LOC)%XY      = TPPROFILER%XY
+  TZPROFILERS(NUMBPROFILER_LOC)%XZ      = TPPROFILER%XZ
+  TZPROFILERS(NUMBPROFILER_LOC)%XLON    = TPPROFILER%XLON
+  TZPROFILERS(NUMBPROFILER_LOC)%XLAT    = TPPROFILER%XLAT
+  TZPROFILERS(NUMBPROFILER_LOC)%NI_M    = TPPROFILER%NI_M
+  TZPROFILERS(NUMBPROFILER_LOC)%NJ_M    = TPPROFILER%NJ_M
+  TZPROFILERS(NUMBPROFILER_LOC)%NI_U    = TPPROFILER%NI_U
+  TZPROFILERS(NUMBPROFILER_LOC)%NJ_V    = TPPROFILER%NJ_V
+  TZPROFILERS(NUMBPROFILER_LOC)%XXMCOEF = TPPROFILER%XXMCOEF
+  TZPROFILERS(NUMBPROFILER_LOC)%XYMCOEF = TPPROFILER%XYMCOEF
+  TZPROFILERS(NUMBPROFILER_LOC)%XXUCOEF = TPPROFILER%XXUCOEF
+  TZPROFILERS(NUMBPROFILER_LOC)%XYVCOEF = TPPROFILER%XYVCOEF
+
+  IF ( ASSOCIATED( TPROFILERS ) ) DEALLOCATE( TPROFILERS ) !Can be done without memory leak because allocatable arrays were
+                                                           !not yet allocated (will be done in PROFILER_ALLOCATE)
+  TPROFILERS => TZPROFILERS
+
+END SUBROUTINE PROFILER_ADD
 
 ! #################################
 SUBROUTINE STATION_ADD( TPSTATION )
@@ -258,7 +440,7 @@ SUBROUTINE STATION_ADD( TPSTATION )
 
   IMPLICIT NONE
 
-  TYPE(TSTATIONDATA), INTENT(IN) :: TPSTATION
+  CLASS(TSTATPROFDATA), INTENT(IN) :: TPSTATION
 
   INTEGER :: JS
   TYPE(TSTATIONDATA), DIMENSION(:), POINTER :: TZSTATIONS
@@ -269,23 +451,42 @@ SUBROUTINE STATION_ADD( TPSTATION )
   DO JS = 1, NUMBSTAT_LOC - 1
     TZSTATIONS(JS) = TSTATIONS(JS)
   END DO
-  TZSTATIONS(NUMBSTAT_LOC) = TPSTATION
+
+  !Copy fields available in TSTATPROFDATA
+  !other fields are not yet set
+  TZSTATIONS(NUMBSTAT_LOC)%CNAME   = TPSTATION%CNAME
+  TZSTATIONS(NUMBSTAT_LOC)%NID     = TPSTATION%NID
+  TZSTATIONS(NUMBSTAT_LOC)%XX      = TPSTATION%XX
+  TZSTATIONS(NUMBSTAT_LOC)%XY      = TPSTATION%XY
+  TZSTATIONS(NUMBSTAT_LOC)%XZ      = TPSTATION%XZ
+  TZSTATIONS(NUMBSTAT_LOC)%XLON    = TPSTATION%XLON
+  TZSTATIONS(NUMBSTAT_LOC)%XLAT    = TPSTATION%XLAT
+  TZSTATIONS(NUMBSTAT_LOC)%NI_M    = TPSTATION%NI_M
+  TZSTATIONS(NUMBSTAT_LOC)%NJ_M    = TPSTATION%NJ_M
+  TZSTATIONS(NUMBSTAT_LOC)%NI_U    = TPSTATION%NI_U
+  TZSTATIONS(NUMBSTAT_LOC)%NJ_V    = TPSTATION%NJ_V
+  TZSTATIONS(NUMBSTAT_LOC)%XXMCOEF = TPSTATION%XXMCOEF
+  TZSTATIONS(NUMBSTAT_LOC)%XYMCOEF = TPSTATION%XYMCOEF
+  TZSTATIONS(NUMBSTAT_LOC)%XXUCOEF = TPSTATION%XXUCOEF
+  TZSTATIONS(NUMBSTAT_LOC)%XYVCOEF = TPSTATION%XYVCOEF
+
   IF ( ASSOCIATED( TSTATIONS ) ) DEALLOCATE( TSTATIONS ) !Can be done without memory leak because allocatable arrays were
                                                          !not yet allocated (will be done in STATION_ALLOCATE)
   TSTATIONS => TZSTATIONS
 
 END SUBROUTINE STATION_ADD
-!
-! ######################################################
-FUNCTION STATION_INTERP_2D( TPSTATION, PA ) RESULT( PB )
-! ######################################################
-  USE MODD_CONF, ONLY: L1D
+
+! ########################################################
+FUNCTION STATPROF_INTERP_2D( TPSTATPROF, PA ) RESULT( PB )
+! ########################################################
+  USE MODD_CONF,       ONLY: L1D
+  USE MODD_PARAMETERS, ONLY: XUNDEF
 
   USE MODE_MSG
 
   IMPLICIT NONE
 
-  TYPE(TSTATIONDATA),   INTENT(IN) :: TPSTATION
+  CLASS(TSTATPROFDATA), INTENT(IN) :: TPSTATPROF
   REAL, DIMENSION(:,:), INTENT(IN) :: PA
   REAL                             :: PB
 
@@ -298,32 +499,34 @@ FUNCTION STATION_INTERP_2D( TPSTATION, PA ) RESULT( PB )
     JI = 2
     JJ = 2
   ELSE
-    JI = TPSTATION%NI_M
-    JJ = TPSTATION%NJ_M
+    JI = TPSTATPROF%NI_M
+    JJ = TPSTATPROF%NJ_M
   END IF
 
   IF (       JI >= 1 .AND. JI < SIZE( PA, 1 ) &
        .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
-    PB = (1.-TPSTATION%XXMCOEF) *  (1.-TPSTATION%XYMCOEF) * PA(JI,JJ)    + &
-         (   TPSTATION%XXMCOEF) *  (1.-TPSTATION%XYMCOEF) * PA(JI+1,JJ)  + &
-         (1.-TPSTATION%XXMCOEF) *  (   TPSTATION%XYMCOEF) * PA(JI,JJ+1)  + &
-         (   TPSTATION%XXMCOEF) *  (   TPSTATION%XYMCOEF) * PA(JI+1,JJ+1)
+    PB = (1.-TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI,   JJ  ) + &
+         (   TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI+1, JJ  ) + &
+         (1.-TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI,   JJ+1) + &
+         (   TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI+1, JJ+1)
   ELSE
-    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATION_INTERP_2D', 'value can not be interpolated' )
+    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_INTERP_2D', 'value can not be interpolated' )
+    PB = XUNDEF
   END IF
 
-END FUNCTION STATION_INTERP_2D
+END FUNCTION STATPROF_INTERP_2D
 
-! ########################################################
-FUNCTION STATION_INTERP_2D_U( TPSTATION, PA ) RESULT( PB )
-! ########################################################
-  USE MODD_CONF, ONLY: L1D
+! ##########################################################
+FUNCTION STATPROF_INTERP_2D_U( TPSTATPROF, PA ) RESULT( PB )
+! ##########################################################
+  USE MODD_CONF,       ONLY: L1D
+  USE MODD_PARAMETERS, ONLY: XUNDEF
 
   USE MODE_MSG
 
   IMPLICIT NONE
 
-  TYPE(TSTATIONDATA),   INTENT(IN) :: TPSTATION
+  CLASS(TSTATPROFDATA), INTENT(IN) :: TPSTATPROF
   REAL, DIMENSION(:,:), INTENT(IN) :: PA
   REAL                             :: PB
 
@@ -336,32 +539,34 @@ FUNCTION STATION_INTERP_2D_U( TPSTATION, PA ) RESULT( PB )
     JI = 2
     JJ = 2
   ELSE
-    JI = TPSTATION%NI_U
-    JJ = TPSTATION%NJ_M
+    JI = TPSTATPROF%NI_U
+    JJ = TPSTATPROF%NJ_M
   END IF
 
   IF (       JI >= 1 .AND. JI < SIZE( PA, 1 ) &
        .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
-    PB = (1.-TPSTATION%XXUCOEF) *  (1.-TPSTATION%XYMCOEF) * PA(JI,JJ)    + &
-         (   TPSTATION%XXUCOEF) *  (1.-TPSTATION%XYMCOEF) * PA(JI+1,JJ)  + &
-         (1.-TPSTATION%XXUCOEF) *  (   TPSTATION%XYMCOEF) * PA(JI,JJ+1)  + &
-         (   TPSTATION%XXUCOEF) *  (   TPSTATION%XYMCOEF) * PA(JI+1,JJ+1)
+    PB = (1.-TPSTATPROF%XXUCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI,   JJ  ) + &
+         (   TPSTATPROF%XXUCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI+1, JJ  ) + &
+         (1.-TPSTATPROF%XXUCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI,   JJ+1) + &
+         (   TPSTATPROF%XXUCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI+1, JJ+1)
   ELSE
-    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATION_INTERP_2D_U', 'value can not be interpolated' )
+    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_INTERP_2D_U', 'value can not be interpolated' )
+    PB = XUNDEF
   END IF
 
-END FUNCTION STATION_INTERP_2D_U
+END FUNCTION STATPROF_INTERP_2D_U
 
-! ########################################################
-FUNCTION STATION_INTERP_2D_V( TPSTATION, PA ) RESULT( PB )
-! ########################################################
-  USE MODD_CONF, ONLY: L1D
+! ##########################################################
+FUNCTION STATPROF_INTERP_2D_V( TPSTATPROF, PA ) RESULT( PB )
+! ##########################################################
+  USE MODD_CONF,       ONLY: L1D
+  USE MODD_PARAMETERS, ONLY: XUNDEF
 
   USE MODE_MSG
 
   IMPLICIT NONE
 
-  TYPE(TSTATIONDATA),   INTENT(IN) :: TPSTATION
+  CLASS(TSTATPROFDATA), INTENT(IN) :: TPSTATPROF
   REAL, DIMENSION(:,:), INTENT(IN) :: PA
   REAL                             :: PB
 
@@ -374,20 +579,149 @@ FUNCTION STATION_INTERP_2D_V( TPSTATION, PA ) RESULT( PB )
     JI = 2
     JJ = 2
   ELSE
-    JI = TPSTATION%NI_M
-    JJ = TPSTATION%NJ_V
+    JI = TPSTATPROF%NI_M
+    JJ = TPSTATPROF%NJ_V
   END IF
 
   IF (       JI >= 1 .AND. JI < SIZE( PA, 1 ) &
        .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
-    PB = (1.-TPSTATION%XXMCOEF) *  (1.-TPSTATION%XYVCOEF) * PA(JI,JJ)    + &
-         (   TPSTATION%XXMCOEF) *  (1.-TPSTATION%XYVCOEF) * PA(JI+1,JJ)  + &
-         (1.-TPSTATION%XXMCOEF) *  (   TPSTATION%XYVCOEF) * PA(JI,JJ+1)  + &
-         (   TPSTATION%XXMCOEF) *  (   TPSTATION%XYVCOEF) * PA(JI+1,JJ+1)
+    PB = (1.-TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYVCOEF) * PA(JI,   JJ  ) + &
+         (   TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYVCOEF) * PA(JI+1, JJ  ) + &
+         (1.-TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYVCOEF) * PA(JI,   JJ+1) + &
+         (   TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYVCOEF) * PA(JI+1, JJ+1)
   ELSE
-    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATION_INTERP_2D_V', 'value can not be interpolated' )
+    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_INTERP_2D_V', 'value can not be interpolated' )
+    PB = XUNDEF
   END IF
 
-END FUNCTION STATION_INTERP_2D_V
+END FUNCTION STATPROF_INTERP_2D_V
+
+! ########################################################
+FUNCTION STATPROF_INTERP_3D( TPSTATPROF, PA ) RESULT( PB )
+! ########################################################
+  USE MODD_CONF,       ONLY: L1D
+  USE MODD_PARAMETERS, ONLY: XUNDEF
+
+  USE MODE_MSG
+
+  IMPLICIT NONE
+
+  CLASS(TSTATPROFDATA),   INTENT(IN) :: TPSTATPROF
+  REAL, DIMENSION(:,:,:), INTENT(IN) :: PA
+  REAL, DIMENSION(SIZE(PA,3))        :: PB
+
+  INTEGER :: JI, JJ, JK
+
+  IF ( SIZE( PA, 1 ) == 2 ) THEN
+    JI = 1
+    JJ = 1
+  ELSE IF ( L1D ) THEN
+    JI = 2
+    JJ = 2
+  ELSE
+    JI = TPSTATPROF%NI_M
+    JJ = TPSTATPROF%NJ_M
+  END IF
+
+  IF (       JI >= 1 .AND. JI < SIZE( PA, 1 ) &
+       .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
+    DO JK = 1, SIZE( PA, 3 )
+      IF ( PA(JI, JJ,   JK) /= XUNDEF .AND. PA(JI+1, JJ,   JK) /= XUNDEF .AND. &
+           PA(JI, JJ+1, JK) /= XUNDEF .AND. PA(JI+1, JJ+1, JK) /= XUNDEF       ) THEN
+        PB(JK) = (1.-TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI,   JJ,   JK) + &
+                 (   TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI+1, JJ,   JK) + &
+                 (1.-TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI,   JJ+1, JK) + &
+                 (   TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI+1, JJ+1, JK)
+      ELSE
+        PB(JK) = XUNDEF
+      END IF
+    END DO
+  ELSE
+    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_INTERP_3D', 'value can not be interpolated' )
+    PB(:) = XUNDEF
+  END IF
+
+END FUNCTION STATPROF_INTERP_3D
+
+! ##########################################################
+FUNCTION STATPROF_INTERP_3D_U( TPSTATPROF, PA ) RESULT( PB )
+! ##########################################################
+  USE MODD_CONF,       ONLY: L1D
+  USE MODD_PARAMETERS, ONLY: XUNDEF
+
+  USE MODE_MSG
+
+  IMPLICIT NONE
+
+  CLASS(TSTATPROFDATA),   INTENT(IN) :: TPSTATPROF
+  REAL, DIMENSION(:,:,:), INTENT(IN) :: PA
+  REAL, DIMENSION(SIZE(PA,3))        :: PB
+
+  INTEGER :: JI, JJ
+
+  IF ( SIZE( PA, 1 ) == 2 ) THEN
+    JI = 1
+    JJ = 1
+  ELSE IF ( L1D ) THEN
+    JI = 2
+    JJ = 2
+  ELSE
+    JI = TPSTATPROF%NI_U
+    JJ = TPSTATPROF%NJ_M
+  END IF
+
+  IF (       JI >= 1 .AND. JI < SIZE( PA, 1 ) &
+       .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
+    PB(:) = (1.-TPSTATPROF%XXUCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI,   JJ,   :) + &
+            (   TPSTATPROF%XXUCOEF) *  (1.-TPSTATPROF%XYMCOEF) * PA(JI+1, JJ,   :) + &
+            (1.-TPSTATPROF%XXUCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI,   JJ+1, :) + &
+            (   TPSTATPROF%XXUCOEF) *  (   TPSTATPROF%XYMCOEF) * PA(JI+1, JJ+1, :)
+  ELSE
+    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_INTERP_3D_U', 'value can not be interpolated' )
+    PB = XUNDEF
+  END IF
+
+END FUNCTION STATPROF_INTERP_3D_U
+
+! ##########################################################
+FUNCTION STATPROF_INTERP_3D_V( TPSTATPROF, PA ) RESULT( PB )
+! ##########################################################
+  USE MODD_CONF,       ONLY: L1D
+  USE MODD_PARAMETERS, ONLY: XUNDEF
+
+  USE MODE_MSG
+
+  IMPLICIT NONE
+
+  CLASS(TSTATPROFDATA),   INTENT(IN) :: TPSTATPROF
+  REAL, DIMENSION(:,:,:), INTENT(IN) :: PA
+  REAL, DIMENSION(SIZE(PA,3))        :: PB
+
+  INTEGER :: JI, JJ
+
+  IF ( SIZE( PA, 1 ) == 2 ) THEN
+    JI = 1
+    JJ = 1
+  ELSE IF ( L1D ) THEN
+    JI = 2
+    JJ = 2
+  ELSE
+    JI = TPSTATPROF%NI_M
+    JJ = TPSTATPROF%NJ_V
+  END IF
+
+  IF (       JI >= 1 .AND. JI < SIZE( PA, 1 ) &
+       .AND. JJ >= 1 .AND. JJ < SIZE( PA, 2 ) ) THEN
+    PB(:) = (1.-TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYVCOEF) * PA(JI,   JJ,   :) + &
+            (   TPSTATPROF%XXMCOEF) *  (1.-TPSTATPROF%XYVCOEF) * PA(JI+1, JJ,   :) + &
+            (1.-TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYVCOEF) * PA(JI,   JJ+1, :) + &
+            (   TPSTATPROF%XXMCOEF) *  (   TPSTATPROF%XYVCOEF) * PA(JI+1, JJ+1, :)
+  ELSE
+    CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'STATPROF_INTERP_3D_V', 'value can not be interpolated' )
+    PB = XUNDEF
+  END IF
+
+END FUNCTION STATPROF_INTERP_3D_V
+
 
 END MODULE MODE_STATPROF_TOOLS
diff --git a/src/MNH/write_profilern.f90 b/src/MNH/write_profilern.f90
index cb9dc945a..e6b744fb9 100644
--- a/src/MNH/write_profilern.f90
+++ b/src/MNH/write_profilern.f90
@@ -19,6 +19,7 @@
 !  P. Wautelet 01/09/2021: fix: correct vertical dimension for ALT and W
 !  P. Wautelet 19/11/2021: bugfix in units for LIMA variables
 !  P. Wautelet 04/02/2022: use TSVLIST to manage metadata of scalar variables
+!  P. Wautelet    04/2022: restructure profilers for better performance, reduce memory usage and correct some problems/bugs
 !-----------------------------------------------------------------
 !      ###########################
 MODULE MODE_WRITE_PROFILER_n
@@ -40,19 +41,31 @@ REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: XWORK6   ! contains temporal serie
 
 contains
 !
-!#####################################
-SUBROUTINE WRITE_PROFILER_n(TPDIAFILE)
-!#####################################
+!#######################################
+SUBROUTINE WRITE_PROFILER_n( TPDIAFILE )
+!#######################################
 !
 !
-!****  *WRITE_PROFILER* - write the balloon and aircraft trajectories and records
-!                      in the diachronic file
+!****  *WRITE_PROFILER* - write the profilers records in the diachronic file
 !
 !*      0. DECLARATIONS
 !          ------------
 !
-USE MODD_IO,         ONLY: TFILEDATA
-use MODD_PROFILER_n, only: NUMBPROFILER
+USE MODD_CONF_n,          ONLY: NRR
+USE MODD_DIM_n,           ONLY: NKMAX
+USE MODD_IO,              ONLY: ISNPROC, ISP, TFILEDATA
+USE MODD_DIAG_IN_RUN,     ONLY: LDIAG_IN_RUN
+USE MODD_MPIF
+USE MODD_NSV,             ONLY: NSV
+USE MODD_PARAMETERS,      ONLY: JPVEXT
+USE MODD_PARAM_n,         ONLY: CCLOUD, CRAD, CTURB
+USE MODD_PRECISION,       ONLY: MNHINT_MPI, MNHREAL_MPI
+USE MODD_PROFILER_n,      only: NUMBPROFILER_LOC, TPROFILERS, tprofilers_time
+USE MODD_RADIATIONS_n,    ONLY: NAER
+USE MODD_TYPE_STATPROF,   ONLY: TPROFILERDATA
+!
+USE MODE_MSG
+USE MODE_STATPROF_TOOLS,  ONLY: PROFILER_ALLOCATE
 !
 IMPLICIT NONE
 !
@@ -65,50 +78,280 @@ TYPE(TFILEDATA),  INTENT(IN) :: TPDIAFILE ! diachronic file to write
 !
 !       0.2  declaration of local variables
 !
-INTEGER :: JI
+INTEGER, PARAMETER :: ITAG = 100
+INTEGER :: IERR
+INTEGER :: IKU
+INTEGER :: JP, JS
+INTEGER :: IDX
+INTEGER :: INUMPROF  ! Total number of profilers (for the current model)
+INTEGER :: IPACKSIZE ! Size of the ZPACK buffer
+INTEGER :: IPOS      ! Position in the ZPACK buffer
+INTEGER :: ISTORE
+INTEGER, DIMENSION(:), ALLOCATABLE :: INPROFPRC    ! Array to store the number of profilers per process (for the current model)
+INTEGER, DIMENSION(:), ALLOCATABLE :: IPROFIDS     ! Intermediate array for MPI communication
+INTEGER, DIMENSION(:), ALLOCATABLE :: IPROFPRCRANK ! Array to store the ranks of the processes where the profilers are
+INTEGER, DIMENSION(:), ALLOCATABLE :: IDS          ! Array to store the profiler number to send
+INTEGER, DIMENSION(:), ALLOCATABLE :: IDISP        ! Array to store the displacements for MPI communications
+REAL,    DIMENSION(:), ALLOCATABLE :: ZPACK        ! Buffer to store raw data of a profiler (used for MPI communication)
+TYPE(TPROFILERDATA) :: TZPROFILER
 !
 !----------------------------------------------------------------------------
 !
-DO JI = 1, NUMBPROFILER
-  CALL PROFILER_DIACHRO_n( TPDIAFILE, JI )
-ENDDO
-!
-!----------------------------------------------------------------------------
+IKU = NKMAX + 2 * JPVEXT
+
+ALLOCATE( INPROFPRC(ISNPROC) )
+ALLOCATE( IDS(NUMBPROFILER_LOC) )
+
+!Gather number of profiler present on each process
+CALL MPI_ALLGATHER( NUMBPROFILER_LOC, 1, MNHINT_MPI, INPROFPRC, 1, MNHINT_MPI, TPDIAFILE%NMPICOMM, IERR )
+
+!Store the identification number of local profilers (these numbers are globals)
+DO JS = 1, NUMBPROFILER_LOC
+  IDS(JS) = TPROFILERS(JS)%NID
+END DO
+
+ALLOCATE( IDISP(ISNPROC) )
+IDISP(1) = 0
+DO JP = 2, ISNPROC
+  IDISP(JP) = IDISP(JP-1) + INPROFPRC(JP-1)
+END DO
+
+INUMPROF = SUM( INPROFPRC(:) )
+ALLOCATE( IPROFIDS(INUMPROF) )
+ALLOCATE( IPROFPRCRANK(INUMPROF) )
+
+!Gather the list of all the profilers of all processes
+CALL MPI_ALLGATHERV( IDS(:), NUMBPROFILER_LOC, MNHINT_MPI, IPROFIDS(:), INPROFPRC(:), &
+                     IDISP(:), MNHINT_MPI, TPDIAFILE%NMPICOMM, IERR )
+
+!Store the rank of each process corresponding to a given profiler
+IDX = 1
+IPROFPRCRANK(:) = -1
+DO JP = 1, ISNPROC
+  DO JS = 1, INPROFPRC(JP)
+    IPROFPRCRANK(IPROFIDS(IDX)) = JP
+    IDX = IDX + 1
+  END DO
+END DO
+
+CALL PROFILER_ALLOCATE( TZPROFILER, SIZE( tprofilers_time%tpdates ) )
+
+!Determine the size of the ZPACK buffer used to transfer profiler data in 1 MPI communication
+IF ( ISNPROC > 1 ) THEN
+  ISTORE = SIZE( TPROFILERS_TIME%TPDATES )
+  IPACKSIZE = 6
+  IPACKSIZE = IPACKSIZE + ISTORE * IKU * ( 16 + NRR + NSV + NAER )
+  IF ( CTURB == 'TKEL') IPACKSIZE = IPACKSIZE + ISTORE * IKU !Tke term
+  IF ( CCLOUD == 'ICE3' .OR. CCLOUD == 'ICE4' ) IPACKSIZE = IPACKSIZE + ISTORE * IKU  !CIZ term
+  IPACKSIZE = IPACKSIZE + 4 * ISTORE
+  IF ( LDIAG_IN_RUN ) THEN
+    IPACKSIZE = IPACKSIZE + ISTORE * 10
+    IF ( CRAD /= 'NONE' )  IPACKSIZE = IPACKSIZE + ISTORE * 4
+    IPACKSIZE = IPACKSIZE + ISTORE * IKU !XTKE_DISS term
+  END IF
+
+  ALLOCATE( ZPACK(IPACKSIZE) )
+END IF
+
+IDX = 1
+
+PROFILER: DO JS = 1, INUMPROF
+  IF ( IPROFPRCRANK(JS) == TPDIAFILE%NMASTER_RANK ) THEN
+    !No communication necessary, the profiler data is already on the writer process
+    IF ( ISP == TPDIAFILE%NMASTER_RANK ) THEN
+      TZPROFILER = TPROFILERS(IDX)
+      IDX = IDX + 1
+    END IF
+  ELSE
+    !The profiler data is not on the writer process
+    IF ( ISP == IPROFPRCRANK(JS) ) THEN
+      ! This process has the data and needs to send it to the writer process
+      IPOS = 1
+      ZPACK(IPOS) = TPROFILERS(IDX)%NID;  IPOS = IPOS + 1
+      ZPACK(IPOS) = TPROFILERS(IDX)%XX;   IPOS = IPOS + 1
+      ZPACK(IPOS) = TPROFILERS(IDX)%XY;   IPOS = IPOS + 1
+      ZPACK(IPOS) = TPROFILERS(IDX)%XZ;   IPOS = IPOS + 1
+      ZPACK(IPOS) = TPROFILERS(IDX)%XLON; IPOS = IPOS + 1
+      ZPACK(IPOS) = TPROFILERS(IDX)%XLAT; IPOS = IPOS + 1
+
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XZON(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XMER(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XFF(:,:),  [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XDD(:,:),  [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XW(:,:),   [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XP(:,:),   [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XZZ(:,:),  [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      IF ( CTURB == 'TKEL') THEN
+        ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XTKE(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      END IF
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XTH(:,:),        [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XTHV(:,:),       [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XVISI(:,:),      [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XVISIKUN(:,:),   [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XCRARE(:,:),     [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XCRARE_ATT(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      IF ( CCLOUD == 'ICE3' .OR. CCLOUD == 'ICE4' ) THEN
+        ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XCIZ(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      END IF
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XLWCZ(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XIWCZ(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XRHOD(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+
+      ZPACK(IPOS:IPOS+ISTORE*IKU*NRR-1)  = RESHAPE( TPROFILERS(IDX)%XR(:,:,:),   [ISTORE*IKU*NRR]  )
+      IPOS = IPOS + ISTORE * IKU * NRR
+      ZPACK(IPOS:IPOS+ISTORE*IKU*NSV-1)  = RESHAPE( TPROFILERS(IDX)%XSV(:,:,:),  [ISTORE*IKU*NSV]  )
+      IPOS = IPOS + ISTORE * IKU * NSV
+      ZPACK(IPOS:IPOS+ISTORE*IKU*NAER-1) = RESHAPE( TPROFILERS(IDX)%XAER(:,:,:), [ISTORE*IKU*NAER] )
+      IPOS = IPOS + ISTORE * IKU * NAER
+
+      ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XIWV(:); IPOS = IPOS + ISTORE
+      ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XZTD(:); IPOS = IPOS + ISTORE
+      ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XZWD(:); IPOS = IPOS + ISTORE
+      ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XZHD(:); IPOS = IPOS + ISTORE
+
+      IF ( LDIAG_IN_RUN ) THEN
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XT2M;    IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XQ2M;    IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XHU2M;   IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XZON10M; IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XMER10M; IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XRN;     IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XH;      IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XLE;     IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XGFLUX;  IPOS = IPOS + ISTORE
+        ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XLEI;    IPOS = IPOS + ISTORE
+        IF ( CRAD /= 'NONE' ) THEN
+          ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XSWD;    IPOS = IPOS + ISTORE
+          ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XSWU;    IPOS = IPOS + ISTORE
+          ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XLWD;    IPOS = IPOS + ISTORE
+          ZPACK(IPOS:IPOS+ISTORE-1) = TPROFILERS(IDX)%XLWU;    IPOS = IPOS + ISTORE
+        END IF
+        ZPACK(IPOS:IPOS+ISTORE*IKU-1) = RESHAPE( TPROFILERS(IDX)%XTKE_DISS(:,:), [ISTORE*IKU] ) ; IPOS = IPOS + ISTORE * IKU
+      END IF
+
+      IF ( IPOS-1 /= IPACKSIZE ) call Print_msg( NVERB_WARNING, 'IO', 'WRITE_PROFILER_n', 'IPOS /= IPACKSIZE (sender side)' )
+
+      CALL MPI_SEND( TPROFILERS(IDX)%CNAME, LEN(TPROFILERS(IDX)%CNAME), MPI_CHARACTER, TPDIAFILE%NMASTER_RANK - 1, &
+                     ITAG, TPDIAFILE%NMPICOMM, IERR )
+      CALL MPI_SEND( ZPACK, IPACKSIZE, MNHREAL_MPI, TPDIAFILE%NMASTER_RANK - 1, ITAG, TPDIAFILE%NMPICOMM, IERR )
+
+      IDX = IDX + 1
+
+    ELSE IF ( ISP == TPDIAFILE%NMASTER_RANK ) THEN
+      ! This process is the writer and will receive the profiler data from its owner
+      CALL MPI_RECV( TZPROFILER%CNAME, LEN(TZPROFILER%CNAME), MPI_CHARACTER, &
+                                                    IPROFPRCRANK(JS) - 1, ITAG, TPDIAFILE%NMPICOMM, MPI_STATUS_IGNORE, IERR )
+      CALL MPI_RECV( ZPACK, IPACKSIZE, MNHREAL_MPI, IPROFPRCRANK(JS) - 1, ITAG, TPDIAFILE%NMPICOMM, MPI_STATUS_IGNORE, IERR )
+
+      IPOS = 1
+      TZPROFILER%NID  = NINT( ZPACK(IPOS) ); IPOS = IPOS + 1
+      TZPROFILER%XX   = ZPACK(IPOS);         IPOS = IPOS + 1
+      TZPROFILER%XY   = ZPACK(IPOS);         IPOS = IPOS + 1
+      TZPROFILER%XZ   = ZPACK(IPOS);         IPOS = IPOS + 1
+      TZPROFILER%XLON = ZPACK(IPOS);         IPOS = IPOS + 1
+      TZPROFILER%XLAT = ZPACK(IPOS);         IPOS = IPOS + 1
+
+      TZPROFILER%XZON(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XMER(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XFF(:,:)  = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XDD(:,:)  = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XW(:,:)   = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XP(:,:)   = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XZZ(:,:)  = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      IF ( CTURB == 'TKEL') THEN
+         TZPROFILER%XTKE(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      END IF
+      TZPROFILER%XTH(:,:)        = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XTHV(:,:)       = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XVISI(:,:)      = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XVISIKUN(:,:)   = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XCRARE(:,:)     = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XCRARE_ATT(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      IF ( CCLOUD == 'ICE3' .OR. CCLOUD == 'ICE4' ) THEN
+        TZPROFILER%XCIZ(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      END IF
+      TZPROFILER%XLWCZ(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XIWCZ(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      TZPROFILER%XRHOD(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+
+      TZPROFILER%XR(:,:,:)   = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU*NRR-1),  [ ISTORE, IKU, NRR ] )
+      IPOS = IPOS + ISTORE * IKU * NRR
+      TZPROFILER%XSV(:,:,:)  = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU*NSV-1),  [ ISTORE, IKU, NSV ] )
+      IPOS = IPOS + ISTORE * IKU * NSV
+      TZPROFILER%XAER(:,:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU*NAER-1), [ ISTORE, IKU, NAER ] )
+      IPOS = IPOS + ISTORE * IKU * NAER
+
+      TZPROFILER%XIWV(:) = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+      TZPROFILER%XZTD(:) = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+      TZPROFILER%XZWD(:) = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+      TZPROFILER%XZHD(:) = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+
+      IF ( LDIAG_IN_RUN ) THEN
+        TZPROFILER%XT2M    = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XQ2M    = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XHU2M   = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XZON10M = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XMER10M = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XRN     = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XH      = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XLE     = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XGFLUX  = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        TZPROFILER%XLEI    = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        IF ( CRAD /= 'NONE' ) THEN
+          TZPROFILER%XSWD = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+          TZPROFILER%XSWU = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+          TZPROFILER%XLWD = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+          TZPROFILER%XLWU = ZPACK(IPOS:IPOS+ISTORE-1) ; IPOS = IPOS + ISTORE
+        END IF
+        TZPROFILER%XTKE_DISS(:,:) = RESHAPE( ZPACK(IPOS:IPOS+ISTORE*IKU-1), [ ISTORE, IKU ] ) ; IPOS = IPOS + ISTORE * IKU
+      END IF
+
+      IF ( IPOS-1 /= IPACKSIZE ) call Print_msg( NVERB_WARNING, 'IO', 'WRITE_PROFILER_n', 'IPOS /= IPACKSIZE (receiver side)' )
+    END IF
+  END IF
+
+  CALL PROFILER_DIACHRO_n( TPDIAFILE, TZPROFILER )
+
+END DO PROFILER
+
+
 END SUBROUTINE WRITE_PROFILER_n
-!----------------------------------------------------------------------------
 
 
-!----------------------------------------------------------------------------
-SUBROUTINE PROFILER_DIACHRO_n( TPDIAFILE, KI )
+! ####################################################
+SUBROUTINE PROFILER_DIACHRO_n( TPDIAFILE, TPPROFILER )
+! ####################################################
 
 use modd_budget,          only: NLVL_CATEGORY, NLVL_SUBCATEGORY, NLVL_GROUP, NLVL_SHAPE, NLVL_TIMEAVG, NLVL_NORM, NLVL_MASK, &
                                 tbudiachrometadata
-USE MODD_DIAG_IN_RUN,     only: LDIAG_IN_RUN
-USE MODD_DUST,            ONLY: LDUST, NMODE_DST
 USE MODD_CH_AEROSOL,      ONLY: LORILAM, JPMODE
+USE MODD_CONF_n,          ONLY: NRR
 USE MODD_CST,             ONLY: XRV
+USE MODD_DIAG_IN_RUN,     only: LDIAG_IN_RUN
+USE MODD_DUST,            ONLY: LDUST, NMODE_DST
+USE MODD_DIM_n,           ONLY: NKMAX
 use modd_field,           only: NMNHDIM_LEVEL, NMNHDIM_LEVEL_W, NMNHDIM_PROFILER_TIME, NMNHDIM_PROFILER_PROC, NMNHDIM_UNUSED, &
                                 tfieldmetadata_base, TYPEREAL
 USE MODD_IO,              ONLY: TFILEDATA
 USE MODD_NSV,             ONLY: tsvlist, nsv, nsv_aer, nsv_aerbeg, nsv_aerend, nsv_dst, nsv_dstbeg, nsv_dstend
-USE MODD_PARAMETERS,      ONLY: XUNDEF
-USE MODD_PARAM_n,         ONLY: CRAD
-USE MODD_PROFILER_n,      ONLY: tprofiler
+USE MODD_PARAMETERS,      ONLY: JPVEXT, XUNDEF
+USE MODD_PARAM_n,         ONLY: CCLOUD, CRAD, CTURB
+USE MODD_PROFILER_n
 USE MODD_RADIATIONS_n,    ONLY: NAER
 USE MODD_SALT,            ONLY: LSALT
+USE MODD_TYPE_STATPROF
 !
 USE MODE_AERO_PSD
 USE MODE_DUST_PSD
 use mode_write_diachro,   only: Write_diachro
 !
-TYPE(TFILEDATA),  INTENT(IN) :: TPDIAFILE ! diachronic file to write
-INTEGER,          INTENT(IN) :: KI
+TYPE(TFILEDATA),     INTENT(IN) :: TPDIAFILE ! diachronic file to write
+TYPE(TPROFILERDATA), INTENT(IN) :: TPPROFILER
 !
 !*      0.2  declaration of local variables for diachro
 !
 character(len=NMNHNAMELGTMAX)                        :: yname
 character(len=NUNITLGTMAX)                           :: yunits
-CHARACTER(LEN=:),                        allocatable :: YGROUP   ! group title
 INTEGER                                              :: IKU
 INTEGER                                              :: IPROC    ! number of variables records
 INTEGER                                              :: JPROC
@@ -116,122 +359,119 @@ integer                                              :: jproc_alt, jproc_w
 INTEGER                                              :: JRR      ! loop counter
 INTEGER                                              :: JSV      ! loop counter
 integer                                              :: ji
-integer                                              :: irr !Number of moist variables
+INTEGER                                              :: ISTORE
 REAL, DIMENSION(:,:,:),                  ALLOCATABLE :: ZRHO
-REAL, DIMENSION(:,:,:), TARGET,          ALLOCATABLE :: ZWORK
-REAL, DIMENSION(:,:,:), POINTER                      :: ZDATA
+REAL, DIMENSION(:,:), TARGET,            ALLOCATABLE :: ZWORK
+REAL, DIMENSION(:,:), POINTER                        :: ZDATA
 REAL, DIMENSION(:,:,:,:),                ALLOCATABLE :: ZSV, ZN0, ZSIG, ZRG
 type(tbudiachrometadata)                             :: tzbudiachro
 type(tfieldmetadata_base), dimension(:), allocatable :: tzfields
 !
 !----------------------------------------------------------------------------
-!
-IF (TPROFILER%X(KI)==XUNDEF) RETURN
-IF (TPROFILER%Y(KI)==XUNDEF) RETURN
 
 ZDATA => Null()
 
-IKU = SIZE(TPROFILER%W,2) !Number of vertical levels
+IKU = NKMAX + 2 * JPVEXT !Number of vertical levels
 !
 !IPROC is too large (not a big problem) due to the separation between vertical profiles and point values
-IPROC = 25 + SIZE(TPROFILER%R,4) + SIZE(TPROFILER%SV,4)
+IPROC = 25 + NRR + NSV
 IF (LDIAG_IN_RUN) IPROC = IPROC + 15
 IF (LORILAM) IPROC = IPROC + JPMODE*3
 IF (LDUST) IPROC = IPROC + NMODE_DST*3
 IF (LDUST .OR. LORILAM .OR. LSALT) IPROC=IPROC+NAER
-IF (SIZE(TPROFILER%TKE  )>0) IPROC = IPROC + 1
-!
-ALLOCATE (XWORK6(1,1,IKU,size(tprofiler%tpdates),1,IPROC))
-ALLOCATE (CCOMMENT(IPROC))
-ALLOCATE (CTITLE  (IPROC))
-ALLOCATE (CUNIT   (IPROC))
-!
-YGROUP = TPROFILER%NAME(KI)
+IF ( CTURB == 'TKEL' ) IPROC = IPROC + 1
+
+ISTORE = SIZE( TPROFILERS_TIME%TPDATES )
+
+ALLOCATE ( XWORK6(1, 1, IKU, ISTORE, 1, IPROC) )
+ALLOCATE ( CCOMMENT(IPROC) )
+ALLOCATE ( CTITLE  (IPROC) )
+ALLOCATE ( CUNIT   (IPROC) )
 !
 !----------------------------------------------------------------------------
 !Treat vertical profiles
 jproc = 0
 
-call Add_profile( 'Th',       'Potential temperature',         'K',      tprofiler%th        )
-call Add_profile( 'Thv',      'Virtual Potential temperature', 'K',      tprofiler%thv       )
-call Add_profile( 'VISI',     'Visibility',                    'km',     tprofiler%visi      )
-call Add_profile( 'VISIKUN',  'Visibility Kunkel',             'km',     tprofiler%visikun   )
-call Add_profile( 'RARE',     'Radar reflectivity',            'dBZ',    tprofiler%crare     )
-call Add_profile( 'RAREatt',  'Radar attenuated reflectivity', 'dBZ',    tprofiler%crare_att )
-call Add_profile( 'P',        'Pressure',                      'Pa',     tprofiler%p         )
-call Add_profile( 'ALT',      'Altitude',                      'm',      tprofiler%zz        )
+call Add_profile( 'Th',       'Potential temperature',         'K',      tpprofiler%xth        )
+call Add_profile( 'Thv',      'Virtual Potential temperature', 'K',      tpprofiler%xthv       )
+call Add_profile( 'VISI',     'Visibility',                    'km',     tpprofiler%xvisi      )
+call Add_profile( 'VISIKUN',  'Visibility Kunkel',             'km',     tpprofiler%xvisikun   )
+call Add_profile( 'RARE',     'Radar reflectivity',            'dBZ',    tpprofiler%xcrare     )
+call Add_profile( 'RAREatt',  'Radar attenuated reflectivity', 'dBZ',    tpprofiler%xcrare_att )
+call Add_profile( 'P',        'Pressure',                      'Pa',     tpprofiler%xp         )
+call Add_profile( 'ALT',      'Altitude',                      'm',      tpprofiler%xzz        )
 !Store position of ALT in the field list. Useful because it is not computed on the same Arakawa-grid points
 jproc_alt = jproc
-call Add_profile( 'ZON_WIND', 'Zonal wind',                    'm s-1',  tprofiler%zon       )
-call Add_profile( 'MER_WIND', 'Meridional wind',               'm s-1',  tprofiler%mer       )
-call Add_profile( 'FF',       'Wind intensity',                'm s-1',  tprofiler%ff        )
-call Add_profile( 'DD',       'Wind direction',                'degree', tprofiler%dd        )
-call Add_profile( 'W',        'Air vertical speed',            'm s-1',  tprofiler%w         )
+call Add_profile( 'ZON_WIND', 'Zonal wind',                    'm s-1',  tpprofiler%xzon       )
+call Add_profile( 'MER_WIND', 'Meridional wind',               'm s-1',  tpprofiler%xmer       )
+call Add_profile( 'FF',       'Wind intensity',                'm s-1',  tpprofiler%xff        )
+call Add_profile( 'DD',       'Wind direction',                'degree', tpprofiler%xdd        )
+call Add_profile( 'W',        'Air vertical speed',            'm s-1',  tpprofiler%xw         )
 !Store position of W in the field list. Useful because it is not computed on the same Arakawa-grid points
 jproc_w = jproc
 
 if ( ldiag_in_run ) &
-  call Add_profile( 'TKE_DISS', 'TKE dissipation rate', 'm2 s-2', tprofiler% tke_diss )
+  call Add_profile( 'TKE_DISS', 'TKE dissipation rate', 'm2 s-2', tpprofiler%xtke_diss )
 
-if ( Size( tprofiler%ciz, 1 ) > 0 ) &
-  call Add_profile( 'CIT',      'Ice concentration',    'kg-3',   tprofiler%ciz )
+if ( ccloud == 'ICE3' .or. ccloud == 'ICE4' ) &
+  call Add_profile( 'CIT',      'Ice concentration',    'kg-3',   tpprofiler%xciz )
 
-irr = Size( tprofiler%r )
-if ( irr >= 1 ) call Add_profile( 'Rv', 'Water vapor mixing ratio',        'kg kg-1', tprofiler%r(:,:,:,1) )
-if ( irr >= 2 ) call Add_profile( 'Rc', 'Liquid cloud water mixing ratio', 'kg kg-1', tprofiler%r(:,:,:,2) )
-if ( irr >= 3 ) call Add_profile( 'Rr', 'Rain water mixing ratio',         'kg kg-1', tprofiler%r(:,:,:,3) )
-if ( irr >= 4 ) call Add_profile( 'Ri', 'Ice cloud water mixing ratio',    'kg kg-1', tprofiler%r(:,:,:,4) )
-if ( irr >= 5 ) call Add_profile( 'Rs', 'Snow mixing ratio',               'kg kg-1', tprofiler%r(:,:,:,5) )
-if ( irr >= 6 ) call Add_profile( 'Rg', 'Graupel mixing ratio',            'kg kg-1', tprofiler%r(:,:,:,6) )
-if ( irr >= 7 ) call Add_profile( 'Rh', 'Hail mixing ratio',               'kg kg-1', tprofiler%r(:,:,:,7) )
+if ( nrr >= 1 ) call Add_profile( 'Rv', 'Water vapor mixing ratio',        'kg kg-1', tpprofiler%xr(:,:,1) )
+if ( nrr >= 2 ) call Add_profile( 'Rc', 'Liquid cloud water mixing ratio', 'kg kg-1', tpprofiler%xr(:,:,2) )
+if ( nrr >= 3 ) call Add_profile( 'Rr', 'Rain water mixing ratio',         'kg kg-1', tpprofiler%xr(:,:,3) )
+if ( nrr >= 4 ) call Add_profile( 'Ri', 'Ice cloud water mixing ratio',    'kg kg-1', tpprofiler%xr(:,:,4) )
+if ( nrr >= 5 ) call Add_profile( 'Rs', 'Snow mixing ratio',               'kg kg-1', tpprofiler%xr(:,:,5) )
+if ( nrr >= 6 ) call Add_profile( 'Rg', 'Graupel mixing ratio',            'kg kg-1', tpprofiler%xr(:,:,6) )
+if ( nrr >= 7 ) call Add_profile( 'Rh', 'Hail mixing ratio',               'kg kg-1', tpprofiler%xr(:,:,7) )
 
-call Add_profile( 'Rhod', 'Density of dry air in moist', 'kg m-3', tprofiler%rhod )
-if ( Size( tprofiler%tke, 1 ) > 0 ) &
-  call Add_profile( 'Tke', 'Turbulent kinetic energy', 'm2 s-2', tprofiler%tke )
+call Add_profile( 'Rhod', 'Density of dry air in moist', 'kg m-3', tpprofiler%xrhod )
+if ( cturb == 'TKEL') &
+  call Add_profile( 'Tke', 'Turbulent kinetic energy', 'm2 s-2', tpprofiler%xtke )
 
-if ( Size( tprofiler%sv, 4 ) > 0  ) then
+if ( nsv > 0  ) then
   ! Scalar variables
-  Allocate( zwork, mold = tprofiler%sv(:,:,:,1) )
+  Allocate( zwork, mold = tpprofiler%xsv(:,:,1) )
   do jsv = 1, nsv
     if ( Trim( tsvlist(jsv)%cunits ) == 'ppv' ) then
       yunits = 'ppb'
-      zwork = tprofiler%sv(:,:,:,jsv) * 1.e9 !*1e9 for conversion ppv->ppb
+      zwork = tpprofiler%xsv(:,:,jsv) * 1.e9 !*1e9 for conversion ppv->ppb
       zdata => zwork
     else
       yunits = Trim( tsvlist(jsv)%cunits )
-      zdata => tprofiler%sv(:,:,:,jsv)
+      zwork = tpprofiler%xsv(:,:,jsv)
+      zdata => zwork
     end if
     call Add_profile( tsvlist(jsv)%cmnhname, '', yunits, zdata )
     zdata => Null()
   end do
   Deallocate( zwork )
 
-  IF ( LORILAM .AND. .NOT.(ANY(TPROFILER%P(:,:,KI) == 0.)) ) THEN
-    ALLOCATE (ZSV (1,iku,size(tprofiler%tpdates),NSV_AER))
-    ALLOCATE (ZRHO(1,iku,size(tprofiler%tpdates)))
-    ALLOCATE (ZN0 (1,iku,size(tprofiler%tpdates),JPMODE))
-    ALLOCATE (ZRG (1,iku,size(tprofiler%tpdates),JPMODE))
-    ALLOCATE (ZSIG(1,iku,size(tprofiler%tpdates),JPMODE))
+  IF ( LORILAM .AND. .NOT.(ANY(TPPROFILER%XP(:,:) == 0.)) ) THEN
+    ALLOCATE (ZSV (1,iku,ISTORE,NSV_AER))
+    ALLOCATE (ZRHO(1,iku,ISTORE))
+    ALLOCATE (ZN0 (1,iku,ISTORE,JPMODE))
+    ALLOCATE (ZRG (1,iku,ISTORE,JPMODE))
+    ALLOCATE (ZSIG(1,iku,ISTORE,JPMODE))
     do ji = 1, iku
-      ZSV(1,ji,:,1:NSV_AER) = TPROFILER%SV(:,ji,KI,NSV_AERBEG:NSV_AEREND)
+      ZSV(1,ji,:,1:NSV_AER) = TPPROFILER%XSV(:,ji,NSV_AERBEG:NSV_AEREND)
     end do
-    IF (SIZE(TPROFILER%R,4) >0) THEN
+    IF ( NRR  > 0) THEN
       ZRHO(1,:,:) = 0.
       do ji = 1, iku
-        DO JRR=1,SIZE(TPROFILER%R,4)
-          ZRHO(1,ji,:) = ZRHO(1,ji,:) + TPROFILER%R(:,ji,KI,JRR)
+        DO JRR = 1, NRR
+          ZRHO(1,ji,:) = ZRHO(1,ji,:) + TPPROFILER%XR(:,ji,JRR)
         ENDDO
-        ZRHO(1,ji,:) = TPROFILER%TH(:,ji,KI) * ( 1. + XRV/XRD*TPROFILER%R(:,ji,KI,1) )  &
+        ZRHO(1,ji,:) = TPPROFILER%XTH(:,ji) * ( 1. + XRV/XRD*TPPROFILER%XR(:,ji,1) )  &
                                              / ( 1. + ZRHO(1,ji,:)                )
       end do
     ELSE
       do ji = 1, iku
-        ZRHO(1,ji,:) = TPROFILER%TH(:,ji,KI)
+        ZRHO(1,ji,:) = TPPROFILER%XTH(:,ji)
       end do
     ENDIF
     do ji = 1, iku
-      ZRHO(1,ji,:) =  TPROFILER%P(:,ji,KI) / &
-                      (XRD *ZRHO(1,ji,:) *((TPROFILER%P(:,ji,KI)/XP00)**(XRD/XCPD)) )
+      ZRHO(1,ji,:) =  TPPROFILER%XP(:,ji) / &
+                      (XRD *ZRHO(1,ji,:) *((TPPROFILER%XP(:,ji)/XP00)**(XRD/XCPD)) )
     end do
     CALL PPP2AERO(ZSV,ZRHO, PSIG3D=ZSIG, PRG3D=ZRG, PN3D=ZN0)
     DO JSV=1,JPMODE
@@ -263,32 +503,32 @@ if ( Size( tprofiler%sv, 4 ) > 0  ) then
     DEALLOCATE (ZSV,ZRHO)
     DEALLOCATE (ZN0,ZRG,ZSIG)
   END IF
-  IF ((LDUST).AND. .NOT.(ANY(TPROFILER%P(:,:,KI) == 0.))) THEN
-    ALLOCATE (ZSV (1,iku,size(tprofiler%tpdates),NSV_DST))
-    ALLOCATE (ZRHO(1,iku,size(tprofiler%tpdates)))
-    ALLOCATE (ZN0 (1,iku,size(tprofiler%tpdates),NMODE_DST))
-    ALLOCATE (ZRG (1,iku,size(tprofiler%tpdates),NMODE_DST))
-    ALLOCATE (ZSIG(1,iku,size(tprofiler%tpdates),NMODE_DST))
+  IF ((LDUST).AND. .NOT.(ANY(TPPROFILER%XP(:,:) == 0.))) THEN
+    ALLOCATE (ZSV (1,iku,ISTORE,NSV_DST))
+    ALLOCATE (ZRHO(1,iku,ISTORE))
+    ALLOCATE (ZN0 (1,iku,ISTORE,NMODE_DST))
+    ALLOCATE (ZRG (1,iku,ISTORE,NMODE_DST))
+    ALLOCATE (ZSIG(1,iku,ISTORE,NMODE_DST))
     do ji = 1, iku
-      ZSV(1,ji,:,1:NSV_DST) = TPROFILER%SV(:,ji,KI,NSV_DSTBEG:NSV_DSTEND)
+      ZSV(1,ji,:,1:NSV_DST) = TPPROFILER%XSV(:,ji,NSV_DSTBEG:NSV_DSTEND)
     end do
-    IF (SIZE(TPROFILER%R,4) >0) THEN
+    IF ( NRR > 0 ) THEN
       ZRHO(1,:,:) = 0.
       do ji = 1, iku
-        DO JRR=1,SIZE(TPROFILER%R,4)
-          ZRHO(1,ji,:) = ZRHO(1,ji,:) + TPROFILER%R(:,ji,KI,JRR)
+        DO JRR = 1, NRR
+          ZRHO(1,ji,:) = ZRHO(1,ji,:) + TPPROFILER%XR(:,ji,JRR)
         ENDDO
-        ZRHO(1,ji,:) = TPROFILER%TH(:,ji,KI) * ( 1. + XRV/XRD*TPROFILER%R(:,ji,KI,1) )  &
+        ZRHO(1,ji,:) = TPPROFILER%XTH(:,ji) * ( 1. + XRV/XRD*TPPROFILER%XR(:,ji,1) )  &
                                              / ( 1. + ZRHO(1,ji,:)                )
       end do
     ELSE
       do ji = 1, iku
-        ZRHO(1,ji,:) = TPROFILER%TH(:,ji,KI)
+        ZRHO(1,ji,:) = TPPROFILER%XTH(:,ji)
       end do
     ENDIF
     do ji = 1, iku
-      ZRHO(1,ji,:) =  TPROFILER%P(:,ji,KI) / &
-                     (XRD *ZRHO(1,ji,:) *((TPROFILER%P(:,ji,KI)/XP00)**(XRD/XCPD)) )
+      ZRHO(1,ji,:) =  TPPROFILER%XP(:,ji) / &
+                     (XRD *ZRHO(1,ji,:) *((TPPROFILER%XP(:,ji)/XP00)**(XRD/XCPD)) )
     end do
     CALL PPP2DUST(ZSV,ZRHO, PSIG3D=ZSIG, PRG3D=ZRG, PN3D=ZN0)
     DO JSV=1,NMODE_DST
@@ -323,7 +563,7 @@ if ( Size( tprofiler%sv, 4 ) > 0  ) then
   if ( ldust .or. lorilam .or. lsalt ) then
     do jsv = 1, naer
       Write( yname, '( a, i1 )' ) 'AEREXT', jsv
-      call Add_profile( yname, 'Aerosol Extinction', '1', tprofiler%aer(:,:,:,jsv) )
+      call Add_profile( yname, 'Aerosol Extinction', '1', tpprofiler%xaer(:,:,jsv) )
     end do
   end if
 end if
@@ -356,12 +596,12 @@ tzbudiachro%clevels  (NLVL_SUBCATEGORY) = ''
 tzbudiachro%ccomments(NLVL_SUBCATEGORY) = ''
 
 tzbudiachro%lleveluse(NLVL_GROUP)       = .true.
-tzbudiachro%clevels  (NLVL_GROUP)       = ygroup
-tzbudiachro%ccomments(NLVL_GROUP)       = 'Data at position of profiler ' // Trim( ygroup )
+tzbudiachro%clevels  (NLVL_GROUP)       = tpprofiler%cname
+tzbudiachro%ccomments(NLVL_GROUP)       = 'Data at position of profiler ' // Trim( tpprofiler%cname )
 
 tzbudiachro%lleveluse(NLVL_SHAPE)       = .true.
 tzbudiachro%clevels  (NLVL_SHAPE)       = 'Vertical_profile'
-tzbudiachro%ccomments(NLVL_SHAPE)       = 'Vertical profiles at position of profiler ' // Trim( ygroup )
+tzbudiachro%ccomments(NLVL_SHAPE)       = 'Vertical profiles at position of profiler ' // Trim( tpprofiler%cname )
 
 tzbudiachro%lleveluse(NLVL_TIMEAVG)     = .false.
 tzbudiachro%clevels  (NLVL_TIMEAVG)     = 'Not_time_averaged'
@@ -397,7 +637,7 @@ tzbudiachro%njh        = 1
 tzbudiachro%nkl        = 1
 tzbudiachro%nkh        = iku
 
-call Write_diachro( tpdiafile, tzbudiachro, tzfields, tprofiler%tpdates, xwork6(:,:,:,:,:,:jproc) )
+call Write_diachro( tpdiafile, tzbudiachro, tzfields, tprofilers_time%tpdates, xwork6(:,:,:,:,:,:jproc) )
 
 Deallocate( tzfields )
 Deallocate( xwork6 )
@@ -405,33 +645,33 @@ Deallocate( xwork6 )
 !----------------------------------------------------------------------------
 !Treat point values
 
-ALLOCATE (XWORK6(1,1,1,size(tprofiler%tpdates),1,IPROC))
+ALLOCATE ( XWORK6(1, 1, 1, ISTORE, 1, IPROC) )
 
 jproc = 0
 
 if ( ldiag_in_run ) then
-  call Add_point( 'T2m',    '2-m temperature',               'K',       tprofiler%t2m    )
-  call Add_point( 'Q2m',    '2-m humidity',                  'kg kg-1', tprofiler%q2m    )
-  call Add_point( 'HU2m',   '2-m relative humidity',         'percent', tprofiler%hu2m   )
-  call Add_point( 'zon10m', '10-m zonal wind',               'm s-1',   tprofiler%zon10m )
-  call Add_point( 'mer10m', '10-m meridian wind',            'm s-1',   tprofiler%mer10m )
-  call Add_point( 'RN',     'Net radiation',                 'W m-2',   tprofiler%rn     )
-  call Add_point( 'H',      'Sensible heat flux',            'W m-2',   tprofiler%h      )
-  call Add_point( 'LE',     'Total Latent heat flux',        'W m-2',   tprofiler%le     )
-  call Add_point( 'G',      'Storage heat flux',             'W m-2',   tprofiler%gflux  )
+  call Add_point( 'T2m',    '2-m temperature',               'K',       tpprofiler%xt2m    )
+  call Add_point( 'Q2m',    '2-m humidity',                  'kg kg-1', tpprofiler%xq2m    )
+  call Add_point( 'HU2m',   '2-m relative humidity',         'percent', tpprofiler%xhu2m   )
+  call Add_point( 'zon10m', '10-m zonal wind',               'm s-1',   tpprofiler%xzon10m )
+  call Add_point( 'mer10m', '10-m meridian wind',            'm s-1',   tpprofiler%xmer10m )
+  call Add_point( 'RN',     'Net radiation',                 'W m-2',   tpprofiler%xrn     )
+  call Add_point( 'H',      'Sensible heat flux',            'W m-2',   tpprofiler%xh      )
+  call Add_point( 'LE',     'Total Latent heat flux',        'W m-2',   tpprofiler%xle     )
+  call Add_point( 'G',      'Storage heat flux',             'W m-2',   tpprofiler%xgflux  )
   if ( crad /= 'NONE' ) then
-    call Add_point( 'SWD',  'Downward short-wave radiation', 'W m-2',   tprofiler%swd    )
-    call Add_point( 'SWU',  'Upward short-wave radiation',   'W m-2',   tprofiler%swu    )
-    call Add_point( 'LWD',  'Downward long-wave radiation',  'W m-2',   tprofiler%lwd    )
-    call Add_point( 'LWU',  'Upward long-wave radiation',    'W m-2',   tprofiler%lwu    )
+    call Add_point( 'SWD',  'Downward short-wave radiation', 'W m-2',   tpprofiler%xswd    )
+    call Add_point( 'SWU',  'Upward short-wave radiation',   'W m-2',   tpprofiler%xswu    )
+    call Add_point( 'LWD',  'Downward long-wave radiation',  'W m-2',   tpprofiler%xlwd    )
+    call Add_point( 'LWU',  'Upward long-wave radiation',    'W m-2',   tpprofiler%xlwu    )
   end if
-  call Add_point( 'LEI',    'Solid Latent heat flux',        'W m-2',   tprofiler%lei    )
+  call Add_point( 'LEI',    'Solid Latent heat flux',        'W m-2',   tpprofiler%xlei    )
 end if
 
-call Add_point( 'IWV', 'Integrated Water Vapour',   'kg m-2', tprofiler%iwv )
-call Add_point( 'ZTD', 'Zenith Tropospheric Delay', 'm',      tprofiler%ztd )
-call Add_point( 'ZWD', 'Zenith Wet Delay',          'm',      tprofiler%zwd )
-call Add_point( 'ZHD', 'Zenith Hydrostatic Delay',  'm',      tprofiler%zhd )
+call Add_point( 'IWV', 'Integrated Water Vapour',   'kg m-2', tpprofiler%xiwv )
+call Add_point( 'ZTD', 'Zenith Tropospheric Delay', 'm',      tpprofiler%xztd )
+call Add_point( 'ZWD', 'Zenith Wet Delay',          'm',      tpprofiler%xzwd )
+call Add_point( 'ZHD', 'Zenith Hydrostatic Delay',  'm',      tpprofiler%xzhd )
 
 Allocate( tzfields( jproc ) )
 
@@ -459,12 +699,12 @@ tzbudiachro%clevels  (NLVL_SUBCATEGORY) = ''
 tzbudiachro%ccomments(NLVL_SUBCATEGORY) = ''
 
 tzbudiachro%lleveluse(NLVL_GROUP)       = .true.
-tzbudiachro%clevels  (NLVL_GROUP)       = ygroup
-tzbudiachro%ccomments(NLVL_GROUP)       = 'Data at position of profiler ' // Trim( ygroup )
+tzbudiachro%clevels  (NLVL_GROUP)       = tpprofiler%cname
+tzbudiachro%ccomments(NLVL_GROUP)       = 'Data at position of profiler ' // Trim( tpprofiler%cname )
 
 tzbudiachro%lleveluse(NLVL_SHAPE)       = .true.
 tzbudiachro%clevels  (NLVL_SHAPE)       = 'Point'
-tzbudiachro%ccomments(NLVL_SHAPE)       = 'Values at position of profiler ' // Trim( ygroup )
+tzbudiachro%ccomments(NLVL_SHAPE)       = 'Values at position of profiler ' // Trim( tpprofiler%cname )
 
 tzbudiachro%lleveluse(NLVL_TIMEAVG)     = .false.
 tzbudiachro%clevels  (NLVL_TIMEAVG)     = 'Not_time_averaged'
@@ -496,7 +736,7 @@ tzbudiachro%njh        = 1
 tzbudiachro%nkl        = 1
 tzbudiachro%nkh        = 1
 
-call Write_diachro( tpdiafile, tzbudiachro, tzfields, tprofiler%tpdates, xwork6(:,:,:,:,:,:jproc) )
+call Write_diachro( tpdiafile, tzbudiachro, tzfields, tprofilers_time%tpdates, xwork6(:,:,:,:,:,:jproc) )
 
 Deallocate( tzfields )
 
@@ -510,13 +750,13 @@ JPROC = JPROC + 1
 CTITLE   (JPROC) = 'LON'
 CUNIT    (JPROC) = 'degree'
 CCOMMENT (JPROC) = 'Longitude'
-XWORK6 (1,1,1,:,1,JPROC) = TPROFILER%LON(KI)
+XWORK6 (1,1,1,:,1,JPROC) = TPPROFILER%XLON
 
 JPROC = JPROC + 1
 CTITLE   (JPROC) = 'LAT'
 CUNIT    (JPROC) = 'degree'
 CCOMMENT (JPROC) = 'Latitude'
-XWORK6 (1,1,1,:,1,JPROC) = TPROFILER%LAT(KI)
+XWORK6 (1,1,1,:,1,JPROC) = TPPROFILER%XLAT
 
 Allocate( tzfields( jproc ) )
 
@@ -535,7 +775,7 @@ tzfields(:)%ndimlist(4) = NMNHDIM_UNUSED
 tzfields(:)%ndimlist(5) = NMNHDIM_UNUSED
 tzfields(:)%ndimlist(6) = NMNHDIM_PROFILER_PROC
 
-call Write_diachro( tpdiafile, tzbudiachro, tzfields, tprofiler%tpdates, xwork6(:,:,:,:,:,:jproc) )
+call Write_diachro( tpdiafile, tzbudiachro, tzfields, tprofilers_time%tpdates, xwork6(:,:,:,:,:,:jproc) )
 
 
 !Necessary because global variables (private inside module)
@@ -544,7 +784,6 @@ Deallocate (ccomment)
 Deallocate (ctitle  )
 Deallocate (cunit   )
 
-
 contains
 
 
@@ -552,23 +791,23 @@ subroutine Add_profile( htitle, hcomment, hunits, pfield )
 
 use mode_msg
 
-character(len=*),       intent(in) :: htitle
-character(len=*),       intent(in) :: hcomment
-character(len=*),       intent(in) :: hunits
-real, dimension(:,:,:), intent(in) :: pfield
+character(len=*),     intent(in) :: htitle
+character(len=*),     intent(in) :: hcomment
+character(len=*),     intent(in) :: hunits
+real, dimension(:,:), intent(in) :: pfield
 
 integer :: jk
 
 jproc = jproc + 1
 
-if ( jproc > iproc ) call Print_msg( NVERB_FATAL, 'IO', 'Add_profile', 'more profiles than expected' )
+if ( jproc > iproc ) call Print_msg( NVERB_FATAL, 'IO', 'Add_profile', 'more processes than expected' )
 
 ctitle(jproc)   = Trim( htitle )
 ccomment(jproc) = Trim( hcomment )
 cunit(jproc)    = Trim( hunits )
 
 do jk = 1, iku
-  xwork6(1, 1, jk, :, 1, jproc) = pfield(:, jk, ki)
+  xwork6(1, 1, jk, :, 1, jproc) = pfield(:, jk)
 end do
 
 end subroutine Add_profile
@@ -578,22 +817,22 @@ subroutine Add_point( htitle, hcomment, hunits, pfield )
 
 use mode_msg
 
-character(len=*),     intent(in) :: htitle
-character(len=*),     intent(in) :: hcomment
-character(len=*),     intent(in) :: hunits
-real, dimension(:,:), intent(in) :: pfield
+character(len=*),   intent(in) :: htitle
+character(len=*),   intent(in) :: hcomment
+character(len=*),   intent(in) :: hunits
+real, dimension(:), intent(in) :: pfield
 
 integer :: jk
 
 jproc = jproc + 1
 
-if ( jproc > iproc ) call Print_msg( NVERB_FATAL, 'IO', 'Add_profile', 'more profiles than expected' )
+if ( jproc > iproc ) call Print_msg( NVERB_FATAL, 'IO', 'Add_point', 'more processes than expected' )
 
 ctitle(jproc)   = Trim( htitle)
 ccomment(jproc) = Trim( hcomment )
 cunit(jproc)    = Trim( hunits )
 
-xwork6(1, 1, 1, :, 1, jproc) = pfield(:, ki)
+xwork6(1, 1, 1, :, 1, jproc) = pfield(:)
 
 end subroutine Add_point
 
diff --git a/src/MNH/write_stationn.f90 b/src/MNH/write_stationn.f90
index 7b35ea4c8..4697b0e87 100644
--- a/src/MNH/write_stationn.f90
+++ b/src/MNH/write_stationn.f90
@@ -3,6 +3,17 @@
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
+! Author:
+!  P. Tulet    15/02/2002
+!
+!  Modifications
+!  P. Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
+!  P. Wautelet 09/10/2020: Write_diachro: use new datatype tpfields
+!  P. Wautelet 03/03/2021: budgets: add tbudiachrometadata type (useful to pass more information to Write_diachro)
+!  P. Wautelet 04/02/2022: use TSVLIST to manage metadata of scalar variables
+!  P. Wautelet    04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
+! --------------------------------------------------------------------------
 !      ###########################
 MODULE MODE_WRITE_STATION_n
 !      ###########################
@@ -28,43 +39,7 @@ SUBROUTINE WRITE_STATION_n( TPDIAFILE )
 ! #####################################
 !
 !
-!!****  *WRITE_STATION* - write the balloon and aircraft trajectories and records
-!!                      in the diachronic file
-!!
-!!    PURPOSE
-!!    -------
-!
-!
-!!**  METHOD
-!!    ------
-!!    
-!!
-!!
-!!
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!    AUTHOR
-!!    ------
-!!      Pierre TULET             * Meteo-France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!     Original 15/02/2002
-!  P. Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
-!  P. Wautelet 09/10/2020: Write_diachro: use new datatype tpfields
-!  P. Wautelet 03/03/2021: budgets: add tbudiachrometadata type (useful to pass more information to Write_diachro)
-!  P. Wautelet 04/02/2022: use TSVLIST to manage metadata of scalar variables
-!  P. Wautelet    04/2022: restructure stations for better performance, reduce memory usage and correct some problems/bugs
-! --------------------------------------------------------------------------
+!****  *WRITE_STATION* - write the stations records in the diachronic file
 !
 !*      0. DECLARATIONS
 !          ------------
@@ -79,6 +54,7 @@ USE MODD_PRECISION,       ONLY: MNHINT_MPI, MNHREAL_MPI
 USE MODD_STATION_n,       only: NUMBSTAT_LOC, TSTATIONS, tstations_time
 USE MODD_TYPE_STATPROF,   ONLY: TSTATIONDATA
 !
+USE MODE_MSG
 USE MODE_STATPROF_TOOLS,  ONLY: STATION_ALLOCATE
 !
 IMPLICIT NONE
@@ -100,7 +76,7 @@ INTEGER :: INUMSTAT  ! Total number of stations (for the current model)
 INTEGER :: IPACKSIZE ! Size of the ZPACK buffer
 INTEGER :: IPOS      ! Position in the ZPACK buffer
 INTEGER :: ISTORE
-INTEGER, DIMENSION(:), ALLOCATABLE :: INSTATPRC    ! Array to store the number of stations per process  (for the current model)
+INTEGER, DIMENSION(:), ALLOCATABLE :: INSTATPRC    ! Array to store the number of stations per process (for the current model)
 INTEGER, DIMENSION(:), ALLOCATABLE :: ISTATIDS     ! Intermediate array for MPI communication
 INTEGER, DIMENSION(:), ALLOCATABLE :: ISTATPRCRANK ! Array to store the ranks of the processes where the stations are
 INTEGER, DIMENSION(:), ALLOCATABLE :: IDS          ! Array to store the station number to send
@@ -222,6 +198,8 @@ STATION: DO JS = 1, INUMSTAT
         ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XSFCO2;    IPOS = IPOS + ISTORE
       END IF
 
+      IF ( IPOS /= IPACKSIZE ) call Print_msg( NVERB_WARNING, 'IO', 'WRITE_STATION_n', 'IPOS /= IPACKSIZE (sender side)' )
+
       CALL MPI_SEND( TSTATIONS(IDX)%CNAME, LEN(TSTATIONS(IDX)%CNAME), MPI_CHARACTER, TPDIAFILE%NMASTER_RANK - 1, &
                      ITAG, TPDIAFILE%NMPICOMM, IERR )
       CALL MPI_SEND( ZPACK, IPACKSIZE, MNHREAL_MPI, TPDIAFILE%NMASTER_RANK - 1, ITAG, TPDIAFILE%NMPICOMM, IERR )
@@ -279,6 +257,8 @@ STATION: DO JS = 1, INUMSTAT
         END IF
         TZSTATION%XSFCO2 =    ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
       END IF
+
+      IF ( IPOS /= IPACKSIZE ) call Print_msg( NVERB_WARNING, 'IO', 'WRITE_STATION_n', 'IPOS /= IPACKSIZE (receiver side)' )
     END IF
   END IF
 
-- 
GitLab