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

Last change on this file since 2182 was 2151, checked in by schwenkel, 8 years ago

Bugfix in lpm_exchange_horiz

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