Home

Context Navigation

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

Last change on this file since 1822 was 1822, checked in by hoffmann, 8 years ago
changes in LPM and bulk cloud microphysics
Property svn:keywords set to `Id`
File size: 39.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	!
[1818]	16	! Copyright 1997-2016 Leibniz Universitaet Hannover
[1036]	17	!--------------------------------------------------------------------------------!
	18	!
[849]	19	! Current revisions:
	20	! ------------------
[1822]	21	! Tails removed. Unused variables removed.
[1784]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: lpm_exchange_horiz.f90 1822 2016-04-07 07:49:42Z hoffmann $
	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	! ------------
[1822]	69	! Exchange of particles 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, &
[1822]	95	ONLY: alloc_factor, deleted_particles, grid_particles, &
	96	ibc_par_lr, ibc_par_ns, min_nr_particle, &
	97	mpi_particle_type, number_of_particles, &
	98	offset_ocean_nzt, offset_ocean_nzt_m1, particles, &
	99	particle_type, prt_count, trlp_count_sum, &
[1359]	100	trlp_count_recv_sum, trnp_count_sum, trnp_count_recv_sum, &
	101	trrp_count_sum, trrp_count_recv_sum, trsp_count_sum, &
[1822]	102	trsp_count_recv_sum, zero_particle
[1320]	103
[849]	104	USE pegrid
	105
	106	IMPLICIT NONE
	107
[1682]	108	INTEGER(iwp), PARAMETER :: NR_2_direction_move = 10000 !<
	109	INTEGER(iwp) :: nr_move_north !<
	110	INTEGER(iwp) :: nr_move_south !<
[1359]	111
	112	TYPE(particle_type), DIMENSION(NR_2_direction_move) :: move_also_north
	113	TYPE(particle_type), DIMENSION(NR_2_direction_move) :: move_also_south
	114
	115	SAVE
	116
	117	PRIVATE
	118	PUBLIC lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
	119
	120	INTERFACE lpm_exchange_horiz
	121	MODULE PROCEDURE lpm_exchange_horiz
	122	END INTERFACE lpm_exchange_horiz
	123
	124	INTERFACE lpm_move_particle
	125	MODULE PROCEDURE lpm_move_particle
	126	END INTERFACE lpm_move_particle
	127
	128	INTERFACE realloc_particles_array
	129	MODULE PROCEDURE realloc_particles_array
	130	END INTERFACE realloc_particles_array
	131
	132	CONTAINS
	133
[1682]	134	!------------------------------------------------------------------------------!
	135	! Description:
	136	! ------------
	137	!> Exchange between subdomains.
	138	!> As soon as one particle has moved beyond the boundary of the domain, it
	139	!> is included in the relevant transfer arrays and marked for subsequent
	140	!> deletion on this PE.
	141	!> First sweep for crossings in x direction. Find out first the number of
	142	!> particles to be transferred and allocate temporary arrays needed to store
	143	!> them.
	144	!> For a one-dimensional decomposition along y, no transfer is necessary,
	145	!> because the particle remains on the PE, but the particle coordinate has to
	146	!> be adjusted.
	147	!------------------------------------------------------------------------------!
[1359]	148	SUBROUTINE lpm_exchange_horiz
	149
	150	IMPLICIT NONE
	151
[1682]	152	INTEGER(iwp) :: i !<
	153	INTEGER(iwp) :: ip !<
	154	INTEGER(iwp) :: j !<
	155	INTEGER(iwp) :: jp !<
	156	INTEGER(iwp) :: kp !<
	157	INTEGER(iwp) :: n !<
	158	INTEGER(iwp) :: trlp_count !<
	159	INTEGER(iwp) :: trlp_count_recv !<
	160	INTEGER(iwp) :: trlpt_count !<
	161	INTEGER(iwp) :: trlpt_count_recv !<
	162	INTEGER(iwp) :: trnp_count !<
	163	INTEGER(iwp) :: trnp_count_recv !<
	164	INTEGER(iwp) :: trnpt_count !<
	165	INTEGER(iwp) :: trnpt_count_recv !<
	166	INTEGER(iwp) :: trrp_count !<
	167	INTEGER(iwp) :: trrp_count_recv !<
	168	INTEGER(iwp) :: trrpt_count !<
	169	INTEGER(iwp) :: trrpt_count_recv !<
	170	INTEGER(iwp) :: trsp_count !<
	171	INTEGER(iwp) :: trsp_count_recv !<
	172	INTEGER(iwp) :: trspt_count !<
	173	INTEGER(iwp) :: trspt_count_recv !<
[849]	174
[1682]	175	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvlp !<
	176	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvnp !<
	177	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvrp !<
	178	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: rvsp !<
	179	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trlp !<
	180	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trnp !<
	181	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trrp !<
	182	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trsp !<
[849]	183
[1359]	184	CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'start' )
	185
[849]	186	#if defined( __parallel )
	187
