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

Last change on this file since 1929 was 1929, checked in by suehring, 5 years ago

several bugfixes in particle model and serial mode

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