Home

Context Navigation

source: palm/trunk/SOURCE/advec_particles.f90 @ 30

Last change on this file since 30 was 16, checked in by raasch, 18 years ago
bugfixes in advec_particles and init_particles
Property svn:keywords set to `Id`
File size: 135.2 KB

Rev	Line
[1]	1	SUBROUTINE advec_particles
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[16]	6	! Bugfix: wrong if-clause from revision 1.32
[1]	7	! TEST: PRINT statements on unit 9 (commented out)
	8	!
	9	! Former revisions:
	10	! -----------------
[3]	11	! $Id: advec_particles.f90 16 2007-02-15 13:16:47Z raasch $
	12	! RCS Log replace by Id keyword, revision history cleaned up
	13	!
[1]	14	! Revision 1.32 2007/02/11 12:48:20 raasch
	15	! Allways the lower level k is used for interpolation
	16	! Bugfix: new particles are released only if end_time_prel > simulated_time
	17	! Bugfix: transfer of particles when x < -0.5*dx (0.0 before), etc.,
	18	! index i,j used instead of cartesian (x,y) coordinate to check for
	19	! transfer because this failed under very rare conditions
	20	! Bugfix: calculation of number of particles with same radius as the current
	21	! particle (cloud droplet code)
	22	!
	23	! Revision 1.31 2006/08/17 09:21:01 raasch
	24	! Two more compilation errors removed from the last revision
	25	!
	26	! Revision 1.30 2006/08/17 09:11:17 raasch
	27	! Two compilation errors removed from the last revision
	28	!
	29	! Revision 1.29 2006/08/04 14:05:01 raasch
	30	! Subgrid scale velocities are (optionally) included for calculating the
	31	! particle advection, new counters trlp_count_sum, etc. for accumulating
	32	! the number of particles exchanged between the subdomains during all
	33	! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z,
	34	! izuf renamed iran, output of particle time series
	35	!
	36	! Revision 1.28 2006/02/23 09:41:21 raasch
	37	! Only a fraction of the particles may have tails, therefore output for unit
	38	! 85 and 90 has changed,
	39	! sorting of particles in gridbox order after each timestep, optimization of
	40	! the droplet growth by collision,
	41	! 0.5 replaced by 0.4999999999 in setting the random fluctuations,
	42	! psl, psr, pdx, etc. are now 1d arrays, i.e. they depend on the particle group,
	43	! improve particle release at PE boundaries, nt_anz renamed
	44	! current_timestep_number, idum (particle_type) renamed tail_id,
	45	! error number argument for handle_netcdf_error,
	46	!
	47	! Revision 1.27 2005/10/20 14:04:15 raasch
	48	! droplet growth by collision completed (still has to be optimized for speed),
	49	! new routine collision_efficiency (see end of this file),
	50	! 2*r replaced by r in the exponential terms (momentum equation for particles),
	51	! number of particles really used is additionally output on the netcdf particle
	52	! data file,
	53	! test print statement removed
	54	!
	55	! Revision 1.26 2005/06/26 19:50:14 raasch
	56	! Cloud droplet features added, radius is used instead of diameter,
	57	! particle type initial data assignment extended,
	58	! output on unit 80 for non-parallel runs
	59	!
	60	! Revision 1.25 2005/05/18 14:53:11 raasch
	61	! Extensions for NetCDF output
	62	!
	63	! Revision 1.24 2005/04/23 08:24:58 raasch
	64	! Revised calculation of output time counters regarding a possible decrease of
	65	! the output time interval in case of restart runs
	66	!
	67	! Revision 1.23 2005/03/26 14:20:05 raasch
	68	! calculation of vertical particle velocity (with inertia) corrected, exp_arg
	69	! had a wrong sign
	70	!
	71	! Revision 1.22 2004/04/30 07:48:57 raasch
	72	! Particle type initial data assignment extended
	73	!
	74	! Revision 1.21 2004/01/28 15:00:48 raasch
	75	! Output of particle infos in subroutine allocate_prt_memory on demand only
	76	!
	77	! Revision 1.20 2003/10/29 08:38:52 raasch
	78	! Module random_function_mod is used, modifications for new particle group
	79	! feature, version number is written on binary file
	80	!
	81	! Revision 1.19 2003/03/16 09:19:32 raasch
	82	! Two underscores (_) are placed in front of all define-strings
	83	!
	84	! Revision 1.18 2003/03/14 13:36:31 raasch
	85	! Error in reflection boundary condition removed: particle speed must also
	86	! be inverted
	87	!
	88	! Revision 1.17 2003/03/04 11:23:20 raasch
	89	! Error in particle inertia part removed (exp_arg must not contain the timestep)
	90	! Particle velocities are also stored in particles in case of zero density ratio
	91	!
	92	! Revision 1.16 2002/09/12 12:55:29 raasch
	93	! Transport of particles with inertia implemented, write density_ratio to
	94	! restart file
	95	!
	96	! Revision 1.15 2002/06/11 12:11:56 raasch
	97	! Declaration of u_int_l, u_int_u, v_int_l, v_int_u added
	98	!
	99	! Revision 1.14 2002/05/02 18:46:34 18:46:34 raasch (Siegfried Raasch)
	100	! Horizontal velocity components for particle advection are now interpolated
	101	! between horizontal grid levels
	102	!
	103	! Revision 1.13 2002/04/24 18:46:11 18:46:11 raasch (Siegfried Raasch)
	104	! Temporary particle-arrays trlp, trnp, trrp and trsp are now allocatable arrays
	105	!
	106	! Revision 1.12 2002/04/16 07:44:07 raasch
	107	! Additional boundary conditions added, particle data for later analysis can be
	108	! written on file
	109	!
	110	! Revision 1.11 2001/11/09 13:38:20 raasch
	111	! Colour information for particle tails is stored in array
	112	! particle_tail_coordinates, adding information to the particle tails
	113	! moved after call of user_particle_attributes
	114	!
	115	! Revision 1.10 2001/08/21 08:19:48 raasch
	116	! Storage of particle tails, precision of vertical advection improved,
	117	! particle attributes (e.g. size and color) can be defined by the user
	118	!
	119	! Revision 1.9 2001/07/12 12:01:27 raasch
	120	! Random start positions of initial particles possible
	121	!
	122	! Revision 1.8 2001/03/29 17:27:10 raasch
	123	! Translation of remaining German identifiers (variables, subroutines, etc.)
	124	!
	125	! Revision 1.7 2001/02/13 10:24:00 raasch
	126	! Particle advection possible in case of galilei transformation
	127	!
	128	! Revision 1.6 2001/01/25 06:51:04 raasch
	129	! Module test_variables removed, writing of particle informations is optional
	130	! now
	131	!
	132	! Revision 1.5 2001/01/02 17:15:29 raasch
	133	! Unit 80 is immediately closed after usage every time. Particle positions
	134	! are written on unit 90 instead of 81 (subroutine write_particles).
	135	!
	136	! Revision 1.4 2000/12/28 12:43:54 raasch
	137	! Packing of particles-array is done explicitly (not by PACK-function, which
	138	! seems to have problems when a large number of PEs is used)
	139	! Missing translations included,
	140	! code is used only optionally (cpp-directives are added)
	141	!
	142	! Revision 1.3 2000/01/20 08:51:37 letzel
	143	! All comments translated into English
	144	!
	145	! Revision 1.2 1999/11/25 16:41:48 raasch
	146	! Indexfehler im nichtparallelen Teil korrigiert
	147	!
	148	! Revision 1.1 1999/11/25 16:16:06 raasch
	149	! Initial revision
	150	!
	151	!
	152	! Description:
	153	! ------------
	154	! Particle advection
	155	!------------------------------------------------------------------------------!
	156	#if defined( __particles )
	157
	158	USE arrays_3d
	159	USE cloud_parameters
	160	USE constants
	161	USE control_parameters
	162	USE cpulog
	163	USE grid_variables
	164	USE indices
	165	USE interfaces
	166	USE netcdf_control
	167	USE particle_attributes
	168	USE pegrid
	169	USE random_function_mod
	170	USE statistics
	171
	172	IMPLICIT NONE
	173
	174	INTEGER :: deleted_particles, deleted_tails, i, ie, ii, inc, is, j, jj, &
	175	js, k, kk, kw, m, n, nc, nn, psi, tlength, trlp_count, &
	176	trlp_count_sum, trlp_count_recv, trlp_count_recv_sum, &
	177	trlpt_count, trlpt_count_recv, trnp_count, trnp_count_sum, &
	178	trnp_count_recv, trnp_count_recv_sum, trnpt_count, &
	179	trnpt_count_recv, trrp_count, trrp_count_sum, trrp_count_recv, &
	180	trrp_count_recv_sum, trrpt_count, trrpt_count_recv, &
	181	trsp_count, trsp_count_sum, trsp_count_recv, &
	182	trsp_count_recv_sum, trspt_count, trspt_count_recv, nd
	183
	184	LOGICAL :: dt_3d_reached, dt_3d_reached_l, prt_position
	185
	186	REAL :: aa, arg, bb, cc, dd, delta_r, dens_ratio, de_dx_int, &
	187	de_dx_int_l, de_dx_int_u, de_dy_int, de_dy_int_l, de_dy_int_u, &
	188	de_dz_int, de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, &
	189	diss_int_u, distance, dt_gap, dt_particle, d_radius, e_a, &
	190	e_int, e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, &
	191	fs_int, gg, lagr_timescale, mean_r, new_r, p_int, pt_int, &
	192	pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, &
	193	ql_int_u, random_gauss, sl_r3, sl_r4, s_r3, s_r4, t_int, &
	194	u_int, u_int_l, u_int_u, vv_int, v_int, v_int_l, v_int_u, &
	195	w_int, w_int_l, w_int_u, x, y
	196
	197	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: de_dx, de_dy, de_dz
	198
	199	REAL, DIMENSION(:,:,:), ALLOCATABLE :: trlpt, trnpt, trrpt, trspt
	200
	201	TYPE(particle_type) :: tmp_particle
	202
	203	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trlp, trnp, trrp, trsp
	204
	205
	206	CALL cpu_log( log_point(25), 'advec_particles', 'start' )
	207
	208	! IF ( number_of_particles /= number_of_tails ) THEN
	209	! WRITE (9,*) '--- advec_particles: #1'
	210	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