[1822]	188	!
	189	!-- Exchange between subdomains.
	190	!-- As soon as one particle has moved beyond the boundary of the domain, it
	191	!-- is included in the relevant transfer arrays and marked for subsequent
	192	!-- deletion on this PE.
	193	!-- First sweep for crossings in x direction. Find out first the number of
	194	!-- particles to be transferred and allocate temporary arrays needed to store
	195	!-- them.
	196	!-- For a one-dimensional decomposition along y, no transfer is necessary,
	197	!-- because the particle remains on the PE, but the particle coordinate has to
	198	!-- be adjusted.
[849]	199	trlp_count = 0
	200	trlpt_count = 0
	201	trrp_count = 0
	202	trrpt_count = 0
	203
	204	trlp_count_recv = 0
	205	trlpt_count_recv = 0
	206	trrp_count_recv = 0
	207	trrpt_count_recv = 0
	208
	209	IF ( pdims(1) /= 1 ) THEN
	210	!
[1359]	211	!-- First calculate the storage necessary for sending and receiving the data.
	212	!-- Compute only first (nxl) and last (nxr) loop iterration.
	213	DO ip = nxl, nxr, nxr - nxl
	214	DO jp = nys, nyn
	215	DO kp = nzb+1, nzt
[849]	216
[1359]	217	number_of_particles = prt_count(kp,jp,ip)
	218	IF ( number_of_particles <= 0 ) CYCLE
	219	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	220	DO n = 1, number_of_particles
	221	IF ( particles(n)%particle_mask ) THEN
	222	i = ( particles(n)%x + 0.5_wp * dx ) * ddx
	223	!
	224	!-- Above calculation does not work for indices less than zero
	225	IF ( particles(n)%x < -0.5_wp * dx ) i = -1
	226
	227	IF ( i < nxl ) THEN
	228	trlp_count = trlp_count + 1
	229	ELSEIF ( i > nxr ) THEN
	230	trrp_count = trrp_count + 1
	231	ENDIF
	232	ENDIF
	233	ENDDO
	234
	235	ENDDO
	236	ENDDO
[849]	237	ENDDO
	238
	239	IF ( trlp_count == 0 ) trlp_count = 1
	240	IF ( trlpt_count == 0 ) trlpt_count = 1
	241	IF ( trrp_count == 0 ) trrp_count = 1
	242	IF ( trrpt_count == 0 ) trrpt_count = 1
	243
	244	ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
	245
[1359]	246	trlp = zero_particle
	247	trrp = zero_particle
[849]	248
	249	trlp_count = 0
	250	trlpt_count = 0
	251	trrp_count = 0
	252	trrpt_count = 0
	253
	254	ENDIF
[1359]	255	!
	256	!-- Compute only first (nxl) and last (nxr) loop iterration
	257	DO ip = nxl, nxr,nxr-nxl
	258	DO jp = nys, nyn
	259	DO kp = nzb+1, nzt
	260	number_of_particles = prt_count(kp,jp,ip)
	261	IF ( number_of_particles <= 0 ) CYCLE
	262	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	263	DO n = 1, number_of_particles
	264	!
	265	!-- Only those particles that have not been marked as 'deleted' may
	266	!-- be moved.
	267	IF ( particles(n)%particle_mask ) THEN
[849]	268
[1359]	269	i = ( particles(n)%x + 0.5_wp * dx ) * ddx
	270	!
	271	!-- Above calculation does not work for indices less than zero
	272	IF ( particles(n)%x < - 0.5_wp * dx ) i = -1
[849]	273
[1359]	274	IF ( i < nxl ) THEN
	275	IF ( i < 0 ) THEN
	276	!
	277	!-- Apply boundary condition along x
	278	IF ( ibc_par_lr == 0 ) THEN
	279	!
	280	!-- Cyclic condition
	281	IF ( pdims(1) == 1 ) THEN
	282	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
	283	particles(n)%origin_x = ( nx + 1 ) * dx + &
	284	particles(n)%origin_x
	285	ELSE
	286	trlp_count = trlp_count + 1
	287	trlp(trlp_count) = particles(n)
	288	trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
	289	trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
	290	( nx + 1 ) * dx
	291	particles(n)%particle_mask = .FALSE.
	292	deleted_particles = deleted_particles + 1
[849]	293
[1359]	294	IF ( trlp(trlp_count)%x >= (nx + 0.5_wp)* dx - 1.0E-12_wp ) THEN
	295	trlp(trlp_count)%x = trlp(trlp_count)%x - 1.0E-10_wp
	296	!++ why is 1 subtracted in next statement???
	297	trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x - 1
	298	ENDIF
[849]	299
[1359]	300	ENDIF
[849]	301
[1359]	302	ELSEIF ( ibc_par_lr == 1 ) THEN
	303	!
	304	!-- Particle absorption
	305	particles(n)%particle_mask = .FALSE.
	306	deleted_particles = deleted_particles + 1
[849]	307
[1359]	308	ELSEIF ( ibc_par_lr == 2 ) THEN
	309	!
	310	!-- Particle reflection
	311	particles(n)%x = -particles(n)%x
	312	particles(n)%speed_x = -particles(n)%speed_x
