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

Last change on this file since 2632 was 2632, checked in by schwenkel, 4 years ago

bugfix in write statement

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