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

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

Bugfixes in CFL check; non-allocated array used in 1D decomposition

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