[849]	313
[1359]	314	ENDIF
	315	ELSE
	316	!
	317	!-- Store particle data in the transfer array, which will be
	318	!-- send to the neighbouring PE
	319	trlp_count = trlp_count + 1
	320	trlp(trlp_count) = particles(n)
	321	particles(n)%particle_mask = .FALSE.
	322	deleted_particles = deleted_particles + 1
[849]	323
[1359]	324	ENDIF
[849]	325
[1359]	326	ELSEIF ( i > nxr ) THEN
	327	IF ( i > nx ) THEN
	328	!
	329	!-- Apply boundary condition along x
	330	IF ( ibc_par_lr == 0 ) THEN
	331	!
	332	!-- Cyclic condition
	333	IF ( pdims(1) == 1 ) THEN
	334	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
	335	particles(n)%origin_x = particles(n)%origin_x - &
	336	( nx + 1 ) * dx
	337	ELSE
	338	trrp_count = trrp_count + 1
	339	trrp(trrp_count) = particles(n)
	340	trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
	341	trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
	342	( nx + 1 ) * dx
	343	particles(n)%particle_mask = .FALSE.
	344	deleted_particles = deleted_particles + 1
[849]	345
[1359]	346	ENDIF
[849]	347
[1359]	348	ELSEIF ( ibc_par_lr == 1 ) THEN
	349	!
	350	!-- Particle absorption
	351	particles(n)%particle_mask = .FALSE.
	352	deleted_particles = deleted_particles + 1
[849]	353
[1359]	354	ELSEIF ( ibc_par_lr == 2 ) THEN
	355	!
	356	!-- Particle reflection
	357	particles(n)%x = 2 * ( nx * dx ) - particles(n)%x
	358	particles(n)%speed_x = -particles(n)%speed_x
[849]	359
[1359]	360	ENDIF
	361	ELSE
	362	!
	363	!-- Store particle data in the transfer array, which will be send
	364	!-- to the neighbouring PE
	365	trrp_count = trrp_count + 1
	366	trrp(trrp_count) = particles(n)
	367	particles(n)%particle_mask = .FALSE.
	368	deleted_particles = deleted_particles + 1
[849]	369
[1359]	370	ENDIF
[849]	371
[1359]	372	ENDIF
	373	ENDIF
	374	ENDDO
	375	ENDDO
	376	ENDDO
[849]	377	ENDDO
	378
	379	!
	380	!-- Send left boundary, receive right boundary (but first exchange how many
	381	!-- and check, if particle storage must be extended)
	382	IF ( pdims(1) /= 1 ) THEN
	383
	384	CALL MPI_SENDRECV( trlp_count, 1, MPI_INTEGER, pleft, 0, &
	385	trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
	386	comm2d, status, ierr )
	387
[1359]	388	ALLOCATE(rvrp(MAX(1,trrp_count_recv)))
[849]	389
[1359]	390	CALL MPI_SENDRECV( trlp(1)%radius, max(1,trlp_count), mpi_particle_type,&
	391	pleft, 1, rvrp(1)%radius, &
	392	max(1,trrp_count_recv), mpi_particle_type, pright, 1,&
[849]	393	comm2d, status, ierr )
	394
[1359]	395	IF ( trrp_count_recv > 0 ) CALL Add_particles_to_gridcell(rvrp(1:trrp_count_recv))
	396
	397	DEALLOCATE(rvrp)
	398
[849]	399	!
	400	!-- Send right boundary, receive left boundary
	401	CALL MPI_SENDRECV( trrp_count, 1, MPI_INTEGER, pright, 0, &
	402	trlp_count_recv, 1, MPI_INTEGER, pleft, 0, &
	403	comm2d, status, ierr )
	404
[1359]	405	ALLOCATE(rvlp(MAX(1,trlp_count_recv)))
[849]	406
[1359]	407	CALL MPI_SENDRECV( trrp(1)%radius, max(1,trrp_count), mpi_particle_type,&
	408	pright, 1, rvlp(1)%radius, &
	409	max(1,trlp_count_recv), mpi_particle_type, pleft, 1, &
[849]	410	comm2d, status, ierr )
	411
[1359]	412	IF ( trlp_count_recv > 0 ) CALL Add_particles_to_gridcell(rvlp(1:trlp_count_recv))
	413
[1822]	414	DEALLOCATE( rvlp )
[1359]	415	DEALLOCATE( trlp, trrp )
[849]	416
	417	ENDIF
	418
	419
	420	!
	421	!-- Check whether particles have crossed the boundaries in y direction. Note
	422	!-- that this case can also apply to particles that have just been received
	423	!-- from the adjacent right or left PE.
	424	!-- Find out first the number of particles to be transferred and allocate
	425	!-- temporary arrays needed to store them.
	426	!-- For a one-dimensional decomposition along x, no transfer is necessary,
	427	!-- because the particle remains on the PE.
