Index: palm/trunk/SOURCE/lpm_exchange_horiz.f90
===================================================================
--- palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 2330)
+++ palm/trunk/SOURCE/lpm_exchange_horiz.f90	(revision 2606)
@@ -25,4 +25,15 @@
 ! -----------------
 ! $Id$
+! Changed particle box locations: center of particle box now coincides 
+! with scalar grid point of same index.
+! lpm_move_particle now moves inner particles that would leave the core 
+! to the respective boundary gridbox. Afterwards they get transferred by 
+! lpm_exchange_horiz. This allows particles to move more than one gridbox 
+! independent of domain composition.  
+! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
+! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
+! lpm_sort -> lpm_sort_timeloop_done
+! 
+! 2330 2017-08-03 14:26:02Z schwenkel
 ! Bugfix: Also for gfortran compilable, function c_sizeof is not used. 
 ! 
@@ -128,6 +139,6 @@
     USE kinds
 
-    USE lpm_pack_arrays_mod,                                                   &
-        ONLY:  lpm_pack_arrays
+    USE lpm_pack_and_sort_mod,                                                   &
+        ONLY:  lpm_pack
 
     USE netcdf_interface,                                                      &
@@ -255,8 +266,8 @@
                 DO  n = 1, number_of_particles
                    IF ( particles(n)%particle_mask )  THEN
-                      i = ( particles(n)%x + 0.5_wp * dx ) * ddx
+                      i = particles(n)%x * ddx
 !
 !--                   Above calculation does not work for indices less than zero
-                      IF ( particles(n)%x < -0.5_wp * dx )  i = -1
+                      IF ( particles(n)%x < 0.0_wp)  i = -1
 
                       IF ( i < nxl )  THEN
@@ -298,8 +309,8 @@
                 IF ( particles(n)%particle_mask )  THEN
 
-                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
+                   i = particles(n)%x * ddx
 !
 !--                Above calculation does not work for indices less than zero
-                   IF ( particles(n)%x < - 0.5_wp * dx )  i = -1
+                   IF ( particles(n)%x < 0.0_wp )  i = -1
 
                    IF ( i <  nxl )  THEN
@@ -323,5 +334,5 @@
                                deleted_particles = deleted_particles + 1
 
