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

Last change on this file since 1093 was 1037, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 36.7 KB
Line 
1 SUBROUTINE lpm_exchange_horiz
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_exchange_horiz.f90 1037 2012-10-22 14:10:22Z raasch $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 851 2012-03-15 14:32:58Z raasch
32! Bugfix: resetting of particle_mask and tail mask moved from end of this
33! routine to lpm
34!
35! 849 2012-03-15 10:35:09Z raasch
36! initial revision (former part of advec_particles)
37!
38!
39! Description:
40! ------------
41! Exchange of particles (and tails) between the subdomains.
42!------------------------------------------------------------------------------!
43
44    USE control_parameters
45    USE cpulog
46    USE grid_variables
47    USE indices
48    USE interfaces
49    USE particle_attributes
50    USE pegrid
51
52    IMPLICIT NONE
53
54    INTEGER ::  i, j, n, nn, tlength, &
55                trlp_count, trlp_count_recv, trlpt_count, trlpt_count_recv, &
56                trnp_count, trnp_count_recv, trnpt_count, trnpt_count_recv, &
57                trrp_count, trrp_count_recv, trrpt_count, trrpt_count_recv, &
58                trsp_count, trsp_count_recv, trspt_count, trspt_count_recv
59
60    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  trlpt, trnpt, trrpt, trspt
61
62    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp, trnp, trrp, trsp
63
64
65#if defined( __parallel )
66
67!
68!-- Exchange between subdomains.
69!-- As soon as one particle has moved beyond the boundary of the domain, it
70!-- is included in the relevant transfer arrays and marked for subsequent
71!-- deletion on this PE.
72!-- First sweep for crossings in x direction. Find out first the number of
73!-- particles to be transferred and allocate temporary arrays needed to store
74!-- them.
75!-- For a one-dimensional decomposition along y, no transfer is necessary,
76!-- because the particle remains on the PE, but the particle coordinate has to
77!-- be adjusted.
78    trlp_count  = 0
79    trlpt_count = 0
80    trrp_count  = 0
81    trrpt_count = 0
82
83    trlp_count_recv   = 0
84    trlpt_count_recv  = 0
85    trrp_count_recv   = 0
86    trrpt_count_recv  = 0
87
88    IF ( pdims(1) /= 1 )  THEN
89!
90!--    First calculate the storage necessary for sending and receiving the data
91       DO  n = 1, number_of_particles
92          i = ( particles(n)%x + 0.5 * dx ) * ddx
93!
94!--       Above calculation does not work for indices less than zero
95          IF ( particles(n)%x < -0.5 * dx )  i = -1
96
97          IF ( i < nxl )  THEN
98             trlp_count = trlp_count + 1
99             IF ( particles(n)%tail_id /= 0 )  trlpt_count = trlpt_count + 1
100          ELSEIF ( i > nxr )  THEN
101             trrp_count = trrp_count + 1
102             IF ( particles(n)%tail_id /= 0 )  trrpt_count = trrpt_count + 1
103          ENDIF
104       ENDDO
105
106       IF ( trlp_count  == 0 )  trlp_count  = 1
107       IF ( trlpt_count == 0 )  trlpt_count = 1
108       IF ( trrp_count  == 0 )  trrp_count  = 1
109       IF ( trrpt_count == 0 )  trrpt_count = 1
110
111       ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
112
113       trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
114                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
115                             0.0, 0, 0, 0, 0 )
116       trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
117                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
118                             0.0, 0, 0, 0, 0 )
119
120       IF ( use_particle_tails )  THEN
121          ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
122                    trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
123          tlength = maximum_number_of_tailpoints * 5
124       ENDIF
125
126       trlp_count  = 0
127       trlpt_count = 0
128       trrp_count  = 0
129       trrpt_count = 0
130
131    ENDIF
132
133    DO  n = 1, number_of_particles
134
135       nn = particles(n)%tail_id
136
137       i = ( particles(n)%x + 0.5 * dx ) * ddx
138!
139!--    Above calculation does not work for indices less than zero
140       IF ( particles(n)%x < - 0.5 * dx )  i = -1
141
142       IF ( i <  nxl )  THEN
143          IF ( i < 0 )  THEN
144!
145!--          Apply boundary condition along x
146             IF ( ibc_par_lr == 0 )  THEN
147!
148!--             Cyclic condition
149                IF ( pdims(1) == 1 )  THEN
150                   particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
151                   particles(n)%origin_x = ( nx + 1 ) * dx + &
152                                           particles(n)%origin_x
153                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
154                      i  = particles(n)%tailpoints
155                      particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
156                                        + particle_tail_coordinates(1:i,1,nn)
157                   ENDIF
158                ELSE
159                   trlp_count = trlp_count + 1
160                   trlp(trlp_count)   = particles(n)
161                   trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
162                   trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
163                                               ( nx + 1 ) * dx
164                   particle_mask(n)  = .FALSE.
165                   deleted_particles = deleted_particles + 1
166
167                   IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.0E-12 ) THEN
168                      trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10
169!++ why is 1 subtracted in next statement???
170                      trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x - 1
171                   ENDIF
172
173                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
174                      trlpt_count = trlpt_count + 1
175                      trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
176                      trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
177                                               trlpt(:,1,trlpt_count)
178                      tail_mask(nn) = .FALSE.
179                      deleted_tails = deleted_tails + 1
180                   ENDIF
181                ENDIF
182
183             ELSEIF ( ibc_par_lr == 1 )  THEN
184!
185!--             Particle absorption
186                particle_mask(n) = .FALSE.
187                deleted_particles = deleted_particles + 1
188                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
189                   tail_mask(nn) = .FALSE.
190                   deleted_tails = deleted_tails + 1
191                ENDIF
192
193             ELSEIF ( ibc_par_lr == 2 )  THEN
194!
195!--             Particle reflection
196                particles(n)%x       = -particles(n)%x
197                particles(n)%speed_x = -particles(n)%speed_x
198
199             ENDIF
200          ELSE
201!
202!--          Store particle data in the transfer array, which will be send
203!--          to the neighbouring PE
204             trlp_count = trlp_count + 1
205             trlp(trlp_count) = particles(n)
206             particle_mask(n) = .FALSE.
207             deleted_particles = deleted_particles + 1
208
209             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
210                trlpt_count = trlpt_count + 1
211                trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
212                tail_mask(nn) = .FALSE.
213                deleted_tails = deleted_tails + 1
214             ENDIF
215         ENDIF
216
217       ELSEIF ( i > nxr )  THEN
218          IF ( i > nx )  THEN
219!
220!--          Apply boundary condition along x
221             IF ( ibc_par_lr == 0 )  THEN
222!
223!--             Cyclic condition
224                IF ( pdims(1) == 1 )  THEN
225                   particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
226                   particles(n)%origin_x = particles(n)%origin_x - &
227                                           ( nx + 1 ) * dx
228                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
229                      i = particles(n)%tailpoints
230                      particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
231                                           + particle_tail_coordinates(1:i,1,nn)
232                   ENDIF
233                ELSE
234                   trrp_count = trrp_count + 1
235                   trrp(trrp_count) = particles(n)
236                   trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
237                   trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
238                                               ( nx + 1 ) * dx
239                   particle_mask(n) = .FALSE.
240                   deleted_particles = deleted_particles + 1
241
242                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
243                      trrpt_count = trrpt_count + 1
244                      trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
245                      trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
246                                               ( nx + 1 ) * dx
247                      tail_mask(nn) = .FALSE.
248                      deleted_tails = deleted_tails + 1
249                   ENDIF
250                ENDIF
251
252             ELSEIF ( ibc_par_lr == 1 )  THEN
253!
254!--             Particle absorption
255                particle_mask(n) = .FALSE.
256                deleted_particles = deleted_particles + 1
257                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
258                   tail_mask(nn) = .FALSE.
259                   deleted_tails = deleted_tails + 1
260                ENDIF
261
262             ELSEIF ( ibc_par_lr == 2 )  THEN
263!
264!--             Particle reflection
265                particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
266                particles(n)%speed_x = -particles(n)%speed_x
267
268             ENDIF
269          ELSE
270!
271!--          Store particle data in the transfer array, which will be send
272!--          to the neighbouring PE
273             trrp_count = trrp_count + 1
274             trrp(trrp_count) = particles(n)
275             particle_mask(n) = .FALSE.
276             deleted_particles = deleted_particles + 1
277
278             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
279                trrpt_count = trrpt_count + 1
280                trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
281                tail_mask(nn) = .FALSE.
282                deleted_tails = deleted_tails + 1
283             ENDIF
284          ENDIF
285
286       ENDIF
287    ENDDO
288
289!
290!-- Send left boundary, receive right boundary (but first exchange how many
291!-- and check, if particle storage must be extended)
292    IF ( pdims(1) /= 1 )  THEN
293
294       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'start' )
295       CALL MPI_SENDRECV( trlp_count,      1, MPI_INTEGER, pleft,  0, &
296                          trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
297                          comm2d, status, ierr )
298
299       IF ( number_of_particles + trrp_count_recv > &
300            maximum_number_of_particles )           &
301       THEN
302          IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
303              message_string = 'maximum_number_of_particles ' //    &
304                               'needs to be increased ' //          &
305                               '&but this is not allowed with ' //  &
306                               'netcdf-data_format < 3'
307             CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 )
308          ELSE
309             CALL lpm_extend_particle_array( trrp_count_recv )
310          ENDIF
311       ENDIF
312
313       CALL MPI_SENDRECV( trlp(1)%age, trlp_count, mpi_particle_type,     &
314                          pleft, 1, particles(number_of_particles+1)%age, &
315                          trrp_count_recv, mpi_particle_type, pright, 1,  &
316                          comm2d, status, ierr )
317
318       IF ( use_particle_tails )  THEN
319
320          CALL MPI_SENDRECV( trlpt_count,      1, MPI_INTEGER, pleft,  0, &
321                             trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
322                             comm2d, status, ierr )
323
324          IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
325          THEN
326             IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
327                message_string = 'maximum_number_of_tails ' //   &
328                                 'needs to be increased ' //     &
329                                 '&but this is not allowed wi'// &
330                                 'th netcdf_data_format < 3'
331                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
332             ELSE
333                CALL lpm_extend_tail_array( trrpt_count_recv )
334             ENDIF
335          ENDIF
336
337          CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,      &
338                             pleft, 1,                                         &
339                             particle_tail_coordinates(1,1,number_of_tails+1), &
340                             trrpt_count_recv*tlength, MPI_REAL, pright, 1,    &
341                             comm2d, status, ierr )
342!
343!--       Update the tail ids for the transferred particles
344          nn = number_of_tails
345          DO  n = number_of_particles+1, number_of_particles+trrp_count_recv
346             IF ( particles(n)%tail_id /= 0 )  THEN
347                nn = nn + 1
348                particles(n)%tail_id = nn
349             ENDIF
350          ENDDO
351
352       ENDIF
353
354       number_of_particles = number_of_particles + trrp_count_recv
355       number_of_tails     = number_of_tails     + trrpt_count_recv
356
357!
358!--    Send right boundary, receive left boundary
359       CALL MPI_SENDRECV( trrp_count,      1, MPI_INTEGER, pright, 0, &
360                          trlp_count_recv, 1, MPI_INTEGER, pleft,  0, &
361                          comm2d, status, ierr )
362
363       IF ( number_of_particles + trlp_count_recv > &
364            maximum_number_of_particles )           &
365       THEN
366          IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
367             message_string = 'maximum_number_of_particles ' //  &
368                              'needs to be increased ' //        &
369                              '&but this is not allowed with '// &
370                              'netcdf_data_format < 3'
371             CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 )
372          ELSE
373             CALL lpm_extend_particle_array( trlp_count_recv )
374          ENDIF
375       ENDIF
376
377       CALL MPI_SENDRECV( trrp(1)%age, trrp_count, mpi_particle_type,      &
378                          pright, 1, particles(number_of_particles+1)%age, &
379                          trlp_count_recv, mpi_particle_type, pleft, 1,    &
380                          comm2d, status, ierr )
381
382       IF ( use_particle_tails )  THEN
383
384          CALL MPI_SENDRECV( trrpt_count,      1, MPI_INTEGER, pright, 0, &
385                             trlpt_count_recv, 1, MPI_INTEGER, pleft,  0, &
386                             comm2d, status, ierr )
387
388          IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
389          THEN
390             IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
391                message_string = 'maximum_number_of_tails ' //   &
392                                 'needs to be increased ' //     &
393                                 '&but this is not allowed wi'// &
394                                 'th netcdf_data_format < 3'
395                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
396             ELSE
397                CALL lpm_extend_tail_array( trlpt_count_recv )
398             ENDIF
399          ENDIF
400
401          CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,      &
402                             pright, 1,                                        &
403                             particle_tail_coordinates(1,1,number_of_tails+1), &
404                             trlpt_count_recv*tlength, MPI_REAL, pleft, 1,     &
405                             comm2d, status, ierr )
406!
407!--       Update the tail ids for the transferred particles
408          nn = number_of_tails
409          DO  n = number_of_particles+1, number_of_particles+trlp_count_recv
410             IF ( particles(n)%tail_id /= 0 )  THEN
411                nn = nn + 1
412                particles(n)%tail_id = nn
413             ENDIF
414          ENDDO
415
416       ENDIF
417
418       number_of_particles = number_of_particles + trlp_count_recv
419       number_of_tails     = number_of_tails     + trlpt_count_recv
420
421       IF ( use_particle_tails )  THEN
422          DEALLOCATE( trlpt, trrpt )
423       ENDIF
424       DEALLOCATE( trlp, trrp ) 
425
426       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'pause' )
427
428    ENDIF
429
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 x, no transfer is necessary,
438!-- because the particle remains on the PE.
439    trsp_count  = 0
440    trspt_count = 0
441    trnp_count  = 0
442    trnpt_count = 0
443
444    trsp_count_recv   = 0
445    trspt_count_recv  = 0
446    trnp_count_recv   = 0
447    trnpt_count_recv  = 0
448
449    IF ( pdims(2) /= 1 )  THEN
450!
451!--    First calculate the storage necessary for sending and receiving the
452!--    data
453       DO  n = 1, number_of_particles
454          IF ( particle_mask(n) )  THEN
455             j = ( particles(n)%y + 0.5 * dy ) * ddy
456!
457!--          Above calculation does not work for indices less than zero
458             IF ( particles(n)%y < -0.5 * dy )  j = -1
459
460             IF ( j < nys )  THEN
461                trsp_count = trsp_count + 1
462                IF ( particles(n)%tail_id /= 0 )  trspt_count = trspt_count+1
463             ELSEIF ( j > nyn )  THEN
464                trnp_count = trnp_count + 1
465                IF ( particles(n)%tail_id /= 0 )  trnpt_count = trnpt_count+1
466             ENDIF
467          ENDIF
468       ENDDO
469
470       IF ( trsp_count  == 0 )  trsp_count  = 1
471       IF ( trspt_count == 0 )  trspt_count = 1
472       IF ( trnp_count  == 0 )  trnp_count  = 1
473       IF ( trnpt_count == 0 )  trnpt_count = 1
474
475       ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
476
477       trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
478                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
479                             0.0, 0, 0, 0, 0 )
480       trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
481                             0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
482                             0.0, 0, 0, 0, 0 )
483
484       IF ( use_particle_tails )  THEN
485          ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
486                    trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
487          tlength = maximum_number_of_tailpoints * 5
488       ENDIF
489
490       trsp_count  = 0
491       trspt_count = 0
492       trnp_count  = 0
493       trnpt_count = 0
494
495    ENDIF
496
497    DO  n = 1, number_of_particles
498
499       nn = particles(n)%tail_id
500!
501!--    Only those particles that have not been marked as 'deleted' may be
502!--    moved.
503       IF ( particle_mask(n) )  THEN
504          j = ( particles(n)%y + 0.5 * dy ) * ddy
505!
506!--       Above calculation does not work for indices less than zero
507          IF ( particles(n)%y < -0.5 * dy )  j = -1
508
509          IF ( j < nys )  THEN
510             IF ( j < 0 )  THEN
511!
512!--             Apply boundary condition along y
513                IF ( ibc_par_ns == 0 )  THEN
514!
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 + &
519                                              particles(n)%origin_y
520                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
521                         i = particles(n)%tailpoints
522                         particle_tail_coordinates(1:i,2,nn) = ( ny+1 ) * dy&
523                                        + particle_tail_coordinates(1:i,2,nn)
524                      ENDIF
525                   ELSE
526                      trsp_count = trsp_count + 1
527                      trsp(trsp_count) = particles(n)
528                      trsp(trsp_count)%y = ( ny + 1 ) * dy + &
529                                           trsp(trsp_count)%y
530                      trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
531                                                  + ( ny + 1 ) * dy
532                      particle_mask(n) = .FALSE.
533                      deleted_particles = deleted_particles + 1
534
535                      IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.0E-12 ) THEN
536                         trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10
537!++ why is 1 subtracted in next statement???
538                         trsp(trsp_count)%origin_y =                        &
539                                                   trsp(trsp_count)%origin_y - 1
540                      ENDIF
541
542                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
543                         trspt_count = trspt_count + 1
544                         trspt(:,:,trspt_count) = &
545                                               particle_tail_coordinates(:,:,nn)
546                         trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
547                                                  trspt(:,2,trspt_count)
548                         tail_mask(nn) = .FALSE.
549                         deleted_tails = deleted_tails + 1
550                      ENDIF
551                   ENDIF
552
553                ELSEIF ( ibc_par_ns == 1 )  THEN
554!
555!--                Particle absorption
556                   particle_mask(n) = .FALSE.
557                   deleted_particles = deleted_particles + 1
558                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
559                      tail_mask(nn) = .FALSE.
560                      deleted_tails = deleted_tails + 1
561                   ENDIF
562
563                ELSEIF ( ibc_par_ns == 2 )  THEN
564!
565!--                Particle reflection
566                   particles(n)%y       = -particles(n)%y
567                   particles(n)%speed_y = -particles(n)%speed_y
568
569                ENDIF
570             ELSE
571!
572!--             Store particle data in the transfer array, which will be send
573!--             to the neighbouring PE
574                trsp_count = trsp_count + 1
575                trsp(trsp_count) = particles(n)
576                particle_mask(n) = .FALSE.
577                deleted_particles = deleted_particles + 1
578
579                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
580                   trspt_count = trspt_count + 1
581                   trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
582                   tail_mask(nn) = .FALSE.
583                   deleted_tails = deleted_tails + 1
584                ENDIF
585             ENDIF
586
587          ELSEIF ( j > nyn )  THEN
588             IF ( j > ny )  THEN
589!
590!--             Apply boundary condition along x
591                IF ( ibc_par_ns == 0 )  THEN
592!
593!--                Cyclic condition
594                   IF ( pdims(2) == 1 )  THEN
595                      particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
596                      particles(n)%origin_y = particles(n)%origin_y - &
597                                              ( ny + 1 ) * dy
598                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
599                         i = particles(n)%tailpoints
600                         particle_tail_coordinates(1:i,2,nn) = - (ny+1) * dy &
601                                           + particle_tail_coordinates(1:i,2,nn)
602                      ENDIF
603                   ELSE
604                      trnp_count = trnp_count + 1
605                      trnp(trnp_count) = particles(n)
606                      trnp(trnp_count)%y = trnp(trnp_count)%y - &
607                                           ( ny + 1 ) * dy
608                      trnp(trnp_count)%origin_y = trnp(trnp_count)%origin_y &
609                                                  - ( ny + 1 ) * dy
610                      particle_mask(n) = .FALSE.
611                      deleted_particles = deleted_particles + 1
612
613                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
614                         trnpt_count = trnpt_count + 1
615                         trnpt(:,:,trnpt_count) = &
616                                               particle_tail_coordinates(:,:,nn)
617                         trnpt(:,2,trnpt_count) = trnpt(:,2,trnpt_count) - &
618                                                  ( ny + 1 ) * dy
619                         tail_mask(nn) = .FALSE.
620                         deleted_tails = deleted_tails + 1
621                      ENDIF
622                   ENDIF
623
624                ELSEIF ( ibc_par_ns == 1 )  THEN
625!
626!--                Particle absorption
627                   particle_mask(n) = .FALSE.
628                   deleted_particles = deleted_particles + 1
629                   IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
630                      tail_mask(nn) = .FALSE.
631                      deleted_tails = deleted_tails + 1
632                   ENDIF
633
634                ELSEIF ( ibc_par_ns == 2 )  THEN
635!
636!--                Particle reflection
637                   particles(n)%y       = 2 * ( ny * dy ) - particles(n)%y
638                   particles(n)%speed_y = -particles(n)%speed_y
639
640                ENDIF
641             ELSE
642!
643!--             Store particle data in the transfer array, which will be send
644!--             to the neighbouring PE
645                trnp_count = trnp_count + 1
646                trnp(trnp_count) = particles(n)
647                particle_mask(n) = .FALSE.
648                deleted_particles = deleted_particles + 1
649
650                IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
651                   trnpt_count = trnpt_count + 1
652                   trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
653                   tail_mask(nn) = .FALSE.
654                   deleted_tails = deleted_tails + 1
655                ENDIF
656             ENDIF
657
658          ENDIF
659       ENDIF
660    ENDDO
661
662!
663!-- Send front boundary, receive back boundary (but first exchange how many
664!-- and check, if particle storage must be extended)
665    IF ( pdims(2) /= 1 )  THEN
666
667       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'continue' )
668       CALL MPI_SENDRECV( trsp_count,      1, MPI_INTEGER, psouth, 0, &
669                          trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
670                          comm2d, status, ierr )
671
672       IF ( number_of_particles + trnp_count_recv > &
673            maximum_number_of_particles )           &
674       THEN
675          IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
676             message_string = 'maximum_number_of_particles ' //  &
677                              'needs to be increased ' //        &
678                              '&but this is not allowed with '// &
679                              'netcdf_data_format < 3'
680             CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 ) 
681          ELSE
682             CALL lpm_extend_particle_array( trnp_count_recv )
683          ENDIF
684       ENDIF
685
686       CALL MPI_SENDRECV( trsp(1)%age, trsp_count, mpi_particle_type,      &
687                          psouth, 1, particles(number_of_particles+1)%age, &
688                          trnp_count_recv, mpi_particle_type, pnorth, 1,   &
689                          comm2d, status, ierr )
690
691       IF ( use_particle_tails )  THEN
692
693          CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
694                             trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
695                             comm2d, status, ierr )
696
697          IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
698          THEN
699             IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
700                message_string = 'maximum_number_of_tails ' //    &
701                                 'needs to be increased ' //      &
702                                 '&but this is not allowed wi' // &
703                                 'th netcdf_data_format < 3'
704                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 ) 
705             ELSE
706                CALL lpm_extend_tail_array( trnpt_count_recv )
707             ENDIF
708          ENDIF
709
710          CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,      &
711                             psouth, 1,                                        &
712                             particle_tail_coordinates(1,1,number_of_tails+1), &
713                             trnpt_count_recv*tlength, MPI_REAL, pnorth, 1,    &
714                             comm2d, status, ierr )
715
716!
717!--       Update the tail ids for the transferred particles
718          nn = number_of_tails
719          DO  n = number_of_particles+1, number_of_particles+trnp_count_recv
720             IF ( particles(n)%tail_id /= 0 )  THEN
721                nn = nn + 1
722                particles(n)%tail_id = nn
723             ENDIF
724          ENDDO
725
726       ENDIF
727
728       number_of_particles = number_of_particles + trnp_count_recv
729       number_of_tails     = number_of_tails     + trnpt_count_recv
730
731!
732!--    Send back boundary, receive front boundary
733       CALL MPI_SENDRECV( trnp_count,      1, MPI_INTEGER, pnorth, 0, &
734                          trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
735                          comm2d, status, ierr )
736
737       IF ( number_of_particles + trsp_count_recv > &
738            maximum_number_of_particles )           &
739       THEN
740          IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
741             message_string = 'maximum_number_of_particles ' //   &
742                              'needs to be increased ' //         &
743                              '&but this is not allowed with ' // &
744                              'netcdf_data_format < 3'
745            CALL message( 'lpm_exch_horiz', 'PA0146', 2, 2, -1, 6, 1 ) 
746          ELSE
747             CALL lpm_extend_particle_array( trsp_count_recv )
748          ENDIF
749       ENDIF
750
751       CALL MPI_SENDRECV( trnp(1)%age, trnp_count, mpi_particle_type,      &
752                          pnorth, 1, particles(number_of_particles+1)%age, &
753                          trsp_count_recv, mpi_particle_type, psouth, 1,   &
754                          comm2d, status, ierr )
755
756       IF ( use_particle_tails )  THEN
757
758          CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
759                             trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
760                             comm2d, status, ierr )
761
762          IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
763          THEN
764             IF ( netcdf_output  .AND.  netcdf_data_format < 3 )  THEN
765                message_string = 'maximum_number_of_tails ' //   &
766                                 'needs to be increased ' //     &
767                                 '&but this is not allowed wi'// &
768                                 'th NetCDF output switched on'
769                CALL message( 'lpm_exch_horiz', 'PA0147', 2, 2, -1, 6, 1 )
770             ELSE
771                CALL lpm_extend_tail_array( trspt_count_recv )
772             ENDIF
773          ENDIF
774
775          CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,      &
776                             pnorth, 1,                                        &
777                             particle_tail_coordinates(1,1,number_of_tails+1), &
778                             trspt_count_recv*tlength, MPI_REAL, psouth, 1,    &
779                             comm2d, status, ierr )
780!
781!--       Update the tail ids for the transferred particles
782          nn = number_of_tails
783          DO  n = number_of_particles+1, number_of_particles+trsp_count_recv
784             IF ( particles(n)%tail_id /= 0 )  THEN
785                nn = nn + 1
786                particles(n)%tail_id = nn
787             ENDIF
788          ENDDO
789
790       ENDIF
791
792       number_of_particles = number_of_particles + trsp_count_recv
793       number_of_tails     = number_of_tails     + trspt_count_recv
794
795       IF ( use_particle_tails )  THEN
796          DEALLOCATE( trspt, trnpt )
797       ENDIF
798       DEALLOCATE( trsp, trnp )
799
800       CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'stop' )
801
802    ENDIF
803
804#else
805
806!
807!-- Apply boundary conditions
808    DO  n = 1, number_of_particles
809
810       nn = particles(n)%tail_id
811
812       IF ( particles(n)%x < -0.5 * dx )  THEN
813
814          IF ( ibc_par_lr == 0 )  THEN
815!
816!--          Cyclic boundary. Relevant coordinate has to be changed.
817             particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
818             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
819                i = particles(n)%tailpoints
820                particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
821                                             particle_tail_coordinates(1:i,1,nn)
822             ENDIF
823          ELSEIF ( ibc_par_lr == 1 )  THEN
824!
825!--          Particle absorption
826             particle_mask(n) = .FALSE.
827             deleted_particles = deleted_particles + 1
828             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
829                tail_mask(nn) = .FALSE.
830                deleted_tails = deleted_tails + 1
831             ENDIF
832          ELSEIF ( ibc_par_lr == 2 )  THEN
833!
834!--          Particle reflection
835             particles(n)%x       = -dx - particles(n)%x
836             particles(n)%speed_x = -particles(n)%speed_x
837          ENDIF
838
839       ELSEIF ( particles(n)%x >= ( nx + 0.5 ) * dx )  THEN
840
841          IF ( ibc_par_lr == 0 )  THEN
842!
843!--          Cyclic boundary. Relevant coordinate has to be changed.
844             particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
845             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
846                i = particles(n)%tailpoints
847                particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
848                                             particle_tail_coordinates(1:i,1,nn)
849             ENDIF
850          ELSEIF ( ibc_par_lr == 1 )  THEN
851!
852!--          Particle absorption
853             particle_mask(n) = .FALSE.
854             deleted_particles = deleted_particles + 1
855             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
856                tail_mask(nn) = .FALSE.
857                deleted_tails = deleted_tails + 1
858             ENDIF
859          ELSEIF ( ibc_par_lr == 2 )  THEN
860!
861!--          Particle reflection
862             particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
863             particles(n)%speed_x = -particles(n)%speed_x
864          ENDIF
865
866       ENDIF
867
868       IF ( particles(n)%y < -0.5 * dy )  THEN
869
870          IF ( ibc_par_ns == 0 )  THEN
871!
872!--          Cyclic boundary. Relevant coordinate has to be changed.
873             particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
874             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
875                i = particles(n)%tailpoints
876                particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
877                                             particle_tail_coordinates(1:i,2,nn)
878             ENDIF
879          ELSEIF ( ibc_par_ns == 1 )  THEN
880!
881!--          Particle absorption
882             particle_mask(n) = .FALSE.
883             deleted_particles = deleted_particles + 1
884             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
885                tail_mask(nn) = .FALSE.
886                deleted_tails = deleted_tails + 1
887             ENDIF
888          ELSEIF ( ibc_par_ns == 2 )  THEN
889!
890!--          Particle reflection
891             particles(n)%y       = -dy - particles(n)%y
892             particles(n)%speed_y = -particles(n)%speed_y
893          ENDIF
894
895       ELSEIF ( particles(n)%y >= ( ny + 0.5 ) * dy )  THEN
896
897          IF ( ibc_par_ns == 0 )  THEN
898!
899!--          Cyclic boundary. Relevant coordinate has to be changed.
900             particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
901             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
902                i = particles(n)%tailpoints
903                particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
904                                             particle_tail_coordinates(1:i,2,nn)
905             ENDIF
906          ELSEIF ( ibc_par_ns == 1 )  THEN
907!
908!--          Particle absorption
909             particle_mask(n) = .FALSE.
910             deleted_particles = deleted_particles + 1
911             IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
912                tail_mask(nn) = .FALSE.
913                deleted_tails = deleted_tails + 1
914             ENDIF
915          ELSEIF ( ibc_par_ns == 2 )  THEN
916!
917!--          Particle reflection
918             particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
919             particles(n)%speed_y = -particles(n)%speed_y
920          ENDIF
921
922       ENDIF
923    ENDDO
924
925#endif
926
927!
928!-- Accumulate the number of particles transferred between the subdomains
929#if defined( __parallel )
930    trlp_count_sum      = trlp_count_sum      + trlp_count
931    trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
932    trrp_count_sum      = trrp_count_sum      + trrp_count
933    trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
934    trsp_count_sum      = trsp_count_sum      + trsp_count
935    trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
936    trnp_count_sum      = trnp_count_sum      + trnp_count
937    trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
938#endif
939
940
941 END SUBROUTINE lpm_exchange_horiz
Note: See TracBrowser for help on using the repository browser.