Index: palm/trunk/SOURCE/lpm_exchange_horiz.f90
===================================================================
--- palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 1874)
+++ palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 1929)
@@ -19,5 +19,13 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfixes: 
+! - reallocation of new particles 
+!   ( did not work for small number of min_nr_particle )
+! - dynamical reallocation of north-south exchange arrays ( particles got lost )
+! - north-south exchange ( nr_move_north, nr_move_south were overwritten by zero )
+! - horizontal particle boundary conditions in serial mode
+!
+! Remove unused variables
+! Descriptions in variable declaration blocks added
 ! 
 ! Former revisions:
@@ -121,6 +129,6 @@
     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
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_north
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_south
 
     SAVE
@@ -161,35 +169,27 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i                                       !<
-    INTEGER(iwp) ::  ip                                      !<
-    INTEGER(iwp) ::  j                                       !<
-    INTEGER(iwp) ::  jp                                      !<
-    INTEGER(iwp) ::  kp                                      !<
-    INTEGER(iwp) ::  n                                       !<
-    INTEGER(iwp) ::  trlp_count                              !<
-    INTEGER(iwp) ::  trlp_count_recv                         !<
-    INTEGER(iwp) ::  trlpt_count                             !<
-    INTEGER(iwp) ::  trlpt_count_recv                        !<
-    INTEGER(iwp) ::  trnp_count                              !<
-    INTEGER(iwp) ::  trnp_count_recv                         !<
-    INTEGER(iwp) ::  trnpt_count                             !<
-    INTEGER(iwp) ::  trnpt_count_recv                        !<
-    INTEGER(iwp) ::  trrp_count                              !<
-    INTEGER(iwp) ::  trrp_count_recv                         !<
-    INTEGER(iwp) ::  trrpt_count                             !<
-    INTEGER(iwp) ::  trrpt_count_recv                        !<
-    INTEGER(iwp) ::  trsp_count                              !<
-    INTEGER(iwp) ::  trsp_count_recv                         !<
-    INTEGER(iwp) ::  trspt_count                             !<
-    INTEGER(iwp) ::  trspt_count_recv                        !<
-
-    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  !<
+    INTEGER(iwp) ::  i                 !< grid index (x) of particle positition
+    INTEGER(iwp) ::  ip                !< index variable along x
+    INTEGER(iwp) ::  j                 !< grid index (y) of particle positition
+    INTEGER(iwp) ::  jp                !< index variable along y
+    INTEGER(iwp) ::  kp                !< index variable along z
+    INTEGER(iwp) ::  n                 !< particle index variable 
+    INTEGER(iwp) ::  trlp_count        !< number of particles send to left PE
+    INTEGER(iwp) ::  trlp_count_recv   !< number of particles receive from right PE
+    INTEGER(iwp) ::  trnp_count        !< number of particles send to north PE
+    INTEGER(iwp) ::  trnp_count_recv   !< number of particles receive from south PE
+    INTEGER(iwp) ::  trrp_count        !< number of particles send to right PE
+    INTEGER(iwp) ::  trrp_count_recv   !< number of particles receive from left PE
+    INTEGER(iwp) ::  trsp_count        !< number of particles send to south PE
+    INTEGER(iwp) ::  trsp_count_recv   !< number of particles receive from north PE
+
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !< particles received from right PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !< particles received from south PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !< particles received from left PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !< particles received from north PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !< particles send to left PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !< particles send to north PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !< particles send to right PE
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !< particles send to south PE
 
     CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
@@ -209,12 +209,8 @@
 !-- be adjusted.
     trlp_count  = 0
-    trlpt_count = 0
     trrp_count  = 0
-    trrpt_count = 0
 
     trlp_count_recv   = 0
-    trlpt_count_recv  = 0
     trrp_count_recv   = 0
-    trrpt_count_recv  = 0
 
     IF ( pdims(1) /= 1 )  THEN
@@ -232,6 +228,6 @@
                    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
+!
+!--                   Above calculation does not work for indices less than zero
                       IF ( particles(n)%x < -0.5_wp * dx )  i = -1
 
@@ -249,7 +245,5 @@
 
        IF ( trlp_count  == 0 )  trlp_count  = 1