[13]	211	! CALL FLUSH_( 9 )
[1]	212	! ENDIF
	213	!
	214	!-- Write particle data on file for later analysis.
	215	!-- This has to be done here (before particles are advected) in order
	216	!-- to allow correct output in case of dt_write_particle_data = dt_prel =
	217	!-- particle_maximum_age. Otherwise (if output is done at the end of this
	218	!-- subroutine), the relevant particles would have been already deleted.
	219	!-- The MOD function allows for changes in the output interval with restart
	220	!-- runs.
	221	!-- Attention: change version number for unit 85 (in routine check_open)
	222	!-- whenever the output format for this unit is changed!
	223	time_write_particle_data = time_write_particle_data + dt_3d
	224	IF ( time_write_particle_data >= dt_write_particle_data ) THEN
	225
	226	CALL cpu_log( log_point_s(40), 'advec_part_io', 'start' )
	227	CALL check_open( 85 )
	228	WRITE ( 85 ) simulated_time, maximum_number_of_particles, &
	229	number_of_particles
	230	WRITE ( 85 ) particles
	231	WRITE ( 85 ) maximum_number_of_tailpoints, maximum_number_of_tails, &
	232	number_of_tails
	233	WRITE ( 85 ) particle_tail_coordinates
	234	CALL close_file( 85 )
	235
	236	IF ( netcdf_output ) CALL output_particles_netcdf
	237
	238	time_write_particle_data = MOD( time_write_particle_data, &
	239	MAX( dt_write_particle_data, dt_3d ) )
	240	CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' )
	241	ENDIF
	242
	243	!
	244	!-- Calculate exponential term used in case of particle inertia for each
	245	!-- of the particle groups
	246	CALL cpu_log( log_point_s(41), 'advec_part_exp', 'start' )
	247	DO m = 1, number_of_particle_groups
	248	IF ( particle_groups(m)%density_ratio /= 0.0 ) THEN
	249	particle_groups(m)%exp_arg = &
	250	4.5 * particle_groups(m)%density_ratio * &
	251	molecular_viscosity / ( particle_groups(m)%radius )**2
	252	particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
	253	dt_3d )
	254	ENDIF
	255	ENDDO
	256	CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' )
	257
	258	! WRITE ( 9, * ) '*** advec_particles: ##0.3'
	259	! CALL FLUSH_( 9 )
	260	! nd = 0
	261	! DO n = 1, number_of_particles
	262	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	263	! ENDDO
	264	! IF ( nd /= deleted_particles ) THEN
	265	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	266	! CALL FLUSH_( 9 )
	267	! CALL MPI_ABORT( comm2d, 9999, ierr )
	268	! ENDIF
	269
	270	!
	271	!-- Particle (droplet) growth by condensation/evaporation and collision
	272	IF ( cloud_droplets ) THEN
	273
	274	!
	275	!-- Reset summation arrays
	276	ql_c = 0.0; ql_v = 0.0; ql_vp = 0.0
	277
	278	!
	279	!-- Particle growth by condensation/evaporation
	280	CALL cpu_log( log_point_s(42), 'advec_part_cond', 'start' )
	281	DO n = 1, number_of_particles
	282	!
	283	!-- Interpolate temperature and humidity.
	284	!-- First determine left, south, and bottom index of the arrays.
	285	i = particles(n)%x * ddx
	286	j = particles(n)%y * ddy
	287	k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if equidistant
	288
	289	x = particles(n)%x - i * dx
	290	y = particles(n)%y - j * dy
	291	aa = x2 + y2
	292	bb = ( dx - x )2 + y2
	293	cc = x2 + ( dy - y )2
	294	dd = ( dx - x )2 + ( dy - y )2
	295	gg = aa + bb + cc + dd
	296
	297	pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) &
	298	+ ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
	299	) / ( 3.0 * gg )
	300
	301	pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) &
	302	+ ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
	303	) / ( 3.0 * gg )
	304
	305	pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
	306	( pt_int_u - pt_int_l )
	307
	308	q_int_l = ( ( gg - aa ) * q(k,j,i) + ( gg - bb ) * q(k,j,i+1) &
	309	+ ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) &
	310	) / ( 3.0 * gg )
	311
	312	q_int_u = ( ( gg-aa ) * q(k+1,j,i) + ( gg-bb ) * q(k+1,j,i+1) &
	313	+ ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) &
	314	) / ( 3.0 * gg )
	315
	316	q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * &
	317	( q_int_u - q_int_l )
	318
	319	ql_int_l = ( ( gg - aa ) * ql(k,j,i) + ( gg - bb ) * ql(k,j,i+1) &
	320	+ ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) &
	321	) / ( 3.0 * gg )
	322
	323	ql_int_u = ( ( gg-aa ) * ql(k+1,j,i) + ( gg-bb ) * ql(k+1,j,i+1) &
	324	+ ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) &
	325	) / ( 3.0 * gg )
	326
	327	ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * &
	328	( ql_int_u - ql_int_l )
	329
	330	!
	331	!-- Calculate real temperature and saturation vapor pressure
	332	p_int = hydro_press(k) + ( particles(n)%z - zu(k) ) / dz * &
	333	( hydro_press(k+1) - hydro_press(k) )
	334	t_int = pt_int * ( p_int / 100000.0 )**0.286
	335
	336	e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) )
	337
	338	!
	339	!-- Current vapor pressure
	340	e_a = q_int * p_int / ( 0.378 * q_int + 0.622 )
	341
	342	!
	343	!-- Change in radius by condensation/evaporation
	344	!-- ATTENTION: this is only an approximation for large radii
	345	arg = particles(n)%radius*2 + 2.0 dt_3d * &
	346	( e_a / e_s - 1.0 ) / &
	347	( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
	348	thermal_conductivity_l + &
	349	rho_l * r_v * t_int / diff_coeff_l / e_s )
	350	IF ( arg < 1.0E-14 ) THEN
	351	new_r = 1.0E-7
	352	ELSE
	353	new_r = SQRT( arg )
	354	ENDIF
	355
	356	delta_r = new_r - particles(n)%radius
	357
	358	! NOTE: this is the correct formula (indipendent of radius).
	359	! nevertheless, it give wrong results for large timesteps
	360	! d_radius = 1.0 / particles(n)%radius
	361	! delta_r = d_radius * ( e_a / e_s - 1.0 - 3.3E-7 / t_int * d_radius + &
	362	! b_cond * d_radius**3 ) / &
	363	! ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
	364	! thermal_conductivity_l + &
	365	! rho_l * r_v * t_int / diff_coeff_l / e_s ) * dt_3d
	366
	367	! new_r = particles(n)%radius + delta_r
	368	! IF ( new_r < 1.0E-7 ) new_r = 1.0E-7
	369
	370	!
	371	!-- Sum up the change in volume of liquid water for the respective grid
	372	!-- volume (this is needed later on for calculating the release of
	373	!-- latent heat)
	374	i = ( particles(n)%x + 0.5 * dx ) * ddx
	375	j = ( particles(n)%y + 0.5 * dy ) * ddy
	376	k = particles(n)%z / dz + 1 ! only exact if equidistant
	377
	378	ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * &
	379	rho_l * 1.33333333 * pi * &
	380	( new_r3 - particles(n)%radius3 ) / &
	381	( rho_surface * dx * dy * dz )
	382	IF ( ql_c(k,j,i) > 100.0 ) THEN
	383	print*,'+++ advec_particles k=',k,' j=',j,' i=',i, &
	384	' ql_c=',ql_c(k,j,i), ' part(',n,')%wf=', &
	385	particles(n)%weight_factor,' delta_r=',delta_r
	386	#if defined( __parallel )
	387	CALL MPI_ABORT( comm2d, 9999, ierr )
	388	#else
	389	STOP
	390	#endif
	391	ENDIF
	392
	393	!
	394	!-- Change the droplet radius
	395	IF ( ( new_r - particles(n)%radius ) < 0.0 .AND. new_r < 0.0 ) &
	396	THEN
	397	print*,'+++ advec_particles #1 k=',k,' j=',j,' i=',i, &
	398	' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, &
	399	' d_radius=',d_radius,' delta_r=',delta_r,&
	400	' particle_radius=',particles(n)%radius
	401	#if defined( __parallel )
	402	CALL MPI_ABORT( comm2d, 9999, ierr )
	403	#else
	404	STOP
	405	#endif
	406	ENDIF
	407	particles(n)%radius = new_r
	408
	409	!
	410	!-- Sum up the total volume of liquid water (needed below for
	411	!-- re-calculating the weighting factors)
	412	ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
	413	particles(n)%radius**3
	414	ENDDO
	415	CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' )
	416
	417	!
	418	!-- Particle growth by collision
	419	CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
	420
	421	DO i = nxl, nxr
	422	DO j = nys, nyn
	423	DO k = nzb+1, nzt
	424	!
	425	!-- Collision requires at least two particles in the box
	426	IF ( prt_count(k,j,i) > 1 ) THEN
	427	!
	428	!-- First, sort particles within the gridbox by their size,
	429	!-- using Shell's method (see Numerical Recipes)
	430	!-- NOTE: In case of using particle tails, the re-sorting of
	431	!-- ---- tails would have to be included here!
	432	psi = prt_start_index(k,j,i) - 1
	433	inc = 1
	434	DO WHILE ( inc <= prt_count(k,j,i) )
	435	inc = 3 * inc + 1
	436	ENDDO
	437
	438	DO WHILE ( inc > 1 )
	439	inc = inc / 3
	440	DO is = inc+1, prt_count(k,j,i)
	441	tmp_particle = particles(psi+is)
	442	js = is
	443	DO WHILE ( particles(psi+js-inc)%radius > &
	444	tmp_particle%radius )
	445	particles(psi+js) = particles(psi+js-inc)
	446	js = js - inc
	447	IF ( js <= inc ) EXIT
	448	ENDDO
	449	particles(psi+js) = tmp_particle
	450	ENDDO
	451	ENDDO
	452
	453	!
	454	!-- Calculate the mean radius of all those particles which
	455	!-- are of smaller or equal size than the current particle
	456	!-- and use this radius for calculating the collision efficiency
	457	psi = prt_start_index(k,j,i)
	458	s_r3 = 0.0
	459	s_r4 = 0.0
	460	DO n = psi, psi+prt_count(k,j,i)-1
	461	!
	462	!-- There may be some particles of size equal to the
	463	!-- current particle but with larger index
	464	sl_r3 = 0.0
	465	sl_r4 = 0.0
	466	DO is = n, psi+prt_count(k,j,i)-2
	467	IF ( particles(is+1)%radius == &
	468	particles(is)%radius ) THEN
	469	sl_r3 = sl_r3 + particles(is+1)%radius**3
	470	sl_r4 = sl_r4 + particles(is+1)%radius**4
	471	ELSE
	472	EXIT
	473	ENDIF
	474	ENDDO
	475
	476	IF ( ( s_r3 + sl_r3 ) > 0.0 ) THEN
	477
	478	mean_r = ( s_r4 + sl_r4 ) / ( s_r3 + sl_r3 )
	479
	480	CALL collision_efficiency( mean_r, &
	481	particles(n)%radius, &
	482	effective_coll_efficiency )
	483
	484	ELSE
	485	effective_coll_efficiency = 0.0
	486	ENDIF
	487
	488	!
	489	!-- Contribution of the current particle to the next one
	490	s_r3 = s_r3 + particles(n)%radius**3
	491	s_r4 = s_r4 + particles(n)%radius**4
	492
	493	IF ( effective_coll_efficiency > 1.0 .OR. &
	494	effective_coll_efficiency < 0.0 ) &
	495	THEN
	496	print*,'+++ advec_particles collision_efficiency ', &
	497	'out of range:', effective_coll_efficiency
	498	#if defined( __parallel )
	499	CALL MPI_ABORT( comm2d, 9999, ierr )
	500	#else
	501	STOP
	502	#endif
	503	ENDIF
	504
	505	!
	506	!-- Interpolation of ...
	507	ii = particles(n)%x * ddx
	508	jj = particles(n)%y * ddy
	509	kk = ( particles(n)%z + 0.5 * dz ) / dz
	510
	511	x = particles(n)%x - ii * dx
	512	y = particles(n)%y - jj * dy
	513	aa = x2 + y2
	514	bb = ( dx - x )2 + y2
	515	cc = x2 + ( dy - y )2
	516	dd = ( dx - x )2 + ( dy - y )2
	517	gg = aa + bb + cc + dd
	518
	519	ql_int_l = ( ( gg-aa ) * ql(kk,jj,ii) + ( gg-bb ) * &
	520	ql(kk,jj,ii+1) &
	521	+ ( gg-cc ) * ql(kk,jj+1,ii) + ( gg-dd ) * &
	522	ql(kk,jj+1,ii+1) &
	523	) / ( 3.0 * gg )
	524
	525	ql_int_u = ( ( gg-aa ) * ql(kk+1,jj,ii) + ( gg-bb ) * &
	526	ql(kk+1,jj,ii+1) &
	527	+ ( gg-cc ) * ql(kk+1,jj+1,ii) + ( gg-dd ) * &
	528	ql(kk+1,jj+1,ii+1) &
	529	) / ( 3.0 * gg )
	530
	531	ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz * &
	532	( ql_int_u - ql_int_l )
	533
	534	!
	535	!-- Interpolate u velocity-component
	536	ii = ( particles(n)%x + 0.5 * dx ) * ddx
	537	jj = particles(n)%y * ddy
	538	kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eq.dist
	539
	540	IF ( ( particles(n)%z - zu(kk) ) > ( 0.5*dz ) ) kk = kk+1
	541
	542	x = particles(n)%x + ( 0.5 - ii ) * dx
	543	y = particles(n)%y - jj * dy
	544	aa = x2 + y2
	545	bb = ( dx - x )2 + y2
	546	cc = x2 + ( dy - y )2
	547	dd = ( dx - x )2 + ( dy - y )2
	548	gg = aa + bb + cc + dd
	549
	550	u_int_l = ( ( gg-aa ) * u(kk,jj,ii) + ( gg-bb ) * &
	551	u(kk,jj,ii+1) &
	552	+ ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) * &
	553	u(kk,jj+1,ii+1) &
	554	) / ( 3.0 * gg ) - u_gtrans
	555	IF ( kk+1 == nzt+1 ) THEN
	556	u_int = u_int_l
	557	ELSE
	558	u_int_u = ( ( gg-aa ) * u(kk+1,jj,ii) + ( gg-bb ) * &
	559	u(kk+1,jj,ii+1) &
	560	+ ( gg-cc ) * u(kk+1,jj+1,ii) + ( gg-dd ) * &
	561	u(kk+1,jj+1,ii+1) &
	562	) / ( 3.0 * gg ) - u_gtrans
	563	u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz * &
	564	( u_int_u - u_int_l )
	565	ENDIF
	566
	567	!
	568	!-- Same procedure for interpolation of the v velocity-compo-
	569	!-- nent (adopt index k from u velocity-component)
	570	ii = particles(n)%x * ddx
	571	jj = ( particles(n)%y + 0.5 * dy ) * ddy
	572
	573	x = particles(n)%x - ii * dx
	574	y = particles(n)%y + ( 0.5 - jj ) * dy
	575	aa = x2 + y2
	576	bb = ( dx - x )2 + y2
	577	cc = x2 + ( dy - y )2
	578	dd = ( dx - x )2 + ( dy - y )2
	579	gg = aa + bb + cc + dd
	580
	581	v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * &
	582	v(kk,jj,ii+1) &
	583	+ ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * &
	584	v(kk,jj+1,ii+1) &
	585	) / ( 3.0 * gg ) - v_gtrans
	586	IF ( kk+1 == nzt+1 ) THEN
	587	v_int = v_int_l
	588	ELSE
	589	v_int_u = ( ( gg-aa ) * v(kk+1,jj,ii) + ( gg-bb ) * &
	590	v(kk+1,jj,ii+1) &
	591	+ ( gg-cc ) * v(kk+1,jj+1,ii) + ( gg-dd ) * &
	592	v(kk+1,jj+1,ii+1) &
	593	) / ( 3.0 * gg ) - v_gtrans
	594	v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz * &
	595	( v_int_u - v_int_l )
	596	ENDIF
	597
	598	!
	599	!-- Same procedure for interpolation of the w velocity-compo-
	600	!-- nent (adopt index i from v velocity-component)
	601	jj = particles(n)%y * ddy
	602	kk = particles(n)%z / dz
	603
	604	x = particles(n)%x - ii * dx
	605	y = particles(n)%y - jj * dy
	606	aa = x2 + y2
	607	bb = ( dx - x )2 + y2
	608	cc = x2 + ( dy - y )2
	609	dd = ( dx - x )2 + ( dy - y )2
	610	gg = aa + bb + cc + dd
	611
	612	w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * &
	613	w(kk,jj,ii+1) &
	614	+ ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * &
	615	w(kk,jj+1,ii+1) &
	616	) / ( 3.0 * gg )
	617	IF ( kk+1 == nzt+1 ) THEN
	618	w_int = w_int_l
	619	ELSE
	620	w_int_u = ( ( gg-aa ) * w(kk+1,jj,ii) + ( gg-bb ) * &
	621	w(kk+1,jj,ii+1) &
	622	+ ( gg-cc ) * w(kk+1,jj+1,ii) + ( gg-dd ) * &
	623	w(kk+1,jj+1,ii+1) &
	624	) / ( 3.0 * gg )
	625	w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz * &
	626	( w_int_u - w_int_l )
	627	ENDIF
	628
	629	!
	630	!-- Change in radius due to collision
	631	delta_r = effective_coll_efficiency * &
	632	ql_int * rho_surface / ( 1.0 - ql_int ) * &
	633	0.25 / rho_l * &
	634	SQRT( ( u_int - particles(n)%speed_x )**2 + &
	635	( v_int - particles(n)%speed_y )**2 + &
	636	( w_int - particles(n)%speed_z )**2 &
	637	) * dt_3d
	638
	639	particles(n)%radius = particles(n)%radius + delta_r
	640
	641	ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3
	642
	643	ENDDO
	644
	645	ENDIF
	646
	647	!
	648	!-- Re-calculate the weighting factor (total liquid water content
	649	!-- must be conserved during collision)
	650	IF ( ql_vp(k,j,i) /= 0.0 ) THEN
	651
	652	ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
	653	!
	654	!-- Re-assign this weighting factor to the particles of the
	655	!-- current gridbox
	656	psi = prt_start_index(k,j,i)
	657	DO n = psi, psi + prt_count(k,j,i)-1
	658	particles(n)%weight_factor = ql_vp(k,j,i)
	659	ENDDO
	660
	661	ENDIF
	662
	663	ENDDO
	664	ENDDO
	665	ENDDO
	666
	667	CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
	668
	669	ENDIF
	670
	671
	672	!
	673	!-- Particle advection.
	674	!-- In case of including the SGS velocities, the LES timestep has probably
	675	!-- to be split into several smaller timesteps because of the Lagrangian
	676	!-- timescale condition. Because the number of timesteps to be carried out is
	677	!-- not known at the beginning, these steps are carried out in an infinite loop
	678	!-- with exit condition.
	679	!
	680	!-- If SGS velocities are used, gradients of the TKE have to be calculated and
	681	!-- boundary conditions have to be set first. Also, horizontally averaged
	682	!-- profiles of the SGS TKE and the resolved-scale velocity variances are
	683	!-- needed.
	684	IF ( use_sgs_for_particles ) THEN
	685
	686	!
	687	!-- TKE gradient along x and y
	688	DO i = nxl, nxr
	689	DO j = nys, nyn
	690	DO k = nzb, nzt+1
	691	de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
	692	de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
	693	ENDDO
	694	ENDDO
	695	ENDDO
	696
	697	!
	698	!-- TKE gradient along z, including bottom and top boundary conditions
	699	DO i = nxl, nxr
	700	DO j = nys, nyn
	701	DO k = nzb+2, nzt-1
	702	de_dz(k,j,i) = 2.0 * sgs_wfw_part * &
	703	( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
	704	ENDDO
	705	de_dz(nzb,j,i) = 0.0
	706	de_dz(nzb+1,j,i) = 2.0 * sgs_wfw_part * &
	707	( e(nzb+2,j,i) - e(nzb+1,j,i) ) / &
	708	( zu(nzb+2) - zu(nzb+1) )
	709	de_dz(nzt,j,i) = 0.0
	710	de_dz(nzt+1,j,i) = 0.0
	711	ENDDO
	712	ENDDO
	713
	714	!
	715	!-- Lateral boundary conditions
	716	CALL exchange_horiz( de_dx, 0, 0 )
	717	CALL exchange_horiz( de_dy, 0, 0 )
	718	CALL exchange_horiz( de_dz, 0, 0 )
	719	CALL exchange_horiz( diss, 0, 0 )
	720
	721	!
	722	!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
	723	!-- velocity variances (they may have been already calculated in routine
	724	!-- flow_statistics).
	725	IF ( .NOT. flow_statistics_called ) THEN
	726	!
	727	!-- First calculate horizontally averaged profiles of the horizontal
	728	!-- velocities.
	729	sums_l(:,1,0) = 0.0
	730	sums_l(:,2,0) = 0.0
	731
	732	DO i = nxl, nxr
	733	DO j = nys, nyn
	734	DO k = nzb_s_outer(j,i), nzt+1
	735	sums_l(k,1,0) = sums_l(k,1,0) + u(k,j,i)
	736	sums_l(k,2,0) = sums_l(k,2,0) + v(k,j,i)
	737	ENDDO
	738	ENDDO
	739	ENDDO
	740
	741	#if defined( __parallel )
	742	!
	743	!-- Compute total sum from local sums
	744	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
	745	MPI_REAL, MPI_SUM, comm2d, ierr )
	746	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
	747	MPI_REAL, MPI_SUM, comm2d, ierr )
	748	#else
	749	sums(:,1) = sums_l(:,1,0)
	750	sums(:,2) = sums_l(:,2,0)
	751	#endif
	752
	753	!
	754	!-- Final values are obtained by division by the total number of grid
	755	!-- points used for the summation.
	756	hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0) ! u
	757	hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0) ! v
	758
	759	!
	760	!-- Now calculate the profiles of SGS TKE and the resolved-scale
	761	!-- velocity variances
	762	sums_l(:,8,0) = 0.0
	763	sums_l(:,30,0) = 0.0
	764	sums_l(:,31,0) = 0.0
	765	sums_l(:,32,0) = 0.0
	766	DO i = nxl, nxr
	767	DO j = nys, nyn
	768	DO k = nzb_s_outer(j,i), nzt+1
	769	sums_l(k,8,0) = sums_l(k,8,0) + e(k,j,i)
	770	sums_l(k,30,0) = sums_l(k,30,0) + &
	771	( u(k,j,i) - hom(k,1,1,0) )**2
	772	sums_l(k,31,0) = sums_l(k,31,0) + &
	773	( v(k,j,i) - hom(k,1,2,0) )**2
	774	sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2
	775	ENDDO
	776	ENDDO
	777	ENDDO
	778
	779	#if defined( __parallel )
	780	!
	781	!-- Compute total sum from local sums
	782	CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
	783	MPI_REAL, MPI_SUM, comm2d, ierr )
	784	CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
	785	MPI_REAL, MPI_SUM, comm2d, ierr )
	786	CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
	787	MPI_REAL, MPI_SUM, comm2d, ierr )
	788	CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
	789	MPI_REAL, MPI_SUM, comm2d, ierr )
	790
	791	#else
	792	sums(:,8) = sums_l(:,8,0)
	793	sums(:,30) = sums_l(:,30,0)
	794	sums(:,31) = sums_l(:,31,0)
	795	sums(:,32) = sums_l(:,32,0)
	796	#endif
	797
	798	!
	799	!-- Final values are obtained by division by the total number of grid
	800	!-- points used for the summation.
	801	hom(:,1,8,0) = sums(:,8) / ngp_2dh_outer(:,0) ! e
	802	hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0) ! u*2
	803	hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0) ! v*2
	804	hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0) ! w*2
	805
	806	ENDIF
	807
	808	ENDIF
	809
	810	!
	811	!-- Initialize variables used for accumulating the number of particles
	812	!-- exchanged between the subdomains during all sub-timesteps (if sgs
	813	!-- velocities are included). These data are output further below on the
	814	!-- particle statistics file.
	815	trlp_count_sum = 0
	816	trlp_count_recv_sum = 0
	817	trrp_count_sum = 0
	818	trrp_count_recv_sum = 0
	819	trsp_count_sum = 0
	820	trsp_count_recv_sum = 0
	821	trnp_count_sum = 0
	822	trnp_count_recv_sum = 0
	823
	824	!
	825	!-- Initialize the variable storing the total time that a particle has advanced
	826	!-- within the timestep procedure
	827	particles(1:number_of_particles)%dt_sum = 0.0
	828
	829	!
	830	!-- Timestep loop.
	831	!-- This loop has to be repeated until the advection time of every particle
	832	!-- (in the total domain!) has reached the LES timestep (dt_3d)
	833	DO
	834
	835	CALL cpu_log( log_point_s(44), 'advec_part_advec', 'start' )
	836
	837	!
	838	!-- Initialize the switch used for the loop exit condition checked at the
	839	!-- end of this loop.
	840	!-- If at least one particle has failed to reach the LES timestep, this
	841	!-- switch will be set false.
	842	dt_3d_reached_l = .TRUE.
	843
	844	!
	845	!-- Initialize variables for the (sub-) timestep, i.e. for marking those
	846	!-- particles to be deleted after the timestep
	847	particle_mask = .TRUE.
	848	deleted_particles = 0
	849	trlp_count_recv = 0
	850	trnp_count_recv = 0
	851	trrp_count_recv = 0
	852	trsp_count_recv = 0
	853	IF ( use_particle_tails ) THEN
	854	tail_mask = .TRUE.
	855	deleted_tails = 0
	856	ENDIF
	857
	858
	859	DO n = 1, number_of_particles
	860	!
	861	!-- Move particles only if the LES timestep has not (approximately) been
	862	!-- reached
	863	IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 ) CYCLE
	864
	865	!
	866	!-- Interpolate u velocity-component, determine left, front, bottom
	867	!-- index of u-array
	868	i = ( particles(n)%x + 0.5 * dx ) * ddx
	869	j = particles(n)%y * ddy
	870	k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if equidistant
	871
	872	!
	873	!-- Interpolation of the velocity components in the xy-plane
	874	x = particles(n)%x + ( 0.5 - i ) * dx
	875	y = particles(n)%y - j * dy
	876	aa = x2 + y2
	877	bb = ( dx - x )2 + y2
	878	cc = x2 + ( dy - y )2
	879	dd = ( dx - x )2 + ( dy - y )2
	880	gg = aa + bb + cc + dd
	881
	882	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
	883	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
	884	) / ( 3.0 * gg ) - u_gtrans
	885	IF ( k+1 == nzt+1 ) THEN
	886	u_int = u_int_l
	887	ELSE
	888	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
	889	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
	890	) / ( 3.0 * gg ) - u_gtrans
	891	u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
	892	( u_int_u - u_int_l )
	893	ENDIF
	894
	895	!
	896	!-- Same procedure for interpolation of the v velocity-component (adopt
	897	!-- index k from u velocity-component)
	898	i = particles(n)%x * ddx
	899	j = ( particles(n)%y + 0.5 * dy ) * ddy
	900
	901	x = particles(n)%x - i * dx
	902	y = particles(n)%y + ( 0.5 - j ) * dy
	903	aa = x2 + y2
	904	bb = ( dx - x )2 + y2
	905	cc = x2 + ( dy - y )2
	906	dd = ( dx - x )2 + ( dy - y )2
	907	gg = aa + bb + cc + dd
	908
	909	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
	910	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
	911	) / ( 3.0 * gg ) - v_gtrans
	912	IF ( k+1 == nzt+1 ) THEN
	913	v_int = v_int_l
	914	ELSE
	915	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
	916	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
	917	) / ( 3.0 * gg ) - v_gtrans
	918	v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
	919	( v_int_u - v_int_l )
	920	ENDIF
	921
	922	!
	923	!-- Same procedure for interpolation of the w velocity-component (adopt
	924	!-- index i from v velocity-component)
	925	IF ( vertical_particle_advection ) THEN
	926	j = particles(n)%y * ddy
	927	k = particles(n)%z / dz
	928
	929	x = particles(n)%x - i * dx
	930	y = particles(n)%y - j * dy
	931	aa = x2 + y2
	932	bb = ( dx - x )2 + y2
	933	cc = x2 + ( dy - y )2
	934	dd = ( dx - x )2 + ( dy - y )2
	935	gg = aa + bb + cc + dd
	936
	937	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
	938	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
	939	) / ( 3.0 * gg )
	940	IF ( k+1 == nzt+1 ) THEN
	941	w_int = w_int_l
	942	ELSE
	943	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
	944	( gg-bb ) * w(k+1,j,i+1) + &
	945	( gg-cc ) * w(k+1,j+1,i) + &
	946	( gg-dd ) * w(k+1,j+1,i+1) &
	947	) / ( 3.0 * gg )
	948	w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
	949	( w_int_u - w_int_l )
	950	ENDIF
	951	ELSE
	952	w_int = 0.0
	953	ENDIF
	954
	955	!
	956	!-- Interpolate and calculate quantities needed for calculating the SGS
	957	!-- velocities
	958	IF ( use_sgs_for_particles ) THEN
	959	!
	960	!-- Interpolate TKE
	961	i = particles(n)%x * ddx
	962	j = particles(n)%y * ddy
	963	k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if eq.dist
	964
	965	x = particles(n)%x - i * dx
	966	y = particles(n)%y - j * dy
	967	aa = x2 + y2
	968	bb = ( dx - x )2 + y2
	969	cc = x2 + ( dy - y )2
	970	dd = ( dx - x )2 + ( dy - y )2
	971	gg = aa + bb + cc + dd
	972
	973	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
	974	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
	975	) / ( 3.0 * gg )
	976
	977	IF ( k+1 == nzt+1 ) THEN
	978	e_int = e_int_l
	979	ELSE
	980	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	981	( gg - bb ) * e(k+1,j,i+1) + &
	982	( gg - cc ) * e(k+1,j+1,i) + &
	983	( gg - dd ) * e(k+1,j+1,i+1) &
	984	) / ( 3.0 * gg )
	985	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
	986	( e_int_u - e_int_l )
	987	ENDIF
	988
	989	!
	990	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and all
	991	!-- position variables from above (TKE))
	992	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	993	( gg - bb ) * de_dx(k,j,i+1) + &
	994	( gg - cc ) * de_dx(k,j+1,i) + &
	995	( gg - dd ) * de_dx(k,j+1,i+1) &
	996	) / ( 3.0 * gg )
	997
	998	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	999	de_dx_int = de_dx_int_l
	1000	ELSE
	1001	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	1002	( gg - bb ) * de_dx(k+1,j,i+1) + &
	1003	( gg - cc ) * de_dx(k+1,j+1,i) + &
	1004	( gg - dd ) * de_dx(k+1,j+1,i+1) &
	1005	) / ( 3.0 * gg )
	1006	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1007	( de_dx_int_u - de_dx_int_l )
	1008	ENDIF
	1009
	1010	!
	1011	!-- Interpolate the TKE gradient along y
	1012	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	1013	( gg - bb ) * de_dy(k,j,i+1) + &
	1014	( gg - cc ) * de_dy(k,j+1,i) + &
	1015	( gg - dd ) * de_dy(k,j+1,i+1) &
	1016	) / ( 3.0 * gg )
	1017	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1018	de_dy_int = de_dy_int_l
	1019	ELSE
	1020	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	1021	( gg - bb ) * de_dy(k+1,j,i+1) + &
	1022	( gg - cc ) * de_dy(k+1,j+1,i) + &
	1023	( gg - dd ) * de_dy(k+1,j+1,i+1) &
	1024	) / ( 3.0 * gg )
	1025	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1026	( de_dy_int_u - de_dy_int_l )
	1027	ENDIF
	1028
	1029	!
	1030	!-- Interpolate the TKE gradient along z
	1031	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	1032	( gg - bb ) * de_dz(k,j,i+1) + &
	1033	( gg - cc ) * de_dz(k,j+1,i) + &
	1034	( gg - dd ) * de_dz(k,j+1,i+1) &
	1035	) / ( 3.0 * gg )
	1036
	1037	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1038	de_dz_int = de_dz_int_l
	1039	ELSE
	1040	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	1041	( gg - bb ) * de_dz(k+1,j,i+1) + &
	1042	( gg - cc ) * de_dz(k+1,j+1,i) + &
	1043	( gg - dd ) * de_dz(k+1,j+1,i+1) &
	1044	) / ( 3.0 * gg )
	1045	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1046	( de_dz_int_u - de_dz_int_l )
	1047	ENDIF
	1048
	1049	!
	1050	!-- Interpolate the dissipation of TKE
	1051	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	1052	( gg - bb ) * diss(k,j,i+1) + &
	1053	( gg - cc ) * diss(k,j+1,i) + &
	1054	( gg - dd ) * diss(k,j+1,i+1) &
	1055	) / ( 3.0 * gg )
	1056
	1057	IF ( k+1 == nzt+1 ) THEN
	1058	diss_int = diss_int_l
	1059	ELSE
	1060	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	1061	( gg - bb ) * diss(k+1,j,i+1) + &
	1062	( gg - cc ) * diss(k+1,j+1,i) + &
	1063	( gg - dd ) * diss(k+1,j+1,i+1) &
	1064	) / ( 3.0 * gg )
	1065	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1066	( diss_int_u - diss_int_l )
	1067	ENDIF
	1068
	1069	!
	1070	!-- Vertically interpolate the horizontally averaged SGS TKE and
	1071	!-- resolved-scale velocity variances and use the interpolated values
	1072	!-- to calculate the coefficient fs, which is a measure of the ratio
	1073	!-- of the subgrid-scale turbulent kinetic energy to the total amount
	1074	!-- of turbulent kinetic energy.
	1075	IF ( k == 0 ) THEN
	1076	e_mean_int = hom(0,1,8,0)
	1077	ELSE
	1078	e_mean_int = hom(k,1,8,0) + &
	1079	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
	1080	( zu(k+1) - zu(k) ) * &
	1081	( particles(n)%z - zu(k) )
	1082	ENDIF
	1083
	1084	kw = particles(n)%z / dz
	1085
	1086	IF ( k == 0 ) THEN
	1087	aa = hom(k+1,1,30,0) * ( particles(n)%z / &
	1088	( 0.5 * ( zu(k+1) - zu(k) ) ) )
	1089	bb = hom(k+1,1,31,0) * ( particles(n)%z / &
	1090	( 0.5 * ( zu(k+1) - zu(k) ) ) )
	1091	cc = hom(kw+1,1,32,0) * ( particles(n)%z / &
	1092	( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
	1093	ELSE
	1094	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
	1095	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
	1096	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
	1097	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
	1098	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
	1099	( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
	1100	ENDIF
	1101
	1102	vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
	1103
	1104	fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
	1105	( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
	1106
	1107	!
	1108	!-- Calculate the Lagrangian timescale according to the suggestion of
	1109	!-- Weil et al. (2004).
	1110	lagr_timescale = ( 4.0 * e_int ) / &
	1111	( 3.0 * fs_int * c_0 * diss_int )
	1112
	1113	!
	1114	!-- Calculate the next particle timestep. dt_gap is the time needed to
	1115	!-- complete the current LES timestep.
	1116	dt_gap = dt_3d - particles(n)%dt_sum
	1117	dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
	1118
	1119	!
	1120	!-- The particle timestep should not be too small in order to prevent
	1121	!-- the number of particle timesteps of getting too large
	1122	IF ( dt_particle < dt_min_part .AND. dt_min_part < dt_gap ) &
	1123	THEN
	1124	dt_particle = dt_min_part
	1125	ENDIF
	1126
	1127	!
	1128	!-- Calculate the SGS velocity components
	1129	IF ( particles(n)%age == 0.0 ) THEN
	1130	!
	1131	!-- For new particles the SGS components are derived from the SGS
	1132	!-- TKE. Limit the Gaussian random number to the interval
	1133	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
	1134	!-- from becoming unrealistically large.
	1135	particles(n)%speed_x_sgs = SQRT( 2.0 * sgs_wfu_part * e_int ) *&
	1136	( random_gauss( iran_part, 5.0 ) &
	1137	- 1.0 )
	1138	particles(n)%speed_y_sgs = SQRT( 2.0 * sgs_wfv_part * e_int ) *&
	1139	( random_gauss( iran_part, 5.0 ) &
	1140	- 1.0 )
	1141	particles(n)%speed_z_sgs = SQRT( 2.0 * sgs_wfw_part * e_int ) *&
	1142	( random_gauss( iran_part, 5.0 ) &
	1143	- 1.0 )
	1144
	1145	ELSE
	1146	!
	1147	!-- For old particles the SGS components are correlated with the
	1148	!-- values from the previous timestep. Random numbers have also to
	1149	!-- be limited (see above).
	1150	particles(n)%speed_x_sgs = particles(n)%speed_x_sgs - &
	1151	fs_int * c_0 * diss_int * particles(n)%speed_x_sgs * &
	1152	dt_particle / ( 4.0 * sgs_wfu_part * e_int ) + &
	1153	( ( 2.0 * sgs_wfu_part * ( e_int - particles(n)%e_m ) / &
	1154	dt_particle ) * particles(n)%speed_x_sgs / &
	1155	( 2.0 * sgs_wfu_part * e_int ) + de_dx_int &
	1156	) * dt_particle / 2.0 + &
	1157	SQRT( fs_int * c_0 * diss_int ) * &
	1158	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	1159	SQRT( dt_particle )
	1160
	1161	particles(n)%speed_y_sgs = particles(n)%speed_y_sgs - &
	1162	fs_int * c_0 * diss_int * particles(n)%speed_y_sgs * &
	1163	dt_particle / ( 4.0 * sgs_wfv_part * e_int ) + &
	1164	( ( 2.0 * sgs_wfv_part * ( e_int - particles(n)%e_m ) / &
	1165	dt_particle ) * particles(n)%speed_y_sgs / &
	1166	( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
	1167	) * dt_particle / 2.0 + &
	1168	SQRT( fs_int * c_0 * diss_int ) * &
	1169	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	1170	SQRT( dt_particle )
	1171
	1172	particles(n)%speed_z_sgs = particles(n)%speed_z_sgs - &
	1173	fs_int * c_0 * diss_int * particles(n)%speed_z_sgs * &
	1174	dt_particle / ( 4.0 * sgs_wfw_part * e_int ) + &
	1175	( ( 2.0 * sgs_wfw_part * ( e_int - particles(n)%e_m ) / &
	1176	dt_particle ) * particles(n)%speed_z_sgs / &
	1177	( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
	1178	) * dt_particle / 2.0 + &
	1179	SQRT( fs_int * c_0 * diss_int ) * &
	1180	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	1181	SQRT( dt_particle )
	1182
	1183	ENDIF
	1184
	1185	u_int = u_int + particles(n)%speed_x_sgs
	1186	v_int = v_int + particles(n)%speed_y_sgs
	1187	w_int = w_int + particles(n)%speed_z_sgs
	1188
	1189	!
	1190	!-- Store the SGS TKE of the current timelevel which is needed for
	1191	!-- for calculating the SGS particle velocities at the next timestep
	1192	particles(n)%e_m = e_int
	1193
	1194	ELSE
	1195	!
	1196	!-- If no SGS velocities are used, only the particle timestep has to
	1197	!-- be set
	1198	dt_particle = dt_3d
	1199
	1200	ENDIF
	1201
	1202	!
	1203	!-- Particle advection
	1204	IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 ) THEN
	1205	!
	1206	!-- Pure passive transport (without particle inertia)
	1207	particles(n)%x = particles(n)%x + u_int * dt_particle
	1208	particles(n)%y = particles(n)%y + v_int * dt_particle
	1209	particles(n)%z = particles(n)%z + w_int * dt_particle
	1210
	1211	particles(n)%speed_x = u_int
	1212	particles(n)%speed_y = v_int
	1213	particles(n)%speed_z = w_int
	1214
	1215	ELSE
	1216	!
	1217	!-- Transport of particles with inertia
	1218	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
	1219	dt_particle
	1220	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
	1221	dt_particle
	1222	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
	1223	dt_particle
	1224
	1225	!
	1226	!-- Update of the particle velocity
	1227	dens_ratio = particle_groups(particles(n)%group)%density_ratio
	1228	IF ( cloud_droplets ) THEN
	1229	exp_arg = 4.5 * dens_ratio * molecular_viscosity / &
	1230	( particles(n)%radius )**2 / &
	1231	( 1.0 + 0.15 * ( 2.0 * particles(n)%radius * &
	1232	SQRT( ( u_int - particles(n)%speed_x )**2 + &
	1233	( v_int - particles(n)%speed_y )**2 + &
	1234	( w_int - particles(n)%speed_z )**2 ) / &
	1235	molecular_viscosity )**0.687 &
	1236	)
	1237	exp_term = EXP( -exp_arg * dt_particle )
	1238	ELSEIF ( use_sgs_for_particles ) THEN
	1239	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1240	exp_term = EXP( -exp_arg * dt_particle )
	1241	ELSE
	1242	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1243	exp_term = particle_groups(particles(n)%group)%exp_term
	1244	ENDIF
	1245	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
	1246	u_int * ( 1.0 - exp_term )
	1247	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
	1248	v_int * ( 1.0 - exp_term )
	1249	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
	1250	( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
	1251	* ( 1.0 - exp_term )
	1252	ENDIF
	1253
	1254	!
	1255	!-- Increment the particle age and the total time that the particle
	1256	!-- has advanced within the particle timestep procedure
	1257	particles(n)%age = particles(n)%age + dt_particle
	1258	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
	1259
	1260	!
	1261	!-- Check whether there is still a particle that has not yet completed
	1262	!-- the total LES timestep
	1263	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 ) THEN
	1264	dt_3d_reached_l = .FALSE.
	1265	ENDIF
	1266
	1267	ENDDO ! advection loop
	1268
	1269	!
	1270	!-- Find out, if all particles on every PE have completed the LES timestep
	1271	!-- and set the switch corespondingly
	1272	#if defined( __parallel )
	1273	CALL MPI_ALLREDUCE( dt_3d_reached_l, dt_3d_reached, 1, MPI_LOGICAL, &
	1274	MPI_LAND, comm2d, ierr )
	1275	#else
	1276	dt_3d_reached = dt_3d_reached_l
	1277	#endif
	1278
	1279	CALL cpu_log( log_point_s(44), 'advec_part_advec', 'stop' )
	1280
	1281	!
	1282	!-- Increment time since last release
	1283	IF ( dt_3d_reached ) time_prel = time_prel + dt_3d
	1284
	1285	! WRITE ( 9, * ) '*** advec_particles: ##0.4'
	1286	! CALL FLUSH_( 9 )
	1287	! nd = 0
	1288	! DO n = 1, number_of_particles
	1289	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	1290	! ENDDO
	1291	! IF ( nd /= deleted_particles ) THEN
	1292	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	1293	! CALL FLUSH_( 9 )
	1294	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1295	! ENDIF
	1296
	1297	!
	1298	!-- If necessary, release new set of particles
	1299	IF ( time_prel >= dt_prel .AND. end_time_prel > simulated_time .AND. &
	1300	dt_3d_reached ) THEN
	1301
	1302	!
	1303	!-- Check, if particle storage must be extended
	1304	IF ( number_of_particles + number_of_initial_particles > &
	1305	maximum_number_of_particles ) THEN
	1306	IF ( netcdf_output ) THEN
	1307	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
	1308	'needs to be increased'
	1309	PRINT*, ' but this is not allowed with ', &
	1310	'NetCDF output switched on'
	1311	#if defined( __parallel )
	1312	CALL MPI_ABORT( comm2d, 9999, ierr )
	1313	#else
	1314	CALL local_stop
	1315	#endif
	1316	ELSE
	1317	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory dt_prel'
	1318	! CALL FLUSH_( 9 )
	1319	CALL allocate_prt_memory( number_of_initial_particles )
	1320	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory dt_prel'
	1321	! CALL FLUSH_( 9 )
	1322	ENDIF
	1323	ENDIF
	1324
	1325	!
	1326	!-- Check, if tail storage must be extended
	1327	IF ( use_particle_tails ) THEN
	1328	IF ( number_of_tails + number_of_initial_tails > &
	1329	maximum_number_of_tails ) THEN
	1330	IF ( netcdf_output ) THEN
	1331	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
	1332	'needs to be increased'
	1333	PRINT*, ' but this is not allowed wi', &
	1334	'th NetCDF output switched on'
	1335	#if defined( __parallel )
	1336	CALL MPI_ABORT( comm2d, 9999, ierr )
	1337	#else
	1338	CALL local_stop
	1339	#endif
	1340	ELSE
	1341	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory dt_prel'
	1342	! CALL FLUSH_( 9 )
	1343	CALL allocate_tail_memory( number_of_initial_tails )
	1344	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory dt_prel'
	1345	! CALL FLUSH_( 9 )
	1346	ENDIF
	1347	ENDIF
	1348	ENDIF
	1349
	1350	!
	1351	!-- The MOD function allows for changes in the output interval with
	1352	!-- restart runs.
	1353	time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
	1354	IF ( number_of_initial_particles /= 0 ) THEN
	1355	is = number_of_particles+1
	1356	ie = number_of_particles+number_of_initial_particles
	1357	particles(is:ie) = initial_particles(1:number_of_initial_particles)
	1358	!
	1359	!-- Add random fluctuation to particle positions. Particles should
	1360	!-- remain in the subdomain.
	1361	IF ( random_start_position ) THEN
	1362	DO n = is, ie
	1363	IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) &
	1364	THEN
	1365	particles(n)%x = particles(n)%x + &
	1366	( random_function( iran ) - 0.5 ) * &
	1367	pdx(particles(n)%group)
	1368	IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN
	1369	particles(n)%x = ( nxl - 0.4999999999 ) * dx
	1370	ELSEIF ( particles(n)%x >= ( nxr + 0.5 ) * dx ) THEN
	1371	particles(n)%x = ( nxr + 0.4999999999 ) * dx
	1372	ENDIF
	1373	ENDIF
	1374	IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) &
	1375	THEN
	1376	particles(n)%y = particles(n)%y + &
	1377	( random_function( iran ) - 0.5 ) * &
	1378	pdy(particles(n)%group)
	1379	IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN
	1380	particles(n)%y = ( nys - 0.4999999999 ) * dy
	1381	ELSEIF ( particles(n)%y >= ( nyn + 0.5 ) * dy ) THEN
	1382	particles(n)%y = ( nyn + 0.4999999999 ) * dy
	1383	ENDIF
	1384	ENDIF
	1385	IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) &
	1386	THEN
	1387	particles(n)%z = particles(n)%z + &
	1388	( random_function( iran ) - 0.5 ) * &
	1389	pdz(particles(n)%group)
	1390	ENDIF
	1391	ENDDO
	1392	ENDIF
	1393
	1394	!
	1395	!-- Set the beginning of the new particle tails and their age
	1396	IF ( use_particle_tails ) THEN
	1397	DO n = is, ie
	1398	!
	1399	!-- New particles which should have a tail, already have got a
	1400	!-- provisional tail id unequal zero (see init_particles)
	1401	IF ( particles(n)%tail_id /= 0 ) THEN
	1402	number_of_tails = number_of_tails + 1
	1403	nn = number_of_tails
	1404	particles(n)%tail_id = nn ! set the final tail id
	1405	particle_tail_coordinates(1,1,nn) = particles(n)%x
	1406	particle_tail_coordinates(1,2,nn) = particles(n)%y
	1407	particle_tail_coordinates(1,3,nn) = particles(n)%z
	1408	particle_tail_coordinates(1,4,nn) = particles(n)%color
	1409	particles(n)%tailpoints = 1
	1410	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
	1411	particle_tail_coordinates(2,1,nn) = particles(n)%x
	1412	particle_tail_coordinates(2,2,nn) = particles(n)%y
	1413	particle_tail_coordinates(2,3,nn) = particles(n)%z
	1414	particle_tail_coordinates(2,4,nn) = particles(n)%color
	1415	particle_tail_coordinates(1:2,5,nn) = 0.0
	1416	particles(n)%tailpoints = 2
	1417	ENDIF
	1418	ENDIF
	1419	ENDDO
	1420	ENDIF
	1421	! WRITE ( 9, * ) '*** advec_particles: after setting the beginning of new tails'
	1422	! CALL FLUSH_( 9 )
	1423
	1424	number_of_particles = number_of_particles + &
	1425	number_of_initial_particles
	1426	ENDIF
	1427
	1428	ENDIF
	1429
	1430	! WRITE ( 9, * ) '*** advec_particles: ##0.5'
	1431	! CALL FLUSH_( 9 )
	1432	! nd = 0
	1433	! DO n = 1, number_of_particles
	1434	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	1435	! ENDDO
	1436	! IF ( nd /= deleted_particles ) THEN
	1437	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	1438	! CALL FLUSH_( 9 )
	1439	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1440	! ENDIF
	1441
	1442	! IF ( number_of_particles /= number_of_tails ) THEN
	1443	! WRITE (9,*) '--- advec_particles: #2'
	1444	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
	1445	! CALL FLUSH_( 9 )
	1446	! ENDIF
	1447	! DO n = 1, number_of_particles
	1448	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	1449	! THEN
	1450	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	1451	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	1452	! CALL FLUSH_( 9 )
	1453	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1454	! ENDIF
	1455	! ENDDO
	1456
	1457	#if defined( __parallel )
	1458	!
	1459	!-- As soon as one particle has moved beyond the boundary of the domain, it
	1460	!-- is included in the relevant transfer arrays and marked for subsequent
	1461	!-- deletion on this PE.
	1462	!-- First run for crossings in x direction. Find out first the number of
	1463	!-- particles to be transferred and allocate temporary arrays needed to store
	1464	!-- them.
	1465	!-- For a one-dimensional decomposition along y, no transfer is necessary,
	1466	!-- because the particle remains on the PE.
	1467	trlp_count = 0
	1468	trlpt_count = 0
	1469	trrp_count = 0
	1470	trrpt_count = 0
	1471	IF ( pdims(1) /= 1 ) THEN
	1472	!
	1473	!-- First calculate the storage necessary for sending and receiving the
	1474	!-- data
	1475	DO n = 1, number_of_particles
	1476	i = ( particles(n)%x + 0.5 * dx ) * ddx
	1477	!
	1478	!-- Above calculation does not work for indices less than zero
	1479	IF ( particles(n)%x < -0.5 * dx ) i = -1
	1480
	1481	IF ( i < nxl ) THEN
	1482	trlp_count = trlp_count + 1
	1483	IF ( particles(n)%tail_id /= 0 ) trlpt_count = trlpt_count + 1
	1484	ELSEIF ( i > nxr ) THEN
	1485	trrp_count = trrp_count + 1
	1486	IF ( particles(n)%tail_id /= 0 ) trrpt_count = trrpt_count + 1
	1487	ENDIF
	1488	ENDDO
	1489	IF ( trlp_count == 0 ) trlp_count = 1
	1490	IF ( trlpt_count == 0 ) trlpt_count = 1
	1491	IF ( trrp_count == 0 ) trrp_count = 1
	1492	IF ( trrpt_count == 0 ) trrpt_count = 1
	1493
	1494	ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
	1495
	1496	trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1497	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1498	0, 0, 0, 0 )
	1499	trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1500	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1501	0, 0, 0, 0 )
	1502
	1503	IF ( use_particle_tails ) THEN
	1504	ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
	1505	trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
	1506	tlength = maximum_number_of_tailpoints * 5
	1507	ENDIF
	1508
	1509	trlp_count = 0
	1510	trlpt_count = 0
	1511	trrp_count = 0
	1512	trrpt_count = 0
	1513
	1514	ENDIF
	1515
	1516	! WRITE ( 9, * ) '*** advec_particles: ##1'
	1517	! CALL FLUSH_( 9 )
	1518	! nd = 0
	1519	! DO n = 1, number_of_particles
	1520	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	1521	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	1522	! THEN
	1523	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	1524	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	1525	! CALL FLUSH_( 9 )
	1526	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1527	! ENDIF
	1528	! ENDDO
	1529	! IF ( nd /= deleted_particles ) THEN
	1530	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	1531	! CALL FLUSH_( 9 )
	1532	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1533	! ENDIF
	1534
	1535	DO n = 1, number_of_particles
	1536
	1537	nn = particles(n)%tail_id
	1538
	1539	i = ( particles(n)%x + 0.5 * dx ) * ddx
	1540	!
	1541	!-- Above calculation does not work for indices less than zero
	1542	IF ( particles(n)%x < - 0.5 * dx ) i = -1
	1543
	1544	IF ( i < nxl ) THEN
	1545	IF ( i < 0 ) THEN
	1546	!
	1547	!-- Apply boundary condition along x
	1548	IF ( ibc_par_lr == 0 ) THEN
	1549	!
	1550	!-- Cyclic condition
	1551	IF ( pdims(1) == 1 ) THEN
	1552	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
	1553	particles(n)%origin_x = ( nx + 1 ) * dx + &
	1554	particles(n)%origin_x
	1555	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1556	i = particles(n)%tailpoints
	1557	particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
	1558	+ particle_tail_coordinates(1:i,1,nn)
	1559	ENDIF
	1560	ELSE
	1561	trlp_count = trlp_count + 1
	1562	trlp(trlp_count) = particles(n)
	1563	trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
	1564	trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
	1565	( nx + 1 ) * dx
	1566	particle_mask(n) = .FALSE.
	1567	deleted_particles = deleted_particles + 1
	1568
	1569	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1570	trlpt_count = trlpt_count + 1
	1571	trlpt(:,:,trlpt_count) = &
	1572	particle_tail_coordinates(:,:,nn)
	1573	trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
	1574	trlpt(:,1,trlpt_count)
	1575	tail_mask(nn) = .FALSE.
	1576	deleted_tails = deleted_tails + 1
	1577	ENDIF
	1578	ENDIF
	1579
	1580	ELSEIF ( ibc_par_lr == 1 ) THEN
	1581	!
	1582	!-- Particle absorption
	1583	particle_mask(n) = .FALSE.
	1584	deleted_particles = deleted_particles + 1
	1585	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1586	tail_mask(nn) = .FALSE.
	1587	deleted_tails = deleted_tails + 1
	1588	ENDIF
	1589
	1590	ELSEIF ( ibc_par_lr == 2 ) THEN
	1591	!
	1592	!-- Particle reflection
	1593	particles(n)%x = -particles(n)%x
	1594	particles(n)%speed_x = -particles(n)%speed_x
	1595
	1596	ENDIF
	1597	ELSE
	1598	!
	1599	!-- Store particle data in the transfer array, which will be send
	1600	!-- to the neighbouring PE
	1601	trlp_count = trlp_count + 1
	1602	trlp(trlp_count) = particles(n)
	1603	particle_mask(n) = .FALSE.
	1604	deleted_particles = deleted_particles + 1
	1605
	1606	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1607	trlpt_count = trlpt_count + 1
	1608	trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
	1609	tail_mask(nn) = .FALSE.
	1610	deleted_tails = deleted_tails + 1
	1611	ENDIF
	1612	ENDIF
	1613
	1614	ELSEIF ( i > nxr ) THEN
	1615	IF ( i > nx ) THEN
	1616	!
	1617	!-- Apply boundary condition along x
	1618	IF ( ibc_par_lr == 0 ) THEN
	1619	!
	1620	!-- Cyclic condition
	1621	IF ( pdims(1) == 1 ) THEN
	1622	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
	1623	particles(n)%origin_x = particles(n)%origin_x - &
	1624	( nx + 1 ) * dx
	1625	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1626	i = particles(n)%tailpoints
	1627	particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
	1628	+ particle_tail_coordinates(1:i,1,nn)
	1629	ENDIF
	1630	ELSE
	1631	trrp_count = trrp_count + 1
	1632	trrp(trrp_count) = particles(n)
	1633	trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
	1634	trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
	1635	( nx + 1 ) * dx
	1636	particle_mask(n) = .FALSE.
	1637	deleted_particles = deleted_particles + 1
	1638
	1639	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1640	trrpt_count = trrpt_count + 1
	1641	trrpt(:,:,trrpt_count) = &
	1642	particle_tail_coordinates(:,:,nn)
	1643	trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
	1644	( nx + 1 ) * dx
	1645	tail_mask(nn) = .FALSE.
	1646	deleted_tails = deleted_tails + 1
	1647	ENDIF
	1648	ENDIF
	1649
	1650	ELSEIF ( ibc_par_lr == 1 ) THEN
	1651	!
	1652	!-- Particle absorption
	1653	particle_mask(n) = .FALSE.
	1654	deleted_particles = deleted_particles + 1
	1655	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1656	tail_mask(nn) = .FALSE.
	1657	deleted_tails = deleted_tails + 1
	1658	ENDIF
	1659
	1660	ELSEIF ( ibc_par_lr == 2 ) THEN
	1661	!
	1662	!-- Particle reflection
	1663	particles(n)%x = 2 * ( nx * dx ) - particles(n)%x
	1664	particles(n)%speed_x = -particles(n)%speed_x
	1665
	1666	ENDIF
	1667	ELSE
	1668	!
	1669	!-- Store particle data in the transfer array, which will be send
	1670	!-- to the neighbouring PE
	1671	trrp_count = trrp_count + 1
	1672	trrp(trrp_count) = particles(n)
	1673	particle_mask(n) = .FALSE.
	1674	deleted_particles = deleted_particles + 1
	1675
	1676	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1677	trrpt_count = trrpt_count + 1
	1678	trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
	1679	tail_mask(nn) = .FALSE.
	1680	deleted_tails = deleted_tails + 1
	1681	ENDIF
	1682	ENDIF
	1683
	1684	ENDIF
	1685	ENDDO
	1686
	1687	! WRITE ( 9, * ) '*** advec_particles: ##2'
	1688	! CALL FLUSH_( 9 )
	1689	! nd = 0
	1690	! DO n = 1, number_of_particles
	1691	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	1692	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	1693	! THEN
	1694	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	1695	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	1696	! CALL FLUSH_( 9 )
	1697	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1698	! ENDIF
	1699	! ENDDO
	1700	! IF ( nd /= deleted_particles ) THEN
	1701	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	1702	! CALL FLUSH_( 9 )
	1703	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1704	! ENDIF
	1705
	1706	!
	1707	!-- Send left boundary, receive right boundary (but first exchange how many
	1708	!-- and check, if particle storage must be extended)
	1709	IF ( pdims(1) /= 1 ) THEN
	1710
	1711	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'start' )
	1712	CALL MPI_SENDRECV( trlp_count, 1, MPI_INTEGER, pleft, 0, &
	1713	trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
	1714	comm2d, status, ierr )
	1715
	1716	IF ( number_of_particles + trrp_count_recv > &
	1717	maximum_number_of_particles ) &
	1718	THEN
	1719	IF ( netcdf_output ) THEN
	1720	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
	1721	'needs to be increased'
	1722	PRINT*, ' but this is not allowed with ', &
	1723	'NetCDF output switched on'
	1724	CALL MPI_ABORT( comm2d, 9999, ierr )
	1725	ELSE
	1726	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trrp'
	1727	! CALL FLUSH_( 9 )
	1728	CALL allocate_prt_memory( trrp_count_recv )
	1729	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trrp'
	1730	! CALL FLUSH_( 9 )
	1731	ENDIF
	1732	ENDIF
	1733
	1734	CALL MPI_SENDRECV( trlp, trlp_count, mpi_particle_type, pleft, 1, &
	1735	particles(number_of_particles+1), trrp_count_recv,&
	1736	mpi_particle_type, pright, 1, &
	1737	comm2d, status, ierr )
	1738
	1739	IF ( use_particle_tails ) THEN
	1740
	1741	CALL MPI_SENDRECV( trlpt_count, 1, MPI_INTEGER, pleft, 0, &
	1742	trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
	1743	comm2d, status, ierr )
	1744
	1745	IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
	1746	THEN
	1747	IF ( netcdf_output ) THEN
	1748	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
	1749	'needs to be increased'
	1750	PRINT*, ' but this is not allowed wi', &
	1751	'th NetCDF output switched on'
	1752	CALL MPI_ABORT( comm2d, 9999, ierr )
	1753	ELSE
	1754	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trrpt'
	1755	! CALL FLUSH_( 9 )
	1756	CALL allocate_tail_memory( trrpt_count_recv )
	1757	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trrpt'
	1758	! CALL FLUSH_( 9 )
	1759	ENDIF
	1760	ENDIF
	1761
	1762	CALL MPI_SENDRECV( trlpt, trlpt_count*tlength, MPI_REAL, pleft, 1,&
	1763	particle_tail_coordinates(1,1,number_of_tails+1), &
	1764	trrpt_count_recv*tlength, MPI_REAL, pright, 1, &
	1765	comm2d, status, ierr )
	1766	!
	1767	!-- Update the tail ids for the transferred particles
	1768	nn = number_of_tails
	1769	DO n = number_of_particles+1, number_of_particles+trrp_count_recv
	1770	IF ( particles(n)%tail_id /= 0 ) THEN
	1771	nn = nn + 1
	1772	particles(n)%tail_id = nn
	1773	ENDIF
	1774	ENDDO
	1775
	1776	ENDIF
	1777
	1778	number_of_particles = number_of_particles + trrp_count_recv
	1779	number_of_tails = number_of_tails + trrpt_count_recv
	1780	! IF ( number_of_particles /= number_of_tails ) THEN
	1781	! WRITE (9,*) '--- advec_particles: #3'
	1782	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
	1783	! CALL FLUSH_( 9 )
	1784	! ENDIF
	1785
	1786	!
	1787	!-- Send right boundary, receive left boundary
	1788	CALL MPI_SENDRECV( trrp_count, 1, MPI_INTEGER, pright, 0, &
	1789	trlp_count_recv, 1, MPI_INTEGER, pleft, 0, &
	1790	comm2d, status, ierr )
	1791
	1792	IF ( number_of_particles + trlp_count_recv > &
	1793	maximum_number_of_particles ) &
	1794	THEN
	1795	IF ( netcdf_output ) THEN
	1796	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
	1797	'needs to be increased'
	1798	PRINT*, ' but this is not allowed with ', &
	1799	'NetCDF output switched on'
	1800	CALL MPI_ABORT( comm2d, 9999, ierr )
	1801	ELSE
	1802	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trlp'
	1803	! CALL FLUSH_( 9 )
	1804	CALL allocate_prt_memory( trlp_count_recv )
	1805	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trlp'
	1806	! CALL FLUSH_( 9 )
	1807	ENDIF
	1808	ENDIF
	1809
	1810	CALL MPI_SENDRECV( trrp, trrp_count, mpi_particle_type, pright, 1, &
	1811	particles(number_of_particles+1), trlp_count_recv,&
	1812	mpi_particle_type, pleft, 1, &
	1813	comm2d, status, ierr )
	1814
	1815	IF ( use_particle_tails ) THEN
	1816
	1817	CALL MPI_SENDRECV( trrpt_count, 1, MPI_INTEGER, pright, 0, &
	1818	trlpt_count_recv, 1, MPI_INTEGER, pleft, 0, &
	1819	comm2d, status, ierr )
	1820
	1821	IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
	1822	THEN
	1823	IF ( netcdf_output ) THEN
	1824	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
	1825	'needs to be increased'
	1826	PRINT*, ' but this is not allowed wi', &
	1827	'th NetCDF output switched on'
	1828	CALL MPI_ABORT( comm2d, 9999, ierr )
	1829	ELSE
	1830	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trlpt'
	1831	! CALL FLUSH_( 9 )
	1832	CALL allocate_tail_memory( trlpt_count_recv )
	1833	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trlpt'
	1834	! CALL FLUSH_( 9 )
	1835	ENDIF
	1836	ENDIF
	1837
	1838	CALL MPI_SENDRECV( trrpt, trrpt_count*tlength, MPI_REAL, pright, &
	1839	1, particle_tail_coordinates(1,1,number_of_tails+1), &
	1840	trlpt_count_recv*tlength, MPI_REAL, pleft, 1, &
	1841	comm2d, status, ierr )
	1842	!
	1843	!-- Update the tail ids for the transferred particles
	1844	nn = number_of_tails
	1845	DO n = number_of_particles+1, number_of_particles+trlp_count_recv
	1846	IF ( particles(n)%tail_id /= 0 ) THEN
	1847	nn = nn + 1
	1848	particles(n)%tail_id = nn
	1849	ENDIF
	1850	ENDDO
	1851
	1852	ENDIF
	1853
	1854	number_of_particles = number_of_particles + trlp_count_recv
	1855	number_of_tails = number_of_tails + trlpt_count_recv
	1856	! IF ( number_of_particles /= number_of_tails ) THEN
	1857	! WRITE (9,*) '--- advec_particles: #4'
	1858	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
	1859	! CALL FLUSH_( 9 )
	1860	! ENDIF
	1861
	1862	IF ( use_particle_tails ) THEN
	1863	DEALLOCATE( trlpt, trrpt )
	1864	ENDIF
	1865	DEALLOCATE( trlp, trrp )
	1866
	1867	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'pause' )
	1868
	1869	ENDIF
	1870
	1871	! WRITE ( 9, * ) '*** advec_particles: ##3'
	1872	! CALL FLUSH_( 9 )
	1873	! nd = 0
	1874	! DO n = 1, number_of_particles
	1875	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	1876	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	1877	! THEN
	1878	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	1879	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	1880	! CALL FLUSH_( 9 )
	1881	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1882	! ENDIF
	1883	! ENDDO
	1884	! IF ( nd /= deleted_particles ) THEN
	1885	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	1886	! CALL FLUSH_( 9 )
	1887	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1888	! ENDIF
	1889
	1890	!
	1891	!-- Check whether particles have crossed the boundaries in y direction. Note
	1892	!-- that this case can also apply to particles that have just been received
	1893	!-- from the adjacent right or left PE.
	1894	!-- Find out first the number of particles to be transferred and allocate
	1895	!-- temporary arrays needed to store them.
	1896	!-- For a one-dimensional decomposition along x, no transfer is necessary,
	1897	!-- because the particle remains on the PE.
	1898	trsp_count = 0
	1899	trspt_count = 0
	1900	trnp_count = 0
	1901	trnpt_count = 0
	1902	IF ( pdims(2) /= 1 ) THEN
	1903	!
	1904	!-- First calculate the storage necessary for sending and receiving the
	1905	!-- data
	1906	DO n = 1, number_of_particles
	1907	IF ( particle_mask(n) ) THEN
	1908	j = ( particles(n)%y + 0.5 * dy ) * ddy
	1909	!
	1910	!-- Above calculation does not work for indices less than zero
	1911	IF ( particles(n)%y < -0.5 * dy ) j = -1
	1912
	1913	IF ( j < nys ) THEN
	1914	trsp_count = trsp_count + 1
	1915	IF ( particles(n)%tail_id /= 0 ) trspt_count = trspt_count+1
	1916	ELSEIF ( j > nyn ) THEN
	1917	trnp_count = trnp_count + 1
	1918	IF ( particles(n)%tail_id /= 0 ) trnpt_count = trnpt_count+1
	1919	ENDIF
	1920	ENDIF
	1921	ENDDO
	1922	IF ( trsp_count == 0 ) trsp_count = 1
	1923	IF ( trspt_count == 0 ) trspt_count = 1
	1924	IF ( trnp_count == 0 ) trnp_count = 1
	1925	IF ( trnpt_count == 0 ) trnpt_count = 1
	1926
	1927	ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
	1928
	1929	trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1930	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1931	0, 0, 0, 0 )
	1932	trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1933	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	1934	0, 0, 0, 0 )
	1935
	1936	IF ( use_particle_tails ) THEN
	1937	ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
	1938	trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
	1939	tlength = maximum_number_of_tailpoints * 5
	1940	ENDIF
	1941
	1942	trsp_count = 0
	1943	trspt_count = 0
	1944	trnp_count = 0
	1945	trnpt_count = 0
	1946
	1947	ENDIF
	1948
	1949	! WRITE ( 9, * ) '*** advec_particles: ##4'
	1950	! CALL FLUSH_( 9 )
	1951	! nd = 0
	1952	! DO n = 1, number_of_particles
	1953	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	1954	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	1955	! THEN
	1956	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	1957	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	1958	! CALL FLUSH_( 9 )
	1959	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1960	! ENDIF
	1961	! ENDDO
	1962	! IF ( nd /= deleted_particles ) THEN
	1963	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	1964	! CALL FLUSH_( 9 )
	1965	! CALL MPI_ABORT( comm2d, 9999, ierr )
	1966	! ENDIF
	1967
	1968	DO n = 1, number_of_particles
	1969
	1970	nn = particles(n)%tail_id
	1971	!
	1972	!-- Only those particles that have not been marked as 'deleted' may be
	1973	!-- moved.
	1974	IF ( particle_mask(n) ) THEN
	1975	j = ( particles(n)%y + 0.5 * dy ) * ddy
	1976	!
	1977	!-- Above calculation does not work for indices less than zero
	1978	IF ( particles(n)%y < -0.5 * dy ) j = -1
	1979
	1980	IF ( j < nys ) THEN
	1981	IF ( j < 0 ) THEN
	1982	!
	1983	!-- Apply boundary condition along y
	1984	IF ( ibc_par_ns == 0 ) THEN
	1985	!
	1986	!-- Cyclic condition
	1987	IF ( pdims(2) == 1 ) THEN
	1988	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
	1989	particles(n)%origin_y = ( ny + 1 ) * dy + &
	1990	particles(n)%origin_y
	1991	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	1992	i = particles(n)%tailpoints
	1993	particle_tail_coordinates(1:i,2,nn) = ( ny+1 ) * dy&
	1994	+ particle_tail_coordinates(1:i,2,nn)
	1995	ENDIF
	1996	ELSE
	1997	trsp_count = trsp_count + 1
	1998	trsp(trsp_count) = particles(n)
	1999	trsp(trsp_count)%y = ( ny + 1 ) * dy + &
	2000	trsp(trsp_count)%y
	2001	trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
	2002	+ ( ny + 1 ) * dy
	2003	particle_mask(n) = .FALSE.
	2004	deleted_particles = deleted_particles + 1
	2005
	2006	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2007	trspt_count = trspt_count + 1
	2008	trspt(:,:,trspt_count) = &
	2009	particle_tail_coordinates(:,:,nn)
	2010	trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
	2011	trspt(:,2,trspt_count)
	2012	tail_mask(nn) = .FALSE.
	2013	deleted_tails = deleted_tails + 1
	2014	ENDIF
	2015	ENDIF
	2016
	2017	ELSEIF ( ibc_par_ns == 1 ) THEN
	2018	!
	2019	!-- Particle absorption
	2020	particle_mask(n) = .FALSE.
	2021	deleted_particles = deleted_particles + 1
	2022	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2023	tail_mask(nn) = .FALSE.
	2024	deleted_tails = deleted_tails + 1
	2025	ENDIF
	2026
	2027	ELSEIF ( ibc_par_ns == 2 ) THEN
	2028	!
	2029	!-- Particle reflection
	2030	particles(n)%y = -particles(n)%y
	2031	particles(n)%speed_y = -particles(n)%speed_y
	2032
	2033	ENDIF
	2034	ELSE
	2035	!
	2036	!-- Store particle data in the transfer array, which will be send
	2037	!-- to the neighbouring PE
	2038	trsp_count = trsp_count + 1
	2039	trsp(trsp_count) = particles(n)
	2040	particle_mask(n) = .FALSE.
	2041	deleted_particles = deleted_particles + 1
	2042
	2043	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2044	trspt_count = trspt_count + 1
	2045	trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
	2046	tail_mask(nn) = .FALSE.
	2047	deleted_tails = deleted_tails + 1
	2048	ENDIF
	2049	ENDIF
	2050
