Home

Context Navigation

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

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