Home

Context Navigation

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

Last change on this file since 2956 was 2809, checked in by schwenkel, 7 years ago
Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
Property svn:keywords set to `Id`
File size: 51.2 KB

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