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

Last change on this file since 2053 was 2001, checked in by knoop, 8 years ago

last commit documented

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