Index: palm/trunk/SOURCE/lpm_exchange_horiz.f90
===================================================================
--- palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 1329)
+++ palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 1359)
@@ -1,3 +1,3 @@
- SUBROUTINE lpm_exchange_horiz
+ MODULE lpm_exchange_horiz_mod
 
 !--------------------------------------------------------------------------------!
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! New particle structure integrated. 
+! Kind definition added to all floating point numbers.
 ! 
 ! Former revisions:
@@ -56,5 +57,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  message_string, netcdf_data_format
+        ONLY:  dz, message_string, netcdf_data_format, simulated_time
 
     USE cpulog,                                                                &
@@ -65,17 +66,22 @@
 
     USE indices,                                                               &
-        ONLY:  nx, nxl, nxr, ny, nyn, nys
+        ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb, nzt
 
     USE kinds
 
+    USE lpm_pack_arrays_mod,                                                   &
+        ONLY:  lpm_pack_arrays
+
     USE particle_attributes,                                                   &
-        ONLY:  deleted_particles, deleted_tails, ibc_par_lr, ibc_par_ns,       &
-               maximum_number_of_particles, maximum_number_of_tails,           &
-               maximum_number_of_tailpoints, mpi_particle_type,                &
-               number_of_tails, number_of_particles, particles, particle_mask, &
-               particle_tail_coordinates, particle_type, tail_mask,            &
-               trlp_count_sum, trlp_count_recv_sum, trnp_count_sum,            &
-               trnp_count_recv_sum, trrp_count_sum, trrp_count_recv_sum,       &
-               trsp_count_sum, trsp_count_recv_sum, use_particle_tails
+        ONLY:  alloc_factor, deleted_particles, deleted_tails, grid_particles, &
+               ibc_par_lr, ibc_par_ns, maximum_number_of_tails,                &
+               maximum_number_of_tailpoints, min_nr_particle,                  &
+               mpi_particle_type, number_of_tails, number_of_particles,        &
+               offset_ocean_nzt, offset_ocean_nzt_m1, particles,               &
+               particle_tail_coordinates, particle_type, prt_count,            &
+               tail_mask, trlp_count_sum,                                      &
+               trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum,       &
+               trrp_count_sum, trrp_count_recv_sum, trsp_count_sum,            &
+               trsp_count_recv_sum, use_particle_tails, zero_particle
 
     USE pegrid
@@ -83,6 +89,40 @@
     IMPLICIT NONE
 
+    INTEGER(iwp), PARAMETER ::  NR_2_direction_move = 10000 !:
+    INTEGER(iwp)            ::  nr_move_north               !:
+    INTEGER(iwp)            ::  nr_move_south               !:
+
+    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_north
+    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_south
+
+    SAVE
+
+    PRIVATE
+    PUBLIC lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
+
+    INTERFACE lpm_exchange_horiz
+       MODULE PROCEDURE lpm_exchange_horiz
+    END INTERFACE lpm_exchange_horiz
+
+    INTERFACE lpm_move_particle
+       MODULE PROCEDURE lpm_move_particle
+    END INTERFACE lpm_move_particle
+
+    INTERFACE realloc_particles_array
+       MODULE PROCEDURE realloc_particles_array
+    END INTERFACE realloc_particles_array
+
+CONTAINS
+
+ SUBROUTINE lpm_exchange_horiz
+
+    IMPLICIT NONE
+
     INTEGER(iwp) ::  i                                       !:
+    INTEGER(iwp) ::  ip                                      !:
     INTEGER(iwp) ::  j                                       !:
+    INTEGER(iwp) ::  jp                                      !:
+    INTEGER(iwp) ::  k                                       !:
+    INTEGER(iwp) ::  kp                                      !:
     INTEGER(iwp) ::  n                                       !:
     INTEGER(iwp) ::  nn                                      !:
@@ -110,9 +150,14 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trspt        !:
 
-
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !:
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !:
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !:
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !:
     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !:
     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !:
     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !:
     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !:
+
+    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
 
 #if defined( __parallel )
@@ -141,18 +186,32 @@
     IF ( pdims(1) /= 1 )  THEN
 !
-!--    First calculate the storage necessary for sending and receiving the data
-       DO  n = 1, number_of_particles
-          i = ( particles(n)%x + 0.5 * dx ) * ddx
-!
-!--       Above calculation does not work for indices less than zero
-          IF ( particles(n)%x < -0.5 * dx )  i = -1
-
-          IF ( i < nxl )  THEN
-             trlp_count = trlp_count + 1
-             IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
-          ELSEIF ( i > nxr )  THEN
-             trrp_count = trrp_count + 1
-             IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
-          ENDIF
+!--    First calculate the storage necessary for sending and receiving the data.
+!--    Compute only first (nxl) and last (nxr) loop iterration.
+       DO  ip = nxl, nxr, nxr - nxl
+          DO  jp = nys, nyn
+             DO  kp = nzb+1, nzt
+
+                number_of_particles = prt_count(kp,jp,ip)
+                IF ( number_of_particles <= 0 )  CYCLE
+                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                DO  n = 1, number_of_particles
+                   IF ( particles(n)%particle_mask )  THEN
+                      i = ( particles(n)%x + 0.5_wp * dx ) * ddx
+   !
+   !--                Above calculation does not work for indices less than zero
+                      IF ( particles(n)%x < -0.5_wp * dx )  i = -1
+
+                      IF ( i < nxl )  THEN
+                         trlp_count = trlp_count + 1
+                         IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
+                      ELSEIF ( i > nxr )  THEN
+                         trrp_count = trrp_count + 1
+                         IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
+                      ENDIF
+                   ENDIF
+                ENDDO
+
+             ENDDO
+          ENDDO
        ENDDO
 
@@ -164,10 +223,6 @@
        ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
 
-       trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0, 0, 0, 0 )
-       trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0, 0, 0, 0 )
+       trlp = zero_particle
+       trrp = zero_particle
 
        IF ( use_particle_tails )  THEN
@@ -183,159 +238,174 @@
 
     ENDIF