-       IF ( trlpt_count == 0 )  trlpt_count = 1
        IF ( trrp_count  == 0 )  trrp_count  = 1
-       IF ( trrpt_count == 0 )  trrpt_count = 1
 
        ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
@@ -259,12 +253,10 @@
 
        trlp_count  = 0
-       trlpt_count = 0
        trrp_count  = 0
-       trrpt_count = 0
 
     ENDIF
 !
 !-- Compute only first (nxl) and last (nxr) loop iterration
-    DO  ip = nxl, nxr,nxr-nxl
+    DO  ip = nxl, nxr, nxr-nxl
        DO  jp = nys, nyn
           DO  kp = nzb+1, nzt
@@ -279,15 +271,15 @@
 
                    i = ( particles(n)%x + 0.5_wp * dx ) * ddx
-   !
-   !--             Above calculation does not work for indices less than zero
+!
+!--                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
+!
+!--                   Apply boundary condition along x
                          IF ( ibc_par_lr == 0 )  THEN
-   !
-   !--                      Cyclic condition
+!
+!--                         Cyclic condition
                             IF ( pdims(1) == 1 )  THEN
                                particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
@@ -312,12 +304,12 @@
 
                          ELSEIF ( ibc_par_lr == 1 )  THEN
-   !
-   !--                      Particle absorption
+!
+!--                         Particle absorption
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
 
                          ELSEIF ( ibc_par_lr == 2 )  THEN
-   !
-   !--                      Particle reflection
+!
+!--                         Particle reflection
                             particles(n)%x       = -particles(n)%x
                             particles(n)%speed_x = -particles(n)%speed_x
@@ -325,7 +317,7 @@
                          ENDIF
                       ELSE
-   !
-   !--                   Store particle data in the transfer array, which will be 
-   !--                   send to the neighbouring PE
+!
+!--                      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)
@@ -337,9 +329,9 @@
                    ELSEIF ( i > nxr )  THEN
                       IF ( i > nx )  THEN
-   !
-   !--                   Apply boundary condition along x
+!
+!--                      Apply boundary condition along x
                          IF ( ibc_par_lr == 0 )  THEN
-   !
-   !--                      Cyclic condition
+!
+!--                         Cyclic condition
                             IF ( pdims(1) == 1 )  THEN
                                particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
@@ -358,12 +350,12 @@
 
                          ELSEIF ( ibc_par_lr == 1 )  THEN
-   !
-   !--                      Particle absorption
+!
+!--                         Particle absorption
                             particles(n)%particle_mask = .FALSE.
                             deleted_particles = deleted_particles + 1
 
                          ELSEIF ( ibc_par_lr == 2 )  THEN
-   !
-   !--                      Particle reflection
+!
+!--                         Particle reflection
                             particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
                             particles(n)%speed_x = -particles(n)%speed_x
@@ -371,7 +363,7 @@
                          ENDIF
                       ELSE
-   !
-   !--                   Store particle data in the transfer array, which will be send
-   !--                   to the neighbouring PE
+!
+!--                      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)
@@ -383,4 +375,5 @@
                    ENDIF
                 ENDIF
+
              ENDDO
           ENDDO
@@ -389,4 +382,12 @@
 
 !
+!-- Allocate arrays required for north-south exchange, as these
+!-- are used directly after particles are exchange along x-direction.
+    ALLOCATE( move_also_north(1:NR_2_direction_move) )
+    ALLOCATE( move_also_south(1:NR_2_direction_move) )
+
+    nr_move_north = 0
+    nr_move_south = 0
+!
 !-- Send left boundary, receive right boundary (but first exchange how many
 !-- and check, if particle storage must be extended)
@@ -427,5 +428,4 @@
 
     ENDIF
-
 
 !
@@ -435,15 +435,11 @@
 !-- Find out first the number of particles to be transferred and allocate
 !-- temporary arrays needed to store them.
-!-- For a one-dimensional decomposition along x, no transfer is necessary,
+!-- For a one-dimensional decomposition along y, no transfer is necessary,
 !-- because the particle remains on the PE.
     trsp_count  = nr_move_south
-    trspt_count = 0
     trnp_count  = nr_move_north
-    trnpt_count = 0
 
     trsp_count_recv   = 0
-    trspt_count_recv  = 0
     trnp_count_recv   = 0
