diff --git a/src/MNH/write_stationn.f90 b/src/MNH/write_stationn.f90
index 8ff4c2fc24bc0fdcb81e9654fb023027168ada25..23d6b40e70cde2fdc11a7b8241da9777f4529660 100644
--- a/src/MNH/write_stationn.f90
+++ b/src/MNH/write_stationn.f90
@@ -7,17 +7,25 @@
 MODULE MODE_WRITE_STATION_n
 !      ###########################
 
+use modd_parameters, only: NCOMMENTLGTMAX, NMNHNAMELGTMAX, NUNITLGTMAX
+
 implicit none
 
 private
 
 public :: WRITE_STATION_n
 
+CHARACTER(LEN=NCOMMENTLGTMAX), DIMENSION(:), ALLOCATABLE :: CCOMMENT ! comment string
+CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:), ALLOCATABLE :: CTITLE   ! title
+CHARACTER(LEN=NUNITLGTMAX),    DIMENSION(:), ALLOCATABLE :: CUNIT    ! physical unit
+
+REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: XWORK6   ! contains temporal serie
+
 contains
 !
-!     ##########################################
-      SUBROUTINE WRITE_STATION_n(TPDIAFILE)
-!     ##########################################
+! #####################################
+SUBROUTINE WRITE_STATION_n( TPDIAFILE )
+! #####################################
 !
 !
 !!****  *WRITE_STATION* - write the balloon and aircraft trajectories and records
@@ -62,16 +70,11 @@ contains
 !          ------------
 !
 USE MODD_ALLSTATION_n,    ONLY: LDIAG_SURFRAD
-use MODD_BUDGET,          ONLY: tbudiachrometadata
-USE MODD_CONF,            ONLY: LCARTESIAN
 USE MODD_CONF_n,          ONLY: NRR
-USE MODD_CST,             ONLY: XRV
 USE MODD_IO,              ONLY: ISNPROC, ISP, TFILEDATA
 USE MODD_MPIF
-USE MODD_NSV,             ONLY: tsvlist, nsv, nsv_aer, nsv_aerbeg, nsv_aerend, &
-                                nsv_dst, nsv_dstbeg, nsv_dstend, nsv_slt, nsv_sltbeg, nsv_sltend
+USE MODD_NSV,             ONLY: nsv
 USE MODD_PARAM_n,         ONLY: CRAD, CSURF, CTURB
-USE MODD_PARAMETERS,      ONLY: XUNDEF
 USE MODD_PRECISION,       ONLY: MNHINT_MPI, MNHREAL_MPI
 USE MODD_STATION_n,       only: NUMBSTAT_LOC, TSTATIONS, tstations_time
 USE MODD_TYPE_STATION,    ONLY: TSTATIONDATA
@@ -151,7 +154,7 @@ IF ( ISNPROC > 1 ) THEN
   IPACKSIZE = IPACKSIZE + ISTORE * ( 5 + NRR + NSV )
   IF ( CTURB == 'TKEL') IPACKSIZE = IPACKSIZE + ISTORE !Tke term
   IF ( CRAD /= 'NONE' ) IPACKSIZE = IPACKSIZE + ISTORE !XTSRAD term
-  IF (LDIAG_SURFRAD) THEN
+  IF ( LDIAG_SURFRAD ) THEN
     IF ( CSURF == 'EXTE' ) IPACKSIZE = IPACKSIZE + ISTORE * 10
     IF ( CRAD /= 'NONE' )  IPACKSIZE = IPACKSIZE + ISTORE * 7
     IPACKSIZE = IPACKSIZE + ISTORE !XSFCO2 term
@@ -194,7 +197,7 @@ STATION: DO JS = 1, INUMSTAT
       IF ( CRAD /= 'NONE' ) THEN
         ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XTSRAD(:); IPOS = IPOS + ISTORE
       END IF
-      IF (LDIAG_SURFRAD) THEN
+      IF ( LDIAG_SURFRAD ) THEN
         IF ( CSURF == 'EXTE') THEN
           ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XT2M;    IPOS = IPOS + ISTORE
           ZPACK(IPOS:IPOS+ISTORE-1) = TSTATIONS(IDX)%XQ2M;    IPOS = IPOS + ISTORE
@@ -252,7 +255,7 @@ STATION: DO JS = 1, INUMSTAT
       IF ( CRAD /= 'NONE' ) THEN
         TZSTATION%XTSRAD(:) = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
       END IF
-      IF (LDIAG_SURFRAD) THEN
+      IF ( LDIAG_SURFRAD ) THEN
         IF ( CSURF == 'EXTE' ) THEN
           TZSTATION%XT2M    = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
           TZSTATION%XQ2M    = ZPACK(IPOS:IPOS+ISTORE-1); IPOS = IPOS + ISTORE
@@ -279,25 +282,28 @@ STATION: DO JS = 1, INUMSTAT
     END IF
   END IF
 
-  CALL STATION_DIACHRO_n( TZSTATION )
+  CALL STATION_DIACHRO_n( TPDIAFILE, TZSTATION )
 
 END DO STATION
-!
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-!
-CONTAINS
-!
-!----------------------------------------------------------------------------
-!----------------------------------------------------------------------------
-!
-!----------------------------------------------------------------------------
-SUBROUTINE STATION_DIACHRO_n( TPSTATION )
 
