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

Last change on this file since 3436 was 3065, checked in by Giersch, 6 years ago

New vertical stretching procedure has been introduced

  • Property svn:keywords set to Id
File size: 51.4 KB
RevLine 
[1873]1!> @file lpm_exchange_horiz.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1036]9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[849]20! Current revisions:
21! ------------------
[1930]22!
[3049]23!
[1930]24! Former revisions:
25! -----------------
26! $Id: lpm_exchange_horiz.f90 3065 2018-06-12 07:03:02Z gronemeier $
[3065]27! dz was replaced by dz(1) to allow for right vertical stretching
28!
29! 3049 2018-05-29 13:52:36Z Giersch
[3049]30! Error messages revised
31!
32! 3046 2018-05-29 08:02:15Z Giersch
[2809]33! Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE.
34! Preprocessor directive(__gfortran) for c_sizeof removed.
35!
36! 2801 2018-02-14 16:01:55Z thiele
[2801]37! Introduce particle transfer in nested models.
38!
39! 2718 2018-01-02 08:49:38Z maronga
[2716]40! Corrected "Former revisions" section
41!
42! 2710 2017-12-19 12:50:11Z suehring
43! Changes from last commit documented
44!
45! 2709 2017-12-19 12:49:40Z suehring
[2710]46! Bugfix in CFL check.
47! Further bugfix, non-allocated array used in case of 1D decomposition.
[2716]48!
49! 2696 2017-12-14 17:12:51Z kanani
50! Change in file header (GPL part)
51!
52! 2632 2017-11-20 15:35:52Z schwenkel
[2632]53! Bugfix in write statement
[2716]54!
[2632]55! 2630 2017-11-20 12:58:20Z schwenkel
[2628]56! Enabled particle advection with grid stretching. Furthermore, the CFL-
57! criterion is checked for every particle at every time step.
58!
59! 2606 2017-11-10 10:36:31Z schwenkel
[2606]60! Changed particle box locations: center of particle box now coincides
61! with scalar grid point of same index.
62! lpm_move_particle now moves inner particles that would leave the core
63! to the respective boundary gridbox. Afterwards they get transferred by
64! lpm_exchange_horiz. This allows particles to move more than one gridbox
65! independent of domain composition. 
66! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
67! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
68! lpm_sort -> lpm_sort_timeloop_done
69!
70! 2330 2017-08-03 14:26:02Z schwenkel
[2330]71! Bugfix: Also for gfortran compilable, function c_sizeof is not used.
72!
73! 2305 2017-07-06 11:18:47Z hoffmann
[2305]74! Improved calculation of particle IDs.
75!
76! 2245 2017-06-02 14:37:10Z schwenkel
[2245]77! Bugfix in Add_particles_to_gridcell:
78! Apply boundary conditions also in y-direction
79!
80! 2151 2017-02-15 10:42:16Z schwenkel
81! Bugfix in lpm_move_particle that prevents
82! false production of additional particles
[1930]83!
[2001]84! 2000 2016-08-20 18:09:15Z knoop
85! Forced header and separation lines into 80 columns
86!
[1937]87! 1936 2016-06-13 13:37:44Z suehring
88! Deallocation of unused memory
89!
[1930]90! 1929 2016-06-09 16:25:25Z suehring
[1929]91! Bugfixes:
92! - reallocation of new particles
93!   ( did not work for small number of min_nr_particle )
94! - dynamical reallocation of north-south exchange arrays ( particles got lost )
95! - north-south exchange ( nr_move_north, nr_move_south were overwritten by zero )
96! - horizontal particle boundary conditions in serial mode
97!
98! Remove unused variables
99! Descriptions in variable declaration blocks added
[1321]100!
[1874]101! 1873 2016-04-18 14:50:06Z maronga
102! Module renamed (removed _mod)
103!
104!
[1851]105! 1850 2016-04-08 13:29:27Z maronga
106! Module renamed
107!
108!
[1823]109! 1822 2016-04-07 07:49:42Z hoffmann
110! Tails removed. Unused variables removed.
111!
[1784]112! 1783 2016-03-06 18:36:17Z raasch
113! new netcdf-module included
114!
[1692]115! 1691 2015-10-26 16:17:44Z maronga
116! Formatting corrections.
117!
[1686]118! 1685 2015-10-08 07:32:13Z raasch
119! bugfix concerning vertical index offset in case of ocean
120!
[1683]121! 1682 2015-10-07 23:56:08Z knoop
122! Code annotations made doxygen readable
123!
[1360]124! 1359 2014-04-11 17:15:14Z hoffmann
125! New particle structure integrated.
126! Kind definition added to all floating point numbers.
127!
[1329]128! 1327 2014-03-21 11:00:16Z raasch
129! -netcdf output queries
130!
[1321]131! 1320 2014-03-20 08:40:49Z raasch
[1320]132! ONLY-attribute added to USE-statements,
133! kind-parameters added to all INTEGER and REAL declaration statements,
134! kinds are defined in new module kinds,
135! comment fields (!:) to be used for variable explanations added to
136! all variable declaration statements
[849]137!
[1319]138! 1318 2014-03-17 13:35:16Z raasch
139! module interfaces removed
140!
[1037]141! 1036 2012-10-22 13:43:42Z raasch
142! code put under GPL (PALM 3.9)
143!
[852]144! 851 2012-03-15 14:32:58Z raasch
145! Bugfix: resetting of particle_mask and tail mask moved from end of this
146! routine to lpm
147!
[850]148! 849 2012-03-15 10:35:09Z raasch
149! initial revision (former part of advec_particles)
[849]150!
[850]151!
[849]152! Description:
153! ------------
[1822]154! Exchange of particles between the subdomains.
[849]155!------------------------------------------------------------------------------!
[1682]156 MODULE lpm_exchange_horiz_mod
157 
[2305]158    USE, INTRINSIC ::  ISO_C_BINDING
[2628]159   
160    USE arrays_3d,                                                             &
161       ONLY:  zw
162   
[1320]163    USE control_parameters,                                                    &
[1783]164        ONLY:  dz, message_string, simulated_time
[1320]165
166    USE cpulog,                                                                &
167        ONLY:  cpu_log, log_point_s
168
169    USE grid_variables,                                                        &
170        ONLY:  ddx, ddy, dx, dy
171
172    USE indices,                                                               &
[1359]173        ONLY:  nx, nxl, nxr, ny, nyn, nys, nzb, nzt
[1320]174
175    USE kinds
176
[2606]177    USE lpm_pack_and_sort_mod,                                                   &
178        ONLY:  lpm_pack
[1359]179
[1783]180    USE netcdf_interface,                                                      &
181        ONLY:  netcdf_data_format
182
[1320]183    USE particle_attributes,                                                   &
[1822]184        ONLY:  alloc_factor, deleted_particles, grid_particles,                &
185               ibc_par_lr, ibc_par_ns, min_nr_particle,                        &
[2305]186               number_of_particles,                                            &
[1822]187               offset_ocean_nzt, offset_ocean_nzt_m1, particles,               &
188               particle_type, prt_count, trlp_count_sum,                       &
[1359]189               trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum,       &
190               trrp_count_sum, trrp_count_recv_sum, trsp_count_sum,            &
[1822]191               trsp_count_recv_sum, zero_particle
[1320]192
[849]193    USE pegrid
194
195    IMPLICIT NONE
196
[1682]197    INTEGER(iwp), PARAMETER ::  NR_2_direction_move = 10000 !<
198    INTEGER(iwp)            ::  nr_move_north               !<
199    INTEGER(iwp)            ::  nr_move_south               !<
[1359]200
[1929]201    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_north
202    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  move_also_south
[1359]203
204    SAVE
205
206    PRIVATE
[1936]207    PUBLIC lpm_exchange_horiz, lpm_move_particle, realloc_particles_array,     &
208           dealloc_particles_array
[1359]209
210    INTERFACE lpm_exchange_horiz
211       MODULE PROCEDURE lpm_exchange_horiz
212    END INTERFACE lpm_exchange_horiz
213
214    INTERFACE lpm_move_particle
215       MODULE PROCEDURE lpm_move_particle
216    END INTERFACE lpm_move_particle
217
218    INTERFACE realloc_particles_array
219       MODULE PROCEDURE realloc_particles_array
220    END INTERFACE realloc_particles_array
221
[1936]222    INTERFACE dealloc_particles_array
223       MODULE PROCEDURE dealloc_particles_array
224    END INTERFACE dealloc_particles_array
[1359]225CONTAINS
226
[1682]227!------------------------------------------------------------------------------!
228! Description:
229! ------------
230!> Exchange between subdomains.
231!> As soon as one particle has moved beyond the boundary of the domain, it
232!> is included in the relevant transfer arrays and marked for subsequent
233!> deletion on this PE.
234!> First sweep for crossings in x direction. Find out first the number of
235!> particles to be transferred and allocate temporary arrays needed to store
236!> them.
237!> For a one-dimensional decomposition along y, no transfer is necessary,
238!> because the particle remains on the PE, but the particle coordinate has to
239!> be adjusted.
240!------------------------------------------------------------------------------!
[1359]241 SUBROUTINE lpm_exchange_horiz
242
243    IMPLICIT NONE
244
[1929]245    INTEGER(iwp) ::  i                 !< grid index (x) of particle positition
246    INTEGER(iwp) ::  ip                !< index variable along x
247    INTEGER(iwp) ::  j                 !< grid index (y) of particle positition
248    INTEGER(iwp) ::  jp                !< index variable along y
249    INTEGER(iwp) ::  kp                !< index variable along z
250    INTEGER(iwp) ::  n                 !< particle index variable
[2305]251    INTEGER(iwp) ::  par_size          !< Particle size in bytes
[1929]252    INTEGER(iwp) ::  trlp_count        !< number of particles send to left PE
253    INTEGER(iwp) ::  trlp_count_recv   !< number of particles receive from right PE
254    INTEGER(iwp) ::  trnp_count        !< number of particles send to north PE
255    INTEGER(iwp) ::  trnp_count_recv   !< number of particles receive from south PE
256    INTEGER(iwp) ::  trrp_count        !< number of particles send to right PE
257    INTEGER(iwp) ::  trrp_count_recv   !< number of particles receive from left PE
258    INTEGER(iwp) ::  trsp_count        !< number of particles send to south PE
259    INTEGER(iwp) ::  trsp_count_recv   !< number of particles receive from north PE
[849]260
[1929]261    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvlp  !< particles received from right PE
262    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvnp  !< particles received from south PE
263    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvrp  !< particles received from left PE
264    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  rvsp  !< particles received from north PE
265    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trlp  !< particles send to left PE
266    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trnp  !< particles send to north PE
267    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trrp  !< particles send to right PE
268    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  trsp  !< particles send to south PE
[849]269
[1359]270    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
271
[849]272#if defined( __parallel )
273
[1822]274!
275!-- Exchange between subdomains.
276!-- As soon as one particle has moved beyond the boundary of the domain, it
277!-- is included in the relevant transfer arrays and marked for subsequent
278!-- deletion on this PE.
279!-- First sweep for crossings in x direction. Find out first the number of
280!-- particles to be transferred and allocate temporary arrays needed to store
281!-- them.
282!-- For a one-dimensional decomposition along y, no transfer is necessary,
283!-- because the particle remains on the PE, but the particle coordinate has to
284!-- be adjusted.
[849]285    trlp_count  = 0
286    trrp_count  = 0
287
288    trlp_count_recv   = 0
289    trrp_count_recv   = 0
290
291    IF ( pdims(1) /= 1 )  THEN
292!
[1359]293!--    First calculate the storage necessary for sending and receiving the data.
294!--    Compute only first (nxl) and last (nxr) loop iterration.
295       DO  ip = nxl, nxr, nxr - nxl
296          DO  jp = nys, nyn
297             DO  kp = nzb+1, nzt
[849]298
[1359]299                number_of_particles = prt_count(kp,jp,ip)
300                IF ( number_of_particles <= 0 )  CYCLE
301                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
302                DO  n = 1, number_of_particles
303                   IF ( particles(n)%particle_mask )  THEN
[2606]304                      i = particles(n)%x * ddx
[1929]305!
306!--                   Above calculation does not work for indices less than zero
[2606]307                      IF ( particles(n)%x < 0.0_wp)  i = -1
[1359]308
309                      IF ( i < nxl )  THEN
310                         trlp_count = trlp_count + 1
311                      ELSEIF ( i > nxr )  THEN
312                         trrp_count = trrp_count + 1
313                      ENDIF
314                   ENDIF
315                ENDDO
316
317             ENDDO
318          ENDDO
[849]319       ENDDO
320
321       IF ( trlp_count  == 0 )  trlp_count  = 1
322       IF ( trrp_count  == 0 )  trrp_count  = 1
323
324       ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
325
[1359]326       trlp = zero_particle
327       trrp = zero_particle
[849]328
329       trlp_count  = 0
330       trrp_count  = 0
331
332    ENDIF
[1359]333!
334!-- Compute only first (nxl) and last (nxr) loop iterration
[1929]335    DO  ip = nxl, nxr, nxr-nxl
[1359]336       DO  jp = nys, nyn
337          DO  kp = nzb+1, nzt
338             number_of_particles = prt_count(kp,jp,ip)
339             IF ( number_of_particles <= 0 ) CYCLE
340             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
341             DO  n = 1, number_of_particles
342!
343!--             Only those particles that have not been marked as 'deleted' may
344!--             be moved.
345                IF ( particles(n)%particle_mask )  THEN
[849]346
[2606]347                   i = particles(n)%x * ddx
[1929]348!
349!--                Above calculation does not work for indices less than zero
[2606]350                   IF ( particles(n)%x < 0.0_wp )  i = -1
[849]351
[1359]352                   IF ( i <  nxl )  THEN
353                      IF ( i < 0 )  THEN
[1929]354!
355!--                   Apply boundary condition along x
[1359]356                         IF ( ibc_par_lr == 0 )  THEN
[1929]357!
358!--                         Cyclic condition
[1359]359                            IF ( pdims(1) == 1 )  THEN
360                               particles(n)%x        = ( nx + 1 ) * dx + particles(n)%x
361                               particles(n)%origin_x = ( nx + 1 ) * dx + &
362                               particles(n)%origin_x
363                            ELSE
364                               trlp_count = trlp_count + 1
365                               trlp(trlp_count)   = particles(n)
366                               trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
367                               trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
368                               ( nx + 1 ) * dx
369                               particles(n)%particle_mask  = .FALSE.
370                               deleted_particles = deleted_particles + 1
[849]371
[2606]372                               IF ( trlp(trlp_count)%x >= (nx + 1)* dx - 1.0E-12_wp )  THEN
[1359]373                                  trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10_wp
374                                  !++ why is 1 subtracted in next statement???
375                                  trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x - 1
376                               ENDIF
[849]377
[1359]378                            ENDIF
[849]379
[1359]380                         ELSEIF ( ibc_par_lr == 1 )  THEN
[1929]381!
382!--                         Particle absorption
[1359]383                            particles(n)%particle_mask = .FALSE.
384                            deleted_particles = deleted_particles + 1
[849]385
[1359]386                         ELSEIF ( ibc_par_lr == 2 )  THEN
[1929]387!
388!--                         Particle reflection
[1359]389                            particles(n)%x       = -particles(n)%x
390                            particles(n)%speed_x = -particles(n)%speed_x
[849]391
[1359]392                         ENDIF
393                      ELSE
[1929]394!
395!--                      Store particle data in the transfer array, which will be
396!--                      send to the neighbouring PE
[1359]397                         trlp_count = trlp_count + 1
398                         trlp(trlp_count) = particles(n)
399                         particles(n)%particle_mask = .FALSE.
400                         deleted_particles = deleted_particles + 1
[849]401
[1359]402                      ENDIF
[849]403
[1359]404                   ELSEIF ( i > nxr )  THEN
405                      IF ( i > nx )  THEN
[1929]406!
407!--                      Apply boundary condition along x
[1359]408                         IF ( ibc_par_lr == 0 )  THEN
[1929]409!
410!--                         Cyclic condition
[1359]411                            IF ( pdims(1) == 1 )  THEN
412                               particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
413                               particles(n)%origin_x = particles(n)%origin_x - &
414                               ( nx + 1 ) * dx
415                            ELSE
416                               trrp_count = trrp_count + 1
417                               trrp(trrp_count) = particles(n)
418                               trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
419                               trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
420                               ( nx + 1 ) * dx
421                               particles(n)%particle_mask = .FALSE.
422                               deleted_particles = deleted_particles + 1
[849]423
[1359]424                            ENDIF
[849]425
[1359]426                         ELSEIF ( ibc_par_lr == 1 )  THEN
[1929]427!
428!--                         Particle absorption
[1359]429                            particles(n)%particle_mask = .FALSE.
430                            deleted_particles = deleted_particles + 1
[849]431
[1359]432                         ELSEIF ( ibc_par_lr == 2 )  THEN
[1929]433!
434!--                         Particle reflection
[1359]435                            particles(n)%x       = 2 * ( nx * dx ) - particles(n)%x
436                            particles(n)%speed_x = -particles(n)%speed_x
[849]437
[1359]438                         ENDIF
439                      ELSE
[1929]440!
441!--                      Store particle data in the transfer array, which will be send
442!--                      to the neighbouring PE
[1359]443                         trrp_count = trrp_count + 1
444                         trrp(trrp_count) = particles(n)
445                         particles(n)%particle_mask = .FALSE.
446                         deleted_particles = deleted_particles + 1
[849]447
[1359]448                      ENDIF
[849]449
[1359]450                   ENDIF
451                ENDIF
[1929]452
[1359]453             ENDDO
454          ENDDO
455       ENDDO
[849]456    ENDDO
457
458!
[2809]459!-- STORAGE_SIZE returns the storage size of argument A in bits. However , it
460!-- is needed in bytes. The function C_SIZEOF which produces this value directly
461!-- causes problems with gfortran. For this reason the use of C_SIZEOF is avoided
462    par_size = STORAGE_SIZE(trlp(1))/8
[2709]463
[2330]464
465!
[1929]466!-- Allocate arrays required for north-south exchange, as these
467!-- are used directly after particles are exchange along x-direction.
468    ALLOCATE( move_also_north(1:NR_2_direction_move) )
469    ALLOCATE( move_also_south(1:NR_2_direction_move) )
470
471    nr_move_north = 0
472    nr_move_south = 0
473!
[849]474!-- Send left boundary, receive right boundary (but first exchange how many
475!-- and check, if particle storage must be extended)
476    IF ( pdims(1) /= 1 )  THEN
477
478       CALL MPI_SENDRECV( trlp_count,      1, MPI_INTEGER, pleft,  0, &
479                          trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
480                          comm2d, status, ierr )
481
[1359]482       ALLOCATE(rvrp(MAX(1,trrp_count_recv)))
[2330]483
[2305]484       CALL MPI_SENDRECV( trlp, max(1,trlp_count)*par_size, MPI_BYTE,&
485                          pleft, 1, rvrp,                            &
486                          max(1,trrp_count_recv)*par_size, MPI_BYTE, pright, 1,&
[849]487                          comm2d, status, ierr )
488
[1359]489       IF ( trrp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvrp(1:trrp_count_recv))
490
491       DEALLOCATE(rvrp)
492
[849]493!
494!--    Send right boundary, receive left boundary
495       CALL MPI_SENDRECV( trrp_count,      1, MPI_INTEGER, pright, 0, &
496                          trlp_count_recv, 1, MPI_INTEGER, pleft,  0, &
497                          comm2d, status, ierr )
498
[1359]499       ALLOCATE(rvlp(MAX(1,trlp_count_recv)))
[2305]500!     
501!--    This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit
502!--    variables in structure particle_type (due to the calculation of par_size)
503       CALL MPI_SENDRECV( trrp, max(1,trrp_count)*par_size, MPI_BYTE,&
504                          pright, 1, rvlp,                           &
505                          max(1,trlp_count_recv)*par_size, MPI_BYTE, pleft, 1, &
[849]506                          comm2d, status, ierr )
507
[1359]508       IF ( trlp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvlp(1:trlp_count_recv))
509
[1822]510       DEALLOCATE( rvlp )
[1359]511       DEALLOCATE( trlp, trrp )
[849]512
513    ENDIF
514
515!
516!-- Check whether particles have crossed the boundaries in y direction. Note
517!-- that this case can also apply to particles that have just been received
518!-- from the adjacent right or left PE.
519!-- Find out first the number of particles to be transferred and allocate
520!-- temporary arrays needed to store them.
[1929]521!-- For a one-dimensional decomposition along y, no transfer is necessary,
[849]522!-- because the particle remains on the PE.
[1359]523    trsp_count  = nr_move_south
524    trnp_count  = nr_move_north
[849]525
526    trsp_count_recv   = 0
527    trnp_count_recv   = 0
528
529    IF ( pdims(2) /= 1 )  THEN
530!
531!--    First calculate the storage necessary for sending and receiving the
532!--    data
[1359]533       DO  ip = nxl, nxr
[1929]534          DO  jp = nys, nyn, nyn-nys    !compute only first (nys) and last (nyn) loop iterration
[1359]535             DO  kp = nzb+1, nzt
536                number_of_particles = prt_count(kp,jp,ip)
537                IF ( number_of_particles <= 0 )  CYCLE
538                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
539                DO  n = 1, number_of_particles
540                   IF ( particles(n)%particle_mask )  THEN
[2606]541                      j = particles(n)%y * ddy
[849]542!
[1359]543!--                   Above calculation does not work for indices less than zero
[2606]544                      IF ( particles(n)%y < 0.0_wp)  j = -1
[849]545
[1359]546                      IF ( j < nys )  THEN
547                         trsp_count = trsp_count + 1
548                      ELSEIF ( j > nyn )  THEN
549                         trnp_count = trnp_count + 1
550                      ENDIF
551                   ENDIF
552                ENDDO
553             ENDDO
554          ENDDO
[849]555       ENDDO
556
557       IF ( trsp_count  == 0 )  trsp_count  = 1
558       IF ( trnp_count  == 0 )  trnp_count  = 1
559
560       ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
561
[1359]562       trsp = zero_particle
563       trnp = zero_particle
[849]564
[1359]565       trsp_count  = nr_move_south
566       trnp_count  = nr_move_north
[1929]567       
[1359]568       trsp(1:nr_move_south) = move_also_south(1:nr_move_south)
569       trnp(1:nr_move_north) = move_also_north(1:nr_move_north)
570
[849]571    ENDIF
572
[1359]573    DO  ip = nxl, nxr
574       DO  jp = nys, nyn, nyn-nys ! compute only first (nys) and last (nyn) loop iterration
575          DO  kp = nzb+1, nzt
576             number_of_particles = prt_count(kp,jp,ip)
577             IF ( number_of_particles <= 0 )  CYCLE
578             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
579             DO  n = 1, number_of_particles
[849]580!
[1359]581!--             Only those particles that have not been marked as 'deleted' may
582!--             be moved.
583                IF ( particles(n)%particle_mask )  THEN
[1929]584
[2606]585                   j = particles(n)%y * ddy
[849]586!
[1359]587!--                Above calculation does not work for indices less than zero
[2606]588                   IF ( particles(n)%y < 0.0_wp )  j = -1
[849]589
[1359]590                   IF ( j < nys )  THEN
591                      IF ( j < 0 )  THEN
[849]592!
[1359]593!--                      Apply boundary condition along y
594                         IF ( ibc_par_ns == 0 )  THEN
[849]595!
[1359]596!--                         Cyclic condition
597                            IF ( pdims(2) == 1 )  THEN
598                               particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
599                               particles(n)%origin_y = ( ny + 1 ) * dy + &
[1929]600                                                     particles(n)%origin_y
[1359]601                            ELSE
[1929]602                               trsp_count         = trsp_count + 1
603                               trsp(trsp_count)   = particles(n)
[1359]604                               trsp(trsp_count)%y = ( ny + 1 ) * dy + &
[1929]605                                                 trsp(trsp_count)%y
[1359]606                               trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
[1929]607                                                + ( ny + 1 ) * dy
[1359]608                               particles(n)%particle_mask = .FALSE.
609                               deleted_particles = deleted_particles + 1
[849]610
[2606]611                               IF ( trsp(trsp_count)%y >= (ny+1)* dy - 1.0E-12_wp )  THEN
[1359]612                                  trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10_wp
613                                  !++ why is 1 subtracted in next statement???
614                                  trsp(trsp_count)%origin_y =                        &
[1929]615                                                  trsp(trsp_count)%origin_y - 1
[1359]616                               ENDIF
[849]617
[1359]618                            ENDIF
[849]619
[1359]620                         ELSEIF ( ibc_par_ns == 1 )  THEN
[849]621!
[1359]622!--                         Particle absorption
623                            particles(n)%particle_mask = .FALSE.
[1929]624                            deleted_particles          = deleted_particles + 1
[849]625
[1359]626                         ELSEIF ( ibc_par_ns == 2 )  THEN
[849]627!
[1359]628!--                         Particle reflection
629                            particles(n)%y       = -particles(n)%y
630                            particles(n)%speed_y = -particles(n)%speed_y
[849]631
[1359]632                         ENDIF
633                      ELSE
[849]634!
[1359]635!--                      Store particle data in the transfer array, which will
636!--                      be send to the neighbouring PE
637                         trsp_count = trsp_count + 1
638                         trsp(trsp_count) = particles(n)
639                         particles(n)%particle_mask = .FALSE.
640                         deleted_particles = deleted_particles + 1
[849]641
[1359]642                      ENDIF
[849]643
[1359]644                   ELSEIF ( j > nyn )  THEN
645                      IF ( j > ny )  THEN
[849]646!
[1929]647!--                       Apply boundary condition along y
[1359]648                         IF ( ibc_par_ns == 0 )  THEN
[849]649!
[1359]650!--                         Cyclic condition
651                            IF ( pdims(2) == 1 )  THEN
[1929]652                               particles(n)%y        = particles(n)%y - ( ny + 1 ) * dy
[1359]653                               particles(n)%origin_y =                         &
[1929]654                                          particles(n)%origin_y - ( ny + 1 ) * dy
[1359]655                            ELSE
[1929]656                               trnp_count         = trnp_count + 1
657                               trnp(trnp_count)   = particles(n)
[1359]658                               trnp(trnp_count)%y =                            &
[1929]659                                          trnp(trnp_count)%y - ( ny + 1 ) * dy
[1359]660                               trnp(trnp_count)%origin_y =                     &
[1929]661                                         trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
[1359]662                               particles(n)%particle_mask = .FALSE.
[1929]663                               deleted_particles          = deleted_particles + 1
[1359]664                            ENDIF
[849]665
[1359]666                         ELSEIF ( ibc_par_ns == 1 )  THEN
[849]667!
[1359]668!--                         Particle absorption
669                            particles(n)%particle_mask = .FALSE.
670                            deleted_particles = deleted_particles + 1
[849]671
[1359]672                         ELSEIF ( ibc_par_ns == 2 )  THEN
[849]673!
[1359]674!--                         Particle reflection
675                            particles(n)%y       = 2 * ( ny * dy ) - particles(n)%y
676                            particles(n)%speed_y = -particles(n)%speed_y
[849]677
[1359]678                         ENDIF
679                      ELSE
[849]680!
[1359]681!--                      Store particle data in the transfer array, which will
682!--                      be send to the neighbouring PE
683                         trnp_count = trnp_count + 1
684                         trnp(trnp_count) = particles(n)
685                         particles(n)%particle_mask = .FALSE.
686                         deleted_particles = deleted_particles + 1
[849]687
[1359]688                      ENDIF
689
690                   ENDIF
[849]691                ENDIF
[1359]692             ENDDO
693          ENDDO
694       ENDDO
[849]695    ENDDO
696
697!
698!-- Send front boundary, receive back boundary (but first exchange how many
699!-- and check, if particle storage must be extended)
700    IF ( pdims(2) /= 1 )  THEN
701
702       CALL MPI_SENDRECV( trsp_count,      1, MPI_INTEGER, psouth, 0, &
703                          trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
704                          comm2d, status, ierr )
705
[1359]706       ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
[2305]707!     
708!--    This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit
709!--    variables in structure particle_type (due to the calculation of par_size)
710       CALL MPI_SENDRECV( trsp, trsp_count*par_size, MPI_BYTE,      &
711                          psouth, 1, rvnp,                             &
712                          trnp_count_recv*par_size, MPI_BYTE, pnorth, 1,   &
[849]713                          comm2d, status, ierr )
714
[1359]715       IF ( trnp_count_recv  > 0 )  CALL Add_particles_to_gridcell(rvnp(1:trnp_count_recv))
716
717       DEALLOCATE(rvnp)
718
[849]719!
720!--    Send back boundary, receive front boundary
721       CALL MPI_SENDRECV( trnp_count,      1, MPI_INTEGER, pnorth, 0, &
722                          trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
723                          comm2d, status, ierr )
724
[1359]725       ALLOCATE(rvsp(MAX(1,trsp_count_recv)))
[2305]726!     
727!--    This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit
728!--    variables in structure particle_type (due to the calculation of par_size)
729       CALL MPI_SENDRECV( trnp, trnp_count*par_size, MPI_BYTE,      &
730                          pnorth, 1, rvsp,                          &
731                          trsp_count_recv*par_size, MPI_BYTE, psouth, 1,   &
[849]732                          comm2d, status, ierr )
733
[1359]734       IF ( trsp_count_recv > 0 )  CALL Add_particles_to_gridcell(rvsp(1:trsp_count_recv))
735
736       DEALLOCATE(rvsp)
737
[849]738       number_of_particles = number_of_particles + trsp_count_recv
739
740       DEALLOCATE( trsp, trnp )
741
742    ENDIF
743
[1929]744    DEALLOCATE( move_also_north )
745    DEALLOCATE( move_also_south )
746
[849]747#else
748
[1929]749    DO  ip = nxl, nxr, nxr-nxl
750       DO  jp = nys, nyn
751          DO  kp = nzb+1, nzt
752             number_of_particles = prt_count(kp,jp,ip)
753             IF ( number_of_particles <= 0 )  CYCLE
754             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
755             DO  n = 1, number_of_particles
[849]756!
[1929]757!--             Apply boundary conditions
[849]758
[2606]759                IF ( particles(n)%x < 0.0_wp )  THEN
[849]760
[1929]761                   IF ( ibc_par_lr == 0 )  THEN
[849]762!
[1929]763!--                   Cyclic boundary. Relevant coordinate has to be changed.
764                      particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
[2606]765                      particles(n)%origin_x = ( nx + 1 ) * dx + &
766                               particles(n)%origin_x
[1929]767                   ELSEIF ( ibc_par_lr == 1 )  THEN
[849]768!
[1929]769!--                   Particle absorption
770                      particles(n)%particle_mask = .FALSE.
771                      deleted_particles = deleted_particles + 1
[1822]772
[1929]773                   ELSEIF ( ibc_par_lr == 2 )  THEN
[849]774!
[1929]775!--                   Particle reflection
776                      particles(n)%x       = -dx - particles(n)%x
777                      particles(n)%speed_x = -particles(n)%speed_x
778                   ENDIF
[849]779
[2606]780                ELSEIF ( particles(n)%x >= ( nx + 1) * dx )  THEN
[849]781
[1929]782                   IF ( ibc_par_lr == 0 )  THEN
[849]783!
[1929]784!--                   Cyclic boundary. Relevant coordinate has to be changed.
785                      particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
[2606]786                      particles(n)%origin_x = particles(n)%origin_x - &
787                               ( nx + 1 ) * dx
[1822]788
[1929]789                   ELSEIF ( ibc_par_lr == 1 )  THEN
[849]790!
[1929]791!--                   Particle absorption
792                      particles(n)%particle_mask = .FALSE.
793                      deleted_particles = deleted_particles + 1
[1822]794
[1929]795                   ELSEIF ( ibc_par_lr == 2 )  THEN
[849]796!
[1929]797!--                   Particle reflection
798                      particles(n)%x       = ( nx + 1 ) * dx - particles(n)%x
799                      particles(n)%speed_x = -particles(n)%speed_x
800                   ENDIF
[849]801
[1929]802                ENDIF
803             ENDDO
804          ENDDO
805       ENDDO
806    ENDDO
[849]807
[1929]808    DO  ip = nxl, nxr
809       DO  jp = nys, nyn, nyn-nys
810          DO  kp = nzb+1, nzt
811             number_of_particles = prt_count(kp,jp,ip)
812             IF ( number_of_particles <= 0 )  CYCLE
813             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
814             DO  n = 1, number_of_particles
[849]815
[2606]816                IF ( particles(n)%y < 0.0_wp)  THEN
[1929]817
818                   IF ( ibc_par_ns == 0 )  THEN
[849]819!
[1929]820!--                   Cyclic boundary. Relevant coordinate has to be changed.
821                      particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
[2606]822                      particles(n)%origin_y = ( ny + 1 ) * dy + &
823                           particles(n)%origin_y
[1822]824
[1929]825                   ELSEIF ( ibc_par_ns == 1 )  THEN
[849]826!
[1929]827!--                   Particle absorption
828                      particles(n)%particle_mask = .FALSE.
829                      deleted_particles = deleted_particles + 1
[1822]830
[1929]831                   ELSEIF ( ibc_par_ns == 2 )  THEN
[849]832!
[1929]833!--                   Particle reflection
834                      particles(n)%y       = -dy - particles(n)%y
835                      particles(n)%speed_y = -particles(n)%speed_y
836                   ENDIF
[849]837
[2606]838                ELSEIF ( particles(n)%y >= ( ny + 1) * dy )  THEN
[849]839
[1929]840                   IF ( ibc_par_ns == 0 )  THEN
[849]841!
[1929]842!--                   Cyclic boundary. Relevant coordinate has to be changed.
843                      particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
[2606]844                      particles(n)%origin_y = particles(n)%origin_y - &
845                                ( ny + 1 ) * dy
[1822]846
[1929]847                   ELSEIF ( ibc_par_ns == 1 )  THEN
[849]848!
[1929]849!--                   Particle absorption
850                      particles(n)%particle_mask = .FALSE.
851                      deleted_particles = deleted_particles + 1
[1822]852
[1929]853                   ELSEIF ( ibc_par_ns == 2 )  THEN
[849]854!
[1929]855!--                   Particle reflection
856                      particles(n)%y       = ( ny + 1 ) * dy - particles(n)%y
857                      particles(n)%speed_y = -particles(n)%speed_y
858                   ENDIF
[849]859
[1929]860                ENDIF
861
862             ENDDO
863          ENDDO
864       ENDDO
[849]865    ENDDO
866#endif
867
868!
869!-- Accumulate the number of particles transferred between the subdomains
870#if defined( __parallel )
871    trlp_count_sum      = trlp_count_sum      + trlp_count
872    trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
873    trrp_count_sum      = trrp_count_sum      + trrp_count
874    trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
875    trsp_count_sum      = trsp_count_sum      + trsp_count
876    trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
877    trnp_count_sum      = trnp_count_sum      + trnp_count
878    trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
879#endif
880
[1359]881    CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'stop' )
[849]882
883 END SUBROUTINE lpm_exchange_horiz
[1359]884
[1682]885!------------------------------------------------------------------------------!
886! Description:
887! ------------
888!> If a particle moves from one processor to another, this subroutine moves
889!> the corresponding elements from the particle arrays of the old grid cells
890!> to the particle arrays of the new grid cells.
891!------------------------------------------------------------------------------!
[1359]892 SUBROUTINE Add_particles_to_gridcell (particle_array)
893
894    IMPLICIT NONE
895
[1929]896    INTEGER(iwp)        ::  ip        !< grid index (x) of particle
897    INTEGER(iwp)        ::  jp        !< grid index (x) of particle
898    INTEGER(iwp)        ::  kp        !< grid index (x) of particle
899    INTEGER(iwp)        ::  n         !< index variable of particle
900    INTEGER(iwp)        ::  pindex    !< dummy argument for new number of particles per grid box
[1359]901
[1682]902    LOGICAL             ::  pack_done !<
[1359]903
[1929]904    TYPE(particle_type), DIMENSION(:), INTENT(IN)  ::  particle_array !< new particles in a grid box
905    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  temp_ns        !< temporary particle array for reallocation
[1359]906
907    pack_done     = .FALSE.
908
[1929]909    DO n = 1, SIZE(particle_array)
[1359]910
[1929]911       IF ( .NOT. particle_array(n)%particle_mask )  CYCLE
912
[2606]913       ip = particle_array(n)%x * ddx
914       jp = particle_array(n)%y * ddy
[3065]915       kp = particle_array(n)%z / dz(1) + 1 + offset_ocean_nzt
[2628]916!
917!--    In case of grid stretching a particle might be above or below the 
918!--    previously calculated particle grid box (indices).     
919       DO WHILE( zw(kp) < particle_array(n)%z )
920          kp = kp + 1
921       ENDDO
[1359]922
[2628]923       DO WHILE( zw(kp-1) > particle_array(n)%z )
924          kp = kp - 1
925       ENDDO
926
[1359]927       IF ( ip >= nxl  .AND.  ip <= nxr  .AND.  jp >= nys  .AND.  jp <= nyn    &
928            .AND.  kp >= nzb+1  .AND.  kp <= nzt)  THEN ! particle stays on processor
929          number_of_particles = prt_count(kp,jp,ip)
930          particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
931
932          pindex = prt_count(kp,jp,ip)+1
933          IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
934             IF ( pack_done )  THEN
935                CALL realloc_particles_array (ip,jp,kp)
936             ELSE
[2606]937                CALL lpm_pack
[1359]938                prt_count(kp,jp,ip) = number_of_particles
939                pindex = prt_count(kp,jp,ip)+1
940                IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) )  THEN
941                   CALL realloc_particles_array (ip,jp,kp)
942                ENDIF
943                pack_done = .TRUE.
944             ENDIF
945          ENDIF
946          grid_particles(kp,jp,ip)%particles(pindex) = particle_array(n)
947          prt_count(kp,jp,ip) = pindex
948       ELSE
[1929]949          IF ( jp <= nys - 1 )  THEN
[1359]950             nr_move_south = nr_move_south+1
[1929]951!
952!--          Before particle information is swapped to exchange-array, check
953!--          if enough memory is allocated. If required, reallocate exchange
954!--          array.
955             IF ( nr_move_south > SIZE(move_also_south) )  THEN
956!
957!--             At first, allocate further temporary array to swap particle
958!--             information.
959                ALLOCATE( temp_ns(SIZE(move_also_south)+NR_2_direction_move) )
960                temp_ns(1:nr_move_south-1) = move_also_south(1:nr_move_south-1)
961                DEALLOCATE( move_also_south )
962                ALLOCATE( move_also_south(SIZE(temp_ns)) )
963                move_also_south(1:nr_move_south-1) = temp_ns(1:nr_move_south-1)
964                DEALLOCATE( temp_ns )
965
966             ENDIF
967
[1359]968             move_also_south(nr_move_south) = particle_array(n)
[1929]969
[1359]970             IF ( jp == -1 )  THEN
[2245]971!
972!--             Apply boundary condition along y
973                IF ( ibc_par_ns == 0 )  THEN
974                   move_also_south(nr_move_south)%y =                          &
975                      move_also_south(nr_move_south)%y + ( ny + 1 ) * dy
976                   move_also_south(nr_move_south)%origin_y =                   &
977                      move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
978                ELSEIF ( ibc_par_ns == 1 )  THEN
979!
980!--                Particle absorption
981                   move_also_south(nr_move_south)%particle_mask = .FALSE.
982                   deleted_particles = deleted_particles + 1
983
984                ELSEIF ( ibc_par_ns == 2 )  THEN
985!
986!--                Particle reflection
987                   move_also_south(nr_move_south)%y       =                    &
988                      -move_also_south(nr_move_south)%y
989                   move_also_south(nr_move_south)%speed_y =                    &
990                      -move_also_south(nr_move_south)%speed_y
991
992                ENDIF
[1359]993             ENDIF
[1929]994          ELSEIF ( jp >= nyn+1 )  THEN
[1359]995             nr_move_north = nr_move_north+1
[1929]996!
997!--          Before particle information is swapped to exchange-array, check
998!--          if enough memory is allocated. If required, reallocate exchange
999!--          array.
1000             IF ( nr_move_north > SIZE(move_also_north) )  THEN
1001!
1002!--             At first, allocate further temporary array to swap particle
1003!--             information.
1004                ALLOCATE( temp_ns(SIZE(move_also_north)+NR_2_direction_move) )
1005                temp_ns(1:nr_move_north-1) = move_also_south(1:nr_move_north-1)
1006                DEALLOCATE( move_also_north )
1007                ALLOCATE( move_also_north(SIZE(temp_ns)) )
1008                move_also_north(1:nr_move_north-1) = temp_ns(1:nr_move_north-1)
1009                DEALLOCATE( temp_ns )
1010
1011             ENDIF
1012
[1359]1013             move_also_north(nr_move_north) = particle_array(n)
1014             IF ( jp == ny+1 )  THEN
[2245]1015!
1016!--             Apply boundary condition along y
1017                IF ( ibc_par_ns == 0 )  THEN
1018
1019                   move_also_north(nr_move_north)%y =                          &
1020                      move_also_north(nr_move_north)%y - ( ny + 1 ) * dy
1021                   move_also_north(nr_move_north)%origin_y =                   &
1022                      move_also_north(nr_move_north)%origin_y - ( ny + 1 ) * dy
1023                ELSEIF ( ibc_par_ns == 1 )  THEN
1024!
1025!--                Particle absorption
1026                   move_also_north(nr_move_north)%particle_mask = .FALSE.
1027                   deleted_particles = deleted_particles + 1
1028
1029                ELSEIF ( ibc_par_ns == 2 )  THEN
1030!
1031!--                Particle reflection
1032                   move_also_north(nr_move_north)%y       =                    &
1033                      -move_also_north(nr_move_north)%y
1034                   move_also_north(nr_move_north)%speed_y =                    &
1035                      -move_also_north(nr_move_north)%speed_y
1036
1037                ENDIF
[1359]1038             ENDIF
1039          ELSE
1040             WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn
1041          ENDIF
1042       ENDIF
1043    ENDDO
1044
1045    RETURN
1046
1047 END SUBROUTINE Add_particles_to_gridcell
1048
1049
1050
1051
[1682]1052!------------------------------------------------------------------------------!
1053! Description:
1054! ------------
1055!> If a particle moves from one grid cell to another (on the current
1056!> processor!), this subroutine moves the corresponding element from the
1057!> particle array of the old grid cell to the particle array of the new grid
1058!> cell.
1059!------------------------------------------------------------------------------!
[1359]1060 SUBROUTINE lpm_move_particle
[2628]1061 
[1359]1062    IMPLICIT NONE
1063
[1929]1064    INTEGER(iwp)        ::  i           !< grid index (x) of particle position
1065    INTEGER(iwp)        ::  ip          !< index variable along x
1066    INTEGER(iwp)        ::  j           !< grid index (y) of particle position
1067    INTEGER(iwp)        ::  jp          !< index variable along y
1068    INTEGER(iwp)        ::  k           !< grid index (z) of particle position
1069    INTEGER(iwp)        ::  kp          !< index variable along z
1070    INTEGER(iwp)        ::  n           !< index variable for particle array
[2606]1071    INTEGER(iwp)        ::  np_before_move !< number of particles per grid box before moving
[1929]1072    INTEGER(iwp)        ::  pindex      !< dummy argument for number of new particle per grid box
[1359]1073
[2606]1074    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles_before_move !< particles before moving
[1359]1075
1076    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
[2628]1077    CALL lpm_check_cfl
[1359]1078    DO  ip = nxl, nxr
1079       DO  jp = nys, nyn
1080          DO  kp = nzb+1, nzt
1081
[2606]1082             np_before_move = prt_count(kp,jp,ip)
1083             IF ( np_before_move <= 0 )  CYCLE
1084             particles_before_move => grid_particles(kp,jp,ip)%particles(1:np_before_move)
[1359]1085             
[2606]1086             DO  n = 1, np_before_move
1087                i = particles_before_move(n)%x * ddx
1088                j = particles_before_move(n)%y * ddy
[2628]1089                k = kp
1090!
1091!--             Find correct vertical particle grid box (necessary in case of grid stretching)
1092!--             Due to the CFL limitations only the neighbouring grid boxes are considered.
1093                IF( zw(k)   < particles_before_move(n)%z ) k = k + 1
1094                IF( zw(k-1) > particles_before_move(n)%z ) k = k - 1 
1095               
[2606]1096!--             For lpm_exchange_horiz to work properly particles need to be moved to the outermost gridboxes
1097!--             of the respective processor. If the particle index is inside the processor the following lines
1098!--             will not change the index
1099                i = MIN ( i , nxr )
1100                i = MAX ( i , nxl )
1101                j = MIN ( j , nyn )
1102                j = MAX ( j , nys )
[2628]1103               
[2606]1104                k = MIN ( k , nzt )
1105                k = MAX ( k , nzb+1 )
[2628]1106               
[1359]1107!
1108!--             Check, if particle has moved to another grid cell.
1109                IF ( i /= ip  .OR.  j /= jp  .OR.  k /= kp )  THEN
[2606]1110!!
1111!--                If the particle stays on the same processor, the particle
1112!--                will be added to the particle array of the new processor.
1113                   number_of_particles = prt_count(k,j,i)
1114                   particles => grid_particles(k,j,i)%particles(1:number_of_particles)
[1359]1115
[2606]1116                   pindex = prt_count(k,j,i)+1
1117                   IF (  pindex > SIZE(grid_particles(k,j,i)%particles)  )     &
1118                   THEN
1119                      CALL realloc_particles_array(i,j,k)
1120                   ENDIF
[1359]1121
[2606]1122                   grid_particles(k,j,i)%particles(pindex) = particles_before_move(n)
1123                   prt_count(k,j,i) = pindex
[1359]1124
[2606]1125                   particles_before_move(n)%particle_mask = .FALSE.
[1359]1126                ENDIF
1127             ENDDO
1128
1129          ENDDO
1130       ENDDO
1131    ENDDO
1132
1133    CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'stop' )
1134
1135    RETURN
1136
1137 END SUBROUTINE lpm_move_particle
[2628]1138 
1139!------------------------------------------------------------------------------!
1140! Description:
1141! ------------
1142!> Check CFL-criterion for each particle. If one particle violated the
1143!> criterion the particle will be deleted and a warning message is given.
1144!------------------------------------------------------------------------------!
1145 SUBROUTINE lpm_check_cfl 
1146       
1147    IMPLICIT NONE
1148   
[2709]1149    INTEGER(iwp)  ::  i !< running index, x-direction
1150    INTEGER(iwp)  ::  j !< running index, y-direction
1151    INTEGER(iwp)  ::  k !< running index, z-direction
1152    INTEGER(iwp)  ::  n !< running index, number of particles
1153
[2628]1154    DO  i = nxl, nxr
1155       DO  j = nys, nyn
1156          DO  k = nzb+1, nzt
1157             number_of_particles = prt_count(k,j,i)
1158             IF ( number_of_particles <= 0 )  CYCLE
1159             particles => grid_particles(k,j,i)%particles(1:number_of_particles)         
1160             DO n = 1, number_of_particles
[2709]1161!
1162!--             Note, check for CFL does not work at first particle timestep
1163!--             when both, age and age_m are zero.
1164                IF ( particles(n)%age - particles(n)%age_m > 0.0_wp )  THEN 
1165                   IF(ABS(particles(n)%speed_x) >                              &
1166                      (dx/(particles(n)%age-particles(n)%age_m))  .OR.         &
1167                      ABS(particles(n)%speed_y) >                              & 
1168                      (dy/(particles(n)%age-particles(n)%age_m))  .OR.         &
1169                      ABS(particles(n)%speed_z) >                              &
1170                      ((zw(k)-zw(k-1))/(particles(n)%age-particles(n)%age_m))) &
1171                   THEN
1172                      WRITE( message_string, * )                               &
[3046]1173                      'Particle violated CFL-criterion: &particle with id ',   &
1174                      particles(n)%id, ' will be deleted!'   
[2801]1175                      CALL message( 'lpm_check_cfl', 'PA0475', 0, 1, -1, 6, 0 )
[2709]1176                      particles(n)%particle_mask= .FALSE.
1177                   ENDIF
[2628]1178                ENDIF
1179             ENDDO
1180          ENDDO
1181       ENDDO
1182    ENDDO   
[1359]1183
[2628]1184 END SUBROUTINE lpm_check_cfl
1185 
[1682]1186!------------------------------------------------------------------------------!
1187! Description:
1188! ------------
[1929]1189!> If the allocated memory for the particle array do not suffice to add arriving
1190!> particles from neighbour grid cells, this subrouting reallocates the
1191!> particle array to assure enough memory is available.
[1682]1192!------------------------------------------------------------------------------!
[1359]1193 SUBROUTINE realloc_particles_array (i,j,k,size_in)
1194
1195    IMPLICIT NONE
1196
[1682]1197    INTEGER(iwp), INTENT(in)                       ::  i              !<
1198    INTEGER(iwp), INTENT(in)                       ::  j              !<
1199    INTEGER(iwp), INTENT(in)                       ::  k              !<
[1822]1200    INTEGER(iwp), INTENT(in), OPTIONAL             ::  size_in        !<
[1359]1201
[1682]1202    INTEGER(iwp)                                   :: old_size        !<
1203    INTEGER(iwp)                                   :: new_size        !<
1204    TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<
1205    TYPE(particle_type), DIMENSION(500)            :: tmp_particles_s !<
[1359]1206
1207    old_size = SIZE(grid_particles(k,j,i)%particles)
1208
1209    IF ( PRESENT(size_in) )   THEN
1210       new_size = size_in
1211    ELSE
[1822]1212       new_size = old_size * ( 1.0_wp + alloc_factor / 100.0_wp )
[1359]1213    ENDIF
1214
[1929]1215    new_size = MAX( new_size, min_nr_particle, old_size + 1 )
[1359]1216
1217    IF ( old_size <= 500 )  THEN
1218
1219       tmp_particles_s(1:old_size) = grid_particles(k,j,i)%particles(1:old_size)
1220
1221       DEALLOCATE(grid_particles(k,j,i)%particles)
1222       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1223
1224       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_s(1:old_size)
1225       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1226
1227    ELSE
1228
1229       ALLOCATE(tmp_particles_d(new_size))
1230       tmp_particles_d(1:old_size) = grid_particles(k,j,i)%particles
1231
1232       DEALLOCATE(grid_particles(k,j,i)%particles)
1233       ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1234
1235       grid_particles(k,j,i)%particles(1:old_size)          = tmp_particles_d(1:old_size)
1236       grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
1237
1238       DEALLOCATE(tmp_particles_d)
1239
1240    ENDIF
1241    particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1242
1243    RETURN
1244 END SUBROUTINE realloc_particles_array
1245
[1936]1246
1247
1248
1249
1250 SUBROUTINE dealloc_particles_array
1251
1252    IMPLICIT NONE
1253
1254    INTEGER(iwp) ::  i
1255    INTEGER(iwp) ::  j
1256    INTEGER(iwp) ::  k
1257    INTEGER(iwp) :: old_size        !<
1258    INTEGER(iwp) :: new_size        !<
1259
1260    LOGICAL                                        :: dealloc 
1261
1262    TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<
1263    TYPE(particle_type), DIMENSION(500)            :: tmp_particles_s !<
1264
1265    DO  i = nxl, nxr
1266       DO  j = nys, nyn
1267          DO  k = nzb+1, nzt
1268!
1269!--          Determine number of active particles
1270             number_of_particles = prt_count(k,j,i)
1271!
1272!--          Determine allocated memory size
1273             old_size = SIZE( grid_particles(k,j,i)%particles )
1274!
1275!--          Check for large unused memory
1276             dealloc = ( ( number_of_particles < min_nr_particle .AND.         &
1277                           old_size            > min_nr_particle )  .OR.       &
1278                         ( number_of_particles > min_nr_particle .AND.         &
1279                           old_size - number_of_particles *                    &
1280                              ( 1.0_wp + 0.01_wp * alloc_factor ) > 0.0_wp ) )
1281                         
1282
1283             IF ( dealloc )  THEN
1284                IF ( number_of_particles < min_nr_particle )  THEN
1285                   new_size = min_nr_particle
1286                ELSE
1287                   new_size = INT( number_of_particles * ( 1.0_wp + 0.01_wp * alloc_factor ) )
1288                ENDIF
1289
1290                IF ( number_of_particles <= 500 )  THEN
1291
1292                   tmp_particles_s(1:number_of_particles) = grid_particles(k,j,i)%particles(1:number_of_particles)
1293
1294                   DEALLOCATE(grid_particles(k,j,i)%particles)
1295                   ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1296
1297                   grid_particles(k,j,i)%particles(1:number_of_particles)          = tmp_particles_s(1:number_of_particles)
1298                   grid_particles(k,j,i)%particles(number_of_particles+1:new_size) = zero_particle
1299
1300                ELSE
1301
1302                   ALLOCATE(tmp_particles_d(number_of_particles))
1303                   tmp_particles_d(1:number_of_particles) = grid_particles(k,j,i)%particles(1:number_of_particles)
1304
1305                   DEALLOCATE(grid_particles(k,j,i)%particles)
1306                   ALLOCATE(grid_particles(k,j,i)%particles(new_size))
1307
1308                   grid_particles(k,j,i)%particles(1:number_of_particles)          = tmp_particles_d(1:number_of_particles)
1309                   grid_particles(k,j,i)%particles(number_of_particles+1:new_size) = zero_particle
1310
1311                   DEALLOCATE(tmp_particles_d)
1312
1313                ENDIF
1314
1315             ENDIF
1316          ENDDO
1317       ENDDO
1318    ENDDO
1319
1320 END SUBROUTINE dealloc_particles_array
1321
1322
[1359]1323END MODULE lpm_exchange_horiz_mod
Note: See TracBrowser for help on using the repository browser.