-    trnpt_count_recv  = 0
 
     IF ( pdims(2) /= 1 )  THEN
@@ -452,5 +448,5 @@
 !--    data
        DO  ip = nxl, nxr
-          DO  jp = nys, nyn, nyn-nys                                 !compute only first (nys) and last (nyn) loop iterration
+          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)
@@ -476,7 +472,5 @@
 
        IF ( trsp_count  == 0 )  trsp_count  = 1
-       IF ( trspt_count == 0 )  trspt_count = 1
        IF ( trnp_count  == 0 )  trnp_count  = 1
-       IF ( trnpt_count == 0 )  trnpt_count = 1
 
        ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
@@ -486,8 +480,6 @@
 
        trsp_count  = nr_move_south
-       trspt_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)
@@ -506,4 +498,5 @@
 !--             be moved.
                 IF ( particles(n)%particle_mask )  THEN
+
                    j = ( particles(n)%y + 0.5_wp * dy ) * ddy
 !
@@ -521,12 +514,12 @@
                                particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
                                particles(n)%origin_y = ( ny + 1 ) * dy + &
-                               particles(n)%origin_y
+                                                     particles(n)%origin_y
                             ELSE
-                               trsp_count = trsp_count + 1
-                               trsp(trsp_count) = particles(n)
+                               trsp_count         = trsp_count + 1
+                               trsp(trsp_count)   = particles(n)
                                trsp(trsp_count)%y = ( ny + 1 ) * dy + &
-                               trsp(trsp_count)%y
+                                                 trsp(trsp_count)%y
                                trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
-                               + ( ny + 1 ) * dy
+                                                + ( ny + 1 ) * dy
                                particles(n)%particle_mask = .FALSE.
                                deleted_particles = deleted_particles + 1
@@ -536,5 +529,5 @@
                                   !++ why is 1 subtracted in next statement???
                                   trsp(trsp_count)%origin_y =                        &
-                                  trsp(trsp_count)%origin_y - 1
+                                                  trsp(trsp_count)%origin_y - 1
                                ENDIF
 
@@ -545,5 +538,5 @@
 !--                         Particle absorption
                             particles(n)%particle_mask = .FALSE.
-                            deleted_particles = deleted_particles + 1
+                            deleted_particles          = deleted_particles + 1
 
                          ELSEIF ( ibc_par_ns == 2 )  THEN
@@ -568,21 +561,21 @@
                       IF ( j > ny )  THEN
 !
-!--                       Apply boundary condition along x
+!--                       Apply boundary condition along y
                          IF ( ibc_par_ns == 0 )  THEN
 !
 !--                         Cyclic condition
                             IF ( pdims(2) == 1 )  THEN
-                               particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
+                               particles(n)%y        = particles(n)%y - ( ny + 1 ) * dy
                                particles(n)%origin_y =                         &
-                                  particles(n)%origin_y - ( ny + 1 ) * dy
+                                          particles(n)%origin_y - ( ny + 1 ) * dy
                             ELSE
-                               trnp_count = trnp_count + 1
-                               trnp(trnp_count) = particles(n)
+                               trnp_count         = trnp_count + 1
+                               trnp(trnp_count)   = particles(n)
                                trnp(trnp_count)%y =                            &
-                                  trnp(trnp_count)%y - ( ny + 1 ) * dy
+                                          trnp(trnp_count)%y - ( ny + 1 ) * dy
                                trnp(trnp_count)%origin_y =                     &
-                                  trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
+                                         trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
                                particles(n)%particle_mask = .FALSE.
-                               deleted_particles = deleted_particles + 1
+                               deleted_particles          = deleted_particles + 1
                             ENDIF
 
@@ -628,5 +621,5 @@
 
        ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
-
+ 
        CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type,      &
                           psouth, 1, rvnp(1)%radius,                             &
@@ -661,95 +654,119 @@
     ENDIF
 
+    DEALLOCATE( move_also_north )
+    DEALLOCATE( move_also_south )
+
 #else
 