-use modd_budget, only: NLVL_CATEGORY, NLVL_SUBCATEGORY, NLVL_GROUP, NLVL_SHAPE, NLVL_TIMEAVG, NLVL_NORM, NLVL_MASK
+END SUBROUTINE WRITE_STATION_n
+
+! ##################################################
+SUBROUTINE STATION_DIACHRO_n( TPDIAFILE, TPSTATION )
+! ##################################################
+
+USE MODD_ALLSTATION_n,    ONLY: LDIAG_SURFRAD
+use modd_budget, only: NLVL_CATEGORY, NLVL_SUBCATEGORY, NLVL_GROUP, NLVL_SHAPE, NLVL_TIMEAVG, NLVL_NORM, NLVL_MASK, &
+                       tbudiachrometadata
+USE MODD_CONF,            ONLY: LCARTESIAN
+USE MODD_CST,             ONLY: XRV
 use modd_field,  only: NMNHDIM_STATION_TIME, NMNHDIM_STATION_PROC, NMNHDIM_UNUSED, &
                        tfieldmetadata_base, TYPEREAL
-use modd_parameters,   only: NCOMMENTLGTMAX, NMNHNAMELGTMAX, NSTATIONNAMELGTMAX, NUNITLGTMAX
+USE MODD_IO,              ONLY:   TFILEDATA
+USE MODD_NSV,             ONLY: nsv, nsv_aer, nsv_aerbeg, nsv_aerend, &
+                                nsv_dst, nsv_dstbeg, nsv_dstend, nsv_slt, nsv_sltbeg, nsv_sltend, &
+                                tsvlist
+USE MODD_PARAM_n,         ONLY: CRAD, CSURF, CTURB
 use modd_station_n,    only: tstations_time
 use modd_type_station, only: tstationdata
 
@@ -306,25 +312,23 @@ USE MODE_DUST_PSD
 USE MODE_SALT_PSD
 use MODE_WRITE_DIACHRO,   ONLY: Write_diachro
 
-TYPE(TSTATIONDATA),   INTENT(IN)       :: TPSTATION
+TYPE(TFILEDATA),    INTENT(IN) :: TPDIAFILE ! diachronic file to write
+TYPE(TSTATIONDATA), INTENT(IN) :: TPSTATION
 !
 !*      0.2  declaration of local variables for diachro
 !
-REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ZWORK6 ! contains temporal series
-REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ZW6    ! contains temporal series to write
+REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: XWORK6 ! contains temporal series
 REAL, DIMENSION(:,:,:,:),     ALLOCATABLE :: ZSV, ZN0, ZSIG, ZRG
-REAL, DIMENSION(:,:,:,:,:),     ALLOCATABLE :: ZPTOTA
+REAL, DIMENSION(:,:,:,:,:),   ALLOCATABLE :: ZPTOTA
 REAL, DIMENSION(:,:,:),       ALLOCATABLE :: ZRHO
 !
-INTEGER, DIMENSION(:),            ALLOCATABLE :: IGRID    ! grid indicator
-CHARACTER(LEN=NSTATIONNAMELGTMAX)             :: YGROUP   ! group title
-CHARACTER(LEN=NCOMMENTLGTMAX), DIMENSION(:), ALLOCATABLE :: YCOMMENT ! comment string
-CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:), ALLOCATABLE :: YTITLE   ! title
-CHARACTER(LEN=NUNITLGTMAX),    DIMENSION(:), ALLOCATABLE :: YUNIT    ! physical unit
+CHARACTER(LEN=NCOMMENTLGTMAX)     :: YCOMMENT ! comment string
+CHARACTER(LEN=NMNHNAMELGTMAX)     :: YTITLE   ! title
 !
 !!! do not forget to increment the IPROC value if you add diagnostic !!!
 INTEGER :: IPROC    ! number of variables records
 !!! do not forget to increment the JPROC value if you add diagnostic !!!
+INTEGER :: ISTORE
 INTEGER :: JPROC    ! loop counter
 INTEGER :: JRR      ! loop counter
 INTEGER :: JSV      ! loop counter
@@ -339,271 +343,106 @@ IF ( CTURB == 'TKEL' ) IPROC = IPROC + 1
 IF (LDIAG_SURFRAD) THEN
   IF(CSURF=="EXTE") IPROC = IPROC + 10
   IF(CRAD/="NONE")  IPROC = IPROC + 7
+  IPROC = IPROC + 1 ! XSFCO2 term
 END IF
 IF (LORILAM) IPROC = IPROC + JPMODE*(3+NSOA+NCARB+NSP)
 IF (LDUST) IPROC = IPROC + NMODE_DST*3
 IF (LSALT) IPROC = IPROC + NMODE_SLT*3
 IF ( CRAD /= 'NONE' )  IPROC = IPROC + 1
-IPROC = IPROC + 1 ! XSFCO2 term
 !
-ALLOCATE (ZWORK6(1,1,1,SIZE(tstations_time%tpdates),1,IPROC))
-ALLOCATE (YCOMMENT(IPROC))
-ALLOCATE (YTITLE  (IPROC))
-ALLOCATE (YUNIT   (IPROC))
-ALLOCATE (IGRID   (IPROC))
+ISTORE = SIZE( TSTATIONS_TIME%TPDATES )
+
+ALLOCATE( XWORK6(1, 1, 1, ISTORE, 1, IPROC) )
+ALLOCATE( CCOMMENT(IPROC) )
+ALLOCATE( CTITLE  (IPROC) )
+ALLOCATE( CUNIT   (IPROC) )
 !
-IGRID  = 1
-YGROUP = TPSTATION%CNAME
 JPROC = 0
 !
 !----------------------------------------------------------------------------
 !