[1359]	428	trsp_count = nr_move_south
[849]	429	trspt_count = 0
[1359]	430	trnp_count = nr_move_north
[849]	431	trnpt_count = 0
	432
	433	trsp_count_recv = 0
	434	trspt_count_recv = 0
	435	trnp_count_recv = 0
	436	trnpt_count_recv = 0
	437
	438	IF ( pdims(2) /= 1 ) THEN
	439	!
	440	!-- First calculate the storage necessary for sending and receiving the
	441	!-- data
[1359]	442	DO ip = nxl, nxr
	443	DO jp = nys, nyn, nyn-nys !compute only first (nys) and last (nyn) loop iterration
	444	DO kp = nzb+1, nzt
	445	number_of_particles = prt_count(kp,jp,ip)
	446	IF ( number_of_particles <= 0 ) CYCLE
	447	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	448	DO n = 1, number_of_particles
	449	IF ( particles(n)%particle_mask ) THEN
	450	j = ( particles(n)%y + 0.5_wp * dy ) * ddy
[849]	451	!
[1359]	452	!-- Above calculation does not work for indices less than zero
	453	IF ( particles(n)%y < -0.5_wp * dy ) j = -1
[849]	454
[1359]	455	IF ( j < nys ) THEN
	456	trsp_count = trsp_count + 1
	457	ELSEIF ( j > nyn ) THEN
	458	trnp_count = trnp_count + 1
	459	ENDIF
	460	ENDIF
	461	ENDDO
	462	ENDDO
	463	ENDDO
[849]	464	ENDDO
	465
	466	IF ( trsp_count == 0 ) trsp_count = 1
	467	IF ( trspt_count == 0 ) trspt_count = 1
	468	IF ( trnp_count == 0 ) trnp_count = 1
	469	IF ( trnpt_count == 0 ) trnpt_count = 1
	470
	471	ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
	472
[1359]	473	trsp = zero_particle
	474	trnp = zero_particle
[849]	475
[1359]	476	trsp_count = nr_move_south
[849]	477	trspt_count = 0
[1359]	478	trnp_count = nr_move_north
[849]	479	trnpt_count = 0
	480
[1359]	481	trsp(1:nr_move_south) = move_also_south(1:nr_move_south)
	482	trnp(1:nr_move_north) = move_also_north(1:nr_move_north)
	483
[849]	484	ENDIF
	485
[1359]	486	DO ip = nxl, nxr
	487	DO jp = nys, nyn, nyn-nys ! compute only first (nys) and last (nyn) loop iterration
	488	DO kp = nzb+1, nzt
	489	number_of_particles = prt_count(kp,jp,ip)
	490	IF ( number_of_particles <= 0 ) CYCLE
	491	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	492	DO n = 1, number_of_particles
[849]	493	!
[1359]	494	!-- Only those particles that have not been marked as 'deleted' may
	495	!-- be moved.
	496	IF ( particles(n)%particle_mask ) THEN
	497	j = ( particles(n)%y + 0.5_wp * dy ) * ddy
[849]	498	!
[1359]	499	!-- Above calculation does not work for indices less than zero
	500	IF ( particles(n)%y < -0.5_wp * dy ) j = -1
[849]	501
[1359]	502	IF ( j < nys ) THEN
	503	IF ( j < 0 ) THEN
[849]	504	!
[1359]	505	!-- Apply boundary condition along y
	506	IF ( ibc_par_ns == 0 ) THEN
[849]	507	!
[1359]	508	!-- Cyclic condition
	509	IF ( pdims(2) == 1 ) THEN
	510	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
	511	particles(n)%origin_y = ( ny + 1 ) * dy + &
	512	particles(n)%origin_y
	513	ELSE
	514	trsp_count = trsp_count + 1
	515	trsp(trsp_count) = particles(n)
	516	trsp(trsp_count)%y = ( ny + 1 ) * dy + &
	517	trsp(trsp_count)%y
	518	trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
	519	+ ( ny + 1 ) * dy
	520	particles(n)%particle_mask = .FALSE.
	521	deleted_particles = deleted_particles + 1
[849]	522
[1359]	523	IF ( trsp(trsp_count)%y >= (ny+0.5_wp)* dy - 1.0E-12_wp ) THEN
	524	trsp(trsp_count)%y = trsp(trsp_count)%y - 1.0E-10_wp
	525	!++ why is 1 subtracted in next statement???
	526	trsp(trsp_count)%origin_y = &
	527	trsp(trsp_count)%origin_y - 1
	528	ENDIF
[849]	529
[1359]	530	ENDIF
[849]	531
[1359]	532	ELSEIF ( ibc_par_ns == 1 ) THEN
[849]	533	!
[1359]	534	!-- Particle absorption
	535	particles(n)%particle_mask = .FALSE.
	536	deleted_particles = deleted_particles + 1
[849]	537
[1359]	538	ELSEIF ( ibc_par_ns == 2 ) THEN
[849]	539	!
[1359]	540	!-- Particle reflection
	541	particles(n)%y = -particles(n)%y
	542	particles(n)%speed_y = -particles(n)%speed_y
[849]	543
[1359]	544	ENDIF
	545	ELSE