-
-    DO  n = 1, number_of_particles
-
-       nn = particles(n)%tail_id
-
-       i = ( particles(n)%x + 0.5 * dx ) * ddx
-!
-!--    Above calculation does not work for indices less than zero
-       IF ( particles(n)%x < - 0.5 * dx )  i = -1
-
-       IF ( i <  nxl )  THEN
-          IF ( i < 0 )  THEN
-!
-!--          Apply boundary condition along x
-             IF ( ibc_par_lr == 0 )  THEN
-!
-!--             Cyclic condition
-                IF ( pdims(1) == 1 )  THEN
-                   particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
-                   particles(n)%origin_x = ( nx + 1 ) * dx + &
-                                           particles(n)%origin_x
-                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                      i  = particles(n)%tailpoints
-                      particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
-                                        + particle_tail_coordinates(1:i,1,nn)
-                   ENDIF
-                ELSE
-                   trlp_count = trlp_count + 1
-                   trlp(trlp_count)   = particles(n)
-                   trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
-                   trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
-                                               ( nx + 1 ) * dx
-                   particle_mask(n)  = .FALSE.
-                   deleted_particles = deleted_particles + 1
-
-                   IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.0E-12 ) THEN
-                      trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10
-!++ why is 1 subtracted in next statement???
-                      trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x - 1
-                   ENDIF
-
-                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                      trlpt_count = trlpt_count + 1
-                      trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
-                      trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
-                                               trlpt(:,1,trlpt_count)
-                      tail_mask(nn) = .FALSE.
-                      deleted_tails = deleted_tails + 1
+!
+!-- Compute only first (nxl) and last (nxr) loop iterration
+    DO  ip = nxl, nxr,nxr-nxl
+       DO  jp = nys, nyn
+          DO  kp = nzb+1, nzt
+             number_of_particles = prt_count(kp,jp,ip)
+             IF ( number_of_particles <= 0 ) CYCLE
+             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+             DO  n = 1, number_of_particles
+
+                nn = particles(n)%tail_id
+!
+!--             Only those particles that have not been marked as 'deleted' may
+!--             be moved.
+                IF ( particles(n)%particle_mask )  THEN
+
+                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
+   !
+   !--             Above calculation does not work for indices less than zero
+                   IF ( particles(n)%x < - 0.5_wp * dx )  i = -1
+
+                   IF ( i <  nxl )  THEN
+                      IF ( i < 0 )  THEN
+   !
+   !--                   Apply boundary condition along x
+                         IF ( ibc_par_lr == 0 )  THEN
+   !
+   !--                      Cyclic condition
+                            IF ( pdims(1) == 1 )  THEN
+                               particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
+                               particles(n)%origin_x = ( nx + 1 ) * dx + &
+                               particles(n)%origin_x
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  i  = particles(n)%tailpoints
+                                  particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
+                                  + particle_tail_coordinates(1:i,1,nn)
+                               ENDIF
+                            ELSE
+                               trlp_count = trlp_count + 1
+                               trlp(trlp_count)   = particles(n)
+                               trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
+                               trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
+                               ( nx + 1 ) * dx
+                               particles(n)%particle_mask  = .FALSE.
+                               deleted_particles = deleted_particles + 1
+
+                               IF ( trlp(trlp_count)%x >= (nx + 0.5_wp)* dx - 1.0E-12_wp )  THEN
+                                  trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10_wp
+                                  !++ why is 1 subtracted in next statement???
+                                  trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x - 1
+                               ENDIF
+
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  trlpt_count = trlpt_count + 1
+                                  trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
+                                  trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
+                                  trlpt(:,1,trlpt_count)
+                                  tail_mask(nn) = .FALSE.
+                                  deleted_tails = deleted_tails + 1
+                               ENDIF
+                            ENDIF
+
+                         ELSEIF ( ibc_par_lr == 1 )  THEN
+   !
+   !--                      Particle absorption
+                            particles(n)%particle_mask = .FALSE.
+                            deleted_particles = deleted_particles + 1
+                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                               tail_mask(nn) = .FALSE.
+                               deleted_tails = deleted_tails + 1
+                            ENDIF
+
+                         ELSEIF ( ibc_par_lr == 2 )  THEN
+   !
+   !--                      Particle reflection
+                            particles(n)%x       = -particles(n)%x
+                            particles(n)%speed_x = -particles(n)%speed_x
+
+                         ENDIF
+                      ELSE
+   !
+   !--                   Store particle data in the transfer array, which will be 
+   !--                   send to the neighbouring PE
+                         trlp_count = trlp_count + 1
+                         trlp(trlp_count) = particles(n)
+                         particles(n)%particle_mask = .FALSE.
+                         deleted_particles = deleted_particles + 1
+
+                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                            trlpt_count = trlpt_count + 1
+                            trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
+                            tail_mask(nn) = .FALSE.
+                            deleted_tails = deleted_tails + 1
+                         ENDIF
+                      ENDIF
+
+                   ELSEIF ( i > nxr )  THEN
+                      IF ( i > nx )  THEN
+   !
+   !--                   Apply boundary condition along x
+                         IF ( ibc_par_lr == 0 )  THEN
+   !
+   !--                      Cyclic condition
+                            IF ( pdims(1) == 1 )  THEN
+                               particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
+                               particles(n)%origin_x = particles(n)%origin_x - &
+                               ( nx + 1 ) * dx
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  i = particles(n)%tailpoints
+                                  particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
+                                  + particle_tail_coordinates(1:i,1,nn)
+                               ENDIF
+                            ELSE
+                               trrp_count = trrp_count + 1
+                               trrp(trrp_count) = particles(n)
+                               trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
+                               trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
+                               ( nx + 1 ) * dx
+                               particles(n)%particle_mask = .FALSE.
+                               deleted_particles = deleted_particles + 1
+
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  trrpt_count = trrpt_count + 1
+                                  trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
+                                  trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
+                                  ( nx + 1 ) * dx
+                                  tail_mask(nn) = .FALSE.
+                                  deleted_tails = deleted_tails + 1
+                               ENDIF
+                            ENDIF
+
+                         ELSEIF ( ibc_par_lr == 1 )  THEN
+   !
+   !--                      Particle absorption
+                            particles(n)%particle_mask = .FALSE.
+                            deleted_particles = deleted_particles + 1
+                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                               tail_mask(nn) = .FALSE.
+                               deleted_tails = deleted_tails + 1
+                            ENDIF
+
+                         ELSEIF ( ibc_par_lr == 2 )  THEN
+   !
+   !--                      Particle reflection
+                            particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
+                            particles(n)%speed_x = -particles(n)%speed_x
+
+                         ENDIF
+                      ELSE
+   !
+   !--                   Store particle data in the transfer array, which will be send
+   !--                   to the neighbouring PE
+                         trrp_count = trrp_count + 1
+                         trrp(trrp_count) = particles(n)
+                         particles(n)%particle_mask = .FALSE.
+                         deleted_particles = deleted_particles + 1
+
+                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                            trrpt_count = trrpt_count + 1
+                            trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
+                            tail_mask(nn) = .FALSE.
+                            deleted_tails = deleted_tails + 1
+                         ENDIF
+                      ENDIF
+
                    ENDIF
                 ENDIF