-JPROC = JPROC + 1
-YTITLE   (JPROC) = 'ZS'
-YUNIT    (JPROC) = 'm'
-YCOMMENT (JPROC) = 'Orography'
-ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XZS
-!
-JPROC = JPROC + 1
-YTITLE   (JPROC) = 'P'
-YUNIT    (JPROC) = 'Pa'
-YCOMMENT (JPROC) = 'Pressure'
-ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XP(:)
-!
-!JPROC = JPROC + 1
-!YTITLE   (JPROC) = 'Z'
-!YUNIT    (JPROC) = 'm'
-!YCOMMENT (JPROC) = 'Z Pos'
-!ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XZ
-!
-IF (LCARTESIAN) THEN
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'X'
-  YUNIT    (JPROC) = 'm'
-  YCOMMENT (JPROC) = 'X Pos'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XX
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'Y'
-  YUNIT    (JPROC) = 'm'
-  YCOMMENT (JPROC) = 'Y Pos'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XY
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'U'
-  YUNIT    (JPROC) = 'm s-1'
-  YCOMMENT (JPROC) = 'Axial velocity'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XZON(:)
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'V'
-  YUNIT    (JPROC) = 'm s-1'
-  YCOMMENT (JPROC) = 'Transversal velocity'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XMER(:)
-ELSE
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'LON'
-  YUNIT    (JPROC) = 'degree'
-  YCOMMENT (JPROC) = 'Longitude'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XLON
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'LAT'
-  YUNIT    (JPROC) = 'degree'
-  YCOMMENT (JPROC) = 'Latitude'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XLAT
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'ZON_WIND'
-  YUNIT    (JPROC) = 'm s-1'
-  YCOMMENT (JPROC) = 'Zonal wind'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XZON(:)
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'MER_WIND'
-  YUNIT    (JPROC) = 'm s-1'
-  YCOMMENT (JPROC) = 'Meridional wind'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XMER(:)
-ENDIF
-!
-JPROC = JPROC + 1
-YTITLE   (JPROC) = 'W'
-YUNIT    (JPROC) = 'm s-1'
-YCOMMENT (JPROC) = 'Air vertical speed'
-ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XW(:)
-!
-JPROC = JPROC + 1
-YTITLE   (JPROC) = 'Th'
-YUNIT    (JPROC) = 'K'
-YCOMMENT (JPROC) = 'Potential temperature'
-ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XTH(:)
-!
-IF (LDIAG_SURFRAD) THEN
-  IF (CSURF=="EXTE") THEN
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'T2m'
-    YUNIT    (JPROC) = 'K'
-    YCOMMENT (JPROC) = '2-m temperature'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XT2M(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'Q2m'
-    YUNIT    (JPROC) = 'kg kg-1'
-    YCOMMENT (JPROC) = '2-m humidity'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XQ2M(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'HU2m'
-    YUNIT    (JPROC) = 'percent'
-    YCOMMENT (JPROC) = '2-m relative humidity'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XHU2M(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'zon10m'
-    YUNIT    (JPROC) = 'm s-1'
-    YCOMMENT (JPROC) = '10-m zonal wind'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XZON10M(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'mer10m'
-    YUNIT    (JPROC) = 'm s-1'
-    YCOMMENT (JPROC) = '10-m meridian wind'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XMER10M(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'RN'
-    YUNIT    (JPROC) = 'W m-2'
-    YCOMMENT (JPROC) = 'Net radiation'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XRN(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'H'
-    YUNIT    (JPROC) = 'W m-2'
-    YCOMMENT (JPROC) = 'Sensible heat flux'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XH(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'LE'
-    YUNIT    (JPROC) = 'W m-2'
-    YCOMMENT (JPROC) = 'Total Latent heat flux'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XLE(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'G'
-    YUNIT    (JPROC) = 'W m-2'
-    YCOMMENT (JPROC) = 'Storage heat flux'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XGFLUX(:)
-    !
-    JPROC = JPROC + 1
-    YTITLE   (JPROC) = 'LEI'
-    YUNIT    (JPROC) = 'W m-2'
-    YCOMMENT (JPROC) = 'Solid Latent heat flux'
-    ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XLEI(:)
-  END IF
- IF (CRAD /= 'NONE') THEN
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'SWD'
-  YUNIT    (JPROC) = 'W m-2'
-  YCOMMENT (JPROC) = 'Downward short-wave radiation'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XSWD(:)
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'SWU'
-  YUNIT    (JPROC) = 'W m-2'
-  YCOMMENT (JPROC) = 'Upward short-wave radiation'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XSWU(:)
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'LWD'
-  YUNIT    (JPROC) = 'W m-2'
-  YCOMMENT (JPROC) = 'Downward long-wave radiation'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XLWD(:)
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'LWU'
-  YUNIT    (JPROC) = 'W m-2'
-  YCOMMENT (JPROC) = 'Upward long-wave radiation'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XLWU(:)
-  JPROC = JPROC + 1
-  !
-  YTITLE   (JPROC) = 'SWDIR'
-  YUNIT    (JPROC) = 'W m-2'
-  YCOMMENT (JPROC) = 'Downward direct short-wave radiation'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XSWDIR(:)
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'SWDIFF'
-  YUNIT    (JPROC) = 'W m-2'
-  YCOMMENT (JPROC) = 'Downward diffuse short-wave radiation'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XSWDIFF(:)
-  !
-  JPROC = JPROC + 1
-  YTITLE   (JPROC) = 'DSTAOD'
-  YUNIT    (JPROC) = 'm'
-  YCOMMENT (JPROC) = 'Dust aerosol optical depth'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XDSTAOD(:)
-  !
- END IF
-ENDIF
-!
-DO JRR=1,SIZE(TPSTATION%XR,2)
-  JPROC = JPROC+1
-  YUNIT    (JPROC) = 'kg kg-1'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XR(:,JRR)
-  IF (JRR==1) THEN
-    YTITLE   (JPROC) = 'Rv'
-    YCOMMENT (JPROC) = 'Water vapor mixing ratio'
-  ELSE IF (JRR==2) THEN
-    YTITLE   (JPROC) = 'Rc'
-    YCOMMENT (JPROC) = 'Liquid cloud water mixing ratio'
-  ELSE IF (JRR==3) THEN
-    YTITLE   (JPROC) = 'Rr'
-    YCOMMENT (JPROC) = 'Rain water mixing ratio'
-  ELSE IF (JRR==4) THEN
-    YTITLE   (JPROC) = 'Ri'
-    YCOMMENT (JPROC) = 'Ice cloud water mixing ratio'
-  ELSE IF (JRR==5) THEN
-    YTITLE   (JPROC) = 'Rs'
-    YCOMMENT (JPROC) = 'Snow mixing ratio'
-  ELSE IF (JRR==6) THEN
-    YTITLE   (JPROC) = 'Rg'
-    YCOMMENT (JPROC) = 'Graupel mixing ratio'
-  ELSE IF (JRR==7) THEN
-    YTITLE   (JPROC) = 'Rh'
-    YCOMMENT (JPROC) = 'Hail mixing ratio'
-  END IF
-END DO
-!
-IF ( CTURB == 'TKEL' ) THEN
-  JPROC = JPROC+1
-  YTITLE   (JPROC) = 'Tke'
-  YUNIT    (JPROC) = 'm2 s-2'
-  YCOMMENT (JPROC) = 'Turbulent kinetic energy'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XTKE(:)
-END IF
-!
-IF (SIZE(TPSTATION%XSV,2)>=1) THEN
+call Add_point( 'ZS', 'Orography', 'm',  SPREAD( tpstation%xzs, 1, istore ) )
+call Add_point( 'P',  'Pressure',  'Pa', tpstation%xp(:) )
+! call Add_point( 'Z', 'Z Pos', 'm', SPREAD( tpstation%xz, 1, istore ) )
+
+if ( lcartesian ) then
+  call Add_point( 'X', 'X Pos', 'm', SPREAD( tpstation%xx, 1, istore ) )
+  call Add_point( 'Y', 'Y Pos', 'm', SPREAD( tpstation%xy, 1, istore ) )
+  call Add_point( 'U', 'Axial velocity',       'm s-1', tpstation%xzon(:) )
+  call Add_point( 'V', 'Transversal velocity', 'm s-1', tpstation%xmer(:) )
+else
+  call Add_point( 'LON', 'Longitude', 'degree', SPREAD( tpstation%xlon, 1, istore ) )
+  call Add_point( 'LAT', 'Latitude',  'degree', SPREAD( tpstation%xlat, 1, istore ) )
+  call Add_point( 'ZON_WIND', 'Zonal wind',      'm s-1', tpstation%xzon(:) )
+  call Add_point( 'MER_WIND', 'Meridional wind', 'm s-1', tpstation%xmer(:) )
+end if
+
+call Add_point( 'W',  'Air vertical speed',    'm s-1', tpstation%xw(:)  )
+call Add_point( 'Th', 'Potential temperature', 'K',     tpstation%xth(:) )
+
+if ( ldiag_surfrad ) then
+  if ( csurf == "EXTE" ) then
+    call Add_point( 'T2m',    '2-m temperature',        'K',       tpstation%xt2m(:)    )
+    call Add_point( 'Q2m',    '2-m humidity',           'kg kg-1', tpstation%xq2m(:)    )
+    call Add_point( 'HU2m',   '2-m relative humidity',  'percent', tpstation%xhu2m(:)   )
+    call Add_point( 'zon10m', '10-m zonal wind',        'm s-1',   tpstation%xzon10m(:) )
+    call Add_point( 'mer10m', '10-m meridian wind',     'm s-1',   tpstation%xmer10m(:) )
+    call Add_point( 'RN',     'Net radiation',          'W m-2',   tpstation%xrn(:)     )
+    call Add_point( 'H',      'Sensible heat flux',     'W m-2',   tpstation%xh(:)      )
+    call Add_point( 'LE',     'Total Latent heat flux', 'W m-2',   tpstation%xle(:)     )
+    call Add_point( 'G',      'Storage heat flux',      'W m-2',   tpstation%xgflux(:)  )
+    call Add_point( 'LEI',    'Solid Latent heat flux', 'W m-2',   tpstation%xlei(:)    )
+  end if
+  if ( crad /= 'NONE' ) then
+    call Add_point( 'SWD',    'Downward short-wave radiation',         'W m-2', tpstation%xswd(:)    )
+    call Add_point( 'SWU',    'Upward short-wave radiation',           'W m-2', tpstation%xswu(:)    )
+    call Add_point( 'LWD',    'Downward long-wave radiation',          'W m-2', tpstation%xlwd(:)    )
+    call Add_point( 'LWU',    'Upward long-wave radiation',            'W m-2', tpstation%xlwu(:)    )
+    call Add_point( 'SWDIR',  'Downward direct short-wave radiation',  'W m-2', tpstation%xswdir(:)  )
+    call Add_point( 'SWDIFF', 'Downward diffuse short-wave radiation', 'W m-2', tpstation%xswdiff(:) )
+    call Add_point( 'DSTAOD', 'Dust aerosol optical depth',            'm',     tpstation%xdstaod(:) )
+  end if
+end if
+
+do jrr = 1, SIZE( tpstation%xr, 2 )
+  select case( jrr )
+    case (1)
+      call Add_point( 'Rv', 'Water vapor mixing ratio',        'kg kg-1', tpstation%xr(:,jrr) )
+    case (2)
+      call Add_point( 'Rc', 'Liquid cloud water mixing ratio', 'kg kg-1', tpstation%xr(:,jrr) )
+    case (3)
+      call Add_point( 'Rr', 'Rain water mixing ratio',         'kg kg-1', tpstation%xr(:,jrr) )
+    case (4)
+      call Add_point( 'Ri', 'Ice cloud water mixing ratio',    'kg kg-1', tpstation%xr(:,jrr) )
+    case (5)
+      call Add_point( 'Rs', 'Snow mixing ratio',               'kg kg-1', tpstation%xr(:,jrr) )
+    case (6)
+      call Add_point( 'Rg', 'Graupel mixing ratio',            'kg kg-1', tpstation%xr(:,jrr) )
+    case (7)
+      call Add_point( 'Rh', 'Hail mixing ratio',               'kg kg-1', tpstation%xr(:,jrr) )
+  end select
+end do
+
+if ( cturb == 'TKEL' ) call Add_point( 'Tke', 'Turbulent kinetic energy', 'm2 s-2', tpstation%xtke(:) )
+
+if ( nsv > 0 ) then
   ! Scalar variables
   DO JSV = 1, NSV
-    JPROC = JPROC + 1
-    YTITLE(JPROC)   = TRIM( TSVLIST(JSV)%CMNHNAME )
-    YCOMMENT(JPROC) = ''
     IF ( TRIM( TSVLIST(JSV)%CUNITS ) == 'ppv' ) THEN
-      YUNIT(JPROC)  = 'ppb'
-      ZWORK6(1,1,1,:,1,JPROC) = TPSTATION%XSV(:,JSV) * 1.e9 !*1e9 for conversion ppv->ppb
+      !*1e9 for conversion ppv->ppb
+      call Add_point( TRIM( TSVLIST(JSV)%CMNHNAME ), '', 'ppb', TPSTATION%XSV(:,JSV) * 1.e9 )
     ELSE
-      YUNIT(JPROC)  = TRIM( TSVLIST(JSV)%CUNITS )
-      ZWORK6(1,1,1,:,1,JPROC) = TPSTATION%XSV(:,JSV)
+      call Add_point( TRIM( TSVLIST(JSV)%CMNHNAME ), '', TSVLIST(JSV)%CUNITS, TPSTATION%XSV(:,JSV) )
     END IF
   END DO
 
   IF ((LORILAM).AND. .NOT.(ANY(TPSTATION%XP(:) == 0.))) THEN
-    ALLOCATE (ZSV(1,1,SIZE(tstations_time%tpdates),NSV_AER))
-    ALLOCATE (ZRHO(1,1,SIZE(tstations_time%tpdates)))
-    ALLOCATE (ZN0(1,1,SIZE(tstations_time%tpdates),JPMODE))
-    ALLOCATE (ZRG(1,1,SIZE(tstations_time%tpdates),JPMODE))
-    ALLOCATE (ZSIG(1,1,SIZE(tstations_time%tpdates),JPMODE))
-    ALLOCATE (ZPTOTA(1,1,SIZE(tstations_time%tpdates),NSP+NCARB+NSOA,JPMODE))
+    ALLOCATE (ZSV(1,1,ISTORE,NSV_AER))
+    ALLOCATE (ZRHO(1,1,ISTORE))
+    ALLOCATE (ZN0(1,1,ISTORE,JPMODE))
+    ALLOCATE (ZRG(1,1,ISTORE,JPMODE))
+    ALLOCATE (ZSIG(1,1,ISTORE,JPMODE))
+    ALLOCATE (ZPTOTA(1,1,ISTORE,NSP+NCARB+NSOA,JPMODE))
     ZSV(1,1,:,1:NSV_AER) = TPSTATION%XSV(:,NSV_AERBEG:NSV_AEREND)
     IF (SIZE(TPSTATION%XR,2) >0) THEN
       ZRHO(1,1,:) = 0.
@@ -622,127 +461,101 @@ IF (SIZE(TPSTATION%XSV,2)>=1) THEN
 
     DO JSV=1,JPMODE
       ! mean radius
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A6,I1)')'AERRGA',JSV
-      YUNIT    (JPROC) = 'um'
-      WRITE(YCOMMENT(JPROC),'(A18,I1)')'RG (nb) AERO MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZRG(1,1,:,JSV)
-      ! standard deviation
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A7,I1)')'AERSIGA',JSV
-      YUNIT    (JPROC) = '  '
-      WRITE(YCOMMENT(JPROC),'(A16,I1)')'SIGMA AERO MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZSIG(1,1,:,JSV)
-      ! particles number
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A6,I1)')'AERN0A',JSV
-      YUNIT    (JPROC) = 'm-3'
-      WRITE(YCOMMENT(JPROC),'(A13,I1)')'N0 AERO MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZN0(1,1,:,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MOC  ',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS OC   AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_OC,JSV)
-
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MBC  ',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS BC   AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_BC,JSV)
+      WRITE(YTITLE,'(A6,I1)')'AERRGA',JSV
+      WRITE(YCOMMENT,'(A18,I1)')'RG (nb) AERO MODE ',JSV
+      call Add_point( ytitle, ycomment, 'um', ZRG(1,1,:,JSV) )
 
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MDST  ',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS DST   AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_DST,JSV)
-
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSO4 ',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SO4  AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SO4,JSV)
-
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MNO3 ',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS NO3  AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_NO3,JSV)
+      ! standard deviation
+      WRITE(YTITLE,'(A7,I1)')'AERSIGA',JSV
+      WRITE(YCOMMENT,'(A16,I1)')'SIGMA AERO MODE ',JSV
+      call Add_point( ytitle, ycomment, '',ZSIG(1,1,:,JSV) )
 
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MH2O ',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS H2O  AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_H2O,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MNH3 ',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS NH3  AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_NH3,JSV)
-      JPROC = JPROC+1
-      IF (NSOA == 10) THEN
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA1',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA1 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA1,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA2',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA2 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA2,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA3',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA3 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA3,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA4',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA4 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA4,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA5',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA5 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA5,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA6',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA6 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA6,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA7',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA7 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA7,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA8',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA8 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA8,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A5,I1)')'MSOA9',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A23,I1)')'MASS SOA9 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA9,JSV)
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A6,I1)')'MSOA10',JSV
-      YUNIT    (JPROC) = 'ug m-3'
-      WRITE(YCOMMENT,'(A24,I1)')'MASS SOA10 AEROSOL MODE ',JSV
-      ZWORK6(1,1,1,:,1,JPROC)=ZPTOTA(1,1,:,JP_AER_SOA10,JSV)
+      ! particles number
+      WRITE(YTITLE,'(A6,I1)')'AERN0A',JSV
+      WRITE(YCOMMENT,'(A13,I1)')'N0 AERO MODE ',JSV
+      call Add_point( ytitle, ycomment, 'm-3', ZN0(1,1,:,JSV) )
+
+      WRITE(YTITLE,'(A5,I1)')'MOC  ',JSV
+      WRITE(CCOMMENT,'(A23,I1)')'MASS OC   AEROSOL MODE ',JSV
+      call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_OC,JSV) )
+
+      WRITE(YTITLE,'(A5,I1)')'MBC  ',JSV
+      WRITE(CCOMMENT,'(A23,I1)')'MASS BC   AEROSOL MODE ',JSV
+      call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_BC,JSV) )
+
+      WRITE(YTITLE,'(A5,I1)')'MDST  ',JSV
+      WRITE(CCOMMENT,'(A23,I1)')'MASS DST   AEROSOL MODE ',JSV
+      call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_DST,JSV) )
+
+      WRITE(YTITLE,'(A5,I1)')'MSO4 ',JSV
+      WRITE(CCOMMENT,'(A23,I1)')'MASS SO4  AEROSOL MODE ',JSV
+      call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SO4,JSV) )
+
+      WRITE(YTITLE,'(A5,I1)')'MNO3 ',JSV
+      WRITE(CCOMMENT,'(A23,I1)')'MASS NO3  AEROSOL MODE ',JSV
+      call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_NO3,JSV) )
+
+      WRITE(YTITLE,'(A5,I1)')'MH2O ',JSV
+      WRITE(CCOMMENT,'(A23,I1)')'MASS H2O  AEROSOL MODE ',JSV
+      call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_H2O,JSV) )
+
+      WRITE(YTITLE,'(A5,I1)')'MNH3 ',JSV
+      WRITE(CCOMMENT,'(A23,I1)')'MASS NH3  AEROSOL MODE ',JSV
+      call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_NH3,JSV) )
+
+      IF ( NSOA == 10 ) THEN
+        WRITE(YTITLE,'(A5,I1)')'MSOA1',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA1 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA1,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA2',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA2 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA2,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA3',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA3 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA3,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA4',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA4 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA4,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA5',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA5 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA5,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA6',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA6 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA6,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA7',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA7 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA7,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA8',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA8 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA8,JSV) )
+
+        WRITE(YTITLE,'(A5,I1)')'MSOA9',JSV
+        WRITE(CCOMMENT,'(A23,I1)')'MASS SOA9 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA9,JSV) )
+
+        WRITE(YTITLE,'(A6,I1)')'MSOA10',JSV
+        WRITE(CCOMMENT,'(A24,I1)')'MASS SOA10 AEROSOL MODE ',JSV
+        call Add_point( ytitle, ycomment, 'ug m-3', ZPTOTA(1,1,:,JP_AER_SOA10,JSV) )
       END IF
