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

Last change on this file since 1949 was 1937, checked in by suehring, 8 years ago

last commit documented

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