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

Last change on this file since 3051 was 3049, checked in by Giersch, 6 years ago

Revision history corrected

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