Home

Context Navigation

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

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