Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4627)
+++ palm/trunk/SOURCE/Makefile	(revision 4628)
@@ -20,9 +20,12 @@
 # Current revisions:
 # ------------------
-# 
-# 
+#
+#
 # Former revisions:
 # -----------------
 # $Id$
+# dependencies for particle data MPI-I/O added
+# 
+# 4564 2020-06-12 14:03:36Z raasch
 # Vertical nesting method of Huq et al. (2019) removed
 # 
@@ -849,4 +852,5 @@
 
 mod_particle_attributes.o: \
+	modules.o \
 	mod_kinds.o
 nesting_offl_mod.o: \
@@ -1090,7 +1094,8 @@
 	virtual_measurement_mod.o
 restart_data_mpi_io_mod.o: \
-	modules.o \
-	mod_kinds.o \
-	exchange_horiz_mod.o \
+	exchange_horiz_mod.o \
+	modules.o \
+	mod_kinds.o \
+	mod_particle_attributes.o \
 	posix_interface_mod.o \
 	shared_memory_io_mod.o
Index: palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4627)
+++ palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4628)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! extensions required for MPI-I/O of particle data to restart files
+!
+! 4616 2020-07-21 10:09:46Z schwenkel
 ! Bugfix in case of strechting: k-calculation limited lower bound of 1
 ! 
@@ -189,9 +192,10 @@
                child_domain,                                                   &
                cloud_droplets, constant_flux_layer, current_timestep_number,   &
-               dt_3d, dt_3d_reached, first_call_lpm, humidity,                 &
+               dt_3d, dt_3d_reached, debug_output, first_call_lpm, humidity,   &
                dt_3d_reached_l, dt_dopts, dz, initializing_actions,            &
                intermediate_timestep_count, intermediate_timestep_count_max,   &
                message_string, molecular_viscosity, ocean_mode,                &
-               particle_maximum_age, restart_data_format_output,               &
+               particle_maximum_age, restart_data_format_input,                &
+               restart_data_format_output,                                     &
                simulated_time, topography, dopts_time_count,                   &
                time_since_reference_point, rho_surface, u_gtrans, v_gtrans,    &
@@ -240,6 +244,16 @@
                id_random_array
 
-    USE restart_data_mpi_io_mod,                                                                   &
-        ONLY:  rd_mpi_io_check_array, rrd_mpi_io, wrd_mpi_io
+    USE restart_data_mpi_io_mod,                                               &
+        ONLY:  rd_mpi_io_check_array,                                          &
+               rd_mpi_io_check_open,                                           &
+               rd_mpi_io_close,                                                &
+               rd_mpi_io_open,                                                 &
+               rd_mpi_io_particle_filetypes,                                   &
+               rrd_mpi_io,                                                     &
+               rrd_mpi_io_global_array,                                        &
+               rrd_mpi_io_particles,                                           &
+               wrd_mpi_io,                                                     &
+               wrd_mpi_io_global_array,                                        &
+               wrd_mpi_io_particles
 
     USE statistics,                                                            &
@@ -251,4 +265,8 @@
                surf_lsm_h,                                                     &
                surf_usm_h
+
+!-- Next lines are in preparation for the output of particle data
+!    USE data_output_particle_mod,                                              &
+!        ONLY:  dop_active, dop_init, dop_output_tseries
 
 #if defined( __parallel )  &&  !defined( __mpifh )
@@ -594,5 +612,5 @@
        write_particle_statistics
 
-       NAMELIST /particle_parameters/ &
+    NAMELIST /particle_parameters/ &
        aero_species, &
        aero_type, &
@@ -619,10 +637,15 @@
        na, &
        number_concentration, &
+       number_of_output_particles, &
        number_of_particle_groups, &
        number_particles_per_gridbox, &
+       oversize, &
        particles_per_point, &
        particle_advection_start, &
        particle_advection_interpolation, &
        particle_maximum_age, &
+       part_output, &
+       part_inc, &
+       part_percent, &
        pdx, &
        pdy, &
@@ -648,4 +671,5 @@
        splitting_mode, &
        step_dealloc, &
+       unlimited_dimension, &
        use_sgs_for_particles, &
        vertical_particle_advection, &
@@ -1215,4 +1239,5 @@
 !-- Initialize collision kernels
     IF ( collision_kernel /= 'none' )  CALL lpm_init_kernels
+
 !
 !-- For the first model run of a possible job chain initialize the
@@ -1241,5 +1266,5 @@
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
-                                      0, 0, 0_idp, .FALSE., -1 )
+                                      0, 0, 0_idp, .FALSE., -1, -1 )
 
        particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
@@ -1304,4 +1329,7 @@
     IF ( nested_run )  CALL pmcp_g_init
 #endif
+
+!-- next line is in preparation for particle data output
+!    CALL dop_init
 
 !
@@ -2995,6 +3023,8 @@
  SUBROUTINE lpm_rrd_local_particles
 
-    CHARACTER (LEN=10) ::  particle_binary_version    !<
-    CHARACTER (LEN=10) ::  version_on_file            !<
+    CHARACTER(LEN=10) ::  particle_binary_version    !<
+    CHARACTER(LEN=10) ::  version_on_file            !<
+
+    CHARACTER(LEN=20) ::  save_restart_data_format_input  !<
 
     INTEGER(iwp) ::  alloc_size !<
@@ -3003,89 +3033,157 @@
     INTEGER(iwp) ::  kp         !<
 
+    INTEGER(idp), ALLOCATABLE, DIMENSION(:,:,:) ::  prt_global_index !<
+
+    LOGICAL ::  save_debug_output  !<
+
     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  tmp_particles !<
 