-
-             ELSEIF ( ibc_par_lr == 1 )  THEN
-!
-!--             Particle absorption
-                particle_mask(n) = .FALSE.
-                deleted_particles = deleted_particles + 1
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   tail_mask(nn) = .FALSE.
-                   deleted_tails = deleted_tails + 1
-                ENDIF
-
-             ELSEIF ( ibc_par_lr == 2 )  THEN
-!
-!--             Particle reflection
-                particles(n)%x       = -particles(n)%x
-                particles(n)%speed_x = -particles(n)%speed_x
-
-             ENDIF
-          ELSE
-!
-!--          Store particle data in the transfer array, which will be send
-!--          to the neighbouring PE
-             trlp_count = trlp_count + 1
-             trlp(trlp_count) = particles(n)
-             particle_mask(n) = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                trlpt_count = trlpt_count + 1
-                trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
-                tail_mask(nn) = .FALSE.
-                deleted_tails = deleted_tails + 1
-             ENDIF
-         ENDIF
-
-       ELSEIF ( i > nxr )  THEN
-          IF ( i > nx )  THEN
-!
-!--          Apply boundary condition along x
-             IF ( ibc_par_lr == 0 )  THEN
-!
-!--             Cyclic condition
-                IF ( pdims(1) == 1 )  THEN
-                   particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
-                   particles(n)%origin_x = particles(n)%origin_x - &
-                                           ( nx + 1 ) * dx
-                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                      i = particles(n)%tailpoints
-                      particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
-                                           + particle_tail_coordinates(1:i,1,nn)
-                   ENDIF
-                ELSE
-                   trrp_count = trrp_count + 1
-                   trrp(trrp_count) = particles(n)
-                   trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
-                   trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
-                                               ( nx + 1 ) * dx
-                   particle_mask(n) = .FALSE.
-                   deleted_particles = deleted_particles + 1
-
-                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                      trrpt_count = trrpt_count + 1
-                      trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
-                      trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
-                                               ( nx + 1 ) * dx
-                      tail_mask(nn) = .FALSE.
-                      deleted_tails = deleted_tails + 1
-                   ENDIF
-                ENDIF
-
-             ELSEIF ( ibc_par_lr == 1 )  THEN
-!
-!--             Particle absorption
-                particle_mask(n) = .FALSE.
-                deleted_particles = deleted_particles + 1
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   tail_mask(nn) = .FALSE.
-                   deleted_tails = deleted_tails + 1
-                ENDIF
-
-             ELSEIF ( ibc_par_lr == 2 )  THEN
-!
-!--             Particle reflection
-                particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
-                particles(n)%speed_x = -particles(n)%speed_x
-
-             ENDIF
-          ELSE
-!
-!--          Store particle data in the transfer array, which will be send
-!--          to the neighbouring PE
-             trrp_count = trrp_count + 1
-             trrp(trrp_count) = particles(n)
-             particle_mask(n) = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                trrpt_count = trrpt_count + 1
-                trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
-                tail_mask(nn) = .FALSE.
-                deleted_tails = deleted_tails + 1
-             ENDIF
-          ENDIF
-
-       ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
     ENDDO
 
@@ -345,27 +415,18 @@
     IF ( pdims(1) /= 1 )  THEN
 
-       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'start' )
        CALL MPI_SENDRECV( trlp_count,      1, MPI_INTEGER, pleft,  0, &
                           trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
                           comm2d, status, ierr )
 
-       IF ( number_of_particles + trrp_count_recv > &
-            maximum_number_of_particles )           &
-       THEN
-          IF ( netcdf_data_format < 3 )  THEN
-              message_string = 'maximum_number_of_particles ' //    &
-                               'needs to be increased ' //          &
-                               '&but this is not allowed with ' //  &
-                               'netcdf-data_format < 3'
-             CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 )
-          ELSE
-             CALL lpm_extend_particle_array( trrp_count_recv )
-          ENDIF
-       ENDIF
-
-       CALL MPI_SENDRECV( trlp(1)%age, trlp_count, mpi_particle_type,     &
-                          pleft, 1, particles(number_of_particles+1)%age, &
-                          trrp_count_recv, mpi_particle_type, pright, 1,  &
+       ALLOCATE(rvrp(MAX(1,trrp_count_recv)))
+
+       CALL MPI_SENDRECV( trlp(1)%radius, max(1,trlp_count), mpi_particle_type,&
+                          pleft, 1, rvrp(1)%radius,                            &
+                          max(1,trrp_count_recv), mpi_particle_type, pright, 1,&
                           comm2d, status, ierr )
+
+       IF ( trrp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvrp(1:trrp_count_recv))
+
+       DEALLOCATE(rvrp)
 
        IF ( use_particle_tails )  THEN
@@ -405,7 +466,4 @@
        ENDIF
 
-       number_of_particles = number_of_particles + trrp_count_recv
-       number_of_tails     = number_of_tails     + trrpt_count_recv
-
 !
 !--    Send right boundary, receive left boundary
@@ -414,22 +472,14 @@
                           comm2d, status, ierr )
 
