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

Last change on this file since 1818 was 1818, checked in by maronga, 6 years ago

last commit documented / copyright update

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