-!
-!-- Read particle data from previous model run.
-!-- First open the input unit.
-    IF ( myid_char == '' )  THEN
-       OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN'//myid_char,                  &
-                  FORM='UNFORMATTED' )
-    ELSE
-       OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN/'//myid_char,                 &
-                  FORM='UNFORMATTED' )
-    ENDIF
-
-!
-!-- First compare the version numbers
-    READ ( 90 )  version_on_file
-    particle_binary_version = '4.0'
-    IF ( TRIM( version_on_file ) /= TRIM( particle_binary_version ) )  THEN
-       message_string = 'version mismatch concerning data from prior ' //      &
-                        'run &version on file = "' //                          &
-                                      TRIM( version_on_file ) //               &
-                        '&version in program = "' //                           &
-                                      TRIM( particle_binary_version ) // '"'
-       CALL message( 'lpm_read_restart_file', 'PA0214', 1, 2, 0, 6, 0 )
-    ENDIF
-
-!
-!-- If less particles are stored on the restart file than prescribed by
-!-- 1, the remainder is initialized by zero_particle to avoid
-!-- errors.
-    zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
-                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
-                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
-                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
-                                   0, 0, 0_idp, .FALSE., -1 )
-!
-!-- Read some particle parameters and the size of the particle arrays,
-!-- allocate them and read their contents.
-    READ ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,                     &
-                 last_particle_release_time, number_of_particle_groups,        &
-                 particle_groups, time_write_particle_data
-
-    ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
-              grid_particles(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-
-    READ ( 90 )  prt_count
-
-    DO  ip = nxl, nxr
-       DO  jp = nys, nyn
-          DO  kp = nzb+1, nzt
-
-             number_of_particles = prt_count(kp,jp,ip)
-             IF ( number_of_particles > 0 )  THEN
-                alloc_size = MAX( INT( number_of_particles *                   &
-                             ( 1.0_wp + alloc_factor / 100.0_wp ) ),           &
-                             1 )
-             ELSE
-                alloc_size = 1
-             ENDIF
-
-             ALLOCATE( grid_particles(kp,jp,ip)%particles(1:alloc_size) )
-
-             IF ( number_of_particles > 0 )  THEN
-                ALLOCATE( tmp_particles(1:number_of_particles) )
-                READ ( 90 )  tmp_particles
-                grid_particles(kp,jp,ip)%particles(1:number_of_particles) = tmp_particles
-                DEALLOCATE( tmp_particles )
-                IF ( number_of_particles < alloc_size )  THEN
-                   grid_particles(kp,jp,ip)%particles(number_of_particles+1:alloc_size) &
-                      = zero_particle
+    IF ( TRIM( restart_data_format_input ) == 'fortran_binary' )  THEN
+
+!
+!--    Read particle data from previous model run.
+!--    First open the input unit.
+       IF ( myid_char == '' )  THEN
+          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN'//myid_char,                  &
+                     FORM='UNFORMATTED' )
+       ELSE
+          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN/'//myid_char,                 &
+                     FORM='UNFORMATTED' )
+       ENDIF
+
+!
+!--    First compare the version numbers
+       READ ( 90 )  version_on_file
+       particle_binary_version = '4.0'
+       IF ( TRIM( version_on_file ) /= TRIM( particle_binary_version ) )  THEN
+          message_string = 'version mismatch concerning data from prior ' //      &
+                           'run &version on file = "' //                          &
+                                         TRIM( version_on_file ) //               &
+                           '&version in program = "' //                           &
+                                         TRIM( particle_binary_version ) // '"'
+          CALL message( 'lpm_read_restart_file', 'PA0214', 1, 2, 0, 6, 0 )
+       ENDIF
+
+!
+!--    If less particles are stored on the restart file than prescribed by
+!--    1, the remainder is initialized by zero_particle to avoid
+!--    errors.
+       zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
+                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
+                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
+                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,     &
+                                      0, 0, 0_idp, .FALSE., -1, -1 )
+!
+!--    Read some particle parameters and the size of the particle arrays,
+!--    allocate them and read their contents.
+       READ ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,                     &
+                    last_particle_release_time, number_of_particle_groups,        &
+                    particle_groups, time_write_particle_data
+
+       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
+                 grid_particles(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+       READ ( 90 )  prt_count
+
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+
+                number_of_particles = prt_count(kp,jp,ip)
+                IF ( number_of_particles > 0 )  THEN
+                   alloc_size = MAX( INT( number_of_particles *                   &
+                                ( 1.0_wp + alloc_factor / 100.0_wp ) ),           &
+                                1 )
+                ELSE
+                   alloc_size = 1
                 ENDIF
-             ELSE
-                grid_particles(kp,jp,ip)%particles(1:alloc_size) = zero_particle
-             ENDIF
-
+
+                ALLOCATE( grid_particles(kp,jp,ip)%particles(1:alloc_size) )
+
+                IF ( number_of_particles > 0 )  THEN
+                   ALLOCATE( tmp_particles(1:number_of_particles) )
+                   READ ( 90 )  tmp_particles
+                   grid_particles(kp,jp,ip)%particles(1:number_of_particles) = tmp_particles
+                   DEALLOCATE( tmp_particles )
+                   IF ( number_of_particles < alloc_size )  THEN
+                      grid_particles(kp,jp,ip)%particles(number_of_particles+1:alloc_size) &
+                         = zero_particle
+                   ENDIF
+                ELSE
+                   grid_particles(kp,jp,ip)%particles(1:alloc_size) = zero_particle
+                ENDIF
+
+             ENDDO
           ENDDO
        ENDDO
-    ENDDO
-
-    CLOSE ( 90 )
+
+       CLOSE ( 90 )
+
+    ELSEIF ( restart_data_format_input(1:3) == 'mpi' )  THEN
+
+       WRITE ( 9, * )  'Here is MPI-IO praticle input ', rd_mpi_io_check_open()
+       FLUSH(9)
+
+       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                        &
+                 grid_particles(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+       ALLOCATE( prt_global_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+!
+!--    Open restart file for read, if not already open, and do not allow usage of
+!--    shared-memory I/O
+       IF ( .NOT. rd_mpi_io_check_open() )  THEN
+          save_restart_data_format_input = restart_data_format_input
+          restart_data_format_input = 'mpi'
+          CALL rd_mpi_io_open( 'READ', 'BININ' )
+          restart_data_format_input = save_restart_data_format_input
+       ENDIF
+
+       CALL  rd_mpi_io_particle_filetypes
+
+       prt_count = 0
+       CALL rrd_mpi_io( 'prt_count', prt_count )
+       CALL rrd_mpi_io( 'prt_global_index', prt_global_index )
+
+!
+!--    Allocate particles arrays
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+
+                number_of_particles = prt_count(kp,jp,ip)
+                IF ( number_of_particles > 0 )  THEN
+                   alloc_size = MAX( INT( number_of_particles *                   &
+                      ( 1.0_wp + alloc_factor / 100.0_wp ) ),           &
+                      1 )
+                ELSE
+                   alloc_size = 1
+                ENDIF
+
+                ALLOCATE( grid_particles(kp,jp,ip)%particles(1:alloc_size) )
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Now read particle data from restart file
+       CALL rrd_mpi_io_particles( 'particles', prt_global_index )
+
+       IF ( .NOT. rd_mpi_io_check_open() )  THEN
+!
+!--       Do not print header a second time to the debug file
+          save_debug_output = debug_output
+          debug_output      = .FALSE.
+          CALL rd_mpi_io_close()
+          debug_output = save_debug_output
+       ENDIF
+
+       DEALLOCATE( prt_global_index )
+
+    ENDIF
 !
 !-- Must be called to sort particles into blocks, which is needed for a fast
 !-- interpolation of the LES fields on the particle position.
     CALL lpm_sort_and_delete
-
 
  END SUBROUTINE lpm_rrd_local_particles
@@ -3256,57 +3354,73 @@
  
     CHARACTER (LEN=10) ::  particle_binary_version   !<
-    CHARACTER (LEN=20) ::  tmp_name                  !< temporary variable
+    CHARACTER (LEN=32) ::  tmp_name                  !< temporary variable
 
     INTEGER(iwp) ::  i                               !< loop index
     INTEGER(iwp) ::  ip                              !<
+    INTEGER(iwp) ::  j                               !< loop index
     INTEGER(iwp) ::  jp                              !<
+    INTEGER(iwp) ::  k                               !< loop index
     INTEGER(iwp) ::  kp                              !< 
-!
-!-- First open the output unit.
-    IF ( myid_char == '' )  THEN
-       OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, &
-                  FORM='UNFORMATTED')
-    ELSE
-       IF ( myid == 0 )  CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' )
+
 #if defined( __parallel )
-!
-!--    Set a barrier in order to allow that thereafter all other processors
-!--    in the directory created by PE0 can open their file
-       CALL MPI_BARRIER( comm2d, ierr )
+    INTEGER      :: ierr                             !<
 #endif
-       OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, &
-                  FORM='UNFORMATTED' )
-    ENDIF
-
-!
-!-- Write the version number of the binary format.
-!-- Attention: After changes to the following output commands the version
-!-- ---------  number of the variable particle_binary_version must be
-!--            changed! Also, the version number and the list of arrays
-!--            to be read in lpm_read_restart_file must be adjusted
-!--            accordingly.
-    particle_binary_version = '4.0'
-    WRITE ( 90 )  particle_binary_version
-
-!
-!-- Write some particle parameters, the size of the particle arrays
-    WRITE ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,                    &
-                  last_particle_release_time, number_of_particle_groups,       &
-                  particle_groups, time_write_particle_data
-
-    WRITE ( 90 )  prt_count
+    INTEGER(iwp) ::  start_index                     !<
+    INTEGER(iwp) ::  start_index_on_pe               !<
+
+    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  nr_particles_global
+    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  nr_particles_local
+
+    INTEGER(idp), ALLOCATABLE, DIMENSION(:,:,:) ::  prt_global_index
+
+
+    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
+!
+!--    First open the output unit.
+       IF ( myid_char == '' )  THEN
+          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, &
+                     FORM='UNFORMATTED')
+       ELSE
+          IF ( myid == 0 )  CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' )
+#if defined( __parallel )
+!
+!--       Set a barrier in order to allow that thereafter all other processors
+!--       in the directory created by PE0 can open their file
+          CALL MPI_BARRIER( comm2d, ierr )
+#endif
+          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, &
+                     FORM='UNFORMATTED' )
+       ENDIF
+
+!
+!--    Write the version number of the binary format.
+!--    Attention: After changes to the following output commands the version
+!--    ---------  number of the variable particle_binary_version must be
+!--               changed! Also, the version number and the list of arrays
+!--               to be read in lpm_read_restart_file must be adjusted
+!--               accordingly.
+       particle_binary_version = '4.0'
+       WRITE ( 90 )  particle_binary_version
+
+!
+!--    Write some particle parameters, the size of the particle arrays
+       WRITE ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,                    &
+                     last_particle_release_time, number_of_particle_groups,       &
+                     particle_groups, time_write_particle_data
+
+       WRITE ( 90 )  prt_count
           
-    DO  ip = nxl, nxr
-       DO  jp = nys, nyn
-          DO  kp = nzb+1, nzt
-             number_of_particles = prt_count(kp,jp,ip)
-             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
-             IF ( number_of_particles <= 0 )  CYCLE
-             WRITE ( 90 )  particles
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+                number_of_particles = prt_count(kp,jp,ip)
+                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                IF ( number_of_particles <= 0 )  CYCLE
+                WRITE ( 90 )  particles
+             ENDDO
           ENDDO
        ENDDO
-    ENDDO
-
-    CLOSE ( 90 )
+
+       CLOSE ( 90 )
 
 #if defined( __parallel )
@@ -3314,5 +3428,4 @@
 #endif
 
-    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
 
        IF ( ALLOCATED( seq_random_array_particles ) )  THEN
@@ -3321,5 +3434,7 @@
        ENDIF
 
+
     ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
+
 
        IF ( ALLOCATED( seq_random_array_particles ) )  THEN
@@ -3330,4 +3445,51 @@
        ENDIF
 
