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

Last change on this file since 1822 was 1822, checked in by hoffmann, 7 years ago

changes in LPM and bulk cloud microphysics

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