-       IF ( number_of_particles + trlp_count_recv > &
-            maximum_number_of_particles )           &
-       THEN
-          IF ( netcdf_data_format < 3 )  THEN
-             message_string = 'maximum_number_of_particles ' //  &
-                              'needs to be increased ' //        &
-                              '&but this is not allowed with '// &
-                              'netcdf_data_format < 3'
-             CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 )
-          ELSE
-             CALL lpm_extend_particle_array( trlp_count_recv )
-          ENDIF
-       ENDIF
-
-       CALL MPI_SENDRECV( trrp(1)%age, trrp_count, mpi_particle_type,      &
-                          pright, 1, particles(number_of_particles+1)%age, &
-                          trlp_count_recv, mpi_particle_type, pleft, 1,    &
+       ALLOCATE(rvlp(MAX(1,trlp_count_recv)))
+
+       CALL MPI_SENDRECV( trrp(1)%radius, max(1,trrp_count), mpi_particle_type,&
+                          pright, 1, rvlp(1)%radius,                           &
+                          max(1,trlp_count_recv), mpi_particle_type, pleft, 1, &
                           comm2d, status, ierr )
+
+       IF ( trlp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvlp(1:trlp_count_recv))
+
+       DEALLOCATE(rvlp)
 
        IF ( use_particle_tails )  THEN
@@ -469,13 +519,11 @@
        ENDIF
 
-       number_of_particles = number_of_particles + trlp_count_recv
-       number_of_tails     = number_of_tails     + trlpt_count_recv
+!       number_of_particles = number_of_particles + trlp_count_recv
+!       number_of_tails     = number_of_tails     + trlpt_count_recv
 
        IF ( use_particle_tails )  THEN
           DEALLOCATE( trlpt, trrpt )
        ENDIF
-       DEALLOCATE( trlp, trrp ) 
-
-       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'pause' )
+       DEALLOCATE( trlp, trrp )
 
     ENDIF
@@ -490,7 +538,7 @@
 !-- For a one-dimensional decomposition along x, no transfer is necessary,
 !-- because the particle remains on the PE.
-    trsp_count  = 0
+    trsp_count  = nr_move_south
     trspt_count = 0
-    trnp_count  = 0
+    trnp_count  = nr_move_north
     trnpt_count = 0
 
@@ -504,19 +552,28 @@
 !--    First calculate the storage necessary for sending and receiving the
 !--    data
-       DO  n = 1, number_of_particles
-          IF ( particle_mask(n) )  THEN
-             j = ( particles(n)%y + 0.5 * dy ) * ddy
-!
-!--          Above calculation does not work for indices less than zero
-             IF ( particles(n)%y < -0.5 * dy )  j = -1
-
-             IF ( j < nys )  THEN
-                trsp_count = trsp_count + 1
-                IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count+1
-             ELSEIF ( j > nyn )  THEN
-                trnp_count = trnp_count + 1
-                IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count+1
-             ENDIF
-          ENDIF
+       DO  ip = nxl, nxr
+          DO  jp = nys, nyn, nyn-nys                                 !compute only first (nys) and last (nyn) loop iterration
+             DO  kp = nzb+1, nzt
+                number_of_particles = prt_count(kp,jp,ip)
+                IF ( number_of_particles <= 0 )  CYCLE
+                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                DO  n = 1, number_of_particles
+                   IF ( particles(n)%particle_mask )  THEN
+                      j = ( particles(n)%y + 0.5_wp * dy ) * ddy
+!
+!--                   Above calculation does not work for indices less than zero
+                      IF ( particles(n)%y < -0.5_wp * dy )  j = -1
+
+                      IF ( j < nys )  THEN
+                         trsp_count = trsp_count + 1
+                         IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count + 1
+                      ELSEIF ( j > nyn )  THEN
+                         trnp_count = trnp_count + 1
+                         IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count + 1
+                      ENDIF
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
        ENDDO
 
@@ -528,10 +585,6 @@
        ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
 
-       trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0, 0, 0, 0 )
-       trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
-                             0.0, 0, 0, 0, 0 )
+       trsp = zero_particle
+       trnp = zero_particle
 
        IF ( use_particle_tails )  THEN
@@ -541,174 +594,185 @@
        ENDIF
 
-       trsp_count  = 0
+       trsp_count  = nr_move_south
        trspt_count = 0
-       trnp_count  = 0
+       trnp_count  = nr_move_north
        trnpt_count = 0
 
+       trsp(1:nr_move_south) = move_also_south(1:nr_move_south)
+       trnp(1:nr_move_north) = move_also_north(1:nr_move_north)
+
     ENDIF
 
