Home

Context Navigation

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

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