Home

Context Navigation

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

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