-    DO  n = 1, number_of_particles
-
-       nn = particles(n)%tail_id
-!
-!--    Only those particles that have not been marked as 'deleted' may be 
-!--    moved.
-       IF ( particle_mask(n) )  THEN
-          j = ( particles(n)%y + 0.5 * dy ) * ddy
-!
-!--       Above calculation does not work for indices less than zero
-          IF ( particles(n)%y < -0.5 * dy )  j = -1
-
-          IF ( j < nys )  THEN
-             IF ( j < 0 )  THEN
-!
-!--             Apply boundary condition along y
-                IF ( ibc_par_ns == 0 )  THEN
-!
-!--                Cyclic condition
-                   IF ( pdims(2) == 1 )  THEN
-                      particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
-                      particles(n)%origin_y = ( ny + 1 ) * dy + &
-                                              particles(n)%origin_y
-                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                         i = particles(n)%tailpoints
-                         particle_tail_coordinates(1:i,2,nn) = ( ny+1 ) * dy&
-                                        + particle_tail_coordinates(1:i,2,nn)
+    DO  ip = nxl, nxr
+       DO  jp = nys, nyn, nyn-nys ! compute only first (nys) and last (nyn) loop iterration
+          DO  kp = nzb+1, nzt
+             number_of_particles = prt_count(kp,jp,ip)
+             IF ( number_of_particles <= 0 )  CYCLE
+             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+             DO  n = 1, number_of_particles
+
+                nn = particles(n)%tail_id
+!
+!--             Only those particles that have not been marked as 'deleted' may
+!--             be moved.
+                IF ( particles(n)%particle_mask )  THEN
+                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
+!
+!--                Above calculation does not work for indices less than zero
+                   IF ( particles(n)%y < -0.5_wp * dy )  j = -1
+
+                   IF ( j < nys )  THEN
+                      IF ( j < 0 )  THEN
+!
+!--                      Apply boundary condition along y
+                         IF ( ibc_par_ns == 0 )  THEN
+!
+!--                         Cyclic condition
+                            IF ( pdims(2) == 1 )  THEN
+                               particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
+                               particles(n)%origin_y = ( ny + 1 ) * dy + &
+                               particles(n)%origin_y
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  i = particles(n)%tailpoints
+                                  particle_tail_coordinates(1:i,2,nn) =        &
+                                     ( ny+1 ) * dy + particle_tail_coordinates(1:i,2,nn)
+                               ENDIF
+                            ELSE
+                               trsp_count = trsp_count + 1
+                               trsp(trsp_count) = particles(n)
+                               trsp(trsp_count)%y = ( ny + 1 ) * dy + &
+                               trsp(trsp_count)%y
+                               trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
+                               + ( ny + 1 ) * dy
+                               particles(n)%particle_mask = .FALSE.
+                               deleted_particles = deleted_particles + 1
+
+                               IF ( trsp(trsp_count)%y >= (ny+0.5_wp)* dy - 1.0E-12_wp )  THEN
+                                  trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10_wp
+                                  !++ why is 1 subtracted in next statement???
+                                  trsp(trsp_count)%origin_y =                        &
+                                  trsp(trsp_count)%origin_y - 1
+                               ENDIF
+
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  trspt_count = trspt_count + 1
+                                  trspt(:,:,trspt_count) = &
+                                  particle_tail_coordinates(:,:,nn)
+                                  trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
+                                  trspt(:,2,trspt_count)
+                                  tail_mask(nn) = .FALSE.
+                                  deleted_tails = deleted_tails + 1
+                               ENDIF
+                            ENDIF
+
+                         ELSEIF ( ibc_par_ns == 1 )  THEN
+!
+!--                         Particle absorption
+                            particles(n)%particle_mask = .FALSE.
+                            deleted_particles = deleted_particles + 1
+                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                               tail_mask(nn) = .FALSE.
+                               deleted_tails = deleted_tails + 1
+                            ENDIF
+
+                         ELSEIF ( ibc_par_ns == 2 )  THEN
+!
+!--                         Particle reflection
+                            particles(n)%y       = -particles(n)%y
+                            particles(n)%speed_y = -particles(n)%speed_y
+
+                         ENDIF
+                      ELSE
+!
+!--                      Store particle data in the transfer array, which will 
+!--                      be send to the neighbouring PE
+                         trsp_count = trsp_count + 1
+                         trsp(trsp_count) = particles(n)
+                         particles(n)%particle_mask = .FALSE.
+                         deleted_particles = deleted_particles + 1
+
+                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                            trspt_count = trspt_count + 1
+                            trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
+                            tail_mask(nn) = .FALSE.
+                            deleted_tails = deleted_tails + 1
+                         ENDIF
                       ENDIF
-                   ELSE
-                      trsp_count = trsp_count + 1
-                      trsp(trsp_count) = particles(n)
-                      trsp(trsp_count)%y = ( ny + 1 ) * dy + &
-                                           trsp(trsp_count)%y
-                      trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
-                                                  + ( ny + 1 ) * dy
-                      particle_mask(n) = .FALSE.
-                      deleted_particles = deleted_particles + 1
-
-                      IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.0E-12 ) THEN
-                         trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10
-!++ why is 1 subtracted in next statement???
-                         trsp(trsp_count)%origin_y =                        &
-                                                   trsp(trsp_count)%origin_y - 1
+
+                   ELSEIF ( j > nyn )  THEN
+                      IF ( j > ny )  THEN
+!
+!--                       Apply boundary condition along x
+                         IF ( ibc_par_ns == 0 )  THEN
+!
+!--                         Cyclic condition
+                            IF ( pdims(2) == 1 )  THEN
+                               particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
+                               particles(n)%origin_y =                         &
+                                  particles(n)%origin_y - ( ny + 1 ) * dy
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  i = particles(n)%tailpoints
+                                  particle_tail_coordinates(1:i,2,nn) =        &
+                                     - (ny+1) * dy + particle_tail_coordinates(1:i,2,nn)
+                               ENDIF
+                            ELSE
+                               trnp_count = trnp_count + 1
+                               trnp(trnp_count) = particles(n)
+                               trnp(trnp_count)%y =                            &
+                                  trnp(trnp_count)%y - ( ny + 1 ) * dy
+                               trnp(trnp_count)%origin_y =                     &
+                                  trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
+                               particles(n)%particle_mask = .FALSE.
+                               deleted_particles = deleted_particles + 1
+                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                                  trnpt_count = trnpt_count + 1
+                                  trnpt(:,:,trnpt_count) =                     &
+                                     particle_tail_coordinates(:,:,nn)
+                                  trnpt(:,2,trnpt_count) =                     &
+                                     trnpt(:,2,trnpt_count) - ( ny + 1 ) * dy
+                                  tail_mask(nn) = .FALSE.
+                                  deleted_tails = deleted_tails + 1
+                               ENDIF
+                            ENDIF
+
+                         ELSEIF ( ibc_par_ns == 1 )  THEN
+!
+!--                         Particle absorption
+                            particles(n)%particle_mask = .FALSE.
+                            deleted_particles = deleted_particles + 1
+                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                               tail_mask(nn) = .FALSE.
+                               deleted_tails = deleted_tails + 1
+                            ENDIF
+
+                         ELSEIF ( ibc_par_ns == 2 )  THEN
+!
+!--                         Particle reflection
+                            particles(n)%y       = 2 * ( ny * dy ) - particles(n)%y
+                            particles(n)%speed_y = -particles(n)%speed_y
+
+                         ENDIF
+                      ELSE
+!
+!--                      Store particle data in the transfer array, which will 
+!--                      be send to the neighbouring PE
+                         trnp_count = trnp_count + 1
+                         trnp(trnp_count) = particles(n)
+                         particles(n)%particle_mask = .FALSE.
+                         deleted_particles = deleted_particles + 1
+
+                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
+                            trnpt_count = trnpt_count + 1
+                            trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
+                            tail_mask(nn) = .FALSE.
+                            deleted_tails = deleted_tails + 1
+                         ENDIF
                       ENDIF
 
