Home

Context Navigation

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

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