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

Last change on this file since 2716 was 2716, checked in by kanani, 4 years ago

Correction of "Former revisions" section

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