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

Last change on this file since 1826 was 1823, checked in by hoffmann, 8 years ago

last commit documented

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