+       ALLOCATE( prt_global_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+#if defined( __parallel )
+!--    TODO: needs to be replaced by standard PALM message
+       IF ( TRIM( restart_data_format_output ) == 'mpi_shared_memory' )   THEN
+          WRITE( 9, * )  'mpi_shared_memory is NOT implemented yet for particle IO'
+          FLUSH( 9 )
+          CALL MPI_ABORT( MPI_COMM_WORLD, 1, ierr )
+       ENDIF
+#endif
+
+       CALL rd_mpi_io_particle_filetypes
+
+       nr_particles_local = 0
+       nr_particles_local(myid) = SUM( prt_count )
+
+#if defined( __parallel )
+       CALL MPI_ALLREDUCE( nr_particles_local, nr_particles_global, numprocs, MPI_INTEGER,         &
+                           MPI_SUM, comm2d, ierr )
+#else
+       nr_particles_global = nr_particles_local
+#endif
+
+       start_index_on_pe = 0
+       IF ( myid > 0 )  THEN
+          DO  i = 1, myid
+             start_index_on_pe = start_index_on_pe + nr_particles_global(i-1)
+          ENDDO
+       ENDIF
+
+       CALL wrd_mpi_io( 'prt_count', prt_count )
+
+       start_index = start_index_on_pe
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb, nzt+1
+                prt_global_index(k,j,i) = start_index
+                start_index             = start_index + prt_count(k,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL wrd_mpi_io( 'prt_global_index', prt_global_index )
+       CALL wrd_mpi_io_particles( 'particles', prt_global_index )
+
+       DEALLOCATE( prt_global_index )
+
     ENDIF
 
@@ -3341,4 +3503,11 @@
 !------------------------------------------------------------------------------! 
  SUBROUTINE lpm_wrd_global
+
+#if defined( __parallel )
+    INTEGER :: ierr  !<
+#endif
+
+    REAL(wp), DIMENSION(4,max_number_of_particle_groups) ::  particle_groups_array  !<
+
  
     IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
@@ -3362,4 +3531,35 @@
        CALL wrd_mpi_io( 'interpolation_simple_predictor', interpolation_simple_predictor )
        CALL wrd_mpi_io( 'interpolation_trilinear', interpolation_trilinear )
+!
+!--    Write some global particle parameters
+!--    In case of Fortran binary format, these variables are written to unit 90
+       CALL wrd_mpi_io( 'bc_par_b', bc_par_b )
+       CALL wrd_mpi_io( 'bc_par_lr', bc_par_lr )
+       CALL wrd_mpi_io( 'bc_par_ns', bc_par_ns )
+       CALL wrd_mpi_io( 'bc_par_t', bc_par_t )
+       CALL wrd_mpi_io( 'last_particle_release_time', last_particle_release_time )
+       CALL wrd_mpi_io( 'number_of_particle_groups', number_of_particle_groups )
+       CALL wrd_mpi_io( 'time_write_particle_data', time_write_particle_data )
+
+!
+!--    Write particle_group informations via 2D array to avoid another overlay in
+!--    restart_data_mpi_io_mod.
+!--    TODO: convert the following to a standard PALM message
+       IF( STORAGE_SIZE( particle_groups(1) ) / (wp*8) /= 4 )  THEN
+          WRITE( 9, * )  'size of structure particle_groups_type has changed '
+          FLUSH( 9 )
+#if defined( __parallel )
+          CALL MPI_ABORT( MPI_COMM_WORLD, 1, ierr )
+#else
+          STOP 'error'
+#endif
+       ENDIF
+
+       particle_groups_array(1,:) = particle_groups(:)%density_ratio
+       particle_groups_array(2,:) = particle_groups(:)%radius
+       particle_groups_array(3,:) = particle_groups(:)%exp_arg
+       particle_groups_array(4,:) = particle_groups(:)%exp_term
+
+       CALL wrd_mpi_io_global_array( 'particle_groups', particle_groups_array )
 
     ENDIF
@@ -3412,8 +3612,46 @@
  SUBROUTINE lpm_rrd_global_mpi
 
+#if defined( __parallel )
+    INTEGER    :: ierr    !<
+#endif
+
+    REAL(wp), DIMENSION(4,max_number_of_particle_groups) ::  particle_groups_array  !<
+
+
     CALL rrd_mpi_io( 'curvature_solution_effects', curvature_solution_effects )
     CALL rrd_mpi_io( 'interpolation_simple_corrector', interpolation_simple_corrector )
     CALL rrd_mpi_io( 'interpolation_simple_predictor', interpolation_simple_predictor )
     CALL rrd_mpi_io( 'interpolation_trilinear', interpolation_trilinear )
+!
+!-- Read some particle parameters.
+!-- In case of Fortran binary format, these variables are read from unit 90.
+    CALL rrd_mpi_io( 'bc_par_b', bc_par_b )
+    CALL rrd_mpi_io( 'bc_par_lr', bc_par_lr )
+    CALL rrd_mpi_io( 'bc_par_ns', bc_par_ns )
+    CALL rrd_mpi_io( 'bc_par_t', bc_par_t )
+    CALL rrd_mpi_io( 'last_particle_release_time', last_particle_release_time )
+    CALL rrd_mpi_io( 'number_of_particle_groups', number_of_particle_groups )
+    CALL rrd_mpi_io( 'time_write_particle_data', time_write_particle_data )
+
+!
+!-- Read particle group information via 2d-array to avoid another overlay in
+!-- restart_data_mpi_io_mod.
+!-- TODO: convert the following to a standard PALM message
+    IF ( STORAGE_SIZE( particle_groups(1) ) / (wp*8) /= 4 )  THEN
+       WRITE( 9, * )  'size of structure particle_groups_type has changed '
+       FLUSH( 9 )
+#if defined( __parallel )
+       CALL MPI_ABORT( MPI_COMM_WORLD, 1, ierr )
+#else
+       STOP 'error'
+#endif
+    ENDIF
+
+    CALL rrd_mpi_io_global_array( 'particle_groups', particle_groups_array )
+
+    particle_groups(:)%density_ratio = particle_groups_array(1,:)
+    particle_groups(:)%radius        = particle_groups_array(2,:)
+    particle_groups(:)%exp_arg       = particle_groups_array(3,:)
+    particle_groups(:)%exp_term      = particle_groups_array(4,:)
 
  END SUBROUTINE lpm_rrd_global_mpi
Index: palm/trunk/SOURCE/mod_particle_attributes.f90
===================================================================
--- palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 4627)
+++ palm/trunk/SOURCE/mod_particle_attributes.f90	(revision 4628)
@@ -20,9 +20,12 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! extensions required for MPI-I/O of particle data to restart files
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -47,5 +50,10 @@
     USE, INTRINSIC ::  ISO_C_BINDING
 
+    USE control_parameters,                                                                        &
+        ONLY: varnamelength
+
     USE kinds
+
+    CHARACTER(LEN=varnamelength), DIMENSION(50) ::  part_output = ' '  !< namelist parameter
 
     INTEGER(iwp) ::  dissipation_classes = 10                     !< namelist parameter (see documentation)
@@ -54,6 +62,8 @@
     INTEGER(iwp) ::  ibc_par_ns                                   !< particle north/south boundary condition dummy
     INTEGER(iwp) ::  ibc_par_t                                    !< particle top boundary condition dummy
+    INTEGER(iwp) ::  number_of_output_particles = 0               !< number of output particles
     INTEGER(iwp) ::  number_of_particles = 0                      !< number of particles for each grid box (3d array is saved on prt_count)
     INTEGER(iwp) ::  number_of_particle_groups = 1                !< namelist parameter (see documentation)
+    INTEGER(iwp) ::  part_inc = 1                                 !< increment of particles in output file
 
     INTEGER(iwp), PARAMETER ::  max_number_of_particle_groups = 10 !< maximum allowed number of particle groups
@@ -62,43 +72,40 @@
     
     LOGICAL ::  particle_advection = .FALSE.              !< parameter to steer the advection of particles
+    LOGICAL ::  unlimited_dimension = .TRUE.              !< umlimited dimension for particle output
     LOGICAL ::  use_sgs_for_particles = .FALSE.           !< namelist parameter (see documentation)    
     LOGICAL ::  wang_kernel = .FALSE.                     !< flag for collision kernel
 
-    REAL(wp) ::  alloc_factor = 20.0_wp                    !< namelist parameter (see documentation)
-    REAL(wp) ::  particle_advection_start = 0.0_wp         !< namelist parameter (see documentation)
+    REAL(wp) ::  alloc_factor = 20.0_wp                   !< namelist parameter (see documentation)
+    REAL(wp) ::  oversize = 100.0_wp                      !< reserve spare particles in output file (in % relative to initial number)
+    REAL(wp) ::  particle_advection_start = 0.0_wp        !< namelist parameter (see documentation)
+    REAL(wp) ::  part_percent = 100.0_wp                  !< percentage of particles in output file
 
-!
-!-- WARNING: For compatibility of derived types, the BIND attribute is required, and interoperable C
-!-- datatypes must be used. These type are hard wired here! So changes in working precision (wp, iwp)
-!-- will not affect the particle_type!
-!-- The main reason for introducing the interoperable datatypes was to avoid compiler warnings of
-!-- the gfortran compiler.
-!-- The BIND attribite is required because of C_F_POINTER usage in the pmc particle interface.
-    TYPE, BIND(C) ::  particle_type
-        REAL(C_DOUBLE) ::  aux1          !< auxiliary multi-purpose feature
-        REAL(C_DOUBLE) ::  aux2          !< auxiliary multi-purpose feature
-        REAL(C_DOUBLE) ::  radius        !< radius of particle
-        REAL(C_DOUBLE) ::  age           !< age of particle
-        REAL(C_DOUBLE) ::  age_m         !<
-        REAL(C_DOUBLE) ::  dt_sum        !<
-        REAL(C_DOUBLE) ::  e_m           !< interpolated sgs tke
-        REAL(C_DOUBLE) ::  origin_x      !< origin x-position of particle (changed cyclic bc)
-        REAL(C_DOUBLE) ::  origin_y      !< origin y-position of particle (changed cyclic bc)
-        REAL(C_DOUBLE) ::  origin_z      !< origin z-position of particle (changed cyclic bc)
-        REAL(C_DOUBLE) ::  rvar1         !<
-        REAL(C_DOUBLE) ::  rvar2         !<
-        REAL(C_DOUBLE) ::  rvar3         !<
-        REAL(C_DOUBLE) ::  speed_x       !< speed of particle in x
-        REAL(C_DOUBLE) ::  speed_y       !< speed of particle in y
-        REAL(C_DOUBLE) ::  speed_z       !< speed of particle in z
-        REAL(C_DOUBLE) ::  weight_factor !< weighting factor
-        REAL(C_DOUBLE) ::  x             !< x-position
-        REAL(C_DOUBLE) ::  y             !< y-position
-        REAL(C_DOUBLE) ::  z             !< z-position
-        INTEGER(C_INT) ::  class         !< radius class needed for collision
-        INTEGER(C_INT) ::  group         !< number of particle group
-        INTEGER(C_LONG_LONG) ::  id            !< particle ID (64 bit integer)
-        LOGICAL(C_BOOL) ::  particle_mask !< if this parameter is set to false the particle will be deleted
-        INTEGER(C_INT) ::  block_nr      !< number for sorting (removable?)
+    TYPE, PUBLIC ::  particle_type
+        REAL(wp)     ::  aux1          !< auxiliary multi-purpose feature
+        REAL(wp)     ::  aux2          !< auxiliary multi-purpose feature
+        REAL(wp)     ::  radius        !< radius of particle
+        REAL(wp)     ::  age           !< age of particle
+        REAL(wp)     ::  age_m         !<
+        REAL(wp)     ::  dt_sum        !<
+        REAL(wp)     ::  e_m           !< interpolated sgs tke
+        REAL(wp)     ::  origin_x      !< origin x-position of particle (changed cyclic bc)
+        REAL(wp)     ::  origin_y      !< origin y-position of particle (changed cyclic bc)
+        REAL(wp)     ::  origin_z      !< origin z-position of particle (changed cyclic bc)
+        REAL(wp)     ::  rvar1         !<
+        REAL(wp)     ::  rvar2         !<
+        REAL(wp)     ::  rvar3         !<
+        REAL(wp)     ::  speed_x       !< speed of particle in x
+        REAL(wp)     ::  speed_y       !< speed of particle in y
+        REAL(wp)     ::  speed_z       !< speed of particle in z
+        REAL(wp)     ::  weight_factor !< weighting factor
+        REAL(wp)     ::  x             !< x-position
+        REAL(wp)     ::  y             !< y-position
+        REAL(wp)     ::  z             !< z-position
+        INTEGER(iwp) ::  class         !< radius class needed for collision
+        INTEGER(iwp) ::  group         !< number of particle group
+        INTEGER(idp) ::  id            !< particle ID (64 bit integer)
+        LOGICAL      ::  particle_mask !< if this parameter is set to false the particle will be deleted
+        INTEGER(iwp) ::  block_nr      !< number for sorting (removable?)
+        INTEGER(iwp) ::  particle_nr=-1   !< particle number for particle IO (increment one
     END TYPE particle_type
 
Index: palm/trunk/SOURCE/posix_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/posix_interface_mod.f90	(revision 4627)
+++ palm/trunk/SOURCE/posix_interface_mod.f90	(revision 4628)
@@ -19,9 +19,12 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! extensions required for MPI-I/O of particle data to restart files
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! Initial version (K. Ketelsen)
 !
@@ -126,4 +129,6 @@
        MODULE PROCEDURE posix_read_int_1d
        MODULE PROCEDURE posix_read_int_2d
+       MODULE PROCEDURE posix_read_i4_3d
+       MODULE PROCEDURE posix_read_i8_3d
        MODULE PROCEDURE posix_read_offset_1d
        MODULE PROCEDURE posix_read_real_1d
@@ -137,4 +142,6 @@
        MODULE PROCEDURE posix_write_int_1d
        MODULE PROCEDURE posix_write_int_2d
+       MODULE PROCEDURE posix_write_i4_3d
+       MODULE PROCEDURE posix_write_i8_3d
        MODULE PROCEDURE posix_write_offset_1d
        MODULE PROCEDURE posix_write_real_1d
@@ -232,5 +239,5 @@
 
 
- SUBROUTINE posix_read_int_2d (fid, data, nw)
+ SUBROUTINE posix_read_int_2d( fid, data, nw )
 
     IMPLICIT NONE
@@ -250,4 +257,46 @@
 
  END SUBROUTINE posix_read_int_2d
+
+
+
+ SUBROUTINE posix_read_i4_3d( fid, data, nw )
+
+    IMPLICIT NONE
+
+    INTEGER(KIND=isp), INTENT(IN), TARGET, DIMENSION(:,:,:) ::  data         !<
+    INTEGER, INTENT(IN)                                     ::  fid          !<
+    INTEGER                                                 ::  nr_byte      !<
+    INTEGER, INTENT(IN)                                     ::  nw           !<
+
+    TYPE(C_PTR)                                             ::  buf          !<
+
+
+    nr_byte = nw * isp
+    buf     = C_LOC( data )
+
+    CALL posix_read( fid, buf, nr_byte )
+
+ END SUBROUTINE posix_read_i4_3d
+
+
+
+ SUBROUTINE posix_read_i8_3d( fid, data, nw )
+
+    IMPLICIT NONE
+
+    INTEGER(KIND=idp), INTENT(IN), TARGET, DIMENSION(:,:,:) ::  data         !<
+    INTEGER, INTENT(IN)                                     ::  fid          !<
+    INTEGER                                                 ::  nr_byte      !<
+    INTEGER, INTENT(IN)                                     ::  nw           !<
+
+    TYPE(C_PTR)                                             ::  buf          !<
+
+
+    nr_byte = nw * idp
+    buf     = C_LOC( data )
+
+    CALL posix_read( fid, buf, nr_byte )
+
+ END SUBROUTINE posix_read_i8_3d
 
 
@@ -448,4 +497,48 @@
 
 
+ SUBROUTINE posix_write_i4_3d( fid, data, nw )
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)                                     ::  fid        !<
+    INTEGER                                                 ::  nr_byte    !<
+    INTEGER, INTENT(IN)                                     ::  nw         !<
+
+    INTEGER(KIND=isp), INTENT(IN), TARGET, DIMENSION(:,:,:) ::  data       !<
+
+    TYPE(C_PTR)                                             ::  buf        !<
+
+
+    nr_byte = nw * isp
+    buf     = C_LOC( data )
+
+    CALL posix_write( fid, buf, nr_byte )
+
+ END SUBROUTINE posix_write_i4_3d
+
+
+
+ SUBROUTINE posix_write_i8_3d( fid, data, nw )
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)                                     ::  fid        !<
+    INTEGER                                                 ::  nr_byte    !<
+    INTEGER, INTENT(IN)                                     ::  nw         !<
+
+    INTEGER(KIND=idp), INTENT(IN), TARGET, DIMENSION(:,:,:) ::  data       !<
+
+    TYPE(C_PTR)                                             ::  buf        !<
+
+
+    nr_byte = nw * idp
+    buf     = C_LOC( data )
+
+    CALL posix_write( fid, buf, nr_byte )
+
+ END SUBROUTINE posix_write_i8_3d
+
+
+
  SUBROUTINE posix_write_offset_1d( fid, data, nw )
 
