Home

Context Navigation

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

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