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

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

revised renaming of modules

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