Home

Context Navigation

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

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