Index: palm/trunk/SOURCE/restart_data_mpi_io_mod.f90
===================================================================
--- palm/trunk/SOURCE/restart_data_mpi_io_mod.f90	(revision 4627)
+++ palm/trunk/SOURCE/restart_data_mpi_io_mod.f90	(revision 4628)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! extensions required for MPI-I/O of particle data to restart files
+!
+! 4622 2020-07-23 09:02:23Z raasch
 ! unused variables removed
 !
@@ -65,6 +68,5 @@
 !
 ! 4495 2020-04-13 20:11:20Z raasch
-! Initial version (K. Ketelsen), adjusted to PALM formatting standards (s. Raasch)
-!
+! Initial version (K. Ketelsen), adjusted to PALM formatting standards (S. Raasch)
 !
 !
@@ -127,4 +129,11 @@
     USE kinds
 
+    USE particle_attributes,                                                                       &
+        ONLY:  grid_particles,                                                                     &
+               particles,                                                                          &
+               particle_type,                                                                      &
+               prt_count,                                                                          &
+               zero_particle
+
     USE pegrid,                                                                                    &
         ONLY:  comm1dx,                                                                            &
@@ -161,5 +170,5 @@
 
     INTEGER(iwp)            ::  comm_io          !< Communicator for MPI-IO
-    INTEGER(iwp)            ::  fh               !< MPI-IO file handle
+    INTEGER(iwp)            ::  fh = -1          !< MPI-IO file handle
 #if defined( __parallel )
     INTEGER(iwp)            ::  fhs = -1         !< MPI-IO file handle to open file with comm2d always
@@ -170,4 +179,6 @@
     INTEGER(iwp)            ::  ft_2d            !< MPI filetype 2D array REAL with outer boundaries
     INTEGER(iwp)            ::  ft_3d            !< MPI filetype 3D array REAL with outer boundaries
+    INTEGER(iwp)            ::  ft_3di4 = -1     !< MPI filetype 3D array INTEGER*4
+    INTEGER(iwp)            ::  ft_3di8 = -1     !< MPI filetype 3D array INTEGER*8
     INTEGER(iwp)            ::  ft_3dsoil        !< MPI filetype for 3d-soil array
 #endif
@@ -181,4 +192,6 @@
     INTEGER(iwp)            ::  win_2di          !<
     INTEGER(iwp)            ::  win_2dr          !<
+    INTEGER(iwp)            ::  win_3di4 = -1    !<
+    INTEGER(iwp)            ::  win_3di8 = -1    !<
     INTEGER(iwp)            ::  win_3dr          !<
     INTEGER(iwp)            ::  win_3ds          !<
