Home

Context Navigation

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

Last change on this file since 1634 was 1360, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 53.4 KB

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