Home

Context Navigation

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

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