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

Last change on this file since 2710 was 2710, checked in by suehring, 4 years ago

changes documented

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