Home

Context Navigation

lpm_exchange_horiz_mod.f90 @ 1867

Last change on this file since 1867 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 39.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/SOURCE/lpm_exchange_horiz_mod.f90 @ 1867

Download in other formats: