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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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