-    ENDDO
+    END DO
 
     DEALLOCATE (ZSV,ZRHO)
     DEALLOCATE (ZN0,ZRG,ZSIG)
   END IF
 
   IF ((LDUST).AND. .NOT.(ANY(TPSTATION%XP(:) == 0.))) THEN
-    ALLOCATE (ZSV(1,1,SIZE(tstations_time%tpdates),NSV_DST))
-    ALLOCATE (ZRHO(1,1,SIZE(tstations_time%tpdates)))
-    ALLOCATE (ZN0(1,1,SIZE(tstations_time%tpdates),NMODE_DST))
-    ALLOCATE (ZRG(1,1,SIZE(tstations_time%tpdates),NMODE_DST))
-    ALLOCATE (ZSIG(1,1,SIZE(tstations_time%tpdates),NMODE_DST))
+    ALLOCATE (ZSV(1,1,ISTORE,NSV_DST))
+    ALLOCATE (ZRHO(1,1,ISTORE))
+    ALLOCATE (ZN0(1,1,ISTORE,NMODE_DST))
+    ALLOCATE (ZRG(1,1,ISTORE,NMODE_DST))
+    ALLOCATE (ZSIG(1,1,ISTORE,NMODE_DST))
     ZSV(1,1,:,1:NSV_DST) = TPSTATION%XSV(:,NSV_DSTBEG:NSV_DSTEND)
     IF (SIZE(TPSTATION%XR,2) >0) THEN
       ZRHO(1,1,:) = 0.