-                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                         trspt_count = trspt_count + 1
-                         trspt(:,:,trspt_count) = &
-                                               particle_tail_coordinates(:,:,nn)
-                         trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
-                                                  trspt(:,2,trspt_count)
-                         tail_mask(nn) = .FALSE.
-                         deleted_tails = deleted_tails + 1
-                      ENDIF
                    ENDIF
-
-                ELSEIF ( ibc_par_ns == 1 )  THEN
-!
-!--                Particle absorption
-                   particle_mask(n) = .FALSE.
-                   deleted_particles = deleted_particles + 1
-                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                      tail_mask(nn) = .FALSE.
-                      deleted_tails = deleted_tails + 1
-                   ENDIF
-
-                ELSEIF ( ibc_par_ns == 2 )  THEN
-!
-!--                Particle reflection
-                   particles(n)%y       = -particles(n)%y
-                   particles(n)%speed_y = -particles(n)%speed_y
-
                 ENDIF
-             ELSE
-!
-!--             Store particle data in the transfer array, which will be send
-!--             to the neighbouring PE
-                trsp_count = trsp_count + 1
-                trsp(trsp_count) = particles(n)
-                particle_mask(n) = .FALSE.
-                deleted_particles = deleted_particles + 1
-
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   trspt_count = trspt_count + 1
-                   trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
-                   tail_mask(nn) = .FALSE.
-                   deleted_tails = deleted_tails + 1
-                ENDIF
-             ENDIF
-
-          ELSEIF ( j > nyn )  THEN
-             IF ( j > ny )  THEN
-!
-!--             Apply boundary condition along x
-                IF ( ibc_par_ns == 0 )  THEN
-!
-!--                Cyclic condition
-                   IF ( pdims(2) == 1 )  THEN
-                      particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
-                      particles(n)%origin_y = particles(n)%origin_y - &
-                                              ( ny + 1 ) * dy
-                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                         i = particles(n)%tailpoints
-                         particle_tail_coordinates(1:i,2,nn) = - (ny+1) * dy &
-                                           + particle_tail_coordinates(1:i,2,nn)
-                      ENDIF
-                   ELSE
-                      trnp_count = trnp_count + 1
-                      trnp(trnp_count) = particles(n)
-                      trnp(trnp_count)%y = trnp(trnp_count)%y - &
-                                           ( ny + 1 ) * dy
-                      trnp(trnp_count)%origin_y = trnp(trnp_count)%origin_y &
-                                                  - ( ny + 1 ) * dy
-                      particle_mask(n) = .FALSE.
-                      deleted_particles = deleted_particles + 1
-
-                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                         trnpt_count = trnpt_count + 1
-                         trnpt(:,:,trnpt_count) = &
-                                               particle_tail_coordinates(:,:,nn)
-                         trnpt(:,2,trnpt_count) = trnpt(:,2,trnpt_count) - &
-                                                  ( ny + 1 ) * dy
-                         tail_mask(nn) = .FALSE.
-                         deleted_tails = deleted_tails + 1
-                      ENDIF
-                   ENDIF
-
-                ELSEIF ( ibc_par_ns == 1 )  THEN
-!
-!--                Particle absorption
-                   particle_mask(n) = .FALSE.
-                   deleted_particles = deleted_particles + 1
-                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                      tail_mask(nn) = .FALSE.
-                      deleted_tails = deleted_tails + 1
-                   ENDIF
-
-                ELSEIF ( ibc_par_ns == 2 )  THEN
-!
-!--                Particle reflection
-                   particles(n)%y       = 2 * ( ny * dy ) - particles(n)%y
-                   particles(n)%speed_y = -particles(n)%speed_y
-
-                ENDIF
-             ELSE
-!
-!--             Store particle data in the transfer array, which will be send
-!--             to the neighbouring PE
-                trnp_count = trnp_count + 1
-                trnp(trnp_count) = particles(n)
-                particle_mask(n) = .FALSE.
-                deleted_particles = deleted_particles + 1
-
-                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
-                   trnpt_count = trnpt_count + 1
-                   trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
-                   tail_mask(nn) = .FALSE.
-                   deleted_tails = deleted_tails + 1
-                ENDIF
-             ENDIF
-
-          ENDIF
-       ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
     ENDDO
 
@@ -718,27 +782,18 @@
     IF ( pdims(2) /= 1 )  THEN
 
-       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'continue' )
        CALL MPI_SENDRECV( trsp_count,      1, MPI_INTEGER, psouth, 0, &
                           trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
                           comm2d, status, ierr )
 
-       IF ( number_of_particles + trnp_count_recv > &
-            maximum_number_of_particles )           &
-       THEN
-          IF ( netcdf_data_format < 3 )  THEN
-             message_string = 'maximum_number_of_particles ' //  &
-                              'needs to be increased ' //        &
-                              '&but this is not allowed with '// &
-                              'netcdf_data_format < 3'
-             CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 ) 
-          ELSE
-             CALL lpm_extend_particle_array( trnp_count_recv )
-          ENDIF
-       ENDIF
-
-       CALL MPI_SENDRECV( trsp(1)%age, trsp_count, mpi_particle_type,      &
-                          psouth, 1, particles(number_of_particles+1)%age, &
+       ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
+
+       CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type,      &
+                          psouth, 1, rvnp(1)%radius,                             &
                           trnp_count_recv, mpi_particle_type, pnorth, 1,   &
                           comm2d, status, ierr )
+
+       IF ( trnp_count_recv  > 0 )  CALL Add_particles_to_gridcell(rvnp(1:trnp_count_recv))
+
+       DEALLOCATE(rvnp)
 
        IF ( use_particle_tails )  THEN
