source: palm/trunk/SOURCE/lpm_exchange_horiz.f90 @ 1736

Last change on this file since 1736 was 1692, checked in by maronga, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 54.4 KB
Line 
1!> @file lpm_exchange_horiz.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2015 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: lpm_exchange_horiz.f90 1692 2015-10-26 16:29:17Z raasch $
26!
27! 1691 2015-10-26 16:17:44Z maronga
28! Formatting corrections.
29!
30! 1685 2015-10-08 07:32:13Z raasch
31! bugfix concerning vertical index offset in case of ocean
32!
33! 1682 2015-10-07 23:56:08Z knoop
34! Code annotations made doxygen readable
35!
36! 1359 2014-04-11 17:15:14Z hoffmann
37! New particle structure integrated.
38! Kind definition added to all floating point numbers.
39!
40! 1327 2014-03-21 11:00:16Z raasch
41! -netcdf output queries
42!
43! 1320 2014-03-20 08:40:49Z raasch
44! ONLY-attribute added to USE-statements,
45! kind-parameters added to all INTEGER and REAL declaration statements,
46! kinds are defined in new module kinds,
47! comment fields (!:) to be used for variable explanations added to
48! all variable declaration statements
49!
50! 1318 2014-03-17 13:35:16Z raasch
51! module interfaces removed
52!
53! 1036 2012-10-22 13:43:42Z raasch
54! code put under GPL (PALM 3.9)
55!
56! 851 2012-03-15 14:32:58Z raasch
57! Bugfix: resetting of particle_mask and tail mask moved from end of this
58! routine to lpm
59!
60! 849 2012-03-15 10:35:09Z raasch
61! initial revision (former part of advec_particles)
62!
63!
64! Description:
65! ------------
66!> Exchange of particles (and tails) between the subdomains.
67!------------------------------------------------------------------------------!
68 MODULE lpm_exchange_horiz_mod
69 
70
71    USE control_parameters,                                                    &
72        ONLY:  dz, message_string, netcdf_data_format, simulated_time
73
74    USE cpulog,                                                                &
75        ONLY:  cpu_log, log_point_s
76
77    USE grid_variables,                                                        &
78        ONLY:  ddx, ddy, dx, dy
79
80    USE indices,                                                               &
81        ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb, nzt
82
83    USE kinds
84
85    USE lpm_pack_arrays_mod,                                                   &
86        ONLY:  lpm_pack_arrays
87
88    USE particle_attributes,                                                   &
89        ONLY:  alloc_factor, deleted_particles, deleted_tails, grid_particles, &
90               ibc_par_lr, ibc_par_ns, maximum_number_of_tails,                &
91               maximum_number_of_tailpoints, min_nr_particle,                  &
92               mpi_particle_type, number_of_tails, number_of_particles,        &
93               offset_ocean_nzt, particles,                                    &
94               particle_tail_coordinates, particle_type, prt_count,            &
95               tail_mask, trlp_count_sum,                                      &
96               trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum,       &
97               trrp_count_sum, trrp_count_recv_sum, trsp_count_sum,            &
98               trsp_count_recv_sum, use_particle_tails, zero_particle
99
100    USE pegrid
101
102    IMPLICIT NONE
103
104    INTEGER(iwp), PARAMETER ::  NR_2_direction_move = 10000 !<
105    INTEGER(iwp)            ::  nr_move_north               !<
106    INTEGER(iwp)            ::  nr_move_south               !<
107
108    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_north
109    TYPE(particle_type), DIMENSION(NR_2_direction_move) ::  move_also_south
110
111    SAVE
112
113    PRIVATE
114    PUBLIC lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
115
116    INTERFACE lpm_exchange_horiz
117       MODULE PROCEDURE lpm_exchange_horiz
118    END INTERFACE lpm_exchange_horiz
119
120    INTERFACE lpm_move_particle
121       MODULE PROCEDURE lpm_move_particle
122    END INTERFACE lpm_move_particle
123
124    INTERFACE realloc_particles_array
125       MODULE PROCEDURE realloc_particles_array
126    END INTERFACE realloc_particles_array
127
128CONTAINS
129
130!------------------------------------------------------------------------------!
131! Description:
132! ------------
133!> Exchange between subdomains.
134!> As soon as one particle has moved beyond the boundary of the domain, it
135!> is included in the relevant transfer arrays and marked for subsequent
136!> deletion on this PE.
137!> First sweep for crossings in x direction. Find out first the number of
138!> particles to be transferred and allocate temporary arrays needed to store
139!> them.
140!> For a one-dimensional decomposition along y, no transfer is necessary,
141!> because the particle remains on the PE, but the particle coordinate has to
142!> be adjusted.
143!------------------------------------------------------------------------------!
144 SUBROUTINE lpm_exchange_horiz
145
146    IMPLICIT NONE
147
148    INTEGER(iwp) ::  i                                       !<
149    INTEGER(iwp) ::  ip                                      !<
150    INTEGER(iwp) ::  j                                       !<
151    INTEGER(iwp) ::  jp                                      !<
152    INTEGER(iwp) ::  k                                       !<
153    INTEGER(iwp) ::  kp                                      !<
154    INTEGER(iwp) ::  n                                       !<
155    INTEGER(iwp) ::  nn                                      !<
156    INTEGER(iwp) ::  tlength                                 !<
157    INTEGER(iwp) ::  trlp_count                              !<
158    INTEGER(iwp) ::  trlp_count_recv                         !<
159    INTEGER(iwp) ::  trlpt_count                             !<
160    INTEGER(iwp) ::  trlpt_count_recv                        !<
161    INTEGER(iwp) ::  trnp_count                              !<
162    INTEGER(iwp) ::  trnp_count_recv                         !<
163    INTEGER(iwp) ::  trnpt_count                             !<
164    INTEGER(iwp) ::  trnpt_count_recv                        !<
165    INTEGER(iwp) ::  trrp_count                              !<
166    INTEGER(iwp) ::  trrp_count_recv                         !<
167    INTEGER(iwp) ::  trrpt_count                             !<
168    INTEGER(iwp) ::  trrpt_count_recv                        !<
169    INTEGER(iwp) ::  trsp_count                              !<
170    INTEGER(iwp) ::  trsp_count_recv                         !<
171    INTEGER(iwp) ::  trspt_count                             !<
172    INTEGER(iwp) ::  trspt_count_recv                        !<
173
174    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trlpt        !<
175    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trnpt        !<
176    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trrpt        !<
177    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  trspt        !<
178
179    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !<
180    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !<
181    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !<
182    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !<
183    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !<
184    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !<
185    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !<
186    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !<
187
188    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
189
190#if defined( __parallel )
191
192    trlp_count  = 0
193    trlpt_count = 0
194    trrp_count  = 0
195    trrpt_count = 0
196
197    trlp_count_recv   = 0
198    trlpt_count_recv  = 0
199    trrp_count_recv   = 0
200    trrpt_count_recv  = 0
201
202    IF ( pdims(1) /= 1 )  THEN
203!
204!--    First calculate the storage necessary for sending and receiving the data.
205!--    Compute only first (nxl) and last (nxr) loop iterration.
206       DO  ip = nxl, nxr, nxr - nxl
207          DO  jp = nys, nyn
208             DO  kp = nzb+1, nzt
209
210                number_of_particles = prt_count(kp,jp,ip)
211                IF ( number_of_particles <= 0 )  CYCLE
212                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
213                DO  n = 1, number_of_particles
214                   IF ( particles(n)%particle_mask )  THEN
215                      i = ( particles(n)%x + 0.5_wp * dx ) * ddx
216   !
217   !--                Above calculation does not work for indices less than zero
218                      IF ( particles(n)%x < -0.5_wp * dx )  i = -1
219
220                      IF ( i < nxl )  THEN
221                         trlp_count = trlp_count + 1
222                         IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
223                      ELSEIF ( i > nxr )  THEN
224                         trrp_count = trrp_count + 1
225                         IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
226                      ENDIF
227                   ENDIF
228                ENDDO
229
230             ENDDO
231          ENDDO
232       ENDDO
233
234       IF ( trlp_count  == 0 )  trlp_count  = 1
235       IF ( trlpt_count == 0 )  trlpt_count = 1
236       IF ( trrp_count  == 0 )  trrp_count  = 1
237       IF ( trrpt_count == 0 )  trrpt_count = 1
238
239       ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
240
241       trlp = zero_particle
242       trrp = zero_particle
243
244       IF ( use_particle_tails )  THEN
245          ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
246                    trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
247          tlength = maximum_number_of_tailpoints * 5
248       ENDIF
249
250       trlp_count  = 0
251       trlpt_count = 0
252       trrp_count  = 0
253       trrpt_count = 0
254
255    ENDIF
256!
257!-- Compute only first (nxl) and last (nxr) loop iterration
258    DO  ip = nxl, nxr,nxr-nxl
259       DO  jp = nys, nyn
260          DO  kp = nzb+1, nzt
261             number_of_particles = prt_count(kp,jp,ip)
262             IF ( number_of_particles <= 0 ) CYCLE
263             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
264             DO  n = 1, number_of_particles
265
266                nn = particles(n)%tail_id
267!
268!--             Only those particles that have not been marked as 'deleted' may
269!--             be moved.
270                IF ( particles(n)%particle_mask )  THEN
271
272                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
273   !
274   !--             Above calculation does not work for indices less than zero
275                   IF ( particles(n)%x < - 0.5_wp * dx )  i = -1
276
277                   IF ( i <  nxl )  THEN
278                      IF ( i < 0 )  THEN
279   !
280   !--                   Apply boundary condition along x
281                         IF ( ibc_par_lr == 0 )  THEN
282   !
283   !--                      Cyclic condition
284                            IF ( pdims(1) == 1 )  THEN
285                               particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
286                               particles(n)%origin_x = ( nx + 1 ) * dx + &
287                               particles(n)%origin_x
288                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
289                                  i  = particles(n)%tailpoints
290                                  particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
291                                  + particle_tail_coordinates(1:i,1,nn)
292                               ENDIF
293                            ELSE
294                               trlp_count = trlp_count + 1
295                               trlp(trlp_count)   = particles(n)
296                               trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
297                               trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
298                               ( nx + 1 ) * dx
299                               particles(n)%particle_mask  = .FALSE.
300                               deleted_particles = deleted_particles + 1
301
302                               IF ( trlp(trlp_count)%x >= (nx + 0.5_wp)* dx - 1.0E-12_wp )  THEN
303                                  trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10_wp
304                                  !++ why is 1 subtracted in next statement???
305                                  trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x - 1
306                               ENDIF
307
308                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
309                                  trlpt_count = trlpt_count + 1
310                                  trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
311                                  trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
312                                  trlpt(:,1,trlpt_count)
313                                  tail_mask(nn) = .FALSE.
314                                  deleted_tails = deleted_tails + 1
315                               ENDIF
316                            ENDIF
317
318                         ELSEIF ( ibc_par_lr == 1 )  THEN
319   !
320   !--                      Particle absorption
321                            particles(n)%particle_mask = .FALSE.
322                            deleted_particles = deleted_particles + 1
323                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
324                               tail_mask(nn) = .FALSE.
325                               deleted_tails = deleted_tails + 1
326                            ENDIF
327
328                         ELSEIF ( ibc_par_lr == 2 )  THEN
329   !
330   !--                      Particle reflection
331                            particles(n)%x       = -particles(n)%x
332                            particles(n)%speed_x = -particles(n)%speed_x
333
334                         ENDIF
335                      ELSE
336   !
337   !--                   Store particle data in the transfer array, which will be
338   !--                   send to the neighbouring PE
339                         trlp_count = trlp_count + 1
340                         trlp(trlp_count) = particles(n)
341                         particles(n)%particle_mask = .FALSE.
342                         deleted_particles = deleted_particles + 1
343
344                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
345                            trlpt_count = trlpt_count + 1
346                            trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
347                            tail_mask(nn) = .FALSE.
348                            deleted_tails = deleted_tails + 1
349                         ENDIF
350                      ENDIF
351
352                   ELSEIF ( i > nxr )  THEN
353                      IF ( i > nx )  THEN
354   !
355   !--                   Apply boundary condition along x
356                         IF ( ibc_par_lr == 0 )  THEN
357   !
358   !--                      Cyclic condition
359                            IF ( pdims(1) == 1 )  THEN
360                               particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
361                               particles(n)%origin_x = particles(n)%origin_x - &
362                               ( nx + 1 ) * dx
363                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
364                                  i = particles(n)%tailpoints
365                                  particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
366                                  + particle_tail_coordinates(1:i,1,nn)
367                               ENDIF
368                            ELSE
369                               trrp_count = trrp_count + 1
370                               trrp(trrp_count) = particles(n)
371                               trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
372                               trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
373                               ( nx + 1 ) * dx
374                               particles(n)%particle_mask = .FALSE.
375                               deleted_particles = deleted_particles + 1
376
377                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
378                                  trrpt_count = trrpt_count + 1
379                                  trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
380                                  trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
381                                  ( nx + 1 ) * dx
382                                  tail_mask(nn) = .FALSE.
383                                  deleted_tails = deleted_tails + 1
384                               ENDIF
385                            ENDIF
386
387                         ELSEIF ( ibc_par_lr == 1 )  THEN
388   !
389   !--                      Particle absorption
390                            particles(n)%particle_mask = .FALSE.
391                            deleted_particles = deleted_particles + 1
392                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
393                               tail_mask(nn) = .FALSE.
394                               deleted_tails = deleted_tails + 1
395                            ENDIF
396
397                         ELSEIF ( ibc_par_lr == 2 )  THEN
398   !
399   !--                      Particle reflection
400                            particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
401                            particles(n)%speed_x = -particles(n)%speed_x
402
403                         ENDIF
404                      ELSE
405   !
406   !--                   Store particle data in the transfer array, which will be send
407   !--                   to the neighbouring PE
408                         trrp_count = trrp_count + 1
409                         trrp(trrp_count) = particles(n)
410                         particles(n)%particle_mask = .FALSE.
411                         deleted_particles = deleted_particles + 1
412
413                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
414                            trrpt_count = trrpt_count + 1
415                            trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
416                            tail_mask(nn) = .FALSE.
417                            deleted_tails = deleted_tails + 1
418                         ENDIF
419                      ENDIF
420
421                   ENDIF
422                ENDIF
423             ENDDO
424          ENDDO
425       ENDDO
426    ENDDO
427
428!
429!-- Send left boundary, receive right boundary (but first exchange how many
430!-- and check, if particle storage must be extended)
431    IF ( pdims(1) /= 1 )  THEN
432
433       CALL MPI_SENDRECV( trlp_count,      1, MPI_INTEGER, pleft,  0, &
434                          trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
435                          comm2d, status, ierr )
436
437       ALLOCATE(rvrp(MAX(1,trrp_count_recv)))
438
439       CALL MPI_SENDRECV( trlp(1)%radius, max(1,trlp_count), mpi_particle_type,&
440                          pleft, 1, rvrp(1)%radius,                            &
441                          max(1,trrp_count_recv), mpi_particle_type, pright, 1,&
442                          comm2d, status, ierr )
443
444       IF ( trrp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvrp(1:trrp_count_recv))
445
446       DEALLOCATE(rvrp)
447
448       IF ( use_particle_tails )  THEN
449
450          CALL MPI_SENDRECV( trlpt_count,      1, MPI_INTEGER, pleft,  0, &
451                             trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
452                             comm2d, status, ierr )
453
454          IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
455          THEN
456             IF ( netcdf_data_format < 3 )  THEN
457                message_string = 'maximum_number_of_tails ' //   &
458                                 'needs to be increased ' //     &
459                                 '&but this is not allowed wi'// &
460                                 'th netcdf_data_format < 3'
461                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
462             ELSE
463                CALL lpm_extend_tail_array( trrpt_count_recv )
464             ENDIF
465          ENDIF
466
467          CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,      &
468                             pleft, 1,                                         &
469                             particle_tail_coordinates(1,1,number_of_tails+1), &
470                             trrpt_count_recv*tlength, MPI_REAL, pright, 1,    &
471                             comm2d, status, ierr )
472!
473!--       Update the tail ids for the transferred particles
474          nn = number_of_tails
475          DO  n = number_of_particles+1, number_of_particles+trrp_count_recv
476             IF ( particles(n)%tail_id /= 0 )  THEN
477                nn = nn + 1
478                particles(n)%tail_id = nn
479             ENDIF
480          ENDDO
481
482       ENDIF
483
484!
485!--    Send right boundary, receive left boundary
486       CALL MPI_SENDRECV( trrp_count,      1, MPI_INTEGER, pright, 0, &
487                          trlp_count_recv, 1, MPI_INTEGER, pleft,  0, &
488                          comm2d, status, ierr )
489
490       ALLOCATE(rvlp(MAX(1,trlp_count_recv)))
491
492       CALL MPI_SENDRECV( trrp(1)%radius, max(1,trrp_count), mpi_particle_type,&
493                          pright, 1, rvlp(1)%radius,                           &
494                          max(1,trlp_count_recv), mpi_particle_type, pleft, 1, &
495                          comm2d, status, ierr )
496
497       IF ( trlp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvlp(1:trlp_count_recv))
498
499       DEALLOCATE(rvlp)
500
501       IF ( use_particle_tails )  THEN
502
503          CALL MPI_SENDRECV( trrpt_count,      1, MPI_INTEGER, pright, 0, &
504                             trlpt_count_recv, 1, MPI_INTEGER, pleft,  0, &
505                             comm2d, status, ierr )
506
507          IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
508          THEN
509             IF ( netcdf_data_format < 3 )  THEN
510                message_string = 'maximum_number_of_tails ' //   &
511                                 'needs to be increased ' //     &
512                                 '&but this is not allowed wi'// &
513                                 'th netcdf_data_format < 3'
514                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
515             ELSE
516                CALL lpm_extend_tail_array( trlpt_count_recv )
517             ENDIF
518          ENDIF
519
520          CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,      &
521                             pright, 1,                                        &
522                             particle_tail_coordinates(1,1,number_of_tails+1), &
523                             trlpt_count_recv*tlength, MPI_REAL, pleft, 1,     &
524                             comm2d, status, ierr )
525!
526!--       Update the tail ids for the transferred particles
527          nn = number_of_tails
528          DO  n = number_of_particles+1, number_of_particles+trlp_count_recv
529             IF ( particles(n)%tail_id /= 0 )  THEN
530                nn = nn + 1
531                particles(n)%tail_id = nn
532             ENDIF
533          ENDDO
534
535       ENDIF
536
537!       number_of_particles = number_of_particles + trlp_count_recv
538!       number_of_tails     = number_of_tails     + trlpt_count_recv
539
540       IF ( use_particle_tails )  THEN
541          DEALLOCATE( trlpt, trrpt )
542       ENDIF
543       DEALLOCATE( trlp, trrp )
544
545    ENDIF
546
547
548!
549!-- Check whether particles have crossed the boundaries in y direction. Note
550!-- that this case can also apply to particles that have just been received
551!-- from the adjacent right or left PE.
552!-- Find out first the number of particles to be transferred and allocate
553!-- temporary arrays needed to store them.
554!-- For a one-dimensional decomposition along x, no transfer is necessary,
555!-- because the particle remains on the PE.
556    trsp_count  = nr_move_south
557    trspt_count = 0
558    trnp_count  = nr_move_north
559    trnpt_count = 0
560
561    trsp_count_recv   = 0
562    trspt_count_recv  = 0
563    trnp_count_recv   = 0
564    trnpt_count_recv  = 0
565
566    IF ( pdims(2) /= 1 )  THEN
567!
568!--    First calculate the storage necessary for sending and receiving the
569!--    data
570       DO  ip = nxl, nxr
571          DO  jp = nys, nyn, nyn-nys                                 !compute only first (nys) and last (nyn) loop iterration
572             DO  kp = nzb+1, nzt
573                number_of_particles = prt_count(kp,jp,ip)
574                IF ( number_of_particles <= 0 )  CYCLE
575                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
576                DO  n = 1, number_of_particles
577                   IF ( particles(n)%particle_mask )  THEN
578                      j = ( particles(n)%y + 0.5_wp * dy ) * ddy
579!
580!--                   Above calculation does not work for indices less than zero
581                      IF ( particles(n)%y < -0.5_wp * dy )  j = -1
582
583                      IF ( j < nys )  THEN
584                         trsp_count = trsp_count + 1
585                         IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count + 1
586                      ELSEIF ( j > nyn )  THEN
587                         trnp_count = trnp_count + 1
588                         IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count + 1
589                      ENDIF
590                   ENDIF
591                ENDDO
592             ENDDO
593          ENDDO
594       ENDDO
595
596       IF ( trsp_count  == 0 )  trsp_count  = 1
597       IF ( trspt_count == 0 )  trspt_count = 1
598       IF ( trnp_count  == 0 )  trnp_count  = 1
599       IF ( trnpt_count == 0 )  trnpt_count = 1
600
601       ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
602
603       trsp = zero_particle
604       trnp = zero_particle
605
606       IF ( use_particle_tails )  THEN
607          ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
608                    trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
609          tlength = maximum_number_of_tailpoints * 5
610       ENDIF
611
612       trsp_count  = nr_move_south
613       trspt_count = 0
614       trnp_count  = nr_move_north
615       trnpt_count = 0
616
617       trsp(1:nr_move_south) = move_also_south(1:nr_move_south)
618       trnp(1:nr_move_north) = move_also_north(1:nr_move_north)
619
620    ENDIF
621
622    DO  ip = nxl, nxr
623       DO  jp = nys, nyn, nyn-nys ! compute only first (nys) and last (nyn) loop iterration
624          DO  kp = nzb+1, nzt
625             number_of_particles = prt_count(kp,jp,ip)
626             IF ( number_of_particles <= 0 )  CYCLE
627             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
628             DO  n = 1, number_of_particles
629
630                nn = particles(n)%tail_id
631!
632!--             Only those particles that have not been marked as 'deleted' may
633!--             be moved.
634                IF ( particles(n)%particle_mask )  THEN
635                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
636!
637!--                Above calculation does not work for indices less than zero
638                   IF ( particles(n)%y < -0.5_wp * dy )  j = -1
639
640                   IF ( j < nys )  THEN
641                      IF ( j < 0 )  THEN
642!
643!--                      Apply boundary condition along y
644                         IF ( ibc_par_ns == 0 )  THEN
645!
646!--                         Cyclic condition
647                            IF ( pdims(2) == 1 )  THEN
648                               particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
649                               particles(n)%origin_y = ( ny + 1 ) * dy + &
650                               particles(n)%origin_y
651                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
652                                  i = particles(n)%tailpoints
653                                  particle_tail_coordinates(1:i,2,nn) =        &
654                                     ( ny+1 ) * dy + particle_tail_coordinates(1:i,2,nn)
655                               ENDIF
656                            ELSE
657                               trsp_count = trsp_count + 1
658                               trsp(trsp_count) = particles(n)
659                               trsp(trsp_count)%y = ( ny + 1 ) * dy + &
660                               trsp(trsp_count)%y
661                               trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
662                               + ( ny + 1 ) * dy
663                               particles(n)%particle_mask = .FALSE.
664                               deleted_particles = deleted_particles + 1
665
666                               IF ( trsp(trsp_count)%y >= (ny+0.5_wp)* dy - 1.0E-12_wp )  THEN
667                                  trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10_wp
668                                  !++ why is 1 subtracted in next statement???
669                                  trsp(trsp_count)%origin_y =                        &
670                                  trsp(trsp_count)%origin_y - 1
671                               ENDIF
672
673                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
674                                  trspt_count = trspt_count + 1
675                                  trspt(:,:,trspt_count) = &
676                                  particle_tail_coordinates(:,:,nn)
677                                  trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
678                                  trspt(:,2,trspt_count)
679                                  tail_mask(nn) = .FALSE.
680                                  deleted_tails = deleted_tails + 1
681                               ENDIF
682                            ENDIF
683
684                         ELSEIF ( ibc_par_ns == 1 )  THEN
685!
686!--                         Particle absorption
687                            particles(n)%particle_mask = .FALSE.
688                            deleted_particles = deleted_particles + 1
689                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
690                               tail_mask(nn) = .FALSE.
691                               deleted_tails = deleted_tails + 1
692                            ENDIF
693
694                         ELSEIF ( ibc_par_ns == 2 )  THEN
695!
696!--                         Particle reflection
697                            particles(n)%y       = -particles(n)%y
698                            particles(n)%speed_y = -particles(n)%speed_y
699
700                         ENDIF
701                      ELSE
702!
703!--                      Store particle data in the transfer array, which will
704!--                      be send to the neighbouring PE
705                         trsp_count = trsp_count + 1
706                         trsp(trsp_count) = particles(n)
707                         particles(n)%particle_mask = .FALSE.
708                         deleted_particles = deleted_particles + 1
709
710                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
711                            trspt_count = trspt_count + 1
712                            trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
713                            tail_mask(nn) = .FALSE.
714                            deleted_tails = deleted_tails + 1
715                         ENDIF
716                      ENDIF
717
718                   ELSEIF ( j > nyn )  THEN
719                      IF ( j > ny )  THEN
720!
721!--                       Apply boundary condition along x
722                         IF ( ibc_par_ns == 0 )  THEN
723!
724!--                         Cyclic condition
725                            IF ( pdims(2) == 1 )  THEN
726                               particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
727                               particles(n)%origin_y =                         &
728                                  particles(n)%origin_y - ( ny + 1 ) * dy
729                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
730                                  i = particles(n)%tailpoints
731                                  particle_tail_coordinates(1:i,2,nn) =        &
732                                     - (ny+1) * dy + particle_tail_coordinates(1:i,2,nn)
733                               ENDIF
734                            ELSE
735                               trnp_count = trnp_count + 1
736                               trnp(trnp_count) = particles(n)
737                               trnp(trnp_count)%y =                            &
738                                  trnp(trnp_count)%y - ( ny + 1 ) * dy
739                               trnp(trnp_count)%origin_y =                     &
740                                  trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
741                               particles(n)%particle_mask = .FALSE.
742                               deleted_particles = deleted_particles + 1
743                               IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
744                                  trnpt_count = trnpt_count + 1
745                                  trnpt(:,:,trnpt_count) =                     &
746                                     particle_tail_coordinates(:,:,nn)
747                                  trnpt(:,2,trnpt_count) =                     &
748                                     trnpt(:,2,trnpt_count) - ( ny + 1 ) * dy
749                                  tail_mask(nn) = .FALSE.
750                                  deleted_tails = deleted_tails + 1
751                               ENDIF
752                            ENDIF
753
754                         ELSEIF ( ibc_par_ns == 1 )  THEN
755!
756!--                         Particle absorption
757                            particles(n)%particle_mask = .FALSE.
758                            deleted_particles = deleted_particles + 1
759                            IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
760                               tail_mask(nn) = .FALSE.
761                               deleted_tails = deleted_tails + 1
762                            ENDIF
763
764                         ELSEIF ( ibc_par_ns == 2 )  THEN
765!
766!--                         Particle reflection
767                            particles(n)%y       = 2 * ( ny * dy ) - particles(n)%y
768                            particles(n)%speed_y = -particles(n)%speed_y
769
770                         ENDIF
771                      ELSE
772!
773!--                      Store particle data in the transfer array, which will
774!--                      be send to the neighbouring PE
775                         trnp_count = trnp_count + 1
776                         trnp(trnp_count) = particles(n)
777                         particles(n)%particle_mask = .FALSE.
778                         deleted_particles = deleted_particles + 1
779
780                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
781                            trnpt_count = trnpt_count + 1
782                            trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
783                            tail_mask(nn) = .FALSE.
784                            deleted_tails = deleted_tails + 1
785                         ENDIF
786                      ENDIF
787
788                   ENDIF
789                ENDIF
790             ENDDO
791          ENDDO
792       ENDDO
793    ENDDO
794
795!
796!-- Send front boundary, receive back boundary (but first exchange how many
797!-- and check, if particle storage must be extended)
798    IF ( pdims(2) /= 1 )  THEN
799
800       CALL MPI_SENDRECV( trsp_count,      1, MPI_INTEGER, psouth, 0, &
801                          trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
802                          comm2d, status, ierr )
803
804       ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
805
806       CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type,      &
807                          psouth, 1, rvnp(1)%radius,                             &
808                          trnp_count_recv, mpi_particle_type, pnorth, 1,   &
809                          comm2d, status, ierr )
810
811       IF ( trnp_count_recv  > 0 )  CALL Add_particles_to_gridcell(rvnp(1:trnp_count_recv))
812
813       DEALLOCATE(rvnp)
814
815       IF ( use_particle_tails )  THEN
816
817          CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
818                             trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
819                             comm2d, status, ierr )
820
821          IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
822          THEN
823             IF ( netcdf_data_format < 3 )  THEN
824                message_string = 'maximum_number_of_tails ' //    &
825                                 'needs to be increased ' //      &
826                                 '&but this is not allowed wi' // &
827                                 'th netcdf_data_format < 3'
828                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
829             ELSE
830                CALL lpm_extend_tail_array( trnpt_count_recv )
831             ENDIF
832          ENDIF
833
834          CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,      &
835                             psouth, 1,                                        &
836                             particle_tail_coordinates(1,1,number_of_tails+1), &
837                             trnpt_count_recv*tlength, MPI_REAL, pnorth, 1,    &
838                             comm2d, status, ierr )
839
840!
841!--       Update the tail ids for the transferred particles
842          nn = number_of_tails
843          DO  n = number_of_particles+1, number_of_particles+trnp_count_recv
844             IF ( particles(n)%tail_id /= 0 )  THEN
845                nn = nn + 1
846                particles(n)%tail_id = nn
847             ENDIF
848          ENDDO
849
850       ENDIF
851
852!       number_of_particles = number_of_particles + trnp_count_recv
853!       number_of_tails     = number_of_tails     + trnpt_count_recv
854
855!
856!--    Send back boundary, receive front boundary
857       CALL MPI_SENDRECV( trnp_count,      1, MPI_INTEGER, pnorth, 0, &
858                          trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
859                          comm2d, status, ierr )
860
861       ALLOCATE(rvsp(MAX(1,trsp_count_recv)))
862
863       CALL MPI_SENDRECV( trnp(1)%radius, trnp_count, mpi_particle_type,      &
864                          pnorth, 1, rvsp(1)%radius,                          &
865                          trsp_count_recv, mpi_particle_type, psouth, 1,   &
866                          comm2d, status, ierr )
867
868       IF ( trsp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvsp(1:trsp_count_recv))
869
870       DEALLOCATE(rvsp)
871
872       IF ( use_particle_tails )  THEN
873
874          CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
875                             trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
876                             comm2d, status, ierr )
877
878          IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
879          THEN
880             IF ( netcdf_data_format < 3 )  THEN
881                message_string = 'maximum_number_of_tails ' //   &
882                                 'needs to be increased ' //     &
883                                 '&but this is not allowed wi'// &
884                                 'th NetCDF output switched on'
885                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
886             ELSE
887                CALL lpm_extend_tail_array( trspt_count_recv )
888             ENDIF
889          ENDIF
890
891          CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,      &
892                             pnorth, 1,                                        &
893                             particle_tail_coordinates(1,1,number_of_tails+1), &
894                             trspt_count_recv*tlength, MPI_REAL, psouth, 1,    &
895                             comm2d, status, ierr )
896!
897!--       Update the tail ids for the transferred particles
898          nn = number_of_tails
899          DO  n = number_of_particles+1, number_of_particles+trsp_count_recv
900             IF ( particles(n)%tail_id /= 0 )  THEN
901                nn = nn + 1
902                particles(n)%tail_id = nn
903             ENDIF
904          ENDDO
905
906       ENDIF
907
908       number_of_particles = number_of_particles + trsp_count_recv
909       number_of_tails     = number_of_tails     + trspt_count_recv
910
911       IF ( use_particle_tails )  THEN
912          DEALLOCATE( trspt, trnpt )
913       ENDIF
914       DEALLOCATE( trsp, trnp )
915
916    ENDIF
917
918#else
919
920!
921!-- Apply boundary conditions
922    DO  n = 1, number_of_particles
923
924       nn = particles(n)%tail_id
925
926       IF ( particles(n)%x < -0.5_wp * dx )  THEN
927
928          IF ( ibc_par_lr == 0 )  THEN
929!
930!--          Cyclic boundary. Relevant coordinate has to be changed.
931             particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
932             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
933                i = particles(n)%tailpoints
934                particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
935                                             particle_tail_coordinates(1:i,1,nn)
936             ENDIF
937          ELSEIF ( ibc_par_lr == 1 )  THEN
938!
939!--          Particle absorption
940             particles(n)%particle_mask = .FALSE.
941             deleted_particles = deleted_particles + 1
942             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
943                tail_mask(nn) = .FALSE.
944                deleted_tails = deleted_tails + 1
945             ENDIF
946          ELSEIF ( ibc_par_lr == 2 )  THEN
947!
948!--          Particle reflection
949             particles(n)%x       = -dx - particles(n)%x
950             particles(n)%speed_x = -particles(n)%speed_x
951          ENDIF
952
953       ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx )  THEN
954
955          IF ( ibc_par_lr == 0 )  THEN
956!
957!--          Cyclic boundary. Relevant coordinate has to be changed.
958             particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
959             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
960                i = particles(n)%tailpoints
961                particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
962                                             particle_tail_coordinates(1:i,1,nn)
963             ENDIF
964          ELSEIF ( ibc_par_lr == 1 )  THEN
965!
966!--          Particle absorption
967             particles(n)%particle_mask = .FALSE.
968             deleted_particles = deleted_particles + 1
969             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
970                tail_mask(nn) = .FALSE.
971                deleted_tails = deleted_tails + 1
972             ENDIF
973          ELSEIF ( ibc_par_lr == 2 )  THEN
974!
975!--          Particle reflection
976             particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
977             particles(n)%speed_x = -particles(n)%speed_x
978          ENDIF
979
980       ENDIF
981
982       IF ( particles(n)%y < -0.5_wp * dy )  THEN
983
984          IF ( ibc_par_ns == 0 )  THEN
985!
986!--          Cyclic boundary. Relevant coordinate has to be changed.
987             particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
988             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
989                i = particles(n)%tailpoints
990                particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
991                                             particle_tail_coordinates(1:i,2,nn)
992             ENDIF
993          ELSEIF ( ibc_par_ns == 1 )  THEN
994!
995!--          Particle absorption
996             particles(n)%particle_mask = .FALSE.
997             deleted_particles = deleted_particles + 1
998             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
999                tail_mask(nn) = .FALSE.
1000                deleted_tails = deleted_tails + 1
1001             ENDIF
1002          ELSEIF ( ibc_par_ns == 2 )  THEN
1003!
1004!--          Particle reflection
1005             particles(n)%y       = -dy - particles(n)%y
1006             particles(n)%speed_y = -particles(n)%speed_y
1007          ENDIF
1008
1009       ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy )  THEN
1010
1011          IF ( ibc_par_ns == 0 )  THEN
1012!
1013!--          Cyclic boundary. Relevant coordinate has to be changed.
1014             particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
1015             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
1016                i = particles(n)%tailpoints
1017                particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
1018                                             particle_tail_coordinates(1:i,2,nn)
1019             ENDIF
1020          ELSEIF ( ibc_par_ns == 1 )  THEN
1021!
1022!--          Particle absorption
1023             particles(n)%particle_mask = .FALSE.
1024             deleted_particles = deleted_particles + 1
1025             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
1026                tail_mask(nn) = .FALSE.
1027                deleted_tails = deleted_tails + 1
1028             ENDIF
1029          ELSEIF ( ibc_par_ns == 2 )  THEN
1030!
1031!--          Particle reflection
1032             particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
1033             particles(n)%speed_y = -particles(n)%speed_y
1034          ENDIF
1035
1036       ENDIF
1037    ENDDO
1038
1039#endif
1040
1041!
1042!-- Accumulate the number of particles transferred between the subdomains
1043#if defined( __parallel )
1044    trlp_count_sum      = trlp_count_sum      + trlp_count
1045    trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
1046    trrp_count_sum      = trrp_count_sum      + trrp_count
1047    trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
1048    trsp_count_sum      = trsp_count_sum      + trsp_count
1049    trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
1050    trnp_count_sum      = trnp_count_sum      + trnp_count
1051    trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
1052#endif
1053
1054    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'stop' )
1055
1056 END SUBROUTINE lpm_exchange_horiz
1057
1058!------------------------------------------------------------------------------!
1059! Description:
1060! ------------
1061!> If a particle moves from one processor to another, this subroutine moves
1062!> the corresponding elements from the particle arrays of the old grid cells
1063!> to the particle arrays of the new grid cells.
1064!------------------------------------------------------------------------------!
1065 SUBROUTINE Add_particles_to_gridcell (particle_array)
1066
1067    IMPLICIT NONE
1068
1069    INTEGER(iwp)        ::  ip        !<
1070    INTEGER(iwp)        ::  jp        !<
1071    INTEGER(iwp)        ::  kp        !<
1072    INTEGER(iwp)        ::  n         !<
1073    INTEGER(iwp)        ::  pindex    !<
1074
1075    LOGICAL             ::  pack_done !<
1076
1077    TYPE(particle_type), DIMENSION(:), INTENT(IN)       :: particle_array
1078
1079    pack_done     = .FALSE.
1080
1081    nr_move_north = 0
1082    nr_move_south = 0
1083
1084    DO n = 1, SIZE(particle_array)
1085       ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx
1086       jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy
1087       kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
1088
1089       IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
1090            .AND.  kp >= nzb+1  .AND.  kp <= nzt)  THEN ! particle stays on processor
1091          number_of_particles = prt_count(kp,jp,ip)
1092          particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
1093
1094          pindex = prt_count(kp,jp,ip)+1
1095          IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
1096             IF ( pack_done )  THEN
1097                CALL realloc_particles_array (ip,jp,kp)
1098             ELSE
1099                CALL lpm_pack_arrays
1100                prt_count(kp,jp,ip) = number_of_particles
1101                pindex = prt_count(kp,jp,ip)+1
1102                IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
1103                   CALL realloc_particles_array (ip,jp,kp)
1104                ENDIF
1105                pack_done = .TRUE.
1106             ENDIF
1107          ENDIF
1108          grid_particles(kp,jp,ip)%particles(pindex) = particle_array(n)
1109          prt_count(kp,jp,ip) = pindex
1110       ELSE
1111          IF ( jp == nys - 1 )  THEN
1112             nr_move_south = nr_move_south+1
1113             move_also_south(nr_move_south) = particle_array(n)
1114             IF ( jp == -1 )  THEN
1115                move_also_south(nr_move_south)%y =                             &
1116                   move_also_south(nr_move_south)%y + ( ny + 1 ) * dy
1117                move_also_south(nr_move_south)%origin_y =                      &
1118                   move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
1119             ENDIF
1120          ELSEIF ( jp == nyn+1 )  THEN
1121             nr_move_north = nr_move_north+1
1122             move_also_north(nr_move_north) = particle_array(n)
1123             IF ( jp == ny+1 )  THEN
1124                move_also_north(nr_move_north)%y =                             &
1125                   move_also_north(nr_move_north)%y - ( ny + 1 ) * dy
1126                move_also_north(nr_move_north)%origin_y =                      &
1127                   move_also_north(nr_move_north)%origin_y - ( ny + 1 ) * dy
1128             ENDIF
1129          ELSE
1130             WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn
1131          ENDIF
1132       ENDIF
1133    ENDDO
1134
1135    RETURN
1136
1137 END SUBROUTINE Add_particles_to_gridcell
1138
1139
1140
1141
1142!------------------------------------------------------------------------------!
1143! Description:
1144! ------------
1145!> If a particle moves from one grid cell to another (on the current
1146!> processor!), this subroutine moves the corresponding element from the
1147!> particle array of the old grid cell to the particle array of the new grid
1148!> cell.
1149!------------------------------------------------------------------------------!
1150 SUBROUTINE lpm_move_particle
1151
1152    IMPLICIT NONE
1153
1154    INTEGER(iwp)        ::  i           !<
1155    INTEGER(iwp)        ::  ip          !<
1156    INTEGER(iwp)        ::  j           !<
1157    INTEGER(iwp)        ::  jp          !<
1158    INTEGER(iwp)        ::  k           !<
1159    INTEGER(iwp)        ::  kp          !<
1160    INTEGER(iwp)        ::  n           !<
1161    INTEGER(iwp)        ::  np_old_cell !<
1162    INTEGER(iwp)        ::  n_start     !<
1163    INTEGER(iwp)        ::  pindex      !<
1164
1165    LOGICAL             ::  pack_done   !<
1166
1167    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_old_cell !<
1168
1169    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
1170
1171    DO  ip = nxl, nxr
1172       DO  jp = nys, nyn
1173          DO  kp = nzb+1, nzt
1174
1175             np_old_cell = prt_count(kp,jp,ip)
1176             IF ( np_old_cell <= 0 )  CYCLE
1177             particles_old_cell => grid_particles(kp,jp,ip)%particles(1:np_old_cell)
1178             n_start = -1
1179             
1180             DO  n = 1, np_old_cell
1181                i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
1182                j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
1183                k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
1184!
1185!--             Check, if particle has moved to another grid cell.
1186                IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
1187!
1188!--                The particle has moved to another grid cell. Now check, if
1189!--                particle stays on the same processor.
1190                   IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.      &
1191                        j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
1192!
1193!--                   If the particle stays on the same processor, the particle
1194!--                   will be added to the particle array of the new processor.
1195                      number_of_particles = prt_count(k,j,i)
1196                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1197
1198                      pindex = prt_count(k,j,i)+1
1199                      IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )  &
1200                      THEN
1201                         n_start = n
1202                         EXIT
1203                      ENDIF
1204
1205                      grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
1206                      prt_count(k,j,i) = pindex
1207
1208                      particles_old_cell(n)%particle_mask = .FALSE.
1209                   ENDIF
1210                ENDIF
1211             ENDDO
1212
1213             IF ( n_start >= 0 )  THEN
1214                pack_done = .FALSE.
1215                DO  n = n_start, np_old_cell
1216                   i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
1217                   j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
1218                   k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
1219                   IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
1220!
1221!--                   Particle is in different box
1222                      IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.   &
1223                           j <= nyn  .AND.  k >= nzb+1  .AND.  k <= nzt)  THEN
1224!
1225!--                      Particle stays on processor
1226                         number_of_particles = prt_count(k,j,i)
1227                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1228
1229                         pindex = prt_count(k,j,i)+1
1230                         IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) &
1231                         THEN
1232                            IF ( pack_done )  THEN
1233                               CALL realloc_particles_array(i,j,k)
1234                               pindex = prt_count(k,j,i)+1
1235                            ELSE
1236                               CALL lpm_pack_arrays
1237                               prt_count(k,j,i) = number_of_particles
1238
1239                               pindex = prt_count(k,j,i)+1
1240!
1241!--                            If number of particles in the new grid box
1242!--                            exceeds its allocated memory, the particle array
1243!--                            will be reallocated
1244                               IF ( pindex > SIZE(grid_particles(k,j,i)%particles) )  THEN
1245                                  CALL realloc_particles_array(i,j,k)
1246                                  pindex = prt_count(k,j,i)+1
1247                               ENDIF
1248
1249                               pack_done = .TRUE.
1250                            ENDIF
1251                         ENDIF
1252
1253                         grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
1254                         prt_count(k,j,i) = pindex
1255
1256                         particles_old_cell(n)%particle_mask = .FALSE.
1257                      ENDIF
1258                   ENDIF
1259                ENDDO
1260             ENDIF
1261          ENDDO
1262       ENDDO
1263    ENDDO
1264
1265    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'stop' )
1266
1267    RETURN
1268
1269 END SUBROUTINE lpm_move_particle
1270
1271!------------------------------------------------------------------------------!
1272! Description:
1273! ------------
1274!> @todo Missing subroutine description.
1275!------------------------------------------------------------------------------!
1276 SUBROUTINE realloc_particles_array (i,j,k,size_in)
1277
1278    IMPLICIT NONE
1279
1280    INTEGER(iwp), INTENT(in)                       ::  i              !<
1281    INTEGER(iwp), INTENT(in)                       ::  j              !<
1282    INTEGER(iwp), INTENT(in)                       ::  k              !<
1283    INTEGER(iwp), INTENT(in), optional             ::  size_in        !<
1284
1285    INTEGER(iwp)                                   :: old_size        !<
1286    INTEGER(iwp)                                   :: new_size        !<
1287    TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<
1288    TYPE(particle_type), DIMENSION(500)            :: tmp_particles_s !<
1289
1290
1291    old_size = SIZE(grid_particles(k,j,i)%particles)
1292
1293    IF ( PRESENT(size_in) )   THEN
1294       new_size = size_in
1295    ELSE
1296       new_size = old_size * ( 1.0 + alloc_factor / 100.0 )
1297    ENDIF
1298
1299    new_size = MAX( new_size, min_nr_particle )
1300
1301    IF ( old_size <= 500 )  THEN
1302
1303       tmp_particles_s(1:old_size) = grid_particles(k,j,i)%particles(1:old_size)
1304
1305       DEALLOCATE(grid_particles(k,j,i)%particles)
1306       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1307
1308       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_s(1:old_size)
1309       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1310
1311    ELSE
1312
1313       ALLOCATE(tmp_particles_d(new_size))
1314       tmp_particles_d(1:old_size) = grid_particles(k,j,i)%particles
1315
1316       DEALLOCATE(grid_particles(k,j,i)%particles)
1317       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1318
1319       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_d(1:old_size)
1320       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1321
1322       DEALLOCATE(tmp_particles_d)
1323
1324    ENDIF
1325    particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1326
1327    RETURN
1328 END SUBROUTINE realloc_particles_array
1329
1330END MODULE lpm_exchange_horiz_mod
Note: See TracBrowser for help on using the repository browser.