@@ -760,33 +573,33 @@ IF (SIZE(TPSTATION%XSV,2)>=1) THEN
     DO JSV=1,NMODE_DST
       ! mean radius
       JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A6,I1)')'DSTRGA',JSV
-      YUNIT    (JPROC) = 'um'
-      WRITE(YCOMMENT(JPROC),'(A18,I1)')'RG (nb) DUST MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZRG(1,1,:,JSV)
+      WRITE(CTITLE(JPROC),'(A6,I1)')'DSTRGA',JSV
+      CUNIT    (JPROC) = 'um'
+      WRITE(CCOMMENT(JPROC),'(A18,I1)')'RG (nb) DUST MODE ',JSV
+      XWORK6 (1,1,1,:,1,JPROC) = ZRG(1,1,:,JSV)
       ! standard deviation
       JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A7,I1)')'DSTSIGA',JSV
-      YUNIT    (JPROC) = '  '
-      WRITE(YCOMMENT(JPROC),'(A16,I1)')'SIGMA DUST MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZSIG(1,1,:,JSV)
+      WRITE(CTITLE(JPROC),'(A7,I1)')'DSTSIGA',JSV
+      CUNIT    (JPROC) = '  '
+      WRITE(CCOMMENT(JPROC),'(A16,I1)')'SIGMA DUST MODE ',JSV
+      XWORK6 (1,1,1,:,1,JPROC) = ZSIG(1,1,:,JSV)
       ! particles number
       JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A6,I1)')'DSTN0A',JSV