@@ -779,6 +834,6 @@
        ENDIF
 
-       number_of_particles = number_of_particles + trnp_count_recv
-       number_of_tails     = number_of_tails     + trnpt_count_recv
+!       number_of_particles = number_of_particles + trnp_count_recv
+!       number_of_tails     = number_of_tails     + trnpt_count_recv
 
 !
@@ -788,22 +843,14 @@
                           comm2d, status, ierr )
 
-       IF ( number_of_particles + trsp_count_recv > &
-            maximum_number_of_particles )           &
-       THEN
-          IF ( netcdf_data_format < 3 )  THEN
-             message_string = 'maximum_number_of_particles ' //   &
-                              'needs to be increased ' //         &
-                              '&but this is not allowed with ' // &
-                              'netcdf_data_format < 3'
-            CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 )  
-          ELSE
-             CALL lpm_extend_particle_array( trsp_count_recv )
-          ENDIF
-       ENDIF
-
-       CALL MPI_SENDRECV( trnp(1)%age, trnp_count, mpi_particle_type,      &
-                          pnorth, 1, particles(number_of_particles+1)%age, &
+       ALLOCATE(rvsp(MAX(1,trsp_count_recv)))
+
+       CALL MPI_SENDRECV( trnp(1)%radius, trnp_count, mpi_particle_type,      &
+                          pnorth, 1, rvsp(1)%radius,                          &
                           trsp_count_recv, mpi_particle_type, psouth, 1,   &
                           comm2d, status, ierr )
+
+       IF ( trsp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvsp(1:trsp_count_recv))
+
+       DEALLOCATE(rvsp)
 
        IF ( use_particle_tails )  THEN
@@ -851,6 +898,4 @@
        DEALLOCATE( trsp, trnp )
 
-       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'stop' )
-
     ENDIF
 
@@ -863,5 +908,5 @@
        nn = particles(n)%tail_id
 
-       IF ( particles(n)%x < -0.5 * dx )  THEN
+       IF ( particles(n)%x < -0.5_wp * dx )  THEN
 
           IF ( ibc_par_lr == 0 )  THEN
@@ -877,5 +922,5 @@
 !
 !--          Particle absorption
-             particle_mask(n) = .FALSE.
+             particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
@@ -890,5 +935,5 @@
           ENDIF
 
-       ELSEIF ( particles(n)%x >= ( nx + 0.5 ) * dx )  THEN
+       ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
 
           IF ( ibc_par_lr == 0 )  THEN
@@ -904,5 +949,5 @@
 !
 !--          Particle absorption
-             particle_mask(n) = .FALSE.
+             particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
@@ -919,5 +964,5 @@
        ENDIF
 
-       IF ( particles(n)%y < -0.5 * dy )  THEN
+       IF ( particles(n)%y < -0.5_wp * dy )  THEN
 
           IF ( ibc_par_ns == 0 )  THEN
@@ -933,5 +978,5 @@
 !
 !--          Particle absorption
-             particle_mask(n) = .FALSE.
+             particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
@@ -946,5 +991,5 @@
           ENDIF
 
-       ELSEIF ( particles(n)%y >= ( ny + 0.5 ) * dy )  THEN
+       ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
 
           IF ( ibc_par_ns == 0 )  THEN
@@ -960,5 +1005,5 @@
 !
 !--          Particle absorption
-             particle_mask(n) = .FALSE.
+             particles(n)%particle_mask = .FALSE.
              deleted_particles = deleted_particles + 1
              IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
@@ -991,4 +1036,270 @@
 #endif
 
+    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'stop' )
 
  END SUBROUTINE lpm_exchange_horiz
