Home

Context Navigation

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

Last change on this file since 2186 was 2151, checked in by schwenkel, 8 years ago
Bugfix in lpm_exchange_horiz
Property svn:keywords set to `Id`
File size: 46.4 KB

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