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

Last change on this file since 2606 was 2606, checked in by schwenkel, 6 years ago

Modified particle box location and further changes in particle model

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