Home

Context Navigation

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

Last change on this file since 1989 was 1937, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 46.0 KB

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