+
+ SUBROUTINE Add_particles_to_gridcell (particle_array)
+
+!
+!-- If a particle moves from one processor to another, this subroutine moves 
+!-- the corresponding elements from the particle arrays of the old grid cells 
+!-- to the particle arrays of the new grid cells.
+
+    IMPLICIT NONE
+
+    INTEGER(iwp)        ::  ip        !:
+    INTEGER(iwp)        ::  jp        !:
+    INTEGER(iwp)        ::  kp        !:
+    INTEGER(iwp)        ::  n         !:
+    INTEGER(iwp)        ::  pindex    !:
+
+    LOGICAL             ::  pack_done !:
+
+    TYPE(particle_type), DIMENSION(:), INTENT(IN)       :: particle_array
+
+    pack_done     = .FALSE.
+
+    nr_move_north = 0
+    nr_move_south = 0
+
+    DO n = 1, SIZE(particle_array)
+       ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx
+       jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy
+       kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt_m1
+
+       IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
+            .AND.  kp >= nzb+1  .AND.  kp <= nzt)  THEN ! particle stays on processor
+          number_of_particles = prt_count(kp,jp,ip)
+          particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+
+          pindex = prt_count(kp,jp,ip)+1
+          IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
+             IF ( pack_done )  THEN
+                CALL realloc_particles_array (ip,jp,kp)
+             ELSE
+                CALL lpm_pack_arrays
+                prt_count(kp,jp,ip) = number_of_particles
+                pindex = prt_count(kp,jp,ip)+1
+                IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
+                   CALL realloc_particles_array (ip,jp,kp)
+                ENDIF
+                pack_done = .TRUE.
+             ENDIF
+          ENDIF
+          grid_particles(kp,jp,ip)%particles(pindex) = particle_array(n)
+          prt_count(kp,jp,ip) = pindex
+       ELSE
+          IF ( jp == nys - 1 )  THEN
+             nr_move_south = nr_move_south+1
+             move_also_south(nr_move_south) = particle_array(n)
+             IF ( jp == -1 )  THEN
+                move_also_south(nr_move_south)%y =                             &
+                   move_also_south(nr_move_south)%y + ( ny + 1 ) * dy
+                move_also_south(nr_move_south)%origin_y =                      &
+                   move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
+             ENDIF
+          ELSEIF ( jp == nyn+1 )  THEN
+             nr_move_north = nr_move_north+1
+             move_also_north(nr_move_north) = particle_array(n)
+             IF ( jp == ny+1 )  THEN
+                move_also_north(nr_move_north)%y =                             &
+                   move_also_north(nr_move_north)%y - ( ny + 1 ) * dy
+                move_also_north(nr_move_north)%origin_y =                      &
+                   move_also_north(nr_move_north)%origin_y - ( ny + 1 ) * dy
+             ENDIF
+          ELSE
+             WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn
+          ENDIF
+       ENDIF
+    ENDDO
+
+    RETURN
+
+ END SUBROUTINE Add_particles_to_gridcell
+
+
+
+
+ SUBROUTINE lpm_move_particle
+
+!
+!-- If a particle moves from one grid cell to another (on the current 
+!-- processor!), this subroutine moves the corresponding element from the
+!-- particle array of the old grid cell to the particle array of the new grid 
+!-- cell.
+
+    IMPLICIT NONE
+
+    INTEGER(iwp)        ::  i           !:
+    INTEGER(iwp)        ::  ip          !:
+    INTEGER(iwp)        ::  j           !:
+    INTEGER(iwp)        ::  jp          !:
+    INTEGER(iwp)        ::  k           !:
+    INTEGER(iwp)        ::  kp          !:
+    INTEGER(iwp)        ::  n           !:
+    INTEGER(iwp)        ::  np_old_cell !:
+    INTEGER(iwp)        ::  n_start     !:
+    INTEGER(iwp)        ::  pindex      !:
+
+    LOGICAL             ::  pack_done   !:
+
+    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !:
+
+    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
+
+    DO  ip = nxl, nxr
+       DO  jp = nys, nyn
+          DO  kp = nzb+1, nzt
+
+             np_old_cell = prt_count(kp,jp,ip)
+             IF ( np_old_cell <= 0 )  CYCLE
+             particles_old_cell => grid_particles(kp,jp,ip)%particles(1:np_old_cell)
+             n_start = -1
+             
+             DO  n = 1, np_old_cell
+                i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
+                j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
+                k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
+!
+!--             Check, if particle has moved to another grid cell.
+                IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
+!
+!--                The particle has moved to another grid cell. Now check, if 
+!--                particle stays on the same processor.
+                   IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.      &
+                        j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
+!
+!--                   If the particle stays on the same processor, the particle
+!--                   will be added to the particle array of the new processor.
+                      number_of_particles = prt_count(k,j,i)
+                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+
+                      pindex = prt_count(k,j,i)+1
+                      IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )  &
+                      THEN
+                         n_start = n
+                         EXIT
+                      ENDIF
+
+                      grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
+                      prt_count(k,j,i) = pindex
+
+                      particles_old_cell(n)%particle_mask = .FALSE.
+                   ENDIF
+                ENDIF
+             ENDDO
+
+             IF ( n_start .GE. 0 )  THEN
+                pack_done = .FALSE.
+                DO  n = n_start, np_old_cell
+                   i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
+                   j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
+                   k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
+                   IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
+!
+!--                   Particle is in different box
+                      IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.   &
+                           j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
+!
+!--                      Particle stays on processor
+                         number_of_particles = prt_count(k,j,i)
+                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+
+                         pindex = prt_count(k,j,i)+1
+                         IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) &
+                         THEN
+                            IF ( pack_done )  THEN
+                               CALL realloc_particles_array(i,j,k)
+                               pindex = prt_count(k,j,i)+1
+                            ELSE
+                               CALL lpm_pack_arrays
+                               prt_count(k,j,i) = number_of_particles
+
+                               pindex = prt_count(k,j,i)+1
+!
+!--                            If number of particles in the new grid box 
+!--                            exceeds its allocated memory, the particle array
+!--                            will be reallocated
+                               IF ( pindex > SIZE(grid_particles(k,j,i)%particles) )  THEN
+                                  CALL realloc_particles_array(i,j,k)
+                                  pindex = prt_count(k,j,i)+1
+                               ENDIF
+
+                               pack_done = .TRUE.
+                            ENDIF
+                         ENDIF
+
+                         grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
+                         prt_count(k,j,i) = pindex
+
+                         particles_old_cell(n)%particle_mask = .FALSE.
+                      ENDIF
+                   ENDIF
+                ENDDO
+             ENDIF
+          ENDDO
+       ENDDO
+    ENDDO
+
+    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'stop' )
+
+    RETURN
+
+ END SUBROUTINE lpm_move_particle
+
+ SUBROUTINE realloc_particles_array (i,j,k,size_in)
+
+    IMPLICIT NONE
+
+    INTEGER(iwp), INTENT(in)                       ::  i              !:
+    INTEGER(iwp), INTENT(in)                       ::  j              !:
+    INTEGER(iwp), INTENT(in)                       ::  k              !:
+    INTEGER(iwp), INTENT(in), optional             ::  size_in        !:
+
+    INTEGER(iwp)                                   :: old_size        !:
+    INTEGER(iwp)                                   :: new_size        !:
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !:
+    TYPE(particle_type), DIMENSION(500)            :: tmp_particles_s !:
+
+
+    old_size = SIZE(grid_particles(k,j,i)%particles)
+
+    IF ( PRESENT(size_in) )   THEN
+       new_size = size_in
+    ELSE
+       new_size = old_size * ( 1.0 + alloc_factor / 100.0 )
+    ENDIF
+
+    new_size = MAX( new_size, min_nr_particle )
+
+    IF ( old_size <= 500 )  THEN
+
+       tmp_particles_s(1:old_size) = grid_particles(k,j,i)%particles(1:old_size)
+
+       DEALLOCATE(grid_particles(k,j,i)%particles)
+       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
+
+       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_s(1:old_size)
+       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
+
+    ELSE
+
+       ALLOCATE(tmp_particles_d(new_size))
+       tmp_particles_d(1:old_size) = grid_particles(k,j,i)%particles
+
+       DEALLOCATE(grid_particles(k,j,i)%particles)
+       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
+
+       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_d(1:old_size)
+       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
+
+       DEALLOCATE(tmp_particles_d)
+
+    ENDIF
+    particles => grid_particles(k,j,i)%particles(1:number_of_particles)
+
+    RETURN
+ END SUBROUTINE realloc_particles_array
+
+END MODULE lpm_exchange_horiz_mod