-      YUNIT    (JPROC) = 'm-3'
-      WRITE(YCOMMENT(JPROC),'(A13,I1)')'N0 DUST MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZN0(1,1,:,JSV)
+      WRITE(CTITLE(JPROC),'(A6,I1)')'DSTN0A',JSV
+      CUNIT    (JPROC) = 'm-3'
+      WRITE(CCOMMENT(JPROC),'(A13,I1)')'N0 DUST MODE ',JSV
+      XWORK6 (1,1,1,:,1,JPROC) = ZN0(1,1,:,JSV)
     ENDDO
     DEALLOCATE (ZSV,ZRHO)
     DEALLOCATE (ZN0,ZRG,ZSIG)
   END IF
 
   IF ((LSALT).AND. .NOT.(ANY(TPSTATION%XP(:) == 0.))) THEN
-    ALLOCATE (ZSV(1,1,SIZE(tstations_time%tpdates),NSV_SLT))
-    ALLOCATE (ZRHO(1,1,SIZE(tstations_time%tpdates)))
-    ALLOCATE (ZN0(1,1,SIZE(tstations_time%tpdates),NMODE_SLT))
-    ALLOCATE (ZRG(1,1,SIZE(tstations_time%tpdates),NMODE_SLT))
-    ALLOCATE (ZSIG(1,1,SIZE(tstations_time%tpdates),NMODE_SLT))
+    ALLOCATE (ZSV(1,1,ISTORE,NSV_SLT))
+    ALLOCATE (ZRHO(1,1,ISTORE))
+    ALLOCATE (ZN0(1,1,ISTORE,NMODE_SLT))
+    ALLOCATE (ZRG(1,1,ISTORE,NMODE_SLT))
+    ALLOCATE (ZSIG(1,1,ISTORE,NMODE_SLT))
     ZSV(1,1,:,1:NSV_SLT) = TPSTATION%XSV(:,NSV_SLTBEG:NSV_SLTEND)
     IF (SIZE(TPSTATION%XR,2) >0) THEN
       ZRHO(1,1,:) = 0.
