Home

Context Navigation

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

Last change on this file since 1725 was 1692, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 54.4 KB

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