[849]	546	!
[1359]	547	!-- Store particle data in the transfer array, which will
	548	!-- be send to the neighbouring PE
	549	trsp_count = trsp_count + 1
	550	trsp(trsp_count) = particles(n)
	551	particles(n)%particle_mask = .FALSE.
	552	deleted_particles = deleted_particles + 1
[849]	553
[1359]	554	ENDIF
[849]	555
[1359]	556	ELSEIF ( j > nyn ) THEN
	557	IF ( j > ny ) THEN
[849]	558	!
[1359]	559	!-- Apply boundary condition along x
	560	IF ( ibc_par_ns == 0 ) THEN
[849]	561	!
[1359]	562	!-- Cyclic condition
	563	IF ( pdims(2) == 1 ) THEN
	564	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
	565	particles(n)%origin_y = &
	566	particles(n)%origin_y - ( ny + 1 ) * dy
	567	ELSE
	568	trnp_count = trnp_count + 1
	569	trnp(trnp_count) = particles(n)
	570	trnp(trnp_count)%y = &
	571	trnp(trnp_count)%y - ( ny + 1 ) * dy
	572	trnp(trnp_count)%origin_y = &
	573	trnp(trnp_count)%origin_y - ( ny + 1 ) * dy
	574	particles(n)%particle_mask = .FALSE.
	575	deleted_particles = deleted_particles + 1
	576	ENDIF
[849]	577
[1359]	578	ELSEIF ( ibc_par_ns == 1 ) THEN
[849]	579	!
[1359]	580	!-- Particle absorption
	581	particles(n)%particle_mask = .FALSE.
	582	deleted_particles = deleted_particles + 1
[849]	583
[1359]	584	ELSEIF ( ibc_par_ns == 2 ) THEN
[849]	585	!
[1359]	586	!-- Particle reflection
	587	particles(n)%y = 2 * ( ny * dy ) - particles(n)%y
	588	particles(n)%speed_y = -particles(n)%speed_y
[849]	589
[1359]	590	ENDIF
	591	ELSE
[849]	592	!
[1359]	593	!-- Store particle data in the transfer array, which will
	594	!-- be send to the neighbouring PE
	595	trnp_count = trnp_count + 1
	596	trnp(trnp_count) = particles(n)
	597	particles(n)%particle_mask = .FALSE.
	598	deleted_particles = deleted_particles + 1
[849]	599
[1359]	600	ENDIF
	601
	602	ENDIF
[849]	603	ENDIF
[1359]	604	ENDDO
	605	ENDDO
	606	ENDDO
[849]	607	ENDDO
	608
	609	!
	610	!-- Send front boundary, receive back boundary (but first exchange how many
	611	!-- and check, if particle storage must be extended)
	612	IF ( pdims(2) /= 1 ) THEN
	613
	614	CALL MPI_SENDRECV( trsp_count, 1, MPI_INTEGER, psouth, 0, &
	615	trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
	616	comm2d, status, ierr )
	617
[1359]	618	ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
[849]	619
[1359]	620	CALL MPI_SENDRECV( trsp(1)%radius, trsp_count, mpi_particle_type, &
	621	psouth, 1, rvnp(1)%radius, &
[849]	622	trnp_count_recv, mpi_particle_type, pnorth, 1, &
	623	comm2d, status, ierr )
	624
[1359]	625	IF ( trnp_count_recv > 0 ) CALL Add_particles_to_gridcell(rvnp(1:trnp_count_recv))
	626
	627	DEALLOCATE(rvnp)
	628
[849]	629	!
	630	!-- Send back boundary, receive front boundary
	631	CALL MPI_SENDRECV( trnp_count, 1, MPI_INTEGER, pnorth, 0, &
	632	trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
	633	comm2d, status, ierr )
	634
[1359]	635	ALLOCATE(rvsp(MAX(1,trsp_count_recv)))
[849]	636
[1359]	637	CALL MPI_SENDRECV( trnp(1)%radius, trnp_count, mpi_particle_type, &
	638	pnorth, 1, rvsp(1)%radius, &
[849]	639	trsp_count_recv, mpi_particle_type, psouth, 1, &
	640	comm2d, status, ierr )
	641
[1359]	642	IF ( trsp_count_recv > 0 ) CALL Add_particles_to_gridcell(rvsp(1:trsp_count_recv))
	643
	644	DEALLOCATE(rvsp)
	645
[849]	646	number_of_particles = number_of_particles + trsp_count_recv
	647
	648	DEALLOCATE( trsp, trnp )
	649
	650	ENDIF
	651
	652	#else
	653
	654	!
	655	!-- Apply boundary conditions
	656	DO n = 1, number_of_particles
	657
[1359]	658	IF ( particles(n)%x < -0.5_wp * dx ) THEN
[849]	659
	660	IF ( ibc_par_lr == 0 ) THEN
	661	!
	662	!-- Cyclic boundary. Relevant coordinate has to be changed.
	663	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
[1822]	664
[849]	665	ELSEIF ( ibc_par_lr == 1 ) THEN
	666	!
	667	!-- Particle absorption
