Home

Context Navigation

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

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