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

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

enable particle advection with grid stretching and some formatation changes in lpm

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