Home

Context Navigation

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

Last change on this file since 2064 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 46.1 KB

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