-!
-!-- Apply boundary conditions
-    DO  n = 1, number_of_particles
-
-       IF ( particles(n)%x < -0.5_wp * dx )  THEN
-
-          IF ( ibc_par_lr == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
-
-          ELSEIF ( ibc_par_lr == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_lr == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%x       = -dx - particles(n)%x
-             particles(n)%speed_x = -particles(n)%speed_x
-          ENDIF
-
-       ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
-
-          IF ( ibc_par_lr == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
-
-          ELSEIF ( ibc_par_lr == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_lr == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
-             particles(n)%speed_x = -particles(n)%speed_x
-          ENDIF
-
-       ENDIF
-
-       IF ( particles(n)%y < -0.5_wp * dy )  THEN
-
-          IF ( ibc_par_ns == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
-
-          ELSEIF ( ibc_par_ns == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_ns == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%y       = -dy - particles(n)%y
-             particles(n)%speed_y = -particles(n)%speed_y
-          ENDIF
-
-       ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
-
-          IF ( ibc_par_ns == 0 )  THEN
-!
-!--          Cyclic boundary. Relevant coordinate has to be changed.
-             particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
-
-          ELSEIF ( ibc_par_ns == 1 )  THEN
-!
-!--          Particle absorption
-             particles(n)%particle_mask = .FALSE.
-             deleted_particles = deleted_particles + 1
-
-          ELSEIF ( ibc_par_ns == 2 )  THEN
-!
-!--          Particle reflection
-             particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
-             particles(n)%speed_y = -particles(n)%speed_y
-          ENDIF
-
-       ENDIF
+    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
+!
+!--             Apply boundary conditions
+
+                IF ( particles(n)%x < -0.5_wp * dx )  THEN
+
+                   IF ( ibc_par_lr == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
+
+                   ELSEIF ( ibc_par_lr == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_lr == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%x       = -dx - particles(n)%x
+                      particles(n)%speed_x = -particles(n)%speed_x
+                   ENDIF
+
+                ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
+
+                   IF ( ibc_par_lr == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
+
+                   ELSEIF ( ibc_par_lr == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_lr == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
+                      particles(n)%speed_x = -particles(n)%speed_x
+                   ENDIF
+
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
     ENDDO
 
+    DO  ip = nxl, nxr
+       DO  jp = nys, nyn, nyn-nys
+          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)%y < -0.5_wp * dy )  THEN
+
+                   IF ( ibc_par_ns == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
+
+                   ELSEIF ( ibc_par_ns == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_ns == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%y       = -dy - particles(n)%y
+                      particles(n)%speed_y = -particles(n)%speed_y
+                   ENDIF
+
+                ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
+
+                   IF ( ibc_par_ns == 0 )  THEN
+!
+!--                   Cyclic boundary. Relevant coordinate has to be changed.
+                      particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
+
+                   ELSEIF ( ibc_par_ns == 1 )  THEN
+!
+!--                   Particle absorption
+                      particles(n)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+
+                   ELSEIF ( ibc_par_ns == 2 )  THEN
+!
+!--                   Particle reflection
+                      particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
+                      particles(n)%speed_y = -particles(n)%speed_y
+                   ENDIF
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDDO
 #endif
 
@@ -782,23 +799,24 @@
     IMPLICIT NONE
 
-    INTEGER(iwp)        ::  ip        !<
-    INTEGER(iwp)        ::  jp        !<
-    INTEGER(iwp)        ::  kp        !<
-    INTEGER(iwp)        ::  n         !<
-    INTEGER(iwp)        ::  pindex    !<
+    INTEGER(iwp)        ::  ip        !< grid index (x) of particle
+    INTEGER(iwp)        ::  jp        !< grid index (x) of particle
+    INTEGER(iwp)        ::  kp        !< grid index (x) of particle
+    INTEGER(iwp)        ::  n         !< index variable of particle
+    INTEGER(iwp)        ::  pindex    !< dummy argument for new number of particles per grid box
 
     LOGICAL             ::  pack_done !<
 
-    TYPE(particle_type), DIMENSION(:), INTENT(IN)       :: particle_array
+    TYPE(particle_type), DIMENSION(:), INTENT(IN)  ::  particle_array !< new particles in a grid box
+    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  temp_ns        !< temporary particle array for reallocation
 
     pack_done     = .FALSE.
 
-    nr_move_north = 0
-    nr_move_south = 0
-
     DO n = 1, SIZE(particle_array)
+
+       IF ( .NOT. particle_array(n)%particle_mask )  CYCLE
+
        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
+       kp =   particle_array(n)%z / dz + 1 + offset_ocean_nzt
 
        IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
@@ -824,7 +842,25 @@
           prt_count(kp,jp,ip) = pindex
        ELSE
-          IF ( jp == nys - 1 )  THEN
+          IF ( jp <= nys - 1 )  THEN
              nr_move_south = nr_move_south+1
+!
+!--          Before particle information is swapped to exchange-array, check 
+!--          if enough memory is allocated. If required, reallocate exchange
+!--          array.
+             IF ( nr_move_south > SIZE(move_also_south) )  THEN
+!
+!--             At first, allocate further temporary array to swap particle 
+!--             information.
+                ALLOCATE( temp_ns(SIZE(move_also_south)+NR_2_direction_move) )
+                temp_ns(1:nr_move_south-1) = move_also_south(1:nr_move_south-1)
+                DEALLOCATE( move_also_south )
+                ALLOCATE( move_also_south(SIZE(temp_ns)) )
+                move_also_south(1:nr_move_south-1) = temp_ns(1:nr_move_south-1)
+                DEALLOCATE( temp_ns )
+
+             ENDIF
+
              move_also_south(nr_move_south) = particle_array(n)
+
              IF ( jp == -1 )  THEN
                 move_also_south(nr_move_south)%y =                             &
@@ -833,6 +869,23 @@
                    move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
              ENDIF
-          ELSEIF ( jp == nyn+1 )  THEN
+          ELSEIF ( jp >= nyn+1 )  THEN
              nr_move_north = nr_move_north+1
+!
+!--          Before particle information is swapped to exchange-array, check 
+!--          if enough memory is allocated. If required, reallocate exchange
+!--          array.
+             IF ( nr_move_north > SIZE(move_also_north) )  THEN
+!
+!--             At first, allocate further temporary array to swap particle 
+!--             information.
+                ALLOCATE( temp_ns(SIZE(move_also_north)+NR_2_direction_move) )
+                temp_ns(1:nr_move_north-1) = move_also_south(1:nr_move_north-1)
+                DEALLOCATE( move_also_north )
+                ALLOCATE( move_also_north(SIZE(temp_ns)) )
+                move_also_north(1:nr_move_north-1) = temp_ns(1:nr_move_north-1)
+                DEALLOCATE( temp_ns )
+
+             ENDIF
+
              move_also_north(nr_move_north) = particle_array(n)
              IF ( jp == ny+1 )  THEN
@@ -867,18 +920,18 @@
     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      !<
+    INTEGER(iwp)        ::  i           !< grid index (x) of particle position
+    INTEGER(iwp)        ::  ip          !< index variable along x
+    INTEGER(iwp)        ::  j           !< grid index (y) of particle position
+    INTEGER(iwp)        ::  jp          !< index variable along y
+    INTEGER(iwp)        ::  k           !< grid index (z) of particle position
+    INTEGER(iwp)        ::  kp          !< index variable along z
+    INTEGER(iwp)        ::  n           !< index variable for particle array
+    INTEGER(iwp)        ::  np_old_cell !< number of particles per grid box before moving
+    INTEGER(iwp)        ::  n_start     !< start index
+    INTEGER(iwp)        ::  pindex      !< dummy argument for number of new particle per grid box
 
     LOGICAL             ::  pack_done   !<
 
-    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !<
+    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !< particles before moving
 
     CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
@@ -947,10 +1000,7 @@
                             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 
@@ -959,5 +1009,4 @@
                                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
 
@@ -987,5 +1036,7 @@
 ! Description:
 ! ------------
-!> @todo Missing subroutine description.
+!> If the allocated memory for the particle array do not suffice to add arriving
+!> particles from neighbour grid cells, this subrouting reallocates the 
+!> particle array to assure enough memory is available. 
 !------------------------------------------------------------------------------!
  SUBROUTINE realloc_particles_array (i,j,k,size_in)
@@ -1012,5 +1063,5 @@
     ENDIF
 
-    new_size = MAX( new_size, min_nr_particle )
+    new_size = MAX( new_size, min_nr_particle, old_size + 1 )
 
     IF ( old_size <= 500 )  THEN
