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

Last change on this file since 1685 was 1685, checked in by raasch, 6 years ago

bugfix concerning vertical index calculation for particles in case of ocean runs

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