Home

Context Navigation

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

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