[16]	2051	ELSEIF ( j > nyn ) THEN
[1]	2052	IF ( j > ny ) THEN
	2053	!
	2054	!-- Apply boundary condition along x
	2055	IF ( ibc_par_ns == 0 ) THEN
	2056	!
	2057	!-- Cyclic condition
	2058	IF ( pdims(2) == 1 ) THEN
	2059	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
	2060	particles(n)%origin_y = particles(n)%origin_y - &
	2061	( ny + 1 ) * dy
	2062	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2063	i = particles(n)%tailpoints
	2064	particle_tail_coordinates(1:i,2,nn) = - (ny+1) * dy&
	2065	+ particle_tail_coordinates(1:i,2,nn)
	2066	ENDIF
	2067	ELSE
	2068	trnp_count = trnp_count + 1
	2069	trnp(trnp_count) = particles(n)
	2070	trnp(trnp_count)%y = trnp(trnp_count)%y - &
	2071	( ny + 1 ) * dy
	2072	trnp(trnp_count)%origin_y = trnp(trnp_count)%origin_y &
	2073	- ( ny + 1 ) * dy
	2074	particle_mask(n) = .FALSE.
	2075	deleted_particles = deleted_particles + 1
	2076
	2077	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2078	trnpt_count = trnpt_count + 1
	2079	trnpt(:,:,trnpt_count) = &
	2080	particle_tail_coordinates(:,:,nn)
	2081	trnpt(:,2,trnpt_count) = trnpt(:,2,trnpt_count) - &
	2082	( ny + 1 ) * dy
	2083	tail_mask(nn) = .FALSE.
	2084	deleted_tails = deleted_tails + 1
	2085	ENDIF
	2086	ENDIF
	2087
	2088	ELSEIF ( ibc_par_ns == 1 ) THEN
	2089	!
	2090	!-- Particle absorption
	2091	particle_mask(n) = .FALSE.
	2092	deleted_particles = deleted_particles + 1
	2093	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2094	tail_mask(nn) = .FALSE.
	2095	deleted_tails = deleted_tails + 1
	2096	ENDIF
	2097
	2098	ELSEIF ( ibc_par_ns == 2 ) THEN
	2099	!
	2100	!-- Particle reflection
	2101	particles(n)%y = 2 * ( ny * dy ) - particles(n)%y
	2102	particles(n)%speed_y = -particles(n)%speed_y
	2103
	2104	ENDIF
	2105	ELSE
	2106	!
	2107	!-- Store particle data in the transfer array, which will be send
	2108	!-- to the neighbouring PE
	2109	trnp_count = trnp_count + 1
	2110	trnp(trnp_count) = particles(n)
	2111	particle_mask(n) = .FALSE.
	2112	deleted_particles = deleted_particles + 1
	2113
	2114	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2115	trnpt_count = trnpt_count + 1
	2116	trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
	2117	tail_mask(nn) = .FALSE.
	2118	deleted_tails = deleted_tails + 1
	2119	ENDIF
	2120	ENDIF
	2121
	2122	ENDIF
	2123	ENDIF
	2124	ENDDO
	2125
	2126	! WRITE ( 9, * ) '*** advec_particles: ##5'
	2127	! CALL FLUSH_( 9 )
	2128	! nd = 0
	2129	! DO n = 1, number_of_particles
	2130	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2131	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2132	! THEN
	2133	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2134	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	2135	! CALL FLUSH_( 9 )
	2136	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2137	! ENDIF
	2138	! ENDDO
	2139	! IF ( nd /= deleted_particles ) THEN
	2140	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	2141	! CALL FLUSH_( 9 )
	2142	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2143	! ENDIF
	2144
	2145	!
	2146	!-- Send front boundary, receive back boundary (but first exchange how many
	2147	!-- and check, if particle storage must be extended)
	2148	IF ( pdims(2) /= 1 ) THEN
	2149
	2150	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'continue' )
	2151	CALL MPI_SENDRECV( trsp_count, 1, MPI_INTEGER, psouth, 0, &
	2152	trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
	2153	comm2d, status, ierr )
	2154
	2155	IF ( number_of_particles + trnp_count_recv > &
	2156	maximum_number_of_particles ) &
	2157	THEN
	2158	IF ( netcdf_output ) THEN
	2159	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
	2160	'needs to be increased'
	2161	PRINT*, ' but this is not allowed with ', &
	2162	'NetCDF output switched on'
	2163	CALL MPI_ABORT( comm2d, 9999, ierr )
	2164	ELSE
	2165	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trnp'
	2166	! CALL FLUSH_( 9 )
	2167	CALL allocate_prt_memory( trnp_count_recv )
	2168	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trnp'
	2169	! CALL FLUSH_( 9 )
	2170	ENDIF
	2171	ENDIF
	2172
	2173	CALL MPI_SENDRECV( trsp, trsp_count, mpi_particle_type, psouth, 1, &
	2174	particles(number_of_particles+1), trnp_count_recv,&
	2175	mpi_particle_type, pnorth, 1, &
	2176	comm2d, status, ierr )
	2177
	2178	IF ( use_particle_tails ) THEN
	2179
	2180	CALL MPI_SENDRECV( trspt_count, 1, MPI_INTEGER, pleft, 0, &
	2181	trnpt_count_recv, 1, MPI_INTEGER, pright, 0, &
	2182	comm2d, status, ierr )
	2183
	2184	IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
	2185	THEN
	2186	IF ( netcdf_output ) THEN
	2187	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
	2188	'needs to be increased'
	2189	PRINT*, ' but this is not allowed wi', &
	2190	'th NetCDF output switched on'
	2191	CALL MPI_ABORT( comm2d, 9999, ierr )
	2192	ELSE
	2193	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trnpt'
	2194	! CALL FLUSH_( 9 )
	2195	CALL allocate_tail_memory( trnpt_count_recv )
	2196	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trnpt'
	2197	! CALL FLUSH_( 9 )
	2198	ENDIF
	2199	ENDIF
	2200
	2201	CALL MPI_SENDRECV( trspt, trspt_count*tlength, MPI_REAL, psouth, &
	2202	1, particle_tail_coordinates(1,1,number_of_tails+1), &
	2203	trnpt_count_recv*tlength, MPI_REAL, pnorth, 1, &
	2204	comm2d, status, ierr )
	2205	!
	2206	!-- Update the tail ids for the transferred particles
	2207	nn = number_of_tails
	2208	DO n = number_of_particles+1, number_of_particles+trnp_count_recv
	2209	IF ( particles(n)%tail_id /= 0 ) THEN
	2210	nn = nn + 1
	2211	particles(n)%tail_id = nn
	2212	ENDIF
	2213	ENDDO
	2214
	2215	ENDIF
	2216
	2217	number_of_particles = number_of_particles + trnp_count_recv
	2218	number_of_tails = number_of_tails + trnpt_count_recv
	2219	! IF ( number_of_particles /= number_of_tails ) THEN
	2220	! WRITE (9,*) '--- advec_particles: #5'
	2221	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
	2222	! CALL FLUSH_( 9 )
	2223	! ENDIF
	2224
	2225	!
	2226	!-- Send back boundary, receive front boundary
	2227	CALL MPI_SENDRECV( trnp_count, 1, MPI_INTEGER, pnorth, 0, &
	2228	trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
	2229	comm2d, status, ierr )
	2230
	2231	IF ( number_of_particles + trsp_count_recv > &
	2232	maximum_number_of_particles ) &
	2233	THEN
	2234	IF ( netcdf_output ) THEN
	2235	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
	2236	'needs to be increased'
	2237	PRINT*, ' but this is not allowed with ', &
	2238	'NetCDF output switched on'
	2239	CALL MPI_ABORT( comm2d, 9999, ierr )
	2240	ELSE
	2241	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trsp'
	2242	! CALL FLUSH_( 9 )
	2243	CALL allocate_prt_memory( trsp_count_recv )
	2244	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trsp'
	2245	! CALL FLUSH_( 9 )
	2246	ENDIF
	2247	ENDIF
	2248
	2249	CALL MPI_SENDRECV( trnp, trnp_count, mpi_particle_type, pnorth, 1, &
	2250	particles(number_of_particles+1), trsp_count_recv,&
	2251	mpi_particle_type, psouth, 1, &
	2252	comm2d, status, ierr )
	2253
	2254	IF ( use_particle_tails ) THEN
	2255
	2256	CALL MPI_SENDRECV( trnpt_count, 1, MPI_INTEGER, pright, 0, &
	2257	trspt_count_recv, 1, MPI_INTEGER, pleft, 0, &
	2258	comm2d, status, ierr )
	2259
	2260	IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
	2261	THEN
	2262	IF ( netcdf_output ) THEN
	2263	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
	2264	'needs to be increased'
	2265	PRINT*, ' but this is not allowed wi', &
	2266	'th NetCDF output switched on'
	2267	CALL MPI_ABORT( comm2d, 9999, ierr )
	2268	ELSE
	2269	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trspt'
	2270	! CALL FLUSH_( 9 )
	2271	CALL allocate_tail_memory( trspt_count_recv )
	2272	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trspt'
	2273	! CALL FLUSH_( 9 )
	2274	ENDIF
	2275	ENDIF
	2276
	2277	CALL MPI_SENDRECV( trnpt, trnpt_count*tlength, MPI_REAL, pnorth, &
	2278	1, particle_tail_coordinates(1,1,number_of_tails+1), &
	2279	trspt_count_recv*tlength, MPI_REAL, psouth, 1, &
	2280	comm2d, status, ierr )
	2281	!
	2282	!-- Update the tail ids for the transferred particles
	2283	nn = number_of_tails
	2284	DO n = number_of_particles+1, number_of_particles+trsp_count_recv
	2285	IF ( particles(n)%tail_id /= 0 ) THEN
	2286	nn = nn + 1
	2287	particles(n)%tail_id = nn
	2288	ENDIF
	2289	ENDDO
	2290
	2291	ENDIF
	2292
	2293	number_of_particles = number_of_particles + trsp_count_recv
	2294	number_of_tails = number_of_tails + trspt_count_recv
	2295	! IF ( number_of_particles /= number_of_tails ) THEN
	2296	! WRITE (9,*) '--- advec_particles: #6'
	2297	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
	2298	! CALL FLUSH_( 9 )
	2299	! ENDIF
	2300
	2301	IF ( use_particle_tails ) THEN
	2302	DEALLOCATE( trspt, trnpt )
	2303	ENDIF
	2304	DEALLOCATE( trsp, trnp )
	2305
	2306	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'stop' )
	2307
	2308	ENDIF
	2309
	2310	! WRITE ( 9, * ) '*** advec_particles: ##6'
	2311	! CALL FLUSH_( 9 )
	2312	! nd = 0
	2313	! DO n = 1, number_of_particles
	2314	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2315	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2316	! THEN
	2317	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2318	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	2319	! CALL FLUSH_( 9 )
	2320	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2321	! ENDIF
	2322	! ENDDO
	2323	! IF ( nd /= deleted_particles ) THEN
	2324	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	2325	! CALL FLUSH_( 9 )
	2326	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2327	! ENDIF
	2328
	2329	#else
	2330
	2331	!
	2332	!-- Apply boundary conditions
	2333	DO n = 1, number_of_particles
	2334
	2335	nn = particles(n)%tail_id
	2336
	2337	IF ( particles(n)%x < -0.5 * dx ) THEN
	2338
	2339	IF ( ibc_par_lr == 0 ) THEN
	2340	!
	2341	!-- Cyclic boundary. Relevant coordinate has to be changed.
	2342	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
	2343	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2344	i = particles(n)%tailpoints
	2345	particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
	2346	particle_tail_coordinates(1:i,1,nn)
	2347	ENDIF
	2348	ELSEIF ( ibc_par_lr == 1 ) THEN
	2349	!
	2350	!-- Particle absorption
	2351	particle_mask(n) = .FALSE.
	2352	deleted_particles = deleted_particles + 1
	2353	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2354	tail_mask(nn) = .FALSE.
	2355	deleted_tails = deleted_tails + 1
	2356	ENDIF
	2357	ELSEIF ( ibc_par_lr == 2 ) THEN
	2358	!
	2359	!-- Particle reflection
	2360	particles(n)%x = -dx - particles(n)%x
	2361	particles(n)%speed_x = -particles(n)%speed_x
	2362	ENDIF
	2363
	2364	ELSEIF ( particles(n)%x >= ( nx + 0.5 ) * dx ) THEN
	2365
	2366	IF ( ibc_par_lr == 0 ) THEN
	2367	!
	2368	!-- Cyclic boundary. Relevant coordinate has to be changed.
	2369	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
	2370	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2371	i = particles(n)%tailpoints
	2372	particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
	2373	particle_tail_coordinates(1:i,1,nn)
	2374	ENDIF
	2375	ELSEIF ( ibc_par_lr == 1 ) THEN
	2376	!
	2377	!-- Particle absorption
	2378	particle_mask(n) = .FALSE.
	2379	deleted_particles = deleted_particles + 1
	2380	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2381	tail_mask(nn) = .FALSE.
	2382	deleted_tails = deleted_tails + 1
	2383	ENDIF
	2384	ELSEIF ( ibc_par_lr == 2 ) THEN
	2385	!
	2386	!-- Particle reflection
	2387	particles(n)%x = ( nx + 1 ) * dx - particles(n)%x
	2388	particles(n)%speed_x = -particles(n)%speed_x
	2389	ENDIF
	2390
	2391	ENDIF
	2392
	2393	IF ( particles(n)%y < -0.5 * dy ) THEN
	2394
	2395	IF ( ibc_par_ns == 0 ) THEN
	2396	!
	2397	!-- Cyclic boundary. Relevant coordinate has to be changed.
	2398	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
	2399	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2400	i = particles(n)%tailpoints
	2401	particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
	2402	particle_tail_coordinates(1:i,2,nn)
	2403	ENDIF
	2404	ELSEIF ( ibc_par_ns == 1 ) THEN
	2405	!
	2406	!-- Particle absorption
	2407	particle_mask(n) = .FALSE.
	2408	deleted_particles = deleted_particles + 1
	2409	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2410	tail_mask(nn) = .FALSE.
	2411	deleted_tails = deleted_tails + 1
	2412	ENDIF
	2413	ELSEIF ( ibc_par_ns == 2 ) THEN
	2414	!
	2415	!-- Particle reflection
	2416	particles(n)%y = -dy - particles(n)%y
	2417	particles(n)%speed_y = -particles(n)%speed_y
	2418	ENDIF
	2419
	2420	ELSEIF ( particles(n)%y >= ( ny + 0.5 ) * dy ) THEN
	2421
	2422	IF ( ibc_par_ns == 0 ) THEN
	2423	!
	2424	!-- Cyclic boundary. Relevant coordinate has to be changed.
	2425	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
	2426	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2427	i = particles(n)%tailpoints
	2428	particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
	2429	particle_tail_coordinates(1:i,2,nn)
	2430	ENDIF
	2431	ELSEIF ( ibc_par_ns == 1 ) THEN
	2432	!
	2433	!-- Particle absorption
	2434	particle_mask(n) = .FALSE.
	2435	deleted_particles = deleted_particles + 1
	2436	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2437	tail_mask(nn) = .FALSE.
	2438	deleted_tails = deleted_tails + 1
	2439	ENDIF
	2440	ELSEIF ( ibc_par_ns == 2 ) THEN
	2441	!
	2442	!-- Particle reflection
	2443	particles(n)%y = ( ny + 1 ) * dy - particles(n)%y
	2444	particles(n)%speed_y = -particles(n)%speed_y
	2445	ENDIF
	2446
	2447	ENDIF
	2448	ENDDO
	2449
	2450	#endif
	2451
	2452	!
	2453	!-- Apply boundary conditions to those particles that have crossed the top or
	2454	!-- bottom boundary and delete those particles, which are older than allowed
	2455	DO n = 1, number_of_particles
	2456
	2457	nn = particles(n)%tail_id
	2458
	2459	IF ( particles(n)%age > particle_maximum_age .AND. &
	2460	particle_mask(n) ) &
	2461	THEN
	2462	particle_mask(n) = .FALSE.
	2463	deleted_particles = deleted_particles + 1
	2464	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2465	tail_mask(nn) = .FALSE.
	2466	deleted_tails = deleted_tails + 1
	2467	ENDIF
	2468	ENDIF
	2469
	2470	IF ( particles(n)%z >= zu(nz) .AND. particle_mask(n) ) THEN
	2471	IF ( ibc_par_t == 1 ) THEN
	2472	!
	2473	!-- Particle absorption
	2474	particle_mask(n) = .FALSE.
	2475	deleted_particles = deleted_particles + 1
	2476	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2477	tail_mask(nn) = .FALSE.
	2478	deleted_tails = deleted_tails + 1
	2479	ENDIF
	2480	ELSEIF ( ibc_par_t == 2 ) THEN
	2481	!
	2482	!-- Particle reflection
	2483	particles(n)%z = 2.0 * zu(nz) - particles(n)%z
	2484	particles(n)%speed_z = -particles(n)%speed_z
	2485	IF ( use_sgs_for_particles .AND. &
	2486	particles(n)%speed_z_sgs > 0.0 ) THEN
	2487	particles(n)%speed_z_sgs = -particles(n)%speed_z_sgs
	2488	ENDIF
	2489	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2490	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
	2491	particle_tail_coordinates(1,3,nn)
	2492	ENDIF
	2493	ENDIF
	2494	ENDIF
	2495	IF ( particles(n)%z < 0.0 .AND. particle_mask(n) ) THEN
	2496	IF ( ibc_par_b == 1 ) THEN
	2497	!
	2498	!-- Particle absorption
	2499	particle_mask(n) = .FALSE.
	2500	deleted_particles = deleted_particles + 1
	2501	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2502	tail_mask(nn) = .FALSE.
	2503	deleted_tails = deleted_tails + 1
	2504	ENDIF
	2505	ELSEIF ( ibc_par_b == 2 ) THEN
	2506	!
	2507	!-- Particle reflection
	2508	particles(n)%z = -particles(n)%z
	2509	particles(n)%speed_z = -particles(n)%speed_z
	2510	IF ( use_sgs_for_particles .AND. &
	2511	particles(n)%speed_z_sgs < 0.0 ) THEN
	2512	particles(n)%speed_z_sgs = -particles(n)%speed_z_sgs
	2513	ENDIF
	2514	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2515	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
	2516	particle_tail_coordinates(1,3,nn)
	2517	ENDIF
	2518	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2519	particle_tail_coordinates(1,3,nn) = &
	2520	-particle_tail_coordinates(1,3,nn)
	2521	ENDIF
	2522	ENDIF
	2523	ENDIF
	2524	ENDDO
	2525
	2526	! WRITE ( 9, * ) '*** advec_particles: ##7'
	2527	! CALL FLUSH_( 9 )
	2528	! nd = 0
	2529	! DO n = 1, number_of_particles
	2530	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2531	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2532	! THEN
	2533	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2534	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	2535	! CALL FLUSH_( 9 )
	2536	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2537	! ENDIF
	2538	! ENDDO
	2539	! IF ( nd /= deleted_particles ) THEN
	2540	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
	2541	! CALL FLUSH_( 9 )
	2542	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2543	! ENDIF
	2544
	2545	!
	2546	!-- Pack particles (eliminate those marked for deletion),
	2547	!-- determine new number of particles
	2548	IF ( number_of_particles > 0 .AND. deleted_particles > 0 ) THEN
	2549	nn = 0
	2550	nd = 0
	2551	DO n = 1, number_of_particles
	2552	IF ( particle_mask(n) ) THEN
	2553	nn = nn + 1
	2554	particles(nn) = particles(n)
	2555	ELSE
	2556	nd = nd + 1
	2557	ENDIF
	2558	ENDDO
	2559	! IF ( nd /= deleted_particles ) THEN
	2560	! WRITE (9,) '** advec_part nd=',nd,' deleted_particles=',deleted_particles
	2561	! CALL FLUSH_( 9 )
	2562	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2563	! ENDIF
	2564
	2565	number_of_particles = number_of_particles - deleted_particles
	2566	!
	2567	!-- Pack the tails, store the new tail ids and re-assign it to the
	2568	!-- respective
	2569	!-- particles
	2570	IF ( use_particle_tails ) THEN
	2571	nn = 0
	2572	nd = 0
	2573	DO n = 1, number_of_tails
	2574	IF ( tail_mask(n) ) THEN
	2575	nn = nn + 1
	2576	particle_tail_coordinates(:,:,nn) = &
	2577	particle_tail_coordinates(:,:,n)
	2578	new_tail_id(n) = nn
	2579	ELSE
	2580	nd = nd + 1
	2581	! WRITE (9,*) '+++ n=',n,' (of ',number_of_tails,' #oftails)'
	2582	! WRITE (9,*) ' id=',new_tail_id(n)
	2583	! CALL FLUSH_( 9 )
	2584	ENDIF
	2585	ENDDO
	2586	ENDIF
	2587
	2588	! IF ( nd /= deleted_tails .AND. use_particle_tails ) THEN
	2589	! WRITE (9,) '** advec_part nd=',nd,' deleted_tails=',deleted_tails
	2590	! CALL FLUSH_( 9 )
	2591	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2592	! ENDIF
	2593
	2594	number_of_tails = number_of_tails - deleted_tails
	2595
	2596	! nn = 0
	2597	DO n = 1, number_of_particles
	2598	IF ( particles(n)%tail_id /= 0 ) THEN
	2599	! nn = nn + 1
	2600	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id > number_of_tails ) THEN
	2601	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2602	! WRITE (9,*) ' tail_id=',particles(n)%tail_id
	2603	! WRITE (9,*) ' new_tail_id=', new_tail_id(particles(n)%tail_id), &
	2604	! ' of (',number_of_tails,')'
	2605	! CALL FLUSH_( 9 )
	2606	! ENDIF
	2607	particles(n)%tail_id = new_tail_id(particles(n)%tail_id)
	2608	ENDIF
	2609	ENDDO
	2610
	2611	! IF ( nn /= number_of_tails .AND. use_particle_tails ) THEN
	2612	! WRITE (9,) '** advec_part #of_tails=',number_of_tails,' nn=',nn
	2613	! CALL FLUSH_( 9 )
	2614	! DO n = 1, number_of_particles
	2615	! WRITE (9,*) 'prt# ',n,' tail_id=',particles(n)%tail_id, &
	2616	! ' x=',particles(n)%x, ' y=',particles(n)%y, &
	2617	! ' z=',particles(n)%z
	2618	! ENDDO
	2619	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2620	! ENDIF
	2621
	2622	ENDIF
	2623
	2624	! IF ( number_of_particles /= number_of_tails ) THEN
	2625	! WRITE (9,*) '--- advec_particles: #7'
	2626	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
	2627	! CALL FLUSH_( 9 )
	2628	! ENDIF
	2629	! WRITE ( 9, * ) '*** advec_particles: ##8'
	2630	! CALL FLUSH_( 9 )
	2631	! DO n = 1, number_of_particles
	2632	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2633	! THEN
	2634	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2635	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	2636	! CALL FLUSH_( 9 )
	2637	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2638	! ENDIF
	2639	! ENDDO
	2640
	2641	!
	2642	!-- Sort particles in the sequence the gridboxes are stored in the memory
	2643	CALL sort_particles
	2644
	2645	! WRITE ( 9, * ) '*** advec_particles: ##9'
	2646	! CALL FLUSH_( 9 )
	2647	! DO n = 1, number_of_particles
	2648	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2649	! THEN
	2650	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2651	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	2652	! CALL FLUSH_( 9 )
	2653	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2654	! ENDIF
	2655	! ENDDO
	2656
	2657	!
	2658	!-- Accumulate the number of particles transferred between the subdomains
	2659	#if defined( __parallel )
	2660	trlp_count_sum = trlp_count_sum + trlp_count
	2661	trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
	2662	trrp_count_sum = trrp_count_sum + trrp_count
	2663	trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
	2664	trsp_count_sum = trsp_count_sum + trsp_count
	2665	trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
	2666	trnp_count_sum = trnp_count_sum + trnp_count
	2667	trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
	2668	#endif
	2669
	2670	IF ( dt_3d_reached ) EXIT
	2671
	2672	ENDDO ! timestep loop
	2673
	2674
	2675	!
	2676	!-- Re-evaluate the weighting factors. After advection, particles within a
	2677	!-- grid box may have different weighting factors if some have been advected
	2678	!-- from a neighbouring box. The factors are re-evaluated so that they are
	2679	!-- the same for all particles of one box. This procedure must conserve the
	2680	!-- liquid water content within one box.
	2681	IF ( cloud_droplets ) THEN
	2682
	2683	CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'start' )
	2684
	2685	ql = 0.0; ql_v = 0.0; ql_vp = 0.0
	2686
	2687	!
	2688	!-- Re-calculate the weighting factors and calculate the liquid water content
	2689	DO i = nxl, nxr
	2690	DO j = nys, nyn
	2691	DO k = nzb, nzt+1
	2692
	2693	!
	2694	!-- Calculate the total volume of particles in the boxes (ql_vp) as
	2695	!-- well as the real volume (ql_v, weighting factor has to be
	2696	!-- included)
	2697	psi = prt_start_index(k,j,i)
	2698	DO n = psi, psi+prt_count(k,j,i)-1
	2699	ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3
	2700
	2701	ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
	2702	particles(n)%radius**3
	2703	ENDDO
	2704
	2705	!
	2706	!-- Re-calculate the weighting factors and calculate the liquid
	2707	!-- water content
	2708	IF ( ql_vp(k,j,i) /= 0.0 ) THEN
	2709	ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
	2710	ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * &
	2711	ql_v(k,j,i) / &
	2712	( rho_surface * dx * dy * dz )
	2713	ELSE
	2714	ql(k,j,i) = 0.0
	2715	ENDIF
	2716
	2717	!
	2718	!-- Re-assign the weighting factor to the particles
	2719	DO n = psi, psi+prt_count(k,j,i)-1
	2720	particles(n)%weight_factor = ql_vp(k,j,i)
	2721	ENDDO
	2722
	2723	ENDDO
	2724	ENDDO
	2725	ENDDO
	2726
	2727	CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'stop' )
	2728
	2729	ENDIF
	2730
	2731	!
	2732	!-- Set particle attributes defined by the user
	2733	CALL user_particle_attributes
	2734	! WRITE ( 9, * ) '*** advec_particles: ##10'
	2735	! CALL FLUSH_( 9 )
	2736	! DO n = 1, number_of_particles
	2737	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2738	! THEN
	2739	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2740	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
	2741	! CALL FLUSH_( 9 )
	2742	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2743	! ENDIF
	2744	! ENDDO
	2745
	2746	!
	2747	!-- If necessary, add the actual particle positions to the particle tails
	2748	IF ( use_particle_tails ) THEN
	2749
	2750	distance = 0.0
	2751	DO n = 1, number_of_particles
	2752
	2753	nn = particles(n)%tail_id
	2754
	2755	IF ( nn /= 0 ) THEN
	2756	!
	2757	!-- Calculate the distance between the actual particle position and the
	2758	!-- next tailpoint
	2759	! WRITE ( 9, * ) '*** advec_particles: ##10.1 nn=',nn
	2760	! CALL FLUSH_( 9 )
	2761	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
	2762	distance = ( particle_tail_coordinates(1,1,nn) - &
	2763	particle_tail_coordinates(2,1,nn) )**2 + &
	2764	( particle_tail_coordinates(1,2,nn) - &
	2765	particle_tail_coordinates(2,2,nn) )**2 + &
	2766	( particle_tail_coordinates(1,3,nn) - &
	2767	particle_tail_coordinates(2,3,nn) )**2
	2768	ENDIF
	2769	! WRITE ( 9, * ) '*** advec_particles: ##10.2'
	2770	! CALL FLUSH_( 9 )
	2771	!
	2772	!-- First, increase the index of all existings tailpoints by one
	2773	IF ( distance >= minimum_tailpoint_distance ) THEN
	2774	DO i = particles(n)%tailpoints, 1, -1
	2775	particle_tail_coordinates(i+1,:,nn) = &
	2776	particle_tail_coordinates(i,:,nn)
	2777	ENDDO
	2778	!
	2779	!-- Increase the counter which contains the number of tailpoints.
	2780	!-- This must always be smaller than the given maximum number of
	2781	!-- tailpoints because otherwise the index bounds of
	2782	!-- particle_tail_coordinates would be exceeded
	2783	IF ( particles(n)%tailpoints < maximum_number_of_tailpoints-1 )&
	2784	THEN
	2785	particles(n)%tailpoints = particles(n)%tailpoints + 1
	2786	ENDIF
	2787	ENDIF
	2788	! WRITE ( 9, * ) '*** advec_particles: ##10.3'
	2789	! CALL FLUSH_( 9 )
	2790	!
	2791	!-- In any case, store the new point at the beginning of the tail
	2792	particle_tail_coordinates(1,1,nn) = particles(n)%x
	2793	particle_tail_coordinates(1,2,nn) = particles(n)%y
	2794	particle_tail_coordinates(1,3,nn) = particles(n)%z
	2795	particle_tail_coordinates(1,4,nn) = particles(n)%color
	2796	! WRITE ( 9, * ) '*** advec_particles: ##10.4'
	2797	! CALL FLUSH_( 9 )
	2798	!
	2799	!-- Increase the age of the tailpoints
	2800	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
	2801	particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) = &
	2802	particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) + &
	2803	dt_3d
	2804	!
	2805	!-- Delete the last tailpoint, if it has exceeded its maximum age
	2806	IF ( particle_tail_coordinates(particles(n)%tailpoints,5,nn) > &
	2807	maximum_tailpoint_age ) THEN
	2808	particles(n)%tailpoints = particles(n)%tailpoints - 1
	2809	ENDIF
	2810	ENDIF
	2811	! WRITE ( 9, * ) '*** advec_particles: ##10.5'
	2812	! CALL FLUSH_( 9 )
	2813
	2814	ENDIF
	2815
	2816	ENDDO
	2817
	2818	ENDIF
	2819	! WRITE ( 9, * ) '*** advec_particles: ##11'
	2820	! CALL FLUSH_( 9 )
	2821
	2822	!
	2823	!-- Write particle statistics on file
	2824	IF ( write_particle_statistics ) THEN
	2825	CALL check_open( 80 )
	2826	#if defined( __parallel )
	2827	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
	2828	number_of_particles, pleft, trlp_count_sum, &
	2829	trlp_count_recv_sum, pright, trrp_count_sum, &
	2830	trrp_count_recv_sum, psouth, trsp_count_sum, &
	2831	trsp_count_recv_sum, pnorth, trnp_count_sum, &
	2832	trnp_count_recv_sum, maximum_number_of_particles
	2833	CALL close_file( 80 )
	2834	#else
	2835	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
	2836	number_of_particles, maximum_number_of_particles
	2837	#endif
	2838	ENDIF
	2839
	2840	CALL cpu_log( log_point(25), 'advec_particles', 'stop' )
	2841
	2842	!
	2843	!-- Formats
	2844	8000 FORMAT (I6,1X,F7.2,4X,I6,5X,4(I3,1X,I4,'/',I4,2X),6X,I6)
	2845
	2846	#endif
	2847	END SUBROUTINE advec_particles
	2848
	2849
	2850	SUBROUTINE allocate_prt_memory( number_of_new_particles )
	2851
	2852	!------------------------------------------------------------------------------!
	2853	! Description:
	2854	! ------------
	2855	! Extend particle memory
	2856	!------------------------------------------------------------------------------!
	2857	#if defined( __particles )
	2858
	2859	USE particle_attributes
	2860
	2861	IMPLICIT NONE
	2862
	2863	INTEGER :: new_maximum_number, number_of_new_particles
	2864
	2865	LOGICAL, DIMENSION(:), ALLOCATABLE :: tmp_particle_mask
	2866
	2867	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles
	2868
	2869
	2870	new_maximum_number = maximum_number_of_particles + &
	2871	MAX( 5*number_of_new_particles, number_of_initial_particles )
	2872
	2873	IF ( write_particle_statistics ) THEN
	2874	CALL check_open( 80 )
	2875	WRITE ( 80, '(''*** Request: '', I7, '' new_maximum_number(prt)'')' ) &
	2876	new_maximum_number
	2877	CALL close_file( 80 )
	2878	ENDIF
	2879
	2880	ALLOCATE( tmp_particles(maximum_number_of_particles), &
	2881	tmp_particle_mask(maximum_number_of_particles) )
	2882	tmp_particles = particles
	2883	tmp_particle_mask = particle_mask
	2884
	2885	DEALLOCATE( particles, particle_mask )
	2886	ALLOCATE( particles(new_maximum_number), &
	2887	particle_mask(new_maximum_number) )
	2888	maximum_number_of_particles = new_maximum_number
	2889
	2890	particles(1:number_of_particles) = tmp_particles(1:number_of_particles)
	2891	particle_mask(1:number_of_particles) = &
	2892	tmp_particle_mask(1:number_of_particles)
	2893	particle_mask(number_of_particles+1:maximum_number_of_particles) = .TRUE.
	2894	DEALLOCATE( tmp_particles, tmp_particle_mask )
	2895
	2896	#endif
	2897	END SUBROUTINE allocate_prt_memory
	2898
	2899
	2900	SUBROUTINE allocate_tail_memory( number_of_new_tails )
	2901
	2902	!------------------------------------------------------------------------------!
	2903	! Description:
	2904	! ------------
	2905	! Extend tail memory
	2906	!------------------------------------------------------------------------------!
	2907	#if defined( __particles )
	2908
	2909	USE particle_attributes
	2910
	2911	IMPLICIT NONE
	2912
	2913	INTEGER :: new_maximum_number, number_of_new_tails
	2914
	2915	LOGICAL, DIMENSION(maximum_number_of_tails) :: tmp_tail_mask
	2916
	2917	REAL, DIMENSION(maximum_number_of_tailpoints,5,maximum_number_of_tails) :: &
	2918	tmp_tail
	2919
	2920
	2921	new_maximum_number = maximum_number_of_tails + &
	2922	MAX( 5*number_of_new_tails, number_of_initial_tails )
	2923
	2924	IF ( write_particle_statistics ) THEN
	2925	CALL check_open( 80 )
	2926	WRITE ( 80, '(''*** Request: '', I5, '' new_maximum_number(tails)'')' ) &
	2927	new_maximum_number
	2928	CALL close_file( 80 )
	2929	ENDIF
	2930	WRITE (9,) '** Request: ',new_maximum_number,' new_maximum_number(tails)'
