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

Last change on this file since 2254 was 2245, checked in by schwenkel, 7 years ago

Bugfix in Add_particles_to_gridcell

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