Home

Context Navigation

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

Last change on this file since 1433 was 1360, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 53.4 KB

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