Home

Context Navigation

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

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