@@ -803,59 +616,39 @@ IF (SIZE(TPSTATION%XSV,2)>=1) THEN
     CALL PPP2SALT(ZSV,ZRHO, PSIG3D=ZSIG, PRG3D=ZRG, PN3D=ZN0)
     DO JSV=1,NMODE_SLT
       ! mean radius
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A6,I1)')'SLTRGA',JSV
-      YUNIT    (JPROC) = 'um'
-      WRITE(YCOMMENT(JPROC),'(A18,I1)')'RG (nb) SALT MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZRG(1,1,:,JSV)
+      WRITE(CTITLE(JPROC),'(A6,I1)')'SLTRGA',JSV
+      WRITE(CCOMMENT(JPROC),'(A18,I1)')'RG (nb) SALT MODE ',JSV
+      call Add_point( ytitle, ycomment, 'um', ZRG(1,1,:,JSV) )
+
       ! standard deviation
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A7,I1)')'SLTSIGA',JSV
-      YUNIT    (JPROC) = '  '
-      WRITE(YCOMMENT(JPROC),'(A16,I1)')'SIGMA DUST MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZSIG(1,1,:,JSV)
+      WRITE(CTITLE(JPROC),'(A7,I1)')'SLTSIGA',JSV
+      WRITE(CCOMMENT(JPROC),'(A16,I1)')'SIGMA DUST MODE ',JSV
+      call Add_point( ytitle, ycomment, '',ZSIG(1,1,:,JSV) )
+
       ! particles number
-      JPROC = JPROC+1
-      WRITE(YTITLE(JPROC),'(A6,I1)')'SLTN0A',JSV
-      YUNIT    (JPROC) = 'm-3'
-      WRITE(YCOMMENT(JPROC),'(A13,I1)')'N0 DUST MODE ',JSV
-      ZWORK6 (1,1,1,:,1,JPROC) = ZN0(1,1,:,JSV)
+      WRITE(CTITLE(JPROC),'(A6,I1)')'SLTN0A',JSV
+      WRITE(CCOMMENT(JPROC),'(A13,I1)')'N0 DUST MODE ',JSV
+      call Add_point( ytitle, ycomment, 'm-3', ZN0(1,1,:,JSV) )
     ENDDO
     DEALLOCATE (ZSV,ZRHO)
     DEALLOCATE (ZN0,ZRG,ZSIG)
   END IF
-END IF
+end if
 
-IF ( CRAD /= 'NONE' ) THEN
-  JPROC = JPROC+1
-  YTITLE   (JPROC) = 'Tsrad'
-  YUNIT    (JPROC) = 'K'
-  YCOMMENT (JPROC) = 'Radiative Surface Temperature'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XTSRAD(:)
-END IF
-!
-! IF (ANY(TPSTATION%XSFCO2(:)/=XUNDEF)) THEN
-  JPROC = JPROC+1
-  YTITLE   (JPROC) = 'SFCO2'
-  YUNIT    (JPROC) = 'mg m-2 s-1'
-  YCOMMENT (JPROC) = 'CO2 Surface Flux'
-  ZWORK6 (1,1,1,:,1,JPROC) = TPSTATION%XSFCO2(:)
-! END IF
+if ( crad /= 'NONE' ) call Add_point( 'Tsrad', 'Radiative Surface Temperature', 'K', tpstation%xtsrad(:) )
+
+if ( ldiag_surfrad ) call Add_point( 'SFCO2', 'CO2 Surface Flux', 'mg m-2 s-1', tpstation%xsfco2(:) )
 !
 !----------------------------------------------------------------------------
 !
 !
-ALLOCATE (ZW6(1,1,1,SIZE(tstations_time%tpdates),1,JPROC))
-ZW6 = ZWORK6(:,:,:,:,:,:JPROC)
-DEALLOCATE(ZWORK6)
-!
 allocate( tzfields( jproc ) )
 
-tzfields(:)%cmnhname  = ytitle(1 : jproc)
+tzfields(:)%cmnhname  = ctitle(1 : jproc)
 tzfields(:)%cstdname  = ''
-tzfields(:)%clongname = ytitle(1 : jproc)
-tzfields(:)%cunits    = yunit(1 : jproc)
-tzfields(:)%ccomment  = ycomment(1 : jproc)
+tzfields(:)%clongname = ctitle(1 : jproc)
+tzfields(:)%cunits    = cunit(1 : jproc)
+tzfields(:)%ccomment  = ccomment(1 : jproc)
 tzfields(:)%ngrid     = 0
 tzfields(:)%ntype     = TYPEREAL
 tzfields(:)%ndims     = 2
@@ -875,12 +668,12 @@ tzbudiachro%clevels  (NLVL_SUBCATEGORY) = ''
 tzbudiachro%ccomments(NLVL_SUBCATEGORY) = ''
 
 tzbudiachro%lleveluse(NLVL_GROUP)       = .true.
-tzbudiachro%clevels  (NLVL_GROUP)       = ygroup
-tzbudiachro%ccomments(NLVL_GROUP)       = 'Values at position of station ' // Trim( ygroup )
+tzbudiachro%clevels  (NLVL_GROUP)       = tpstation%cname
+tzbudiachro%ccomments(NLVL_GROUP)       = 'Values at position of station ' // Trim( tpstation%cname )
 
 tzbudiachro%lleveluse(NLVL_SHAPE)       = .false.
 tzbudiachro%clevels  (NLVL_SHAPE)       = 'Point'
-tzbudiachro%ccomments(NLVL_SHAPE)       = 'Values at position of station ' // Trim( ygroup )
+tzbudiachro%ccomments(NLVL_SHAPE)       = 'Values at position of station ' // Trim( tpstation%cname )
 
 tzbudiachro%lleveluse(NLVL_TIMEAVG)     = .false.
 tzbudiachro%clevels  (NLVL_TIMEAVG)     = 'Not_time_averaged'
@@ -912,18 +705,45 @@ tzbudiachro%njh        = 1
 tzbudiachro%nkl        = 1
 tzbudiachro%nkh        = 1
 
-call Write_diachro( tpdiafile, tzbudiachro, tzfields, tstations_time%tpdates, zw6 )
+call Write_diachro( tpdiafile, tzbudiachro, tzfields, tstations_time%tpdates, xwork6(:,:,:,:,:,:jproc) )
 
 deallocate( tzfields )
 
-DEALLOCATE (ZW6)
-DEALLOCATE (YCOMMENT)
-DEALLOCATE (YTITLE  )
-DEALLOCATE (YUNIT   )
-DEALLOCATE (IGRID   )
+!Necessary because global variables (private inside module)
+Deallocate( xwork6  )
+Deallocate (ccomment)
+Deallocate (ctitle  )
+Deallocate (cunit   )
+
 !----------------------------------------------------------------------------
+
+contains
+
+! ######################################################
+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
+
+integer :: jk
+
+jproc = jproc + 1
+
+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(:)
+
+end subroutine Add_point
+
 END SUBROUTINE STATION_DIACHRO_n
-!
-END SUBROUTINE WRITE_STATION_n
 
 END MODULE MODE_WRITE_STATION_n