Home

Context Navigation

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

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