Home

Context Navigation

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

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