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

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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