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

Last change on this file since 4013 was 4013, checked in by raasch, 2 years ago

index bugfix concerning reallocation of particle array

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