Home

Context Navigation

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

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