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

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

Improved calculation of particle IDs.

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