[1359]	668	particles(n)%particle_mask = .FALSE.
[849]	669	deleted_particles = deleted_particles + 1
[1822]	670
[849]	671	ELSEIF ( ibc_par_lr == 2 ) THEN
	672	!
	673	!-- Particle reflection
	674	particles(n)%x = -dx - particles(n)%x
	675	particles(n)%speed_x = -particles(n)%speed_x
	676	ENDIF
	677
[1359]	678	ELSEIF ( particles(n)%x >= ( nx + 0.5_wp ) * dx ) THEN
[849]	679
	680	IF ( ibc_par_lr == 0 ) THEN
	681	!
	682	!-- Cyclic boundary. Relevant coordinate has to be changed.
	683	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
[1822]	684
[849]	685	ELSEIF ( ibc_par_lr == 1 ) THEN
	686	!
	687	!-- Particle absorption
[1359]	688	particles(n)%particle_mask = .FALSE.
[849]	689	deleted_particles = deleted_particles + 1
[1822]	690
[849]	691	ELSEIF ( ibc_par_lr == 2 ) THEN
	692	!
	693	!-- Particle reflection
	694	particles(n)%x = ( nx + 1 ) * dx - particles(n)%x
	695	particles(n)%speed_x = -particles(n)%speed_x
	696	ENDIF
	697
	698	ENDIF
	699
[1359]	700	IF ( particles(n)%y < -0.5_wp * dy ) THEN
[849]	701
	702	IF ( ibc_par_ns == 0 ) THEN
	703	!
	704	!-- Cyclic boundary. Relevant coordinate has to be changed.
	705	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
[1822]	706
[849]	707	ELSEIF ( ibc_par_ns == 1 ) THEN
	708	!
	709	!-- Particle absorption
[1359]	710	particles(n)%particle_mask = .FALSE.
[849]	711	deleted_particles = deleted_particles + 1
[1822]	712
[849]	713	ELSEIF ( ibc_par_ns == 2 ) THEN
	714	!
	715	!-- Particle reflection
	716	particles(n)%y = -dy - particles(n)%y
	717	particles(n)%speed_y = -particles(n)%speed_y
	718	ENDIF
	719
[1359]	720	ELSEIF ( particles(n)%y >= ( ny + 0.5_wp ) * dy ) THEN
[849]	721
	722	IF ( ibc_par_ns == 0 ) THEN
	723	!
	724	!-- Cyclic boundary. Relevant coordinate has to be changed.
	725	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
[1822]	726
[849]	727	ELSEIF ( ibc_par_ns == 1 ) THEN
	728	!
	729	!-- Particle absorption
[1359]	730	particles(n)%particle_mask = .FALSE.
[849]	731	deleted_particles = deleted_particles + 1
[1822]	732
[849]	733	ELSEIF ( ibc_par_ns == 2 ) THEN
	734	!
	735	!-- Particle reflection
	736	particles(n)%y = ( ny + 1 ) * dy - particles(n)%y
	737	particles(n)%speed_y = -particles(n)%speed_y
	738	ENDIF
	739
	740	ENDIF
	741	ENDDO
	742
	743	#endif
	744
	745	!
	746	!-- Accumulate the number of particles transferred between the subdomains
	747	#if defined( __parallel )
	748	trlp_count_sum = trlp_count_sum + trlp_count
	749	trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
	750	trrp_count_sum = trrp_count_sum + trrp_count
	751	trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
	752	trsp_count_sum = trsp_count_sum + trsp_count
	753	trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
	754	trnp_count_sum = trnp_count_sum + trnp_count
	755	trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
	756	#endif
	757
[1359]	758	CALL cpu_log( log_point_s(23), 'lpm_exchange_horiz', 'stop' )
[849]	759
	760	END SUBROUTINE lpm_exchange_horiz
[1359]	761
[1682]	762	!------------------------------------------------------------------------------!
	763	! Description:
	764	! ------------
	765	!> If a particle moves from one processor to another, this subroutine moves
	766	!> the corresponding elements from the particle arrays of the old grid cells
	767	!> to the particle arrays of the new grid cells.
	768	!------------------------------------------------------------------------------!
[1359]	769	SUBROUTINE Add_particles_to_gridcell (particle_array)
	770
	771	IMPLICIT NONE
	772
[1682]	773	INTEGER(iwp) :: ip !<
	774	INTEGER(iwp) :: jp !<
	775	INTEGER(iwp) :: kp !<
	776	INTEGER(iwp) :: n !<
	777	INTEGER(iwp) :: pindex !<
[1359]	778
[1682]	779	LOGICAL :: pack_done !<
[1359]	780
	781	TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_array
	782
	783	pack_done = .FALSE.
	784
	785	nr_move_north = 0
	786	nr_move_south = 0
	787
	788	DO n = 1, SIZE(particle_array)
	789	ip = ( particle_array(n)%x + 0.5_wp * dx ) * ddx
	790	jp = ( particle_array(n)%y + 0.5_wp * dy ) * ddy
