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

Last change on this file since 2801 was 2801, checked in by thiele, 6 years ago

Introduce particle transfer in nested models

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