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

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

Remaining error messages revised, comments extended

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