Home

Context Navigation

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

Last change on this file since 2709 was 2709, checked in by suehring, 6 years ago
Bugfixes in CFL check; non-allocated array used in 1D decomposition
Property svn:keywords set to `Id`
File size: 51.3 KB

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