Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |