Home

Context Navigation

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

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