[15]	2931	! CALL FLUSH_( 9 )
[1]	2932
	2933	tmp_tail(:,:,1:number_of_tails) = &
	2934	particle_tail_coordinates(:,:,1:number_of_tails)
	2935	tmp_tail_mask(1:number_of_tails) = tail_mask(1:number_of_tails)
	2936
	2937	DEALLOCATE( new_tail_id, particle_tail_coordinates, tail_mask )
	2938	ALLOCATE( new_tail_id(new_maximum_number), &
	2939	particle_tail_coordinates(maximum_number_of_tailpoints,5, &
	2940	new_maximum_number), &
	2941	tail_mask(new_maximum_number) )
	2942	maximum_number_of_tails = new_maximum_number
	2943
	2944	particle_tail_coordinates = 0.0
	2945	particle_tail_coordinates(:,:,1:number_of_tails) = &
	2946	tmp_tail(:,:,1:number_of_tails)
	2947	tail_mask(1:number_of_tails) = tmp_tail_mask(1:number_of_tails)
	2948	tail_mask(number_of_tails+1:maximum_number_of_tails) = .TRUE.
	2949
	2950	#endif
	2951	END SUBROUTINE allocate_tail_memory
	2952
	2953
	2954	SUBROUTINE output_particles_netcdf
	2955	#if defined( __netcdf )
	2956
	2957	USE control_parameters
	2958	USE netcdf_control
	2959	USE particle_attributes
	2960
	2961	IMPLICIT NONE
	2962
	2963
	2964	CALL check_open( 108 )
	2965
	2966	!
	2967	!-- Update the NetCDF time axis
	2968	prt_time_count = prt_time_count + 1
	2969
	2970	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, (/ simulated_time /), &
	2971	start = (/ prt_time_count /), count = (/ 1 /) )
	2972	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 1 )
	2973
	2974	!
	2975	!-- Output the real number of particles used
	2976	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, &
	2977	(/ number_of_particles /), &
	2978	start = (/ prt_time_count /), count = (/ 1 /) )
	2979	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 2 )
	2980
	2981	!
	2982	!-- Output all particle attributes
	2983	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age, &
	2984	start = (/ 1, prt_time_count /), &
	2985	count = (/ maximum_number_of_particles /) )
	2986	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 3 )
	2987
	2988	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize, &
	2989	start = (/ 1, prt_time_count /), &
	2990	count = (/ maximum_number_of_particles /) )
	2991	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 4 )
	2992
	2993	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x, &
	2994	start = (/ 1, prt_time_count /), &
	2995	count = (/ maximum_number_of_particles /) )
	2996	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 5 )
	2997
	2998	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y, &
	2999	start = (/ 1, prt_time_count /), &
	3000	count = (/ maximum_number_of_particles /) )
	3001	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 6 )
	3002
	3003	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z, &
	3004	start = (/ 1, prt_time_count /), &
	3005	count = (/ maximum_number_of_particles /) )
	3006	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 7 )
	3007
	3008	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius, &
	3009	start = (/ 1, prt_time_count /), &
	3010	count = (/ maximum_number_of_particles /) )
	3011	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 8 )
	3012
	3013	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x, &
	3014	start = (/ 1, prt_time_count /), &
	3015	count = (/ maximum_number_of_particles /) )
	3016	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 9 )
	3017
	3018	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y, &
	3019	start = (/ 1, prt_time_count /), &
	3020	count = (/ maximum_number_of_particles /) )
	3021	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 10 )
	3022
	3023	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z, &
	3024	start = (/ 1, prt_time_count /), &
	3025	count = (/ maximum_number_of_particles /) )
	3026	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 11 )
	3027
	3028	nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10),particles%weight_factor,&
	3029	start = (/ 1, prt_time_count /), &
	3030	count = (/ maximum_number_of_particles /) )
	3031	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 12 )
	3032
	3033	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x, &
	3034	start = (/ 1, prt_time_count /), &
	3035	count = (/ maximum_number_of_particles /) )
	3036	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 13 )
	3037
	3038	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y, &
	3039	start = (/ 1, prt_time_count /), &
	3040	count = (/ maximum_number_of_particles /) )
	3041	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 14 )
	3042
	3043	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z, &
	3044	start = (/ 1, prt_time_count /), &
	3045	count = (/ maximum_number_of_particles /) )
	3046	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 15 )
	3047
	3048	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%color, &
	3049	start = (/ 1, prt_time_count /), &
	3050	count = (/ maximum_number_of_particles /) )
	3051	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 16 )
	3052
	3053	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group, &
	3054	start = (/ 1, prt_time_count /), &
	3055	count = (/ maximum_number_of_particles /) )
	3056	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 17 )
	3057
	3058	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16), particles%tailpoints, &
	3059	start = (/ 1, prt_time_count /), &
	3060	count = (/ maximum_number_of_particles /) )
	3061	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 18 )
	3062
	3063	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id, &
	3064	start = (/ 1, prt_time_count /), &
	3065	count = (/ maximum_number_of_particles /) )
	3066	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 19 )
	3067
	3068	#endif
	3069	END SUBROUTINE output_particles_netcdf
	3070
	3071
	3072	SUBROUTINE write_particles
	3073
	3074	!------------------------------------------------------------------------------!
	3075	! Description:
	3076	! ------------
	3077	! Write particle data on restart file
	3078	!------------------------------------------------------------------------------!
	3079	#if defined( __particles )
	3080
	3081	USE control_parameters
	3082	USE particle_attributes
	3083	USE pegrid
	3084
	3085	IMPLICIT NONE
	3086
	3087	CHARACTER (LEN=10) :: particle_binary_version
	3088
	3089	!
	3090	!-- First open the output unit.
	3091	IF ( myid_char == '' ) THEN
	3092	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, &
	3093	FORM='UNFORMATTED')
	3094	ELSE
	3095	IF ( myid == 0 ) CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' )
	3096	#if defined( __parallel )
	3097	!
	3098	!-- Set a barrier in order to allow that thereafter all other processors
	3099	!-- in the directory created by PE0 can open their file
	3100	CALL MPI_BARRIER( comm2d, ierr )
	3101	#endif
	3102	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, &
	3103	FORM='UNFORMATTED' )
	3104	ENDIF
	3105
	3106	!
	3107	!-- Write the version number of the binary format.
	3108	!-- Attention: After changes to the following output commands the version
	3109	!-- --------- number of the variable particle_binary_version must be changed!
	3110	!-- Also, the version number and the list of arrays to be read in
	3111	!-- init_particles must be adjusted accordingly.
	3112	particle_binary_version = '3.0'
	3113	WRITE ( 90 ) particle_binary_version
	3114
	3115	!
	3116	!-- Write some particle parameters, the size of the particle arrays as well as
	3117	!-- other dvrp-plot variables.
	3118	WRITE ( 90 ) bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, &
	3119	maximum_number_of_particles, maximum_number_of_tailpoints, &
	3120	maximum_number_of_tails, number_of_initial_particles, &
	3121	number_of_particles, number_of_particle_groups, &
	3122	number_of_tails, particle_groups, time_prel, &
	3123	time_write_particle_data, uniform_particles
	3124
	3125	IF ( number_of_initial_particles /= 0 ) WRITE ( 90 ) initial_particles
	3126
	3127	WRITE ( 90 ) prt_count, prt_start_index
	3128	WRITE ( 90 ) particles
	3129
	3130	IF ( use_particle_tails ) THEN
	3131	WRITE ( 90 ) particle_tail_coordinates
	3132	ENDIF
	3133
	3134	CLOSE ( 90 )
	3135
	3136	#endif
	3137	END SUBROUTINE write_particles
	3138
	3139
	3140	SUBROUTINE collision_efficiency( mean_r, r, e)
	3141	!------------------------------------------------------------------------------!
	3142	! Description:
	3143	! ------------
	3144	! Interpolate collision efficiency from table
	3145	!------------------------------------------------------------------------------!
	3146	#if defined( __particles )
	3147
	3148	IMPLICIT NONE
	3149
	3150	INTEGER :: i, j, k
	3151
	3152	LOGICAL, SAVE :: first = .TRUE.
	3153
	3154	REAL :: aa, bb, cc, dd, dx, dy, e, gg, mean_r, mean_rm, r, rm, &
	3155	x, y
	3156
	3157	REAL, DIMENSION(1:9), SAVE :: collected_r = 0.0
	3158	REAL, DIMENSION(1:19), SAVE :: collector_r = 0.0
	3159	REAL, DIMENSION(1:9,1:19), SAVE :: ef = 0.0
	3160
	3161	mean_rm = mean_r * 1.0E06
	3162	rm = r * 1.0E06
	3163
	3164	IF ( first ) THEN
	3165	collected_r = (/ 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0 /)
	3166	collector_r = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0, 150.0,&
	3167	200.0, 300.0, 400.0, 500.0, 600.0, 1000.0, 1400.0, &
	3168	1800.0, 2400.0, 3000.0 /)
	3169	ef(:,1) = (/0.017, 0.027, 0.037, 0.052, 0.052, 0.052, 0.052, 0.0, 0.0 /)
	3170	ef(:,2) = (/0.001, 0.016, 0.027, 0.060, 0.12, 0.17, 0.17, 0.17, 0.0 /)
	3171	ef(:,3) = (/0.001, 0.001, 0.02, 0.13, 0.28, 0.37, 0.54, 0.55, 0.47/)
	3172	ef(:,4) = (/0.001, 0.001, 0.02, 0.23, 0.4, 0.55, 0.7, 0.75, 0.75/)
	3173	ef(:,5) = (/0.01, 0.01, 0.03, 0.3, 0.4, 0.58, 0.73, 0.75, 0.79/)
	3174	ef(:,6) = (/0.01, 0.01, 0.13, 0.38, 0.57, 0.68, 0.80, 0.86, 0.91/)
	3175	ef(:,7) = (/0.01, 0.085, 0.23, 0.52, 0.68, 0.76, 0.86, 0.92, 0.95/)
	3176	ef(:,8) = (/0.01, 0.14, 0.32, 0.60, 0.73, 0.81, 0.90, 0.94, 0.96/)
	3177	ef(:,9) = (/0.025, 0.25, 0.43, 0.66, 0.78, 0.83, 0.92, 0.95, 0.96/)
	3178	ef(:,10)= (/0.039, 0.3, 0.46, 0.69, 0.81, 0.87, 0.93, 0.95, 0.96/)
	3179	ef(:,11)= (/0.095, 0.33, 0.51, 0.72, 0.82, 0.87, 0.93, 0.96, 0.97/)
	3180	ef(:,12)= (/0.098, 0.36, 0.51, 0.73, 0.83, 0.88, 0.93, 0.96, 0.97/)
	3181	ef(:,13)= (/0.1, 0.36, 0.52, 0.74, 0.83, 0.88, 0.93, 0.96, 0.97/)
	3182	ef(:,14)= (/0.17, 0.4, 0.54, 0.72, 0.83, 0.88, 0.94, 0.98, 1.0 /)
	3183	ef(:,15)= (/0.15, 0.37, 0.52, 0.74, 0.82, 0.88, 0.94, 0.98, 1.0 /)
	3184	ef(:,16)= (/0.11, 0.34, 0.49, 0.71, 0.83, 0.88, 0.94, 0.95, 1.0 /)
	3185	ef(:,17)= (/0.08, 0.29, 0.45, 0.68, 0.8, 0.86, 0.96, 0.94, 1.0 /)
	3186	ef(:,18)= (/0.04, 0.22, 0.39, 0.62, 0.75, 0.83, 0.92, 0.96, 1.0 /)
	3187	ef(:,19)= (/0.02, 0.16, 0.33, 0.55, 0.71, 0.81, 0.90, 0.94, 1.0 /)
	3188	ENDIF
	3189
	3190	DO k = 1, 8
	3191	IF ( collected_r(k) <= mean_rm ) i = k
	3192	ENDDO
	3193
	3194	DO k = 1, 18
	3195	IF ( collector_r(k) <= rm ) j = k
	3196	ENDDO
	3197
	3198	IF ( rm < 10.0 ) THEN
	3199	e = 0.0
	3200	ELSEIF ( mean_rm < 2.0 ) THEN
	3201	e = 0.001
	3202	ELSEIF ( mean_rm >= 25.0 ) THEN
	3203	IF( j <= 3 ) e = 0.55
	3204	IF( j == 4 ) e = 0.8
	3205	IF( j == 5 ) e = 0.9
	3206	IF( j >=6 ) e = 1.0
	3207	ELSEIF ( rm >= 3000.0 ) THEN
	3208	e = 1.0
	3209	ELSE
	3210	x = mean_rm - collected_r(i)
	3211	y = rm - collected_r(j)
	3212	dx = collected_r(i+1) - collected_r(i)
	3213	dy = collector_r(j+1) - collector_r(j)
	3214	aa = x2 + y2
	3215	bb = ( dx - x )2 + y2
	3216	cc = x2 + ( dy - y )2
	3217	dd = ( dx - x )2 + ( dy - y )2
	3218	gg = aa + bb + cc + dd
	3219
	3220	e = ( (gg-aa)ef(i,j) + (gg-bb)ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
	3221	(gg-dd)ef(i+1,j+1) ) / (3.0gg)
	3222	ENDIF
	3223
	3224	#endif
	3225	END SUBROUTINE collision_efficiency
	3226
	3227
	3228
	3229	SUBROUTINE sort_particles
	3230
	3231	!------------------------------------------------------------------------------!
	3232	! Description:
	3233	! ------------
	3234	! Sort particles in the sequence the grid boxes are stored in memory
	3235	!------------------------------------------------------------------------------!
	3236	#if defined( __particles )
	3237
	3238	USE arrays_3d
	3239	USE control_parameters
	3240	USE cpulog
	3241	USE grid_variables
	3242	USE indices
	3243	USE interfaces
	3244	USE particle_attributes
	3245
	3246	IMPLICIT NONE
	3247
	3248	INTEGER :: i, ilow, j, k, n
	3249
	3250	TYPE(particle_type), DIMENSION(1:number_of_particles) :: particles_temp
	3251
	3252
	3253	CALL cpu_log( log_point_s(47), 'sort_particles', 'start' )
	3254
	3255	!
	3256	!-- Initialize the array used for counting and indexing the particles
	3257	prt_count = 0
	3258
	3259	!
	3260	!-- Count the particles per gridbox
	3261	DO n = 1, number_of_particles
	3262
	3263	i = ( particles(n)%x + 0.5 * dx ) * ddx
	3264	j = ( particles(n)%y + 0.5 * dy ) * ddy
	3265	k = particles(n)%z / dz + 1 ! only exact if equidistant
	3266
	3267	prt_count(k,j,i) = prt_count(k,j,i) + 1
	3268
	3269	IF ( i < nxl .OR. i > nxr .OR. j < nys .OR. j > nyn .OR. k < nzb+1 .OR. &
	3270	k > nzt ) THEN
	3271	PRINT*, '+++ sort_particles: particle out of range: i=', i, ' j=', &
	3272	j, ' k=', k
	3273	PRINT*, ' nxl=', nxl, ' nxr=', nxr, &
	3274	' nys=', nys, ' nyn=', nyn, &
	3275	' nzb=', nzb, ' nzt=', nzt
	3276	ENDIF
	3277
	3278	ENDDO
	3279
	3280	!
	3281	!-- Calculate the lower indices of those ranges of the particles-array
	3282	!-- containing particles which belong to the same gridpox i,j,k
	3283	ilow = 1
	3284	DO i = nxl, nxr
	3285	DO j = nys, nyn
	3286	DO k = nzb+1, nzt
	3287	prt_start_index(k,j,i) = ilow
	3288	ilow = ilow + prt_count(k,j,i)
	3289	ENDDO
	3290	ENDDO
	3291	ENDDO
	3292
	3293	!
	3294	!-- Sorting the particles
	3295	DO n = 1, number_of_particles
	3296
	3297	i = ( particles(n)%x + 0.5 * dx ) * ddx
	3298	j = ( particles(n)%y + 0.5 * dy ) * ddy
	3299	k = particles(n)%z / dz + 1 ! only exact if equidistant
	3300
	3301	particles_temp(prt_start_index(k,j,i)) = particles(n)
	3302
	3303	prt_start_index(k,j,i) = prt_start_index(k,j,i) + 1
	3304
	3305	ENDDO
	3306
	3307	particles(1:number_of_particles) = particles_temp
	3308
	3309	!
	3310	!-- Reset the index array to the actual start position
	3311	DO i = nxl, nxr
	3312	DO j = nys, nyn
	3313	DO k = nzb+1, nzt
	3314	prt_start_index(k,j,i) = prt_start_index(k,j,i) - prt_count(k,j,i)
	3315	ENDDO
	3316	ENDDO
	3317	ENDDO
	3318
	3319	CALL cpu_log( log_point_s(47), 'sort_particles', 'stop' )
	3320
	3321	#endif
	3322	END SUBROUTINE sort_particles

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |