Home

Context Navigation

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

Last change on this file since 2501 was 2330, checked in by schwenkel, 7 years ago
Bugfix for gfortran
Property svn:keywords set to `Id`
File size: 49.2 KB

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