Home

Context Navigation

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

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