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

Last change on this file since 1932 was 1930, checked in by suehring, 8 years ago

last commit documented

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