source: palm/trunk/SOURCE/lpm_exchange_horiz_mod.f90 @ 1851

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

last commit documented

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