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

Last change on this file since 2330 was 2330, checked in by schwenkel, 4 years ago

Bugfix for gfortran

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