Home

Context Navigation

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

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