-                               IF ( trlp(trlp_count)%x >= (nx + 0.5_wp)* dx - 1.0E-12_wp )  THEN
+                               IF ( trlp(trlp_count)%x >= (nx + 1)* 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???
@@ -501,8 +512,8 @@
                 DO  n = 1, number_of_particles
                    IF ( particles(n)%particle_mask )  THEN
-                      j = ( particles(n)%y + 0.5_wp * dy ) * ddy
+                      j = particles(n)%y * ddy
 !
 !--                   Above calculation does not work for indices less than zero
-                      IF ( particles(n)%y < -0.5_wp * dy )  j = -1
+                      IF ( particles(n)%y < 0.0_wp)  j = -1
 
                       IF ( j < nys )  THEN
@@ -545,8 +556,8 @@
                 IF ( particles(n)%particle_mask )  THEN
 
-                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
+                   j = particles(n)%y * ddy
 !
 !--                Above calculation does not work for indices less than zero
-                   IF ( particles(n)%y < -0.5_wp * dy )  j = -1
+                   IF ( particles(n)%y < 0.0_wp )  j = -1
 
                    IF ( j < nys )  THEN
@@ -571,5 +582,5 @@
                                deleted_particles = deleted_particles + 1
 
-                               IF ( trsp(trsp_count)%y >= (ny+0.5_wp)* dy - 1.0E-12_wp )  THEN
+                               IF ( trsp(trsp_count)%y >= (ny+1)* 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???
@@ -719,5 +730,5 @@
 !--             Apply boundary conditions
 
-                IF ( particles(n)%x < -0.5_wp * dx )  THEN
+                IF ( particles(n)%x < 0.0_wp )  THEN
 
                    IF ( ibc_par_lr == 0 )  THEN
@@ -725,5 +736,6 @@
 !--                   Cyclic boundary. Relevant coordinate has to be changed.
                       particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
-
+                      particles(n)%origin_x = ( nx + 1 ) * dx + &
+                               particles(n)%origin_x
                    ELSEIF ( ibc_par_lr == 1 )  THEN
 !
@@ -739,5 +751,5 @@
                    ENDIF
 
-                ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
+                ELSEIF ( particles(n)%x >= ( nx + 1) * dx )  THEN
 
                    IF ( ibc_par_lr == 0 )  THEN
@@ -745,4 +757,6 @@
 !--                   Cyclic boundary. Relevant coordinate has to be changed.
                       particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
+                      particles(n)%origin_x = particles(n)%origin_x - &
+                               ( nx + 1 ) * dx
 
                    ELSEIF ( ibc_par_lr == 1 )  THEN
@@ -773,5 +787,5 @@
              DO  n = 1, number_of_particles
 
-                IF ( particles(n)%y < -0.5_wp * dy )  THEN
+                IF ( particles(n)%y < 0.0_wp)  THEN
 
                    IF ( ibc_par_ns == 0 )  THEN
@@ -779,4 +793,6 @@
 !--                   Cyclic boundary. Relevant coordinate has to be changed.
                       particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
+                      particles(n)%origin_y = ( ny + 1 ) * dy + &
+                           particles(n)%origin_y
 
                    ELSEIF ( ibc_par_ns == 1 )  THEN
@@ -793,5 +809,5 @@
                    ENDIF
 
-                ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
+                ELSEIF ( particles(n)%y >= ( ny + 1) * dy )  THEN
 
                    IF ( ibc_par_ns == 0 )  THEN
@@ -799,4 +815,6 @@
 !--                   Cyclic boundary. Relevant coordinate has to be changed.
                       particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
+                      particles(n)%origin_y = particles(n)%origin_y - &
+                                ( ny + 1 ) * dy
 
                    ELSEIF ( ibc_par_ns == 1 )  THEN
@@ -866,7 +884,7 @@
        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
+       ip = particle_array(n)%x * ddx
+       jp = particle_array(n)%y * ddy
+       kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
 
        IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
@@ -880,5 +898,5 @@
                 CALL realloc_particles_array (ip,jp,kp)
              ELSE
-                CALL lpm_pack_arrays
+                CALL lpm_pack
                 prt_count(kp,jp,ip) = number_of_particles
                 pindex = prt_count(kp,jp,ip)+1
@@ -1014,11 +1032,8 @@
     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)        ::  np_before_move !< number of particles per grid box before moving
     INTEGER(iwp)        ::  pindex      !< dummy argument for number of new particle per grid box
 
-    LOGICAL             ::  pack_done   !<
-
-    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !< particles before moving
+    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_before_move !< particles before moving
 
     CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
@@ -1028,90 +1043,44 @@
           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
+             np_before_move = prt_count(kp,jp,ip)
+             IF ( np_before_move <= 0 )  CYCLE
+             particles_before_move => grid_particles(kp,jp,ip)%particles(1:np_before_move)
              
-             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
+             DO  n = 1, np_before_move
+                i = particles_before_move(n)%x * ddx
+                j = particles_before_move(n)%y * ddy
+                k = particles_before_move(n)%z / dz + 1 + offset_ocean_nzt
+
+!--             For lpm_exchange_horiz to work properly particles need to be moved to the outermost gridboxes
+!--             of the respective processor. If the particle index is inside the processor the following lines
+!--             will not change the index
+                i = MIN ( i , nxr )
+                i = MAX ( i , nxl )
+                j = MIN ( j , nyn )
+                j = MAX ( j , nys )
+                k = MIN ( k , nzt )
+                k = MAX ( k , nzb+1 )
 !
 !--             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.
+!!
+!--                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
+                      CALL realloc_particles_array(i,j,k)
                    ENDIF
+
+                   grid_particles(k,j,i)%particles(pindex) = particles_before_move(n)
+                   prt_count(k,j,i) = pindex
+
+                   particles_before_move(n)%particle_mask = .FALSE.
                 ENDIF
              ENDDO
 
-             IF ( n_start >= 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)
-                            ELSE
-                               CALL lpm_pack_arrays
-                               prt_count(k,j,i) = number_of_particles
-!
-!--                            Set pindex to next index of the modified
-!--                            particle-array (due to lpm_pack_arrays)
-                               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)
-                               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
