Home

Context Navigation

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

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