Home

Context Navigation

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

Last change on this file since 2834 was 2809, checked in by schwenkel, 7 years ago
Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
Property svn:keywords set to `Id`
File size: 51.2 KB

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