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

Last change on this file since 1874 was 1874, checked in by maronga, 5 years ago

last commit documented

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