@@ -190,5 +203,5 @@
     INTEGER(KIND=rd_offset_kind) ::  header_position  !<
 
-    INTEGER(iwp), DIMENSION(:,:), POINTER, CONTIGUOUS ::  array_2di  !<
+    INTEGER(iwp), DIMENSION(:,:), POINTER, CONTIGUOUS   ::  array_2di   !<
 
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_end_index     !<
@@ -196,4 +209,6 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_start_index   !<
 
+    INTEGER(isp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  array_3di4  !<
+    INTEGER(idp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  array_3di8  !<
 
     LOGICAL ::  all_pes_write                 !< all PEs have data to write
@@ -279,4 +294,10 @@
     TYPE(domain_decomposition_grid_features) ::  prerun_grid   !< grid variables for the prerun
 
+!
+!-- MPI_INTEGER8 is not standard MPI, but is supported on most MPI distibutions
+!-- If not suppported, a workaround could be enabled with the following preprocessor directive
+!#if defined( __NO_INTEGER8)
+!    INTEGER ::  MPI_INTEGER8  !< MPI data type INTEGER8
+!#endif
 
     SAVE
@@ -296,4 +317,8 @@
     END INTERFACE rd_mpi_io_close
 
+    INTERFACE rd_mpi_io_check_open
+       MODULE PROCEDURE rd_mpi_io_check_open
+    END INTERFACE rd_mpi_io_check_open
+
     INTERFACE rd_mpi_io_open
        MODULE PROCEDURE rd_mpi_io_open
@@ -304,4 +329,6 @@
        MODULE PROCEDURE rrd_mpi_io_int
        MODULE PROCEDURE rrd_mpi_io_int_2d
+       MODULE PROCEDURE rrd_mpi_io_int4_3d
+       MODULE PROCEDURE rrd_mpi_io_int8_3d
        MODULE PROCEDURE rrd_mpi_io_logical
        MODULE PROCEDURE rrd_mpi_io_real
@@ -324,4 +351,12 @@
     END INTERFACE rrd_mpi_io_surface
 
+    INTERFACE rrd_mpi_io_particles
+       MODULE PROCEDURE rrd_mpi_io_particles
+    END INTERFACE rrd_mpi_io_particles
+
+    INTERFACE rd_mpi_io_particle_filetypes
+       MODULE PROCEDURE rd_mpi_io_particle_filetypes
+    END INTERFACE rd_mpi_io_particle_filetypes
+
     INTERFACE rd_mpi_io_surface_filetypes
        MODULE PROCEDURE rd_mpi_io_surface_filetypes
@@ -332,4 +367,6 @@
        MODULE PROCEDURE wrd_mpi_io_int
        MODULE PROCEDURE wrd_mpi_io_int_2d
+       MODULE PROCEDURE wrd_mpi_io_int4_3d
+       MODULE PROCEDURE wrd_mpi_io_int8_3d
        MODULE PROCEDURE wrd_mpi_io_logical
        MODULE PROCEDURE wrd_mpi_io_real
@@ -347,4 +384,8 @@
     END INTERFACE wrd_mpi_io_global_array
 
+    INTERFACE wrd_mpi_io_particles
+       MODULE PROCEDURE wrd_mpi_io_particles
+    END INTERFACE wrd_mpi_io_particles
+
     INTERFACE wrd_mpi_io_surface
        MODULE PROCEDURE wrd_mpi_io_surface
@@ -353,12 +394,16 @@
 
     PUBLIC  rd_mpi_io_check_array,                                                                 &
+            rd_mpi_io_check_open,                                                                  &
             rd_mpi_io_close,                                                                       &
             rd_mpi_io_open,                                                                        &
+            rd_mpi_io_particle_filetypes,                                                          &
+            rd_mpi_io_surface_filetypes,                                                           &
             rrd_mpi_io,                                                                            &
             rrd_mpi_io_global_array,                                                               &
+            rrd_mpi_io_particles,                                                                  &
             rrd_mpi_io_surface,                                                                    &
-            rd_mpi_io_surface_filetypes,                                                           &
             wrd_mpi_io,                                                                            &
             wrd_mpi_io_global_array,                                                               &
+            wrd_mpi_io_particles,                                                                  &
             wrd_mpi_io_surface
 
@@ -394,6 +439,4 @@
     TYPE(C_PTR)                   ::  buf_ptr  !<
 #endif
-
-!    write(9,*) 'Here is rd_mpi_io_open',nx,nx_on_file,ny,ny_on_file,TRIM(action)   !kk may become Debug Output
 
     offset = 0
@@ -1415,4 +1458,142 @@
 ! Description:
 ! ------------
+!> Read 3d-INTEGER*4 array with MPI-IO
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE rrd_mpi_io_int4_3d( name, data )
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
+
+    INTEGER(iwp)                       ::  i       !<
+
+#if defined( __parallel )
+    INTEGER, DIMENSION(rd_status_size) ::  status  !<
+#endif
+
+    LOGICAL                            ::  found   !<
+
+    INTEGER(isp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data  !<
+
+
+    found = .FALSE.
+    data  = -1.0
+
+    DO  i = 1, tgh%nr_arrays
+       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
+          array_position = array_offset(i)
+          found = .TRUE.
+          EXIT
+       ENDIF
+    ENDDO
+
+    IF ( found )  THEN
+
+#if defined( __parallel )
+       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited # of cores is inactive
+       IF( sm_io%iam_io_pe )  THEN
+          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_3di4, 'native',              &
+                                  MPI_INFO_NULL, ierr )
+          CALL MPI_FILE_READ_ALL( fh, array_3di4, SIZE( array_3di4 ), MPI_INTEGER, status, ierr )
+       ENDIF
+       CALL sm_io%sm_node_barrier()
+#else
+       CALL posix_lseek( fh, array_position )
+       CALL posix_read(fh, array_3di4, SIZE( array_3di4 ) )
+#endif
+       IF ( include_total_domain_boundaries )  THEN
+          DO  i = iog%nxl, iog%nxr
+             data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3di4(:,i,iog%nys:iog%nyn)
+          ENDDO
+       ELSE
+          DO  i = nxl, nxr
+             data(:,nys:nyn,i) = array_3di4(:,i,nys:nyn)
+          ENDDO
+       ENDIF
+
+    ELSE
+
+       message_string = '3d-INTEGER*4 array "' // TRIM( name ) // '" not found in restart file'
+       CALL message( 'rrd_mpi_io_int4_3d', 'PA0722', 3, 2, 0, 6, 0 )
+
+    ENDIF
+
+
+ END SUBROUTINE rrd_mpi_io_int4_3d
+
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Read 3d-INTEGER*8 array with MPI-IO
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE rrd_mpi_io_int8_3d( name, data )
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
+
+    INTEGER(iwp)                       ::  i       !<
+
+#if defined( __parallel )
+    INTEGER, DIMENSION(rd_status_size) ::  status  !<
+#endif
+
+    LOGICAL                            ::  found   !<
+
+    INTEGER(idp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data  !<
+
+
+    found = .FALSE.
+    data  = -1.0
+
+    DO  i = 1, tgh%nr_arrays
+       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
+          array_position = array_offset(i)
+          found = .TRUE.
+          EXIT
+       ENDIF
+    ENDDO
+
+    IF ( found )  THEN
+
+#if defined( __parallel )
+          CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited # of cores is inactive
+          IF( sm_io%iam_io_pe )  THEN
+             CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER8, ft_3di8, 'native', MPI_INFO_NULL, &
+                                     ierr )
+             CALL MPI_FILE_READ_ALL( fh, array_3di8, SIZE( array_3di8 ), MPI_INTEGER8, status, ierr )
+          ENDIF
+          CALL sm_io%sm_node_barrier()
+#else
+          CALL posix_lseek( fh, array_position )
+          CALL posix_read(fh, array_3di8, SIZE( array_3di8 ) )
+#endif
+          IF ( include_total_domain_boundaries )  THEN
+             DO  i = iog%nxl, iog%nxr
+                data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3di8(:,i,iog%nys:iog%nyn)
+             ENDDO
+          ELSE
+             DO  i = nxl, nxr
+                data(:,nys:nyn,i) = array_3di8(:,i,nys:nyn)
+             ENDDO
+          ENDIF
+
+    ELSE
+
+       message_string = '3d-INTEGER*8 array "' // TRIM( name ) // '" not found in restart file'
+       CALL message( 'rrd_mpi_io_int8_3d', 'PA0722', 3, 2, 0, 6, 0 )
+
+    ENDIF
+
+
+ END SUBROUTINE rrd_mpi_io_int8_3d
+
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Read 2d-REAL array with MPI-IO
 !--------------------------------------------------------------------------------------------------!
@@ -1949,4 +2130,140 @@
 ! Description:
 ! ------------
+!> Write 3d-INTEGER*4 array with MPI-IO
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE wrd_mpi_io_int4_3d( name, data )
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
+
+    INTEGER(iwp)                       ::  i       !<
+#if defined( __parallel )
+    INTEGER, DIMENSION(rd_status_size) ::  status  !<
+#endif
+    INTEGER(isp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data !<
+
+
+    IF ( header_array_index == max_nr_arrays )  THEN
+       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
+    ENDIF
+
+    array_names(header_array_index)  = name
+    array_offset(header_array_index) = array_position
+    header_array_index = header_array_index + 1
+
+    IF ( include_total_domain_boundaries )  THEN
+!
+!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
+!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
+!--    index order of the array in the same way, i.e. the first dimension should be along x and the
+!--    second along y. For this reason, the original PALM data need to be swaped.
+       DO  i = iog%nxl, iog%nxr
+          array_3di4(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
+       ENDDO
+
+    ELSE
+!
+!--    Prepare output of 3d-REAL-array without ghost layers
+       DO  i = nxl, nxr
+           array_3di4(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
+       ENDDO
+
+    ENDIF
+#if defined( __parallel )
+    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
+    IF ( sm_io%iam_io_pe )  THEN
+       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_3di4, 'native', MPI_INFO_NULL, ierr )
+       CALL MPI_FILE_WRITE_ALL( fh, array_3di4, SIZE( array_3di4 ), MPI_INTEGER, status, ierr )
+    ENDIF
+    CALL sm_io%sm_node_barrier()
+#else
+    CALL posix_lseek( fh, array_position )
+    CALL posix_write( fh, array_3di4, SIZE( array_3di4 ) )
+#endif
+!
+!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
+!-- Maybe a compiler problem.
+    array_position = array_position + INT(     (nz+2), KIND = rd_offset_kind ) *                   &
+                                      INT( (iog%ny+1), KIND = rd_offset_kind ) *                   &
+                                      INT( (iog%nx+1), KIND = rd_offset_kind ) * isp
+
+    write(9,*) 'array_position int4_3d ',trim(name),' ',array_position
+
+ END SUBROUTINE wrd_mpi_io_int4_3d
+
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write 3d-INTEGER*8 array with MPI-IO
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE wrd_mpi_io_int8_3d( name, data )
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
+
+    INTEGER(iwp)                       ::  i       !<
+#if defined( __parallel )
+    INTEGER, DIMENSION(rd_status_size) ::  status  !<
+#endif
+    INTEGER(idp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data !<
+
+
+    IF ( header_array_index == max_nr_arrays )  THEN
+       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
+    ENDIF
+
+    array_names(header_array_index)  = name
+    array_offset(header_array_index) = array_position
+    header_array_index = header_array_index + 1
+
+    IF ( include_total_domain_boundaries )  THEN
+!
+!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
+!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
+!--    index order of the array in the same way, i.e. the first dimension should be along x and the
+!--    second along y. For this reason, the original PALM data need to be swaped.
+       DO  i = iog%nxl, iog%nxr
+          array_3di8(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
+       ENDDO
+
+    ELSE
+!
+!--    Prepare output of 3d-REAL-array without ghost layers
+       DO  i = nxl, nxr
+           array_3di8(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
+       ENDDO
+
+    ENDIF
+#if defined( __parallel )
+    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
+    IF ( sm_io%iam_io_pe )  THEN
+       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER8, ft_3di8, 'native', MPI_INFO_NULL, ierr )
+       CALL MPI_FILE_WRITE_ALL( fh, array_3di8, SIZE( array_3di8 ), MPI_INTEGER8, status, ierr )
+    ENDIF
+    CALL sm_io%sm_node_barrier()
+#else
+    CALL posix_lseek( fh, array_position )
+    CALL posix_write( fh, array_3di8, SIZE( array_3di8 ) )
+#endif
+!
+!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
+!-- Maybe a compiler problem.
+    array_position = array_position + INT(     (nz+2), KIND = rd_offset_kind ) *                   &
+                                      INT( (iog%ny+1), KIND = rd_offset_kind ) *                   &
+                                      INT( (iog%nx+1), KIND = rd_offset_kind ) * dp
+
+    write(9,*) 'array_position int8_3d ',trim(name),' ',array_position
+
+ END SUBROUTINE wrd_mpi_io_int8_3d
+
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Write 3d-REAL array with MPI-IO
 !--------------------------------------------------------------------------------------------------!
@@ -2007,4 +2324,6 @@
                                       INT( (iog%ny+1), KIND = rd_offset_kind ) *                   &
                                       INT( (iog%nx+1), KIND = rd_offset_kind ) * wp
+
+    write(9,*) 'array_position real3d ',trim(name),' ',array_position
 
  END SUBROUTINE wrd_mpi_io_real_3d
@@ -2576,4 +2895,123 @@
 ! Description:
 ! ------------
+!> Read particle data array with MPI-IO.
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE rrd_mpi_io_particles( name, prt_global_index )
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN)                       ::  name            !<
+    INTEGER(idp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  prt_global_index      !<
+
+    INTEGER(iwp)                       ::  array_size      !<
+    INTEGER(iwp)                       ::  byte_column     !<
+    INTEGER(iwp)                       ::  i               !<
+    INTEGER(iwp)                       ::  ind             !<
+    INTEGER(iwp)                       ::  j               !<
+    INTEGER(iwp)                       ::  k               !<
+    INTEGER(iwp)                       ::  n               !<
+    INTEGER(iwp)                       ::  particle_size   !<
+
+    INTEGER(KIND=rd_offset_kind)       ::  disp            !<
+    INTEGER(KIND=rd_offset_kind)       ::  offset          !<
+    INTEGER(KIND=rd_offset_kind)       ::  prt_nr_bytes    !<
+
+    LOGICAL                            ::  found           !<
+
+    REAL(dp)                           :: rr               !< there is no data type INTEGER*8 in MPI
+    REAL(dp)                           :: rs               !< use REAL*8 to compute max offset
+
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET :: prt_data   !<
+
+#if defined( __parallel )
+    INTEGER, DIMENSION(rd_status_size) ::  status          !<
+#else
+    TYPE(C_PTR)                        ::  buf
+#endif
+
+    found = .FALSE.
+
+    DO  i = 1, tgh%nr_arrays
+       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
+          array_position = array_offset(i)
+          found = .TRUE.
+          EXIT
+       ENDIF
+    ENDDO
+
+    IF ( found )  THEN
+
+       offset = 0
+
+       particle_size = STORAGE_SIZE(zero_particle) / 8  ! 8 here means number of bits/byte, NOT wp
+
+       array_size = 0
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             array_size = MAX( array_size, SUM(prt_count(:,j,i)) )
+          ENDDO
+       ENDDO
+
+       write(9,*) 'particle_size_read ',particle_size,array_size,array_position,sum(prt_global_index)
+
+       ALLOCATE( prt_data(MAX(array_size,1)) )
+
+!
+!--    Write columns of particle
+#if defined( __parallel )
+       CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
+#endif
+       prt_nr_bytes = 0
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             disp         = array_position + prt_global_index(nzb,j,i) * particle_size
+             byte_column  = SUM( prt_count(:,j,i) ) * particle_size
+             prt_nr_bytes = MAX( disp+byte_column, prt_nr_bytes )
+
+#if defined( __parallel )
+             CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
+             IF ( byte_column > 0 )  THEN
+                CALL MPI_FILE_SEEK( fh, disp, MPI_SEEK_SET, ierr )
+                CALL MPI_FILE_READ( fh, prt_data, byte_column, MPI_BYTE, status, ierr )
+             ENDIF
+             CALL sm_io%sm_node_barrier()
+#else
+             buf = C_LOC(prt_data)     ! use C_PTR to avaid another overlay in posix interface
+             CALL posix_lseek( fh, disp )
+             CALL posix_read( fh, buf, byte_column )
+#endif
+             ind = 1
+             DO  k = nzb, nzt+1
+                DO  n = 1, prt_count(k,j,i)
+                   grid_particles(k,j,i)%particles(n) = prt_data(ind)
+                   ind = ind+1
+                ENDDO
+             ENDDO
+
+          ENDDO
+       ENDDO
+
+#if defined( __parallel )
+       rs = prt_nr_bytes
+       CALL MPI_ALLREDUCE( rs, rr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, comm2d, ierr )
+       prt_nr_bytes = rr
+#else
+       rr = rs
+#endif
+       array_position = prt_nr_bytes
+
+       write(9,*) 'array_position after particle read ',array_position,prt_nr_bytes,rs
+
+       DEALLOCATE( prt_data )
+
+    ENDIF
+
+ END SUBROUTINE rrd_mpi_io_particles
+
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Read 1d-REAL surface data array with MPI-IO.
 !--------------------------------------------------------------------------------------------------!
@@ -2599,9 +3037,11 @@
 #if defined( __parallel )
     INTEGER, DIMENSION(rd_status_size)  ::  status   !<
+#else
+    TYPE(C_PTR)                         ::  buf      !<
 #endif
 
     LOGICAL                             ::  found    !<
 
-    REAL(wp), INTENT(OUT), DIMENSION(:) ::  data     !<
+    REAL(wp), INTENT(OUT), DIMENSION(:), TARGET ::  data  !<
 
 
@@ -2677,6 +3117,9 @@
                       ierr )
 #else
+!
+!--                Use C_PTR here, because posix read does not work with indexed array
+                   buf = C_LOC( data(m_start_index(j_f,i_f)) )
                    CALL posix_lseek( fh, disp_f )
-                   CALL posix_read( fh, data(m_start_index(j_f,i_f):), nr_bytes_f )
+                   CALL posix_read( fh, buf, nr_bytes_f )
 #endif
                    disp_f     = disp
@@ -2697,12 +3140,4 @@
 
     ENDIF
-
-!    IF ( lo_first_index == 1 )  THEN
-!       IF ( debug_level >= 2 .AND. nr_val > 0 )  WRITE(9,*)  'r_surf_1 ', TRIM( name ), ' ', nr_val, SUM( data(1:nr_val) )
-!    ELSE
-!       IF ( debug_level >= 2 .AND. nr_val > 0 )  WRITE(9,*)  'r_surf_next ', TRIM( name ), ' ', &
-!                                                             lo_first_index,nr_val, SUM( data(1:nr_val) )
-!    ENDIF
-
 
  CONTAINS
@@ -2928,4 +3363,114 @@
 ! Description:
 ! ------------
+!> Write particle data with MPI-IO.
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE wrd_mpi_io_particles( name, prt_global_index )
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN)                       ::  name            !<
+    INTEGER(idp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  prt_global_index      !<
+
+    INTEGER(iwp)                       ::  array_size      !<
+    INTEGER(iwp)                       ::  byte_column     !<
+    INTEGER(iwp)                       ::  i               !<
+    INTEGER(iwp)                       ::  ind             !<
+    INTEGER(iwp)                       ::  j               !<
+    INTEGER(iwp)                       ::  k               !<
+    INTEGER(iwp)                       ::  n               !<
+    INTEGER(iwp)                       ::  particle_size   !<
+
+    INTEGER(KIND=rd_offset_kind)       ::  disp            !<
+    INTEGER(KIND=rd_offset_kind)       ::  offset          !<
+    INTEGER(KIND=rd_offset_kind)       ::  prt_nr_bytes    !<
+
+    REAL(dp)                           :: rs               !< use REAL*8 to compute max offset
+    REAL(dp)                           :: rr               !< there is no data type INTEGER*8 in MPI
+
+
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET ::  prt_data   !<
+
+#if defined( __parallel )
+    INTEGER, DIMENSION(rd_status_size) ::  status          !<
+#else
+    TYPE(C_PTR)                        ::  buf
+#endif
+
+
+    offset = 0
+
+    array_names(header_array_index)  = TRIM(name)
+    array_offset(header_array_index) = array_position
+    header_array_index = header_array_index + 1
+
+    particle_size = STORAGE_SIZE( zero_particle ) / 8
+
+    array_size = 0
+    DO  i = nxl, nxr
+      DO  j = nys, nyn
+         array_size = MAX( array_size, SUM(prt_count(:,j,i)) )
+       ENDDO
+    ENDDO
+
+    ALLOCATE( prt_data(MAX(array_size,1)) )
+
+!
+!-- Write columns of particles.
+!-- The particles of a column are stored sequentially in the first dimension of the particles array.
+!-- Store only the particles of one cell would be possible with this setup, but the I/O portions
+!-- for a maximum of 100 particles are not big enough.
+#if defined( __parallel )
+    CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
+#endif
+    prt_nr_bytes = 0
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
+          disp         = array_position + prt_global_index(nzb,j,i) * particle_size
+          byte_column  = SUM( prt_count(:,j,i) ) * particle_size
+          prt_nr_bytes = MAX( disp+byte_column, prt_nr_bytes )
+
+          ind = 1
+          DO  k = nzb, nzt+1
+             DO  n = 1, prt_count(k,j,i)
+                prt_data(ind) = grid_particles(k,j,i)%particles(n)
+                ind = ind+1
+             ENDDO
+          ENDDO
+
+#if defined( __parallel )
+          CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
+          IF ( byte_column > 0 )  THEN
+             CALL MPI_FILE_SEEK( fh, disp, MPI_SEEK_SET, ierr )
+             CALL MPI_FILE_WRITE( fh, prt_data, byte_column, MPI_BYTE, status, ierr )
+          ENDIF
+          CALL sm_io%sm_node_barrier()
+#else
+          buf = C_LOC(prt_data)  ! use C_PTR to avoid another overlay in posix interface
+          CALL posix_lseek( fh, disp )
+          CALL posix_write( fh, buf, byte_column )
+#endif
+       ENDDO
+    ENDDO
+
+#if defined( __parallel )
+    rs = prt_nr_bytes
+    CALL MPI_ALLREDUCE( rs, rr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, comm2d, ierr )
+    prt_nr_bytes = rr
+#else
+    rr = rs
+#endif
+    array_position = prt_nr_bytes
+
+    write(9,*) 'array_position after particle ',array_position,prt_nr_bytes,rs
+
+    DEALLOCATE( prt_data )
+
+ END SUBROUTINE wrd_mpi_io_particles
+
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Write 1d-REAL surface data array with MPI-IO.
 !--------------------------------------------------------------------------------------------------!
@@ -3043,4 +3588,5 @@
 
 
+
 !--------------------------------------------------------------------------------------------------!
 ! Description:
@@ -3207,7 +3753,19 @@
     ENDIF
 
+    fh = -1
+
     restart_file_size = array_position / ( 1024.0_dp * 1024.0_dp )
 
  END SUBROUTINE rd_mpi_io_close
+
+
+
+ FUNCTION rd_mpi_io_check_open()  RESULT ( isopen )
+
+    LOGICAL ::  isopen
+
+    isopen = ( fh /= -1 )
+
+ END FUNCTION rd_mpi_io_check_open
 
 
@@ -3567,4 +4125,125 @@
 ! Description:
 ! ------------
+!> This subroutine creates file types to access 3d-INTEGER*4 arrays and 3d-INTEGER*8 arrays
+!> distributed in blocks among processes to a single file that contains the global arrays.
+!> These types are only used for particle data.
+!--------------------------------------------------------------------------------------------------!
+  SUBROUTINE rd_mpi_io_particle_filetypes
+
+    IMPLICIT NONE
+
+    INTEGER, DIMENSION(3) ::  dims3   !<
+    INTEGER, DIMENSION(3) ::  lize3   !<
+    INTEGER, DIMENSION(3) ::  start3  !<
+
+    TYPE(domain_decomposition_grid_features) ::  save_io_grid  !< temporary variable to store grid settings
+
+!
+!-- MPI_INTEGER8 is not standard MPI, but is supported on most MPI distibutions
+!-- If not suppported, a workaround could be enabled with the following preprocessor directive
+!#if defined( __NO_INTEGER8)
+!    CALL MPI_TYPE_CONTIGUOUS( 2, MPI_INTEGER, MPI_INTEGER8, ierr )
+!    CALL MPI_TYPE_COMMIT( MPI_INTEGER8, ierr )
+!#endif
+
+    IF ( sm_io%is_sm_active() )  THEN
+       save_io_grid = sm_io%io_grid
+    ENDIF
+
+    IF( .NOT. pe_active_for_read )  RETURN
+
+!
+!-- Decision, if storage with or without ghost layers.
+!-- Please note that the indexing of the global array always starts at 0, even in Fortran.
+!-- Therefore the PE local indices have to be shifted by nbgp in the case with ghost layers.
+    IF ( include_total_domain_boundaries )  THEN
+
+       iog%nxl = nxl + nbgp
+       iog%nxr = nxr + nbgp
+       iog%nys = nys + nbgp
+       iog%nyn = nyn + nbgp
+       iog%nnx = nnx
+       iog%nny = nny
+       iog%nx  = nx + 2 * nbgp
+       iog%ny  = ny + 2 * nbgp
+       IF ( myidx == 0 )  THEN
+          iog%nxl = iog%nxl - nbgp
+          iog%nnx = iog%nnx + nbgp
+       ENDIF
+       IF ( myidx == npex-1  .OR.  npex == -1 )  THEN   ! npex == 1 if -D__parallel not set
+          iog%nxr = iog%nxr + nbgp
+          iog%nnx = iog%nnx + nbgp
+       ENDIF
+       IF ( myidy == 0 )  THEN
+          iog%nys = iog%nys - nbgp
+          iog%nny = iog%nny + nbgp
+       ENDIF
+       IF ( myidy == npey-1  .OR.  npey == -1 )  THEN   ! npey == 1 if -D__parallel not set
+          iog%nyn = iog%nyn + nbgp
+          iog%nny = iog%nny + nbgp
+       ENDIF
+
+       CALL sm_io%sm_adjust_outer_boundary()
+
+    ELSE
+
+       iog%nxl = nxl
+       iog%nxr = nxr
+       iog%nys = nys
+       iog%nyn = nyn
+       iog%nnx = nnx
+       iog%nny = nny
+       iog%nx  = nx
+       iog%ny  = ny
+
+    ENDIF
+
+    IF ( sm_io%is_sm_active() )  THEN
+#if defined( __parallel )
+       CALL sm_io%sm_allocate_shared( array_3di4, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr,&
+                                      sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3di4 )
+       CALL sm_io%sm_allocate_shared( array_3di8, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr,&
+                                      sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3di8 )
+#endif
+    ELSE
+       ALLOCATE( array_3di4(nzb:nzt+1,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
+       ALLOCATE( array_3di8(nzb:nzt+1,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
+
+       sm_io%io_grid = iog
+    ENDIF
+
+!
+!-- Create filetype for 3d INTEGER array
+    dims3(1)  = nz + 2
+    dims3(2)  = iog%nx + 1
+    dims3(3)  = iog%ny + 1
+
+    lize3(1)  = dims3(1)
+    lize3(2)  = sm_io%io_grid%nnx
+    lize3(3)  = sm_io%io_grid%nny
+
+    start3(1) = nzb
+    start3(2) = sm_io%io_grid%nxl
+    start3(3) = sm_io%io_grid%nys
+
+#if defined( __parallel )
+    IF ( sm_io%iam_io_pe )  THEN
+       CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_INTEGER,     &
+                                      ft_3di4, ierr )
+       CALL MPI_TYPE_COMMIT( ft_3di4, ierr )
+
+       CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_INTEGER8,    &
+                                      ft_3di8, ierr )
+       CALL MPI_TYPE_COMMIT( ft_3di8, ierr )
+    ENDIF
+#endif
+
+ END SUBROUTINE rd_mpi_io_particle_filetypes
+
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> This subroutine creates file types to access 3d-soil arrays distributed in blocks among processes
 !> to a single file that contains the global arrays. It is not required for the serial mode.
@@ -3655,8 +4334,24 @@
     ENDIF
 
+!
+!-- Free last particle filetypes
+    IF ( sm_io%iam_io_pe .AND. ft_3di4 /= -1 )  THEN
+       CALL MPI_TYPE_FREE( ft_3di4, ierr )
+       CALL MPI_TYPE_FREE( ft_3di8, ierr )
+    ENDIF
+
+    IF ( sm_io%is_sm_active() .AND.  win_3di4 /= -1 )  THEN
+       CALL sm_io%sm_free_shared( win_3di4 )
+       CALL sm_io%sm_free_shared( win_3di8 )
+    ENDIF
+
     ft_surf  = -1
     win_surf = -1
 #else
-    DEALLOCATE( array_2d, array_2di, array_3d )
+    IF ( ASSOCIATED(array_2d)   )  DEALLOCATE( array_2d )
+    IF ( ASSOCIATED(array_2di)  )  DEALLOCATE( array_2di )
+    IF ( ASSOCIATED(array_3d)   )  DEALLOCATE( array_3d )
+    IF ( ASSOCIATED(array_3di4) )  DEALLOCATE( array_3di4 )
+    IF ( ASSOCIATED(array_3di8) )  DEALLOCATE( array_3di8 )
 #endif
 
Index: palm/trunk/SOURCE/shared_memory_io_mod.f90
===================================================================
--- palm/trunk/SOURCE/shared_memory_io_mod.f90	(revision 4627)
+++ palm/trunk/SOURCE/shared_memory_io_mod.f90	(revision 4628)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! extensions required for MPI-I/O of particle data to restart files
+!
+! 4620 2020-07-22 14:11:16Z raasch
 ! bugfix: variable definition changed
 !
@@ -85,4 +88,6 @@
     USE kinds,                                                                                     &
         ONLY: dp,                                                                                  &
+              idp,                                                                                 &
+              isp,                                                                                 &
               iwp,                                                                                 &
               sp,                                                                                  &
@@ -192,10 +197,13 @@
           PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_3d_64
           PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_3d_32
+          PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_3di_32
+          PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_3di_64
 
           GENERIC, PUBLIC ::  sm_allocate_shared =>                                                &
-                                               sm_allocate_shared_1d_64, sm_allocate_shared_1d_32, &
-                                               sm_allocate_shared_2d_64, sm_allocate_shared_2d_32, &
-                                               sm_allocate_shared_2di,   sm_allocate_shared_3d_64, &
-                                               sm_allocate_shared_3d_32, sm_allocate_shared_1di
+                                              sm_allocate_shared_1d_64,  sm_allocate_shared_1d_32, &
+                                              sm_allocate_shared_2d_64,  sm_allocate_shared_2d_32, &
+                                              sm_allocate_shared_2di,    sm_allocate_shared_3d_64, &
+                                              sm_allocate_shared_3d_32,  sm_allocate_shared_1di,   &
+                                              sm_allocate_shared_3di_32, sm_allocate_shared_3di_64
 #endif
     END TYPE sm_class
@@ -1030,4 +1038,130 @@
  END SUBROUTINE sm_allocate_shared_3d_32
 
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Allocate shared 3d-REAL (32 bit) array on ALL threads
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE sm_allocate_shared_3di_32( this, p3, d1s, d1e, d2s, d2e, d3s, d3e, win )
+
+    IMPLICIT NONE
+
+    CLASS(sm_class), INTENT(inout)      ::  this
+
+    INTEGER                             ::  disp_unit
+    INTEGER, INTENT(IN)                 ::  d1e
+    INTEGER, INTENT(IN)                 ::  d1s
+    INTEGER, INTENT(IN)                 ::  d2e
+    INTEGER, INTENT(IN)                 ::  d2s
+    INTEGER, INTENT(IN)                 ::  d3e
+    INTEGER, INTENT(IN)                 ::  d3s
+    INTEGER, SAVE                       ::  pe_from = 0
+    INTEGER, INTENT(OUT)                ::  win
+
+    INTEGER(KIND=MPI_ADDRESS_KIND)      ::  rem_size
+    INTEGER(KIND=MPI_ADDRESS_KIND)      ::  wsize
+
+    INTEGER, DIMENSION(3)               ::  buf_shape
+
+    INTEGER(isp), DIMENSION(:,:,:), POINTER ::  buf
+    INTEGER(isp), DIMENSION(:,:,:), POINTER ::  p3
+
+    TYPE(C_PTR), SAVE                   ::  base_ptr
+    TYPE(C_PTR), SAVE                   ::  rem_ptr
+
+
+    IF ( this%no_shared_memory_in_this_run )  RETURN
+!
+!-- Allocate shared memory on node rank 0 threads.
+    IF ( this%sh_rank == pe_from )  THEN
+       wsize = ( d3e - d3s + 1 ) * ( d2e - d2s + 1 ) * ( d1e - d1s + 1 )
+    ELSE
+       wsize = 1
+    ENDIF
+
+    wsize = wsize * isp ! Please note, size is always in bytes, independently of the displacement
+                       ! unit
+
+    CALL MPI_WIN_ALLOCATE_SHARED( wsize, isp, MPI_INFO_NULL, this%comm_shared, base_ptr, win, ierr )
+!
+!-- Get C-pointer of the memory located on node-rank pe_from (sh_rank == pe_from)
+    CALL MPI_WIN_SHARED_QUERY( win, pe_from, rem_size, disp_unit, rem_ptr, ierr )
+!
+!-- Convert C- to Fortran-pointer
+    buf_shape(3) = d3e - d3s + 1
+    buf_shape(2) = d2e - d2s + 1
+    buf_shape(1) = d1e - d1s + 1
+    CALL C_F_POINTER( rem_ptr, buf, buf_shape )
+    p3(d1s:,d2s:,d3s:) => buf
+!
+!-- Allocate shared memory in round robin on all PEs of a node.
+    pe_from = MOD( pe_from, this%sh_npes )
+
+ END SUBROUTINE sm_allocate_shared_3di_32
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Allocate shared 3d-REAL (64 bit) array on ALL threads
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE sm_allocate_shared_3di_64( this, p3, d1s, d1e, d2s, d2e, d3s, d3e, win )
+
+    IMPLICIT NONE
+
+    CLASS(sm_class), INTENT(inout)      ::  this         !<
+
+    INTEGER                             ::  disp_unit    !<
+    INTEGER, INTENT(IN)                 ::  d1e          !<
+    INTEGER, INTENT(IN)                 ::  d1s          !<
+    INTEGER, INTENT(IN)                 ::  d2e          !<
+    INTEGER, INTENT(IN)                 ::  d2s          !<
+    INTEGER, INTENT(IN)                 ::  d3e          !<
+    INTEGER, INTENT(IN)                 ::  d3s          !<
+    INTEGER, SAVE                       ::  pe_from = 0  !<
+    INTEGER, INTENT(OUT)                ::  win          !<
+
+    INTEGER(KIND=MPI_ADDRESS_KIND)      ::  rem_size     !<
+    INTEGER(KIND=MPI_ADDRESS_KIND)      ::  wsize        !<
+
+    INTEGER, DIMENSION(3)               ::  buf_shape    !<
+
+    INTEGER(idp), DIMENSION(:,:,:), POINTER ::  buf          !<
+    INTEGER(idp), DIMENSION(:,:,:), POINTER ::  p3           !<
+
+    TYPE(C_PTR), SAVE                   ::  base_ptr     !<
+    TYPE(C_PTR), SAVE                   ::  rem_ptr      !<
+
+
+    IF ( this%no_shared_memory_in_this_run )  RETURN
+!
+!-- Allocate shared memory on node rank 0 threads.
+    IF ( this%sh_rank == pe_from )  THEN
+       wsize = ( d3e - d3s + 1 ) * ( d2e - d2s + 1 ) * ( d1e - d1s + 1 )
+    ELSE
+       wsize = 1
+    ENDIF
+
+    wsize = wsize * idp ! Please note, size is always in bytes, independently of the displacement
+                        ! unit
+
+    CALL MPI_WIN_ALLOCATE_SHARED( wsize, idp, MPI_INFO_NULL, this%comm_shared, base_ptr, win, ierr )
+!
+!-- Get C-pointer of the memory located on node-rank pe_from (sh_rank == pe_from)
+    CALL MPI_WIN_SHARED_QUERY( win, pe_from, rem_size, disp_unit, rem_ptr, ierr )
+!
+!-- Convert C- to Fortran-pointer
+    buf_shape(3) = d3e - d3s + 1
+    buf_shape(2) = d2e - d2s + 1
+    buf_shape(1) = d1e - d1s + 1
+    CALL C_F_POINTER( rem_ptr, buf, buf_shape )
+    p3(d1s:,d2s:,d3s:) => buf
+!
+!-- Allocate shared memory in round robin on all PEs of a node.
+    pe_from = MOD( pe_from, this%sh_npes )
+
+ END SUBROUTINE sm_allocate_shared_3di_64
+
 #endif
 
