Home

Context Navigation

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

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