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

Last change on this file since 1552 was 1360, checked in by hoffmann, 10 years ago

last commit documented

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