Home

Context Navigation

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

Last change on this file since 2709 was 2709, checked in by suehring, 6 years ago
Bugfixes in CFL check; non-allocated array used in 1D decomposition
Property svn:keywords set to `Id`
File size: 51.3 KB

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