[1685]	791	kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt
[1359]	792
	793	IF ( ip >= nxl .AND. ip <= nxr .AND. jp >= nys .AND. jp <= nyn &
	794	.AND. kp >= nzb+1 .AND. kp <= nzt) THEN ! particle stays on processor
	795	number_of_particles = prt_count(kp,jp,ip)
	796	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
	797
	798	pindex = prt_count(kp,jp,ip)+1
	799	IF( pindex > SIZE(grid_particles(kp,jp,ip)%particles) ) THEN
	800	IF ( pack_done ) THEN
	801	CALL realloc_particles_array (ip,jp,kp)
	802	ELSE
	803	CALL lpm_pack_arrays
	804	prt_count(kp,jp,ip) = number_of_particles
	805	pindex = prt_count(kp,jp,ip)+1
	806	IF ( pindex > SIZE(grid_particles(kp,jp,ip)%particles) ) THEN
	807	CALL realloc_particles_array (ip,jp,kp)
	808	ENDIF
	809	pack_done = .TRUE.
	810	ENDIF
	811	ENDIF
	812	grid_particles(kp,jp,ip)%particles(pindex) = particle_array(n)
	813	prt_count(kp,jp,ip) = pindex
	814	ELSE
	815	IF ( jp == nys - 1 ) THEN
	816	nr_move_south = nr_move_south+1
	817	move_also_south(nr_move_south) = particle_array(n)
	818	IF ( jp == -1 ) THEN
	819	move_also_south(nr_move_south)%y = &
	820	move_also_south(nr_move_south)%y + ( ny + 1 ) * dy
	821	move_also_south(nr_move_south)%origin_y = &
	822	move_also_south(nr_move_south)%origin_y + ( ny + 1 ) * dy
	823	ENDIF
	824	ELSEIF ( jp == nyn+1 ) THEN
	825	nr_move_north = nr_move_north+1
	826	move_also_north(nr_move_north) = particle_array(n)
	827	IF ( jp == ny+1 ) THEN
	828	move_also_north(nr_move_north)%y = &
	829	move_also_north(nr_move_north)%y - ( ny + 1 ) * dy
	830	move_also_north(nr_move_north)%origin_y = &
	831	move_also_north(nr_move_north)%origin_y - ( ny + 1 ) * dy
	832	ENDIF
	833	ELSE
	834	WRITE(0,'(a,8i7)') 'particle out of range ',myid,ip,jp,kp,nxl,nxr,nys,nyn
	835	ENDIF
	836	ENDIF
	837	ENDDO
	838
	839	RETURN
	840
	841	END SUBROUTINE Add_particles_to_gridcell
	842
	843
	844
	845
[1682]	846	!------------------------------------------------------------------------------!
	847	! Description:
	848	! ------------
	849	!> If a particle moves from one grid cell to another (on the current
	850	!> processor!), this subroutine moves the corresponding element from the
	851	!> particle array of the old grid cell to the particle array of the new grid
	852	!> cell.
	853	!------------------------------------------------------------------------------!
[1359]	854	SUBROUTINE lpm_move_particle
	855
	856	IMPLICIT NONE
	857
[1682]	858	INTEGER(iwp) :: i !<
	859	INTEGER(iwp) :: ip !<
	860	INTEGER(iwp) :: j !<
	861	INTEGER(iwp) :: jp !<
	862	INTEGER(iwp) :: k !<
	863	INTEGER(iwp) :: kp !<
	864	INTEGER(iwp) :: n !<
	865	INTEGER(iwp) :: np_old_cell !<
	866	INTEGER(iwp) :: n_start !<
	867	INTEGER(iwp) :: pindex !<
[1359]	868
[1682]	869	LOGICAL :: pack_done !<
[1359]	870
[1682]	871	TYPE(particle_type), DIMENSION(:), POINTER :: particles_old_cell !<
[1359]	872
	873	CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'start' )
	874
	875	DO ip = nxl, nxr
	876	DO jp = nys, nyn
	877	DO kp = nzb+1, nzt
	878
	879	np_old_cell = prt_count(kp,jp,ip)
	880	IF ( np_old_cell <= 0 ) CYCLE
	881	particles_old_cell => grid_particles(kp,jp,ip)%particles(1:np_old_cell)
	882	n_start = -1
	883
	884	DO n = 1, np_old_cell
	885	i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
	886	j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
	887	k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
	888	!
	889	!-- Check, if particle has moved to another grid cell.
	890	IF ( i /= ip .OR. j /= jp .OR. k /= kp ) THEN
	891	!
	892	!-- The particle has moved to another grid cell. Now check, if
	893	!-- particle stays on the same processor.
	894	IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. &
	895	j <= nyn .AND. k >= nzb+1 .AND. k <= nzt) THEN
	896	!
	897	!-- If the particle stays on the same processor, the particle
	898	!-- will be added to the particle array of the new processor.
	899	number_of_particles = prt_count(k,j,i)
	900	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	901
	902	pindex = prt_count(k,j,i)+1
	903	IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) &
	904	THEN
	905	n_start = n
	906	EXIT
	907	ENDIF
	908
	909	grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
	910	prt_count(k,j,i) = pindex
	911
	912	particles_old_cell(n)%particle_mask = .FALSE.
	913	ENDIF
	914	ENDIF
	915	ENDDO
	916
[1691]	917	IF ( n_start >= 0 ) THEN
[1359]	918	pack_done = .FALSE.
	919	DO n = n_start, np_old_cell
	920	i = ( particles_old_cell(n)%x + 0.5_wp * dx ) * ddx
	921	j = ( particles_old_cell(n)%y + 0.5_wp * dy ) * ddy
	922	k = particles_old_cell(n)%z / dz + 1 + offset_ocean_nzt
	923	IF ( i /= ip .OR. j /= jp .OR. k /= kp ) THEN
	924	!
	925	!-- Particle is in different box
	926	IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. &
	927	j <= nyn .AND. k >= nzb+1 .AND. k <= nzt) THEN
	928	!
	929	!-- Particle stays on processor
	930	number_of_particles = prt_count(k,j,i)
	931	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	932
	933	pindex = prt_count(k,j,i)+1
	934	IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) &
	935	THEN
	936	IF ( pack_done ) THEN
	937	CALL realloc_particles_array(i,j,k)
	938	pindex = prt_count(k,j,i)+1
	939	ELSE
	940	CALL lpm_pack_arrays
	941	prt_count(k,j,i) = number_of_particles
	942
	943	pindex = prt_count(k,j,i)+1
	944	!
	945	!-- If number of particles in the new grid box
	946	!-- exceeds its allocated memory, the particle array
	947	!-- will be reallocated
	948	IF ( pindex > SIZE(grid_particles(k,j,i)%particles) ) THEN
	949	CALL realloc_particles_array(i,j,k)
	950	pindex = prt_count(k,j,i)+1
	951	ENDIF
	952
	953	pack_done = .TRUE.
	954	ENDIF
	955	ENDIF
	956
	957	grid_particles(k,j,i)%particles(pindex) = particles_old_cell(n)
	958	prt_count(k,j,i) = pindex
	959
	960	particles_old_cell(n)%particle_mask = .FALSE.
	961	ENDIF
	962	ENDIF
	963	ENDDO
	964	ENDIF
	965	ENDDO
	966	ENDDO
	967	ENDDO
	968
	969	CALL cpu_log( log_point_s(41), 'lpm_move_particle', 'stop' )
	970
	971	RETURN
	972
	973	END SUBROUTINE lpm_move_particle
	974
[1682]	975	!------------------------------------------------------------------------------!
	976	! Description:
	977	! ------------
	978	!> @todo Missing subroutine description.
	979	!------------------------------------------------------------------------------!
[1359]	980	SUBROUTINE realloc_particles_array (i,j,k,size_in)
	981
	982	IMPLICIT NONE
	983
[1682]	984	INTEGER(iwp), INTENT(in) :: i !<
	985	INTEGER(iwp), INTENT(in) :: j !<
	986	INTEGER(iwp), INTENT(in) :: k !<
[1822]	987	INTEGER(iwp), INTENT(in), OPTIONAL :: size_in !<
[1359]	988
[1682]	989	INTEGER(iwp) :: old_size !<
	990	INTEGER(iwp) :: new_size !<
	991	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles_d !<
	992	TYPE(particle_type), DIMENSION(500) :: tmp_particles_s !<
[1359]	993
	994
	995	old_size = SIZE(grid_particles(k,j,i)%particles)
	996
	997	IF ( PRESENT(size_in) ) THEN
	998	new_size = size_in
	999	ELSE
[1822]	1000	new_size = old_size * ( 1.0_wp + alloc_factor / 100.0_wp )
[1359]	1001	ENDIF
	1002
	1003	new_size = MAX( new_size, min_nr_particle )
	1004
	1005	IF ( old_size <= 500 ) THEN
	1006
	1007	tmp_particles_s(1:old_size) = grid_particles(k,j,i)%particles(1:old_size)
	1008
	1009	DEALLOCATE(grid_particles(k,j,i)%particles)
	1010	ALLOCATE(grid_particles(k,j,i)%particles(new_size))
	1011
	1012	grid_particles(k,j,i)%particles(1:old_size) = tmp_particles_s(1:old_size)
	1013	grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
	1014
	1015	ELSE
	1016
	1017	ALLOCATE(tmp_particles_d(new_size))
	1018	tmp_particles_d(1:old_size) = grid_particles(k,j,i)%particles
	1019
	1020	DEALLOCATE(grid_particles(k,j,i)%particles)
	1021	ALLOCATE(grid_particles(k,j,i)%particles(new_size))
	1022
	1023	grid_particles(k,j,i)%particles(1:old_size) = tmp_particles_d(1:old_size)
	1024	grid_particles(k,j,i)%particles(old_size+1:new_size) = zero_particle
	1025
	1026	DEALLOCATE(tmp_particles_d)
	1027
	1028	ENDIF
	1029	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	1030
	1031	RETURN
	1032	END SUBROUTINE realloc_particles_array
	1033
	1034	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 |