Home

Context Navigation

advec_particles.f90 @ 825

Last change on this file since 825 was 825, checked in by raasch, 13 years ago

New:
---

droplet growth by condensation may include curvature and solution effects.
Steered by new inipar-parameter curvature_solution_effects.
(advec_particles, check_parameters, header, init_cloud_physics, init_particles, modules, parin, read_var_list, write_var_list)

mean/min/max particle radius added as output quantity. (data_output_ptseries, modules)

Changed:

Initialisation of temporary particle array for resorting removed.
(advec_particles)

particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
(advec_particles, data_output_ptseries, modules, init_particles, particle_boundary_conds)

routine wang_kernel and respective module renamed lpm_collision_kernels.
Package (particle) parameters wang_collision_kernel and turbulence_effects_on_collision
replaced by parameter collision_kernel.
(Makefile, advec_particles, check_parameters, diffusion_e, init_3d_model, modules, package_parin, time_integration, new: lpm_collision_kernels, deleted: wang_kernel)

Errors:

Property svn:keywords set to Id

File size: 184.6 KB

Rev	Line
[1]	1	SUBROUTINE advec_particles
	2
	3	!------------------------------------------------------------------------------!
[230]	4	! Current revisions:
[759]	5	! ------------------
[825]	6	! droplet growth by condensation may include curvature and solution effects,
	7	! initialisation of temporary particle array for resorting removed,
	8	! particle attributes speed_x\|y\|z_sgs renamed rvar1\|2\|3,
	9	! module wang_kernel_mod renamed lpm_collision_kernels_mod,
	10	! wang_collision_kernel renamed wang_kernel
[482]	11	!
	12	! Former revisions:
	13	! -----------------
	14	! $Id: advec_particles.f90 825 2012-02-19 03:03:44Z raasch $
	15	!
[800]	16	! 799 2011-12-21 17:48:03Z franke
	17	! Implementation of Wang collision kernel and corresponding new parameter
	18	! wang_collision_kernel
	19	!
[799]	20	! 792 2011-12-01 raasch
	21	! particle arrays (particles, particles_temp) implemented as pointers in
	22	! order to speed up sorting (see routine sort_particles)
	23	!
[760]	24	! 759 2011-09-15 13:58:31Z raasch
	25	! Splitting of parallel I/O (routine write_particles)
	26	!
[668]	27	! 667 2010-12-23 12:06:00Z suehring/gryschka
	28	! Declaration of de_dx, de_dy, de_dz adapted to additional
	29	! ghost points. Furthermore the calls of exchange_horiz were modified.
	30	!
[623]	31	! 622 2010-12-10 08:08:13Z raasch
	32	! optional barriers included in order to speed up collective operations
	33	!
[520]	34	! 519 2010-03-19 05:30:02Z raasch
	35	! NetCDF4 output format allows size of particle array to be extended
	36	!
[482]	37	! 420 2010-01-13 15:10:53Z franke
[420]	38	! Own weighting factor for every cloud droplet is implemented and condensation
	39	! and collision of cloud droplets are adjusted accordingly. +delta_v, -s_r3,
	40	! -s_r4
	41	! Initialization of variables for the (sub-) timestep is moved to the beginning
	42	! in order to enable deletion of cloud droplets due to collision processes.
	43	! Collision efficiency for large cloud droplets has changed according to table
	44	! of Rogers and Yau. (collision_efficiency)
	45	! Bugfix: calculation of cloud droplet velocity
	46	! Bugfix: transfer of particles at south/left edge when new position
	47	! y>=(ny+0.5)dy-1.e-12 or x>=(nx+0.5)dx-1.e-12, very rare
	48	! Bugfix: calculation of y (collision_efficiency)
[392]	49	!
	50	! 336 2009-06-10 11:19:35Z raasch
[264]	51	! Particle attributes are set with new routine set_particle_attributes.
[336]	52	! Vertical particle advection depends on the particle group.
[263]	53	! Output of NetCDF messages with aid of message handling routine.
[230]	54	! Output of messages replaced by message handling routine
[262]	55	! Bugfix: error in check, if particles moved further than one subdomain length.
	56	! This check must not be applied for newly released particles
[229]	57	! Bugfix: several tail counters are initialized, particle_tail_coordinates is
	58	! only written to file if its third index is > 0
[139]	59	!
[226]	60	! 212 2008-11-11 09:09:24Z raasch
	61	! Bugfix in calculating k index in case of oceans runs (sort_particles)
	62	!
[198]	63	! 150 2008-02-29 08:19:58Z raasch
	64	! Bottom boundary condition and vertical index calculations adjusted for
	65	! ocean runs.
	66	!
[139]	67	! 119 2007-10-17 10:27:13Z raasch
[116]	68	! Sorting of particles is controlled by dt_sort_particles and moved from
	69	! the SGS timestep loop after the end of this loop.
[114]	70	! Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
	71	! numbers along y
[119]	72	! Small bugfixes in the SGS part
[77]	73	!
[110]	74	! 106 2007-08-16 14:30:26Z raasch
	75	! remaining variables iran changed to iran_part
	76	!
[98]	77	! 95 2007-06-02 16:48:38Z raasch
	78	! hydro_press renamed hyp
	79	!
[77]	80	! 75 2007-03-22 09:54:05Z raasch
[59]	81	! Particle reflection at vertical walls implemented in new subroutine
	82	! particle_boundary_conds,
	83	! vertical walls are regarded in the SGS model,
[64]	84	! + user_advec_particles, particles-package is now part of the defaut code,
	85	! array arguments in sendrecv calls have to refer to first element (1) due to
[75]	86	! mpich (mpiI) interface requirements,
	87	! 2nd+3rd argument removed from exchange horiz
[1]	88	!
[39]	89	! 16 2007-02-15 13:16:47Z raasch
	90	! Bugfix: wrong if-clause from revision 1.32
	91	!
	92	! r4 \| raasch \| 2007-02-13 12:33:16 +0100 (Tue, 13 Feb 2007)
[3]	93	! RCS Log replace by Id keyword, revision history cleaned up
	94	!
[1]	95	! Revision 1.32 2007/02/11 12:48:20 raasch
	96	! Allways the lower level k is used for interpolation
	97	! Bugfix: new particles are released only if end_time_prel > simulated_time
	98	! Bugfix: transfer of particles when x < -0.5*dx (0.0 before), etc.,
	99	! index i,j used instead of cartesian (x,y) coordinate to check for
	100	! transfer because this failed under very rare conditions
	101	! Bugfix: calculation of number of particles with same radius as the current
	102	! particle (cloud droplet code)
	103	!
	104	! Revision 1.31 2006/08/17 09:21:01 raasch
	105	! Two more compilation errors removed from the last revision
	106	!
	107	! Revision 1.30 2006/08/17 09:11:17 raasch
	108	! Two compilation errors removed from the last revision
	109	!
	110	! Revision 1.29 2006/08/04 14:05:01 raasch
	111	! Subgrid scale velocities are (optionally) included for calculating the
	112	! particle advection, new counters trlp_count_sum, etc. for accumulating
	113	! the number of particles exchanged between the subdomains during all
	114	! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z,
	115	! izuf renamed iran, output of particle time series
	116	!
	117	! Revision 1.1 1999/11/25 16:16:06 raasch
	118	! Initial revision
	119	!
	120	!
	121	! Description:
	122	! ------------
	123	! Particle advection
	124	!------------------------------------------------------------------------------!
	125
	126	USE arrays_3d
	127	USE cloud_parameters
	128	USE constants
	129	USE control_parameters
	130	USE cpulog
	131	USE grid_variables
	132	USE indices
	133	USE interfaces
[825]	134	USE lpm_collision_kernels_mod
[1]	135	USE netcdf_control
	136	USE particle_attributes
	137	USE pegrid
	138	USE random_function_mod
	139	USE statistics
	140
	141	IMPLICIT NONE
	142
[824]	143	INTEGER :: agp, deleted_particles, deleted_tails, i, ie, ii, inc, &
	144	internal_timestep_count, is, j, &
	145	jj, js, jtry, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi, &
[799]	146	tlength, trlp_count, trlp_count_sum, trlp_count_recv, &
[60]	147	trlp_count_recv_sum, trlpt_count, trlpt_count_recv, &
	148	trnp_count, trnp_count_sum, trnp_count_recv, &
	149	trnp_count_recv_sum, trnpt_count, trnpt_count_recv, &
	150	trrp_count, trrp_count_sum, trrp_count_recv, &
[1]	151	trrp_count_recv_sum, trrpt_count, trrpt_count_recv, &
	152	trsp_count, trsp_count_sum, trsp_count_recv, &
[799]	153	trsp_count_recv_sum, trspt_count, trspt_count_recv
[1]	154
[60]	155	INTEGER :: gp_outside_of_building(1:8)
	156
[824]	157	INTEGER, PARAMETER :: maxtry = 40 ! for Rosenbrock method
	158
[1]	159	LOGICAL :: dt_3d_reached, dt_3d_reached_l, prt_position
	160
[824]	161	REAL :: aa, afactor, arg, bb, cc, dd, ddenom, delta_r, delta_v, &
	162	dens_ratio, de_dt, de_dt_min, de_dx_int, de_dx_int_l, &
	163	de_dx_int_u, de_dy_int, de_dy_int_l, de_dy_int_u, de_dz_int, &
	164	de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, diss_int_u, &
	165	distance, drdt, drdt_ini, drdt_m, dt_ros, dt_ros_last, &
	166	dt_ros_next, dt_ros_sum, dt_ros_sum_ini, d2rdt2, d2rdtdr, &
	167	dt_gap, dt_particle, dt_particle_m, d_sum, e_a, e_int, &
	168	e_int_l, e_int_u, e_mean_int, e_s, err_ros, errmax, exp_arg, &
	169	exp_term, fs_int, gg, g1, g2, g3, g4, integral, &
	170	lagr_timescale, lw_max, mean_r, new_r, p_int, pt_int, &
	171	pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, &
	172	ql_int_u, random_gauss, r_ros, r_ros_ini, sigma, sl_r3, sl_r4, &
	173	t_int, u_int, u_int_l, u_int_u,vv_int, v_int, v_int_l, &
	174	v_int_u, w_int, w_int_l, w_int_u, x, y
	175	!
	176	!-- Parameters for Rosenbrock method
	177	REAL, PARAMETER :: a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0, &
	178	a2x = 1.0, a3x = 3.0/5.0, b1 = 19.0/9.0, b2 = 0.5, &
	179	b3 = 25.0/108.0, b4 = 125.0/108.0, c21 = -8.0, &
	180	c31 = 372.0/25.0, c32 = 12.0/5.0, &
	181	c41 = -112.0/125.0, c42 = -54.0/125.0, &
	182	c43 = -2.0/5.0, c1x = 0.5, c2x= -3.0/2.0, &
	183	c3x = 121.0/50.0, c4x = 29.0/250.0, &
	184	errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0, &
	185	e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, &
	186	pgrow = -0.25, pshrnk = -1.0/3.0, shrnk = 0.5
[1]	187
[60]	188	REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
[799]	189	REAL, DIMENSION(:,:), ALLOCATABLE :: kern
[60]	190
	191	REAL :: location(1:30,1:3)
	192
[667]	193	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: de_dx, de_dy, de_dz
[1]	194
	195	REAL, DIMENSION(:,:,:), ALLOCATABLE :: trlpt, trnpt, trrpt, trspt
	196
	197	TYPE(particle_type) :: tmp_particle
	198
	199	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trlp, trnp, trrp, trsp
	200
	201
	202	CALL cpu_log( log_point(25), 'advec_particles', 'start' )
	203
[825]	204
[1]	205	!
	206	!-- Write particle data on file for later analysis.
	207	!-- This has to be done here (before particles are advected) in order
	208	!-- to allow correct output in case of dt_write_particle_data = dt_prel =
	209	!-- particle_maximum_age. Otherwise (if output is done at the end of this
	210	!-- subroutine), the relevant particles would have been already deleted.
	211	!-- The MOD function allows for changes in the output interval with restart
	212	!-- runs.
	213	!-- Attention: change version number for unit 85 (in routine check_open)
	214	!-- whenever the output format for this unit is changed!
	215	time_write_particle_data = time_write_particle_data + dt_3d
	216	IF ( time_write_particle_data >= dt_write_particle_data ) THEN
	217
	218	CALL cpu_log( log_point_s(40), 'advec_part_io', 'start' )
	219	CALL check_open( 85 )
	220	WRITE ( 85 ) simulated_time, maximum_number_of_particles, &
	221	number_of_particles
	222	WRITE ( 85 ) particles
	223	WRITE ( 85 ) maximum_number_of_tailpoints, maximum_number_of_tails, &
	224	number_of_tails
[229]	225	IF ( maximum_number_of_tails > 0 ) THEN
	226	WRITE ( 85 ) particle_tail_coordinates
	227	ENDIF
[1]	228	CALL close_file( 85 )
	229
	230	IF ( netcdf_output ) CALL output_particles_netcdf
	231
	232	time_write_particle_data = MOD( time_write_particle_data, &
	233	MAX( dt_write_particle_data, dt_3d ) )
	234	CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' )
	235	ENDIF
	236
	237	!
[420]	238	!-- Initialize variables for the (sub-) timestep, i.e. for
	239	!-- marking those particles to be deleted after the timestep
	240	particle_mask = .TRUE.
	241	deleted_particles = 0
	242	trlp_count_recv = 0
	243	trnp_count_recv = 0
	244	trrp_count_recv = 0
	245	trsp_count_recv = 0
	246	trlpt_count_recv = 0
	247	trnpt_count_recv = 0
	248	trrpt_count_recv = 0
	249	trspt_count_recv = 0
	250	IF ( use_particle_tails ) THEN
	251	tail_mask = .TRUE.
	252	ENDIF
	253	deleted_tails = 0
	254
	255	!
[1]	256	!-- Calculate exponential term used in case of particle inertia for each
	257	!-- of the particle groups
	258	CALL cpu_log( log_point_s(41), 'advec_part_exp', 'start' )
	259	DO m = 1, number_of_particle_groups
	260	IF ( particle_groups(m)%density_ratio /= 0.0 ) THEN
	261	particle_groups(m)%exp_arg = &
	262	4.5 * particle_groups(m)%density_ratio * &
	263	molecular_viscosity / ( particle_groups(m)%radius )**2
	264	particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
	265	dt_3d )
	266	ENDIF
	267	ENDDO
	268	CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' )
	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
[420]	277	delta_r=0.0; delta_v=0.0
[1]	278
	279	!
	280	!-- Particle growth by condensation/evaporation
	281	CALL cpu_log( log_point_s(42), 'advec_part_cond', 'start' )
	282	DO n = 1, number_of_particles
	283	!
	284	!-- Interpolate temperature and humidity.
	285	!-- First determine left, south, and bottom index of the arrays.
	286	i = particles(n)%x * ddx
	287	j = particles(n)%y * ddy
[150]	288	k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz &
	289	+ offset_ocean_nzt ! only exact if equidistant
[1]	290
	291	x = particles(n)%x - i * dx
	292	y = particles(n)%y - j * dy
	293	aa = x2 + y2
	294	bb = ( dx - x )2 + y2
	295	cc = x2 + ( dy - y )2
	296	dd = ( dx - x )2 + ( dy - y )2
	297	gg = aa + bb + cc + dd
	298
	299	pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) &
	300	+ ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
	301	) / ( 3.0 * gg )
	302
	303	pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) &
	304	+ ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
	305	) / ( 3.0 * gg )
	306
	307	pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
	308	( pt_int_u - pt_int_l )
	309
	310	q_int_l = ( ( gg - aa ) * q(k,j,i) + ( gg - bb ) * q(k,j,i+1) &
	311	+ ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) &
	312	) / ( 3.0 * gg )
	313
	314	q_int_u = ( ( gg-aa ) * q(k+1,j,i) + ( gg-bb ) * q(k+1,j,i+1) &
	315	+ ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) &
	316	) / ( 3.0 * gg )
	317
	318	q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * &
	319	( q_int_u - q_int_l )
	320
	321	ql_int_l = ( ( gg - aa ) * ql(k,j,i) + ( gg - bb ) * ql(k,j,i+1) &
	322	+ ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) &
	323	) / ( 3.0 * gg )
	324
	325	ql_int_u = ( ( gg-aa ) * ql(k+1,j,i) + ( gg-bb ) * ql(k+1,j,i+1) &
	326	+ ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) &
	327	) / ( 3.0 * gg )
	328
	329	ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * &
	330	( ql_int_u - ql_int_l )
	331
	332	!
	333	!-- Calculate real temperature and saturation vapor pressure
[95]	334	p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) )
[1]	335	t_int = pt_int * ( p_int / 100000.0 )**0.286
	336
	337	e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) )
	338
	339	!
	340	!-- Current vapor pressure
	341	e_a = q_int * p_int / ( 0.378 * q_int + 0.622 )
	342
	343	!
	344	!-- Change in radius by condensation/evaporation
[825]	345	IF ( particles(n)%radius >= 1.0E-6 .OR. &
	346	.NOT. curvature_solution_effects ) THEN
	347	!
	348	!-- Approximation for large radii, where curvature and solution
	349	!-- effects can be neglected
	350	arg = particles(n)%radius*2 + 2.0 dt_3d * &
	351	( e_a / e_s - 1.0 ) / &
	352	( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
	353	thermal_conductivity_l + &
	354	rho_l * r_v * t_int / diff_coeff_l / e_s )
	355	IF ( arg < 1.0E-14 ) THEN
	356	new_r = 1.0E-7
	357	ELSE
	358	new_r = SQRT( arg )
	359	ENDIF
	360	ENDIF
	361
	362	IF ( curvature_solution_effects .AND. &
	363	( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) &
[824]	364	THEN
	365	!
	366	!-- Curvature and solutions effects are included in growth equation.
	367	!-- Change in Radius is calculated with 4th-order Rosenbrock method
	368	!-- for stiff o.d.e's with monitoring local truncation error to adjust
	369	!-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
	370	!-- For larger radii the simple analytic method (see ELSE) gives
	371	!-- almost the same results.
	372	!
	373	!-- Surface tension after (Straka, 2009)
	374	sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
	375
	376	r_ros = particles(n)%radius
	377	dt_ros_sum = 0.0 ! internal integrated time (s)
	378	internal_timestep_count = 0
	379
	380	ddenom = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) + &
	381	( l_v / ( r_v * t_int ) - 1.0 ) * &
	382	rho_l * l_v / ( thermal_conductivity_l * t_int )&
	383	)
	384
	385	afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
	386
	387	IF ( particles(n)%rvar3 == -9999999.9 ) THEN
	388	!
	389	!-- First particle timestep. Derivative has to be calculated.
	390	drdt_m = ddenom / r_ros * ( e_a / e_s - 1.0 - &
	391	afactor / r_ros + &
	392	bfactor / r_ros**3 )
	393	ELSE
	394	!
	395	!-- Take value from last PALM timestep
	396	drdt_m = particles(n)%rvar3
	397	ENDIF
	398	!
	399	!-- Take internal timestep values from the end of last PALM timestep
	400	dt_ros_last = particles(n)%rvar1
	401	dt_ros_next = particles(n)%rvar2
	402	!
	403	!-- Internal timestep must not be larger than PALM timestep
	404	dt_ros = MIN( dt_ros_next, dt_3d )
	405	!
	406	!-- Integrate growth equation in time unless PALM timestep is reached
	407	DO WHILE ( dt_ros_sum < dt_3d )
	408
	409	internal_timestep_count = internal_timestep_count + 1
	410
	411	!
	412	!-- Derivative at starting value
	413	drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
	414	bfactor / r_ros**3 )
	415	drdt_ini = drdt
	416	dt_ros_sum_ini = dt_ros_sum
	417	r_ros_ini = r_ros
	418
	419	!
	420	!-- Calculate time derivative of dr/dt
	421	d2rdt2 = ( drdt - drdt_m ) / dt_ros_last
	422
	423	!
	424	!-- Calculate radial derivative of dr/dt
	425	d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
	426	2.0 * afactor / r_ros**3 - &
	427	4.0 * bfactor / r_ros**5 )
	428	!
	429	!-- Adjust stepsize unless required accuracy is reached
	430	DO jtry = 1, maxtry+1
	431
	432	IF ( jtry == maxtry+1 ) THEN
	433	message_string = 'maxtry > 40 in Rosenbrock method'
	434	CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 )
	435	ENDIF
	436
	437	aa = 1.0 / ( gam * dt_ros ) - d2rdtdr
	438	g1 = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
	439	r_ros = r_ros_ini + a21 * g1
	440	drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - &
	441	afactor / r_ros + &
	442	bfactor / r_ros**3 )
	443
	444	g2 = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
	445	/ aa
	446	r_ros = r_ros_ini + a31 * g1 + a32 * g2
	447	drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - &
	448	afactor / r_ros + &
	449	bfactor / r_ros**3 )
	450
	451	g3 = ( drdt + dt_ros * c3x * d2rdt2 + &
	452	( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
	453	g4 = ( drdt + dt_ros * c4x * d2rdt2 + &
	454	( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
	455	r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
	456
	457	dt_ros_sum = dt_ros_sum_ini + dt_ros
	458
	459	IF ( dt_ros_sum == dt_ros_sum_ini ) THEN
	460	message_string = 'zero stepsize in Rosenbrock method'
	461	CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 )
	462	ENDIF
	463	!
	464	!-- Calculate error
	465	err_ros = e1g1 + e2g2 + e3g3 + e4g4
	466	errmax = 0.0
	467	errmax = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
	468	!
	469	!-- Leave loop if accuracy is sufficient, otherwise try again
	470	!-- with a reduced stepsize
	471	IF ( errmax <= 1.0 ) THEN
	472	EXIT
	473	ELSE
	474	dt_ros = SIGN( &
	475	MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
	476	shrnk * ABS( dt_ros ) ), &
	477	dt_ros &
	478	)
	479	ENDIF
	480
	481	ENDDO ! loop for stepsize adjustment
	482
	483	!
	484	!-- Calculate next internal timestep
	485	IF ( errmax > errcon ) THEN
	486	dt_ros_next = 0.9 * dt_ros * errmax**pgrow
	487	ELSE
	488	dt_ros_next = grow * dt_ros
	489	ENDIF
	490
	491	!
	492	!-- Estimated timestep is reduced if the PALM time step is exceeded
	493	dt_ros_last = dt_ros
	494	IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d ) THEN
	495	dt_ros = dt_3d - dt_ros_sum
	496	ELSE
	497	dt_ros = dt_ros_next
	498	ENDIF
	499
	500	drdt_m = drdt
	501
	502	ENDDO
	503	!
	504	!-- Store derivative and internal timestep values for next PALM step
	505	particles(n)%rvar1 = dt_ros_last
	506	particles(n)%rvar2 = dt_ros_next
	507	particles(n)%rvar3 = drdt_m
	508
	509	new_r = r_ros
	510	!
	511	!-- Radius should not fall below 1E-8 because Rosenbrock method may
	512	!-- lead to errors otherwise
	513	new_r = MAX( new_r, 1.0E-8 )
	514
[1]	515	ENDIF
	516
	517	delta_r = new_r - particles(n)%radius
	518
	519	!
	520	!-- Sum up the change in volume of liquid water for the respective grid
	521	!-- volume (this is needed later on for calculating the release of
	522	!-- latent heat)
	523	i = ( particles(n)%x + 0.5 * dx ) * ddx
	524	j = ( particles(n)%y + 0.5 * dy ) * ddy
[150]	525	k = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
	526	! only exact if equidistant
[1]	527
	528	ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * &
	529	rho_l * 1.33333333 * pi * &
	530	( new_r3 - particles(n)%radius3 ) / &
	531	( rho_surface * dx * dy * dz )
	532	IF ( ql_c(k,j,i) > 100.0 ) THEN
[274]	533	WRITE( message_string, * ) 'k=',k,' j=',j,' i=',i, &
[230]	534	' ql_c=',ql_c(k,j,i), ' &part(',n,')%wf=', &
[1]	535	particles(n)%weight_factor,' delta_r=',delta_r
[276]	536	CALL message( 'advec_particles', 'PA0143', 2, 2, -1, 6, 1 )
[1]	537	ENDIF
	538
	539	!
	540	!-- Change the droplet radius
	541	IF ( ( new_r - particles(n)%radius ) < 0.0 .AND. new_r < 0.0 ) &
	542	THEN
[274]	543	WRITE( message_string, * ) '#1 k=',k,' j=',j,' i=',i, &
	544	' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, &
[824]	545	' &delta_r=',delta_r, &
[1]	546	' particle_radius=',particles(n)%radius
[276]	547	CALL message( 'advec_particles', 'PA0144', 2, 2, -1, 6, 1 )
[1]	548	ENDIF
	549	particles(n)%radius = new_r
	550
	551	!
	552	!-- Sum up the total volume of liquid water (needed below for
	553	!-- re-calculating the weighting factors)
[420]	554	ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
[1]	555	particles(n)%radius**3
	556	ENDDO
	557	CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' )
	558
	559	!
	560	!-- Particle growth by collision
[825]	561	IF ( collision_kernel /= 'none' ) THEN
	562
[1]	563	CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
[420]	564
[1]	565	DO i = nxl, nxr
	566	DO j = nys, nyn
	567	DO k = nzb+1, nzt
	568	!
	569	!-- Collision requires at least two particles in the box
	570	IF ( prt_count(k,j,i) > 1 ) THEN
	571	!
	572	!-- First, sort particles within the gridbox by their size,
	573	!-- using Shell's method (see Numerical Recipes)
	574	!-- NOTE: In case of using particle tails, the re-sorting of
	575	!-- ---- tails would have to be included here!
	576	psi = prt_start_index(k,j,i) - 1
	577	inc = 1
	578	DO WHILE ( inc <= prt_count(k,j,i) )
	579	inc = 3 * inc + 1
	580	ENDDO
	581
	582	DO WHILE ( inc > 1 )
	583	inc = inc / 3
	584	DO is = inc+1, prt_count(k,j,i)
	585	tmp_particle = particles(psi+is)
	586	js = is
[420]	587	DO WHILE ( particles(psi+js-inc)%radius > &
[1]	588	tmp_particle%radius )
	589	particles(psi+js) = particles(psi+js-inc)
	590	js = js - inc
	591	IF ( js <= inc ) EXIT
	592	ENDDO
	593	particles(psi+js) = tmp_particle
	594	ENDDO
	595	ENDDO
	596
[799]	597	psi = prt_start_index(k,j,i)
	598	pse = psi + prt_count(k,j,i)-1
	599
[825]	600	IF ( hall_kernel .OR. wang_kernel ) THEN
[799]	601
[825]	602	ALLOCATE( kern(psi:pse,psi:pse) )
	603
[1]	604	!
[825]	605	!-- Calculate collision kernel for all particles in grid box
	606	CALL colker( i, j, k, kern )
[799]	607	!
[825]	608	!-- Collison kernel is calculated in cm**3/s but needed
	609	!-- in m**3/s
	610	kern = kern * 1.0E-6 ! to be moved to colker
[799]	611
[825]	612	DO n = pse, psi+1, -1
[799]	613
[825]	614	integral = 0.0
	615	lw_max = 0.0
[799]	616
	617	!
[825]	618	!-- Calculate growth of collector particle radius using all
	619	!-- droplets smaller than current droplet
	620	DO is = psi, n-1
[799]	621
[825]	622	integral = integral + &
	623	( particles(is)%radius*3 kern(n,is) *&
	624	particles(is)%weight_factor )
[799]	625	!
[825]	626	!-- Calculate volume of liquid water of the collected
	627	!-- droplets which is the maximum liquid water available
	628	!-- for droplet growth
	629	lw_max = lw_max + ( particles(is)%radius*3 &
	630	particles(is)%weight_factor )
[799]	631
[825]	632	ENDDO
[799]	633
	634	!
[825]	635	!-- Change in radius of collector droplet due to collision
	636	delta_r = 1.0 / ( 3.0 * particles(n)%radius*2 ) &
	637	integral * dt_3d * ddx * ddy / dz
[799]	638
	639	!
[825]	640	!-- Change in volume of collector droplet due to collision
	641	delta_v = particles(n)%weight_factor &
	642	* ( ( particles(n)%radius + delta_r )**3 &
	643	- particles(n)%radius**3 )
[799]	644
[825]	645	IF ( lw_max < delta_v .AND. delta_v > 0.0 ) THEN
	646	!-- replace by message call
	647	PRINT*, 'Particle has grown to much because the', &
	648	' change of volume of particle is larger', &
	649	' than liquid water available!'
[799]	650
[825]	651	ELSEIF ( lw_max == delta_v .AND. delta_v > 0.0 ) THEN
[799]	652
[825]	653	DO is = psi, n-1
[799]	654
[825]	655	particles(is)%weight_factor = 0.0
	656	particle_mask(is) = .FALSE.
	657	deleted_particles = deleted_particles + 1
[799]	658
[825]	659	ENDDO
[799]	660
[825]	661	ELSEIF ( lw_max > delta_v .AND. delta_v > 0.0 ) THEN
[799]	662	!
[825]	663	!-- Calculate new weighting factor of collected droplets
	664	DO is = psi, n-1
[799]	665
[825]	666	particles(is)%weight_factor = &
	667	particles(is)%weight_factor - &
	668	( ( kern(n,is) * particles(is)%weight_factor &
	669	/ integral ) * delta_v )
[799]	670
[825]	671	IF ( particles(is)%weight_factor < 0.0 ) THEN
	672	WRITE( message_string, * ) 'negative ', &
	673	'weighting factor: ', &
[799]	674	particles(is)%weight_factor
[825]	675	CALL message( 'advec_particles', '', 2, 2,&
	676	-1, 6, 1 )
[799]	677
[825]	678	ELSEIF ( particles(is)%weight_factor == 0.0 ) &
	679	THEN
[799]	680
[825]	681	particles(is)%weight_factor = 0.0
	682	particle_mask(is) = .FALSE.
	683	deleted_particles = deleted_particles + 1
[799]	684
[825]	685	ENDIF
[799]	686
[825]	687	ENDDO
[799]	688
[825]	689	ENDIF
[799]	690
[825]	691	particles(n)%radius = particles(n)%radius + delta_r
	692	ql_vp(k,j,i) = ql_vp(k,j,i) + &
	693	particles(n)%weight_factor &
[799]	694	* particles(n)%radius**3
	695
[825]	696	ENDDO
[799]	697
[825]	698	DEALLOCATE( kern )
[799]	699
[825]	700	ELSEIF ( palm_kernel ) THEN
[799]	701	!
[825]	702	!-- PALM collision kernel
	703	!
	704	!-- Calculate the mean radius of all those particles which
	705	!-- are of smaller size than the current particle and
	706	!-- use this radius for calculating the collision efficiency
	707	DO n = psi+prt_count(k,j,i)-1, psi+1, -1
[799]	708
[825]	709	sl_r3 = 0.0
	710	sl_r4 = 0.0
[799]	711
[825]	712	DO is = n-1, psi, -1
	713	IF ( particles(is)%radius < particles(n)%radius ) &
	714	THEN
	715	sl_r3 = sl_r3 + particles(is)%weight_factor &
	716	( particles(is)%radius*3 )
	717	sl_r4 = sl_r4 + particles(is)%weight_factor &
	718	( particles(is)%radius*4 )
	719	ENDIF
[799]	720	ENDDO
	721
[825]	722	IF ( ( sl_r3 ) > 0.0 ) THEN
	723	mean_r = ( sl_r4 ) / ( sl_r3 )
[420]	724
[825]	725	CALL collision_efficiency( mean_r, &
[420]	726	particles(n)%radius, &
[1]	727	effective_coll_efficiency )
	728
[825]	729	ELSE
	730	effective_coll_efficiency = 0.0
	731	ENDIF
[1]	732
[825]	733	IF ( effective_coll_efficiency > 1.0 .OR. &
	734	effective_coll_efficiency < 0.0 ) &
	735	THEN
	736	WRITE( message_string, * ) 'collision_efficien' , &
	737	'cy out of range:' ,effective_coll_efficiency
	738	CALL message( 'advec_particles', 'PA0145', 2, 2, &
	739	-1, 6, 1 )
	740	ENDIF
[1]	741
	742	!
[825]	743	!-- Interpolation of ...
	744	ii = particles(n)%x * ddx
	745	jj = particles(n)%y * ddy
	746	kk = ( particles(n)%z + 0.5 * dz ) / dz
[1]	747
[825]	748	x = particles(n)%x - ii * dx
	749	y = particles(n)%y - jj * dy
	750	aa = x2 + y2
	751	bb = ( dx - x )2 + y2
	752	cc = x2 + ( dy - y )2
	753	dd = ( dx - x )2 + ( dy - y )2
	754	gg = aa + bb + cc + dd
[1]	755
[825]	756	ql_int_l = ( (gg-aa) * ql(kk,jj,ii) + (gg-bb) * &
[1]	757	ql(kk,jj,ii+1) &
[825]	758	+ (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) * &
[1]	759	ql(kk,jj+1,ii+1) &
[825]	760	) / ( 3.0 * gg )
[1]	761
[825]	762	ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii) + (gg-bb) * &
[1]	763	ql(kk+1,jj,ii+1) &
[825]	764	+ (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) * &
[1]	765	ql(kk+1,jj+1,ii+1) &
[825]	766	) / ( 3.0 * gg )
[1]	767
[825]	768	ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
	769	( ql_int_u - ql_int_l )
[1]	770
	771	!
[825]	772	!-- Interpolate u velocity-component
	773	ii = ( particles(n)%x + 0.5 * dx ) * ddx
	774	jj = particles(n)%y * ddy
	775	kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist
[1]	776
[825]	777	IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1
[1]	778
[825]	779	x = particles(n)%x + ( 0.5 - ii ) * dx
	780	y = particles(n)%y - jj * dy
	781	aa = x2 + y2
	782	bb = ( dx - x )2 + y2
	783	cc = x2 + ( dy - y )2
	784	dd = ( dx - x )2 + ( dy - y )2
	785	gg = aa + bb + cc + dd
[1]	786
[825]	787	u_int_l = ( (gg-aa) * u(kk,jj,ii) + (gg-bb) * &
[420]	788	u(kk,jj,ii+1) &
[825]	789	+ (gg-cc) * u(kk,jj+1,ii) + (gg-dd) * &
[420]	790	u(kk,jj+1,ii+1) &
[825]	791	) / ( 3.0 * gg ) - u_gtrans
	792	IF ( kk+1 == nzt+1 ) THEN
	793	u_int = u_int_l
	794	ELSE
	795	u_int_u = ( (gg-aa) * u(kk+1,jj,ii) + (gg-bb) * &
[1]	796	u(kk+1,jj,ii+1) &
[825]	797	+ (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) * &
[1]	798	u(kk+1,jj+1,ii+1) &
[825]	799	) / ( 3.0 * gg ) - u_gtrans
	800	u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
	801	* ( u_int_u - u_int_l )
	802	ENDIF
[1]	803
	804	!
[825]	805	!-- Same procedure for interpolation of the v velocity-com-
	806	!-- ponent (adopt index k from u velocity-component)
	807	ii = particles(n)%x * ddx
	808	jj = ( particles(n)%y + 0.5 * dy ) * ddy
[1]	809
[825]	810	x = particles(n)%x - ii * dx
	811	y = particles(n)%y + ( 0.5 - jj ) * dy
	812	aa = x2 + y2
	813	bb = ( dx - x )2 + y2
	814	cc = x2 + ( dy - y )2
	815	dd = ( dx - x )2 + ( dy - y )2
	816	gg = aa + bb + cc + dd
[1]	817
[825]	818	v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * &
[1]	819	v(kk,jj,ii+1) &
[825]	820	+ ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * &
[1]	821	v(kk,jj+1,ii+1) &
[825]	822	) / ( 3.0 * gg ) - v_gtrans
	823	IF ( kk+1 == nzt+1 ) THEN
	824	v_int = v_int_l
	825	ELSE
	826	v_int_u = ( (gg-aa) * v(kk+1,jj,ii) + (gg-bb) * &
[1]	827	v(kk+1,jj,ii+1) &
[825]	828	+ (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) * &
[1]	829	v(kk+1,jj+1,ii+1) &
[825]	830	) / ( 3.0 * gg ) - v_gtrans
	831	v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
	832	* ( v_int_u - v_int_l )
	833	ENDIF
[1]	834
	835	!
[825]	836	!-- Same procedure for interpolation of the w velocity-com-
	837	!-- ponent (adopt index i from v velocity-component)
	838	jj = particles(n)%y * ddy
	839	kk = particles(n)%z / dz
[1]	840
[825]	841	x = particles(n)%x - ii * dx
	842	y = particles(n)%y - jj * dy
	843	aa = x2 + y2
	844	bb = ( dx - x )2 + y2
	845	cc = x2 + ( dy - y )2
	846	dd = ( dx - x )2 + ( dy - y )2
	847	gg = aa + bb + cc + dd
[1]	848
[825]	849	w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * &
[1]	850	w(kk,jj,ii+1) &
[825]	851	+ ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * &
[1]	852	w(kk,jj+1,ii+1) &
[825]	853	) / ( 3.0 * gg )
	854	IF ( kk+1 == nzt+1 ) THEN
	855	w_int = w_int_l
	856	ELSE
	857	w_int_u = ( (gg-aa) * w(kk+1,jj,ii) + (gg-bb) * &
[1]	858	w(kk+1,jj,ii+1) &
[825]	859	+ (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) * &
[1]	860	w(kk+1,jj+1,ii+1) &
[825]	861	) / ( 3.0 * gg )
	862	w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
	863	* ( w_int_u - w_int_l )
	864	ENDIF
[1]	865
	866	!
[825]	867	!-- Change in radius due to collision
	868	delta_r = effective_coll_efficiency / 3.0 &
	869	* pi * sl_r3 * ddx * ddy / dz &
	870	* SQRT( ( u_int - particles(n)%speed_x )**2 &
	871	+ ( v_int - particles(n)%speed_y )**2 &
	872	+ ( w_int - particles(n)%speed_z )**2 &
	873	) * dt_3d
[420]	874	!
[825]	875	!-- Change in volume due to collision
	876	delta_v = particles(n)%weight_factor &
	877	* ( ( particles(n)%radius + delta_r )**3 &
	878	- particles(n)%radius**3 )
[1]	879
[420]	880	!
[825]	881	!-- Check if collected particles provide enough LWC for
	882	!-- volume change of collector particle
	883	IF ( delta_v >= sl_r3 .AND. sl_r3 > 0.0 ) THEN
[1]	884
[825]	885	delta_r = ( ( sl_r3/particles(n)%weight_factor ) &
	886	+ particles(n)%radius3 )( 1./3. ) &
	887	- particles(n)%radius
[1]	888
[825]	889	DO is = n-1, psi, -1
	890	IF ( particles(is)%radius < &
	891	particles(n)%radius ) THEN
	892	particles(is)%weight_factor = 0.0
	893	particle_mask(is) = .FALSE.
	894	deleted_particles = deleted_particles + 1
	895	ENDIF
	896	ENDDO
[1]	897
[825]	898	ELSE IF ( delta_v < sl_r3 .AND. sl_r3 > 0.0 ) THEN
[1]	899
[825]	900	DO is = n-1, psi, -1
	901	IF ( particles(is)%radius < particles(n)%radius &
	902	.AND. sl_r3 > 0.0 ) THEN
	903	particles(is)%weight_factor = &
[420]	904	( ( particles(is)%weight_factor &
	905	* ( particles(is)%radius**3 ) ) &
	906	- ( delta_v &
	907	* particles(is)%weight_factor &
	908	* ( particles(is)%radius**3 ) &
	909	/ sl_r3 ) ) &
	910	/ ( particles(is)%radius**3 )
[1]	911
[825]	912	IF ( particles(is)%weight_factor < 0.0 ) THEN
	913	WRITE( message_string, * ) 'negative ', &
[420]	914	'weighting factor: ', &
	915	particles(is)%weight_factor
[825]	916	CALL message( 'advec_particles', '', 2, &
	917	2, -1, 6, 1 )
	918	ENDIF
[420]	919	ENDIF
[825]	920	ENDDO
	921	ENDIF
[420]	922
[825]	923	particles(n)%radius = particles(n)%radius + delta_r
	924	ql_vp(k,j,i) = ql_vp(k,j,i) + &
	925	particles(n)%weight_factor * &
	926	( particles(n)%radius**3 )
	927	ENDDO
[1]	928
[825]	929	ENDIF ! collision kernel
	930
[420]	931	ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor &
[825]	932	* particles(psi)%radius**3
[420]	933
[825]	934
[420]	935	ELSE IF ( prt_count(k,j,i) == 1 ) THEN
[825]	936
[420]	937	psi = prt_start_index(k,j,i)
[825]	938	ql_vp(k,j,i) = particles(psi)%weight_factor * &
	939	particles(psi)%radius**3
[1]	940	ENDIF
[825]	941
[420]	942	!
[825]	943	!-- Check if condensation of LWC was conserved during collision
	944	!-- process
	945	IF ( ql_v(k,j,i) /= 0.0 ) THEN
	946	IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001 .OR. &
	947	ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 ) THEN
[420]	948	WRITE( message_string, * ) 'LWC is not conserved during',&
	949	' collision! ', &
	950	'LWC after condensation: ', &
	951	ql_v(k,j,i), &
	952	' LWC after collision: ', &
	953	ql_vp(k,j,i)
[825]	954	CALL message( 'advec_particles', '', 2, 2, -1, 6, &
	955	1 )
[420]	956	ENDIF
	957	ENDIF
[1]	958
	959	ENDDO
	960	ENDDO
	961	ENDDO
	962
[825]	963	ENDIF ! collision handling
[799]	964
[1]	965	CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
	966
[825]	967	ENDIF ! cloud droplet handling
[1]	968
	969
	970	!
	971	!-- Particle advection.
	972	!-- In case of including the SGS velocities, the LES timestep has probably
	973	!-- to be split into several smaller timesteps because of the Lagrangian
	974	!-- timescale condition. Because the number of timesteps to be carried out is
	975	!-- not known at the beginning, these steps are carried out in an infinite loop
	976	!-- with exit condition.
	977	!
	978	!-- If SGS velocities are used, gradients of the TKE have to be calculated and
	979	!-- boundary conditions have to be set first. Also, horizontally averaged
	980	!-- profiles of the SGS TKE and the resolved-scale velocity variances are
	981	!-- needed.
	982	IF ( use_sgs_for_particles ) THEN
	983
	984	!
	985	!-- TKE gradient along x and y
	986	DO i = nxl, nxr
	987	DO j = nys, nyn
	988	DO k = nzb, nzt+1
[57]	989
	990	IF ( k <= nzb_s_inner(j,i-1) .AND. &
	991	k > nzb_s_inner(j,i) .AND. &
	992	k > nzb_s_inner(j,i+1) ) THEN
	993	de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
	994	( e(k,j,i+1) - e(k,j,i) ) * ddx
	995	ELSEIF ( k > nzb_s_inner(j,i-1) .AND. &
	996	k > nzb_s_inner(j,i) .AND. &
	997	k <= nzb_s_inner(j,i+1) ) THEN
	998	de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
	999	( e(k,j,i) - e(k,j,i-1) ) * ddx
	1000	ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j,i+1) ) &
	1001	THEN
	1002	de_dx(k,j,i) = 0.0
	1003	ELSEIF ( k < nzb_s_inner(j,i-1) .AND. k < nzb_s_inner(j,i) ) &
	1004	THEN
	1005	de_dx(k,j,i) = 0.0
	1006	ELSE
	1007	de_dx(k,j,i) = sgs_wfu_part * &
	1008	( e(k,j,i+1) - e(k,j,i-1) ) * ddx
	1009	ENDIF
	1010
	1011	IF ( k <= nzb_s_inner(j-1,i) .AND. &
	1012	k > nzb_s_inner(j,i) .AND. &
	1013	k > nzb_s_inner(j+1,i) ) THEN
	1014	de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
	1015	( e(k,j+1,i) - e(k,j,i) ) * ddy
	1016	ELSEIF ( k > nzb_s_inner(j-1,i) .AND. &
	1017	k > nzb_s_inner(j,i) .AND. &
	1018	k <= nzb_s_inner(j+1,i) ) THEN
	1019	de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
	1020	( e(k,j,i) - e(k,j-1,i) ) * ddy
	1021	ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j+1,i) ) &
	1022	THEN
	1023	de_dy(k,j,i) = 0.0
	1024	ELSEIF ( k < nzb_s_inner(j-1,i) .AND. k < nzb_s_inner(j,i) ) &
	1025	THEN
	1026	de_dy(k,j,i) = 0.0
	1027	ELSE
	1028	de_dy(k,j,i) = sgs_wfv_part * &
	1029	( e(k,j+1,i) - e(k,j-1,i) ) * ddy
	1030	ENDIF
	1031
[1]	1032	ENDDO
	1033	ENDDO
	1034	ENDDO
	1035
	1036	!
	1037	!-- TKE gradient along z, including bottom and top boundary conditions
	1038	DO i = nxl, nxr
	1039	DO j = nys, nyn
[57]	1040
	1041	DO k = nzb_s_inner(j,i)+2, nzt-1
[1]	1042	de_dz(k,j,i) = 2.0 * sgs_wfw_part * &
	1043	( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
	1044	ENDDO
[57]	1045
	1046	k = nzb_s_inner(j,i)
	1047	de_dz(nzb:k,j,i) = 0.0
	1048	de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) &
	1049	/ ( zu(k+2) - zu(k+1) )
[1]	1050	de_dz(nzt,j,i) = 0.0
	1051	de_dz(nzt+1,j,i) = 0.0
	1052	ENDDO
[57]	1053	ENDDO
[1]	1054
	1055	!
	1056	!-- Lateral boundary conditions
[667]	1057	CALL exchange_horiz( de_dx, nbgp )
	1058	CALL exchange_horiz( de_dy, nbgp )
	1059	CALL exchange_horiz( de_dz, nbgp )
	1060	CALL exchange_horiz( diss, nbgp )
[1]	1061
	1062	!
	1063	!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
	1064	!-- velocity variances (they may have been already calculated in routine
	1065	!-- flow_statistics).
	1066	IF ( .NOT. flow_statistics_called ) THEN
	1067	!
	1068	!-- First calculate horizontally averaged profiles of the horizontal
	1069	!-- velocities.
	1070	sums_l(:,1,0) = 0.0
	1071	sums_l(:,2,0) = 0.0
	1072
	1073	DO i = nxl, nxr
	1074	DO j = nys, nyn
	1075	DO k = nzb_s_outer(j,i), nzt+1
	1076	sums_l(k,1,0) = sums_l(k,1,0) + u(k,j,i)
	1077	sums_l(k,2,0) = sums_l(k,2,0) + v(k,j,i)
	1078	ENDDO
	1079	ENDDO
	1080	ENDDO
	1081
	1082	#if defined( __parallel )
	1083	!
	1084	!-- Compute total sum from local sums
[622]	1085	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	1086	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
	1087	MPI_REAL, MPI_SUM, comm2d, ierr )
[622]	1088	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	1089	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
	1090	MPI_REAL, MPI_SUM, comm2d, ierr )
	1091	#else
	1092	sums(:,1) = sums_l(:,1,0)
	1093	sums(:,2) = sums_l(:,2,0)
	1094	#endif
	1095
	1096	!
	1097	!-- Final values are obtained by division by the total number of grid
[420]	1098	!-- points used for the summation.
[1]	1099	hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0) ! u
	1100	hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0) ! v
	1101
	1102	!
	1103	!-- Now calculate the profiles of SGS TKE and the resolved-scale
	1104	!-- velocity variances
	1105	sums_l(:,8,0) = 0.0
	1106	sums_l(:,30,0) = 0.0
	1107	sums_l(:,31,0) = 0.0
	1108	sums_l(:,32,0) = 0.0
	1109	DO i = nxl, nxr
	1110	DO j = nys, nyn
	1111	DO k = nzb_s_outer(j,i), nzt+1
	1112	sums_l(k,8,0) = sums_l(k,8,0) + e(k,j,i)
	1113	sums_l(k,30,0) = sums_l(k,30,0) + &
	1114	( u(k,j,i) - hom(k,1,1,0) )**2
	1115	sums_l(k,31,0) = sums_l(k,31,0) + &
	1116	( v(k,j,i) - hom(k,1,2,0) )**2
	1117	sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2
	1118	ENDDO
	1119	ENDDO
	1120	ENDDO
	1121
	1122	#if defined( __parallel )
	1123	!
	1124	!-- Compute total sum from local sums
[622]	1125	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	1126	CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
	1127	MPI_REAL, MPI_SUM, comm2d, ierr )
[622]	1128	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	1129	CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
	1130	MPI_REAL, MPI_SUM, comm2d, ierr )
[622]	1131	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	1132	CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
	1133	MPI_REAL, MPI_SUM, comm2d, ierr )
[622]	1134	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	1135	CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
	1136	MPI_REAL, MPI_SUM, comm2d, ierr )
[420]	1137
[1]	1138	#else
	1139	sums(:,8) = sums_l(:,8,0)
	1140	sums(:,30) = sums_l(:,30,0)
	1141	sums(:,31) = sums_l(:,31,0)
	1142	sums(:,32) = sums_l(:,32,0)
	1143	#endif
	1144
	1145	!
	1146	!-- Final values are obtained by division by the total number of grid
	1147	!-- points used for the summation.
	1148	hom(:,1,8,0) = sums(:,8) / ngp_2dh_outer(:,0) ! e
	1149	hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0) ! u*2
	1150	hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0) ! v*2
	1151	hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0) ! w*2
	1152
	1153	ENDIF
	1154
[420]	1155	ENDIF
[1]	1156
	1157	!
	1158	!-- Initialize variables used for accumulating the number of particles
	1159	!-- exchanged between the subdomains during all sub-timesteps (if sgs
	1160	!-- velocities are included). These data are output further below on the
	1161	!-- particle statistics file.
	1162	trlp_count_sum = 0
	1163	trlp_count_recv_sum = 0
	1164	trrp_count_sum = 0
	1165	trrp_count_recv_sum = 0
	1166	trsp_count_sum = 0
	1167	trsp_count_recv_sum = 0
	1168	trnp_count_sum = 0
	1169	trnp_count_recv_sum = 0
	1170
	1171	!
	1172	!-- Initialize the variable storing the total time that a particle has advanced
	1173	!-- within the timestep procedure
	1174	particles(1:number_of_particles)%dt_sum = 0.0
	1175
	1176	!
	1177	!-- Timestep loop.
	1178	!-- This loop has to be repeated until the advection time of every particle
	1179	!-- (in the total domain!) has reached the LES timestep (dt_3d)
	1180	DO
	1181
	1182	CALL cpu_log( log_point_s(44), 'advec_part_advec', 'start' )
	1183
	1184	!
	1185	!-- Initialize the switch used for the loop exit condition checked at the
	1186	!-- end of this loop.
	1187	!-- If at least one particle has failed to reach the LES timestep, this
	1188	!-- switch will be set false.
	1189	dt_3d_reached_l = .TRUE.
	1190
	1191	DO n = 1, number_of_particles
	1192	!
	1193	!-- Move particles only if the LES timestep has not (approximately) been
	1194	!-- reached
	1195	IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 ) CYCLE
	1196
	1197	!
	1198	!-- Interpolate u velocity-component, determine left, front, bottom
	1199	!-- index of u-array
	1200	i = ( particles(n)%x + 0.5 * dx ) * ddx
	1201	j = particles(n)%y * ddy
[150]	1202	k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz &
	1203	+ offset_ocean_nzt ! only exact if equidistant
[420]	1204
[1]	1205	!
	1206	!-- Interpolation of the velocity components in the xy-plane
	1207	x = particles(n)%x + ( 0.5 - i ) * dx
	1208	y = particles(n)%y - j * dy
	1209	aa = x2 + y2
	1210	bb = ( dx - x )2 + y2
	1211	cc = x2 + ( dy - y )2
	1212	dd = ( dx - x )2 + ( dy - y )2
	1213	gg = aa + bb + cc + dd
	1214
	1215	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
	1216	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
	1217	) / ( 3.0 * gg ) - u_gtrans
	1218	IF ( k+1 == nzt+1 ) THEN
	1219	u_int = u_int_l
	1220	ELSE
	1221	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
	1222	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
	1223	) / ( 3.0 * gg ) - u_gtrans
	1224	u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1225	( u_int_u - u_int_l )
	1226	ENDIF
	1227
	1228	!
	1229	!-- Same procedure for interpolation of the v velocity-component (adopt
	1230	!-- index k from u velocity-component)
	1231	i = particles(n)%x * ddx
	1232	j = ( particles(n)%y + 0.5 * dy ) * ddy
	1233
	1234	x = particles(n)%x - i * dx
	1235	y = particles(n)%y + ( 0.5 - j ) * dy
	1236	aa = x2 + y2
	1237	bb = ( dx - x )2 + y2
	1238	cc = x2 + ( dy - y )2
	1239	dd = ( dx - x )2 + ( dy - y )2
	1240	gg = aa + bb + cc + dd
	1241
	1242	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
	1243	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
	1244	) / ( 3.0 * gg ) - v_gtrans
	1245	IF ( k+1 == nzt+1 ) THEN
	1246	v_int = v_int_l
	1247	ELSE
	1248	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
	1249	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
	1250	) / ( 3.0 * gg ) - v_gtrans
	1251	v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1252	( v_int_u - v_int_l )
	1253	ENDIF
	1254
	1255	!
	1256	!-- Same procedure for interpolation of the w velocity-component (adopt
	1257	!-- index i from v velocity-component)
[336]	1258	IF ( vertical_particle_advection(particles(n)%group) ) THEN
[1]	1259	j = particles(n)%y * ddy
[150]	1260	k = particles(n)%z / dz + offset_ocean_nzt_m1
[1]	1261
	1262	x = particles(n)%x - i * dx
	1263	y = particles(n)%y - j * dy
	1264	aa = x2 + y2
	1265	bb = ( dx - x )2 + y2
	1266	cc = x2 + ( dy - y )2
	1267	dd = ( dx - x )2 + ( dy - y )2
	1268	gg = aa + bb + cc + dd
	1269
	1270	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
	1271	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
	1272	) / ( 3.0 * gg )
	1273	IF ( k+1 == nzt+1 ) THEN
	1274	w_int = w_int_l
	1275	ELSE
	1276	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
	1277	( gg-bb ) * w(k+1,j,i+1) + &
	1278	( gg-cc ) * w(k+1,j+1,i) + &
	1279	( gg-dd ) * w(k+1,j+1,i+1) &
	1280	) / ( 3.0 * gg )
	1281	w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
	1282	( w_int_u - w_int_l )
	1283	ENDIF
	1284	ELSE
	1285	w_int = 0.0
	1286	ENDIF
	1287
	1288	!
	1289	!-- Interpolate and calculate quantities needed for calculating the SGS
	1290	!-- velocities
	1291	IF ( use_sgs_for_particles ) THEN
	1292	!
	1293	!-- Interpolate TKE
	1294	i = particles(n)%x * ddx
	1295	j = particles(n)%y * ddy
[150]	1296	k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz &
	1297	+ offset_ocean_nzt ! only exact if eq.dist
[1]	1298
[420]	1299	IF ( topography == 'flat' ) THEN
[1]	1300
[57]	1301	x = particles(n)%x - i * dx
	1302	y = particles(n)%y - j * dy
	1303	aa = x2 + y2
	1304	bb = ( dx - x )2 + y2
	1305	cc = x2 + ( dy - y )2
	1306	dd = ( dx - x )2 + ( dy - y )2
	1307	gg = aa + bb + cc + dd
[1]	1308
[57]	1309	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
	1310	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
[1]	1311	) / ( 3.0 * gg )
	1312
[57]	1313	IF ( k+1 == nzt+1 ) THEN
	1314	e_int = e_int_l
	1315	ELSE
	1316	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	1317	( gg - bb ) * e(k+1,j,i+1) + &
	1318	( gg - cc ) * e(k+1,j+1,i) + &
	1319	( gg - dd ) * e(k+1,j+1,i+1) &
	1320	) / ( 3.0 * gg )
	1321	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1322	( e_int_u - e_int_l )
	1323	ENDIF
	1324
[1]	1325	!
[57]	1326	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
	1327	!-- all position variables from above (TKE))
	1328	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	1329	( gg - bb ) * de_dx(k,j,i+1) + &
	1330	( gg - cc ) * de_dx(k,j+1,i) + &
	1331	( gg - dd ) * de_dx(k,j+1,i+1) &
[1]	1332	) / ( 3.0 * gg )
	1333
[57]	1334	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1335	de_dx_int = de_dx_int_l
	1336	ELSE
	1337	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	1338	( gg - bb ) * de_dx(k+1,j,i+1) + &
	1339	( gg - cc ) * de_dx(k+1,j+1,i) + &
	1340	( gg - dd ) * de_dx(k+1,j+1,i+1) &
	1341	) / ( 3.0 * gg )
	1342	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1343	( de_dx_int_u - de_dx_int_l )
	1344	ENDIF
	1345
[1]	1346	!
[57]	1347	!-- Interpolate the TKE gradient along y
	1348	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	1349	( gg - bb ) * de_dy(k,j,i+1) + &
	1350	( gg - cc ) * de_dy(k,j+1,i) + &
	1351	( gg - dd ) * de_dy(k,j+1,i+1) &
[1]	1352	) / ( 3.0 * gg )
[57]	1353	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1354	de_dy_int = de_dy_int_l
	1355	ELSE
	1356	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	1357	( gg - bb ) * de_dy(k+1,j,i+1) + &
	1358	( gg - cc ) * de_dy(k+1,j+1,i) + &
	1359	( gg - dd ) * de_dy(k+1,j+1,i+1) &
	1360	) / ( 3.0 * gg )
	1361	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1362	( de_dy_int_u - de_dy_int_l )
	1363	ENDIF
[1]	1364
	1365	!
[57]	1366	!-- Interpolate the TKE gradient along z
	1367	IF ( particles(n)%z < 0.5 * dz ) THEN
	1368	de_dz_int = 0.0
	1369	ELSE
	1370	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	1371	( gg - bb ) * de_dz(k,j,i+1) + &
	1372	( gg - cc ) * de_dz(k,j+1,i) + &
	1373	( gg - dd ) * de_dz(k,j+1,i+1) &
	1374	) / ( 3.0 * gg )
[1]	1375
[57]	1376	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1377	de_dz_int = de_dz_int_l
	1378	ELSE
	1379	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	1380	( gg - bb ) * de_dz(k+1,j,i+1) + &
	1381	( gg - cc ) * de_dz(k+1,j+1,i) + &
	1382	( gg - dd ) * de_dz(k+1,j+1,i+1) &
	1383	) / ( 3.0 * gg )
	1384	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1385	( de_dz_int_u - de_dz_int_l )
	1386	ENDIF
	1387	ENDIF
[1]	1388
	1389	!
[57]	1390	!-- Interpolate the dissipation of TKE
	1391	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	1392	( gg - bb ) * diss(k,j,i+1) + &
	1393	( gg - cc ) * diss(k,j+1,i) + &
	1394	( gg - dd ) * diss(k,j+1,i+1) &
	1395	) / ( 3.0 * gg )
[1]	1396
[57]	1397	IF ( k+1 == nzt+1 ) THEN
	1398	diss_int = diss_int_l
	1399	ELSE
	1400	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	1401	( gg - bb ) * diss(k+1,j,i+1) + &
	1402	( gg - cc ) * diss(k+1,j+1,i) + &
	1403	( gg - dd ) * diss(k+1,j+1,i+1) &
	1404	) / ( 3.0 * gg )
	1405	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1406	( diss_int_u - diss_int_l )
	1407	ENDIF
	1408
[1]	1409	ELSE
	1410
	1411	!
[57]	1412	!-- In case that there are buildings it has to be determined
	1413	!-- how many of the gridpoints defining the particle box are
	1414	!-- situated within a building
	1415	!-- gp_outside_of_building(1): i,j,k
	1416	!-- gp_outside_of_building(2): i,j+1,k
	1417	!-- gp_outside_of_building(3): i,j,k+1
	1418	!-- gp_outside_of_building(4): i,j+1,k+1
	1419	!-- gp_outside_of_building(5): i+1,j,k
	1420	!-- gp_outside_of_building(6): i+1,j+1,k
	1421	!-- gp_outside_of_building(7): i+1,j,k+1
	1422	!-- gp_outside_of_building(8): i+1,j+1,k+1
[420]	1423
[57]	1424	gp_outside_of_building = 0
	1425	location = 0.0
	1426	num_gp = 0
	1427
	1428	IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	1429	num_gp = num_gp + 1
	1430	gp_outside_of_building(1) = 1
	1431	location(num_gp,1) = i * dx
[119]	1432	location(num_gp,2) = j * dy
[57]	1433	location(num_gp,3) = k * dz - 0.5 * dz
	1434	ei(num_gp) = e(k,j,i)
	1435	dissi(num_gp) = diss(k,j,i)
	1436	de_dxi(num_gp) = de_dx(k,j,i)
	1437	de_dyi(num_gp) = de_dy(k,j,i)
	1438	de_dzi(num_gp) = de_dz(k,j,i)
	1439	ENDIF
	1440
	1441	IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
	1442	THEN
	1443	num_gp = num_gp + 1
	1444	gp_outside_of_building(2) = 1
	1445	location(num_gp,1) = i * dx
[119]	1446	location(num_gp,2) = (j+1) * dy
[57]	1447	location(num_gp,3) = k * dz - 0.5 * dz
	1448	ei(num_gp) = e(k,j+1,i)
[420]	1449	dissi(num_gp) = diss(k,j+1,i)
[57]	1450	de_dxi(num_gp) = de_dx(k,j+1,i)
	1451	de_dyi(num_gp) = de_dy(k,j+1,i)
	1452	de_dzi(num_gp) = de_dz(k,j+1,i)
	1453	ENDIF
	1454
	1455	IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	1456	num_gp = num_gp + 1
	1457	gp_outside_of_building(3) = 1
	1458	location(num_gp,1) = i * dx
	1459	location(num_gp,2) = j * dy
	1460	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	1461	ei(num_gp) = e(k+1,j,i)
[420]	1462	dissi(num_gp) = diss(k+1,j,i)
[57]	1463	de_dxi(num_gp) = de_dx(k+1,j,i)
	1464	de_dyi(num_gp) = de_dy(k+1,j,i)
	1465	de_dzi(num_gp) = de_dz(k+1,j,i)
	1466	ENDIF
	1467
	1468	IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
	1469	THEN
	1470	num_gp = num_gp + 1
	1471	gp_outside_of_building(4) = 1
	1472	location(num_gp,1) = i * dx
	1473	location(num_gp,2) = (j+1) * dy
	1474	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	1475	ei(num_gp) = e(k+1,j+1,i)
[420]	1476	dissi(num_gp) = diss(k+1,j+1,i)
[57]	1477	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	1478	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	1479	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	1480	ENDIF
[799]	1481
[57]	1482	IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
	1483	THEN
	1484	num_gp = num_gp + 1
	1485	gp_outside_of_building(5) = 1
	1486	location(num_gp,1) = (i+1) * dx
	1487	location(num_gp,2) = j * dy
	1488	location(num_gp,3) = k * dz - 0.5 * dz
	1489	ei(num_gp) = e(k,j,i+1)
[420]	1490	dissi(num_gp) = diss(k,j,i+1)
[57]	1491	de_dxi(num_gp) = de_dx(k,j,i+1)
	1492	de_dyi(num_gp) = de_dy(k,j,i+1)
	1493	de_dzi(num_gp) = de_dz(k,j,i+1)
	1494	ENDIF
	1495
	1496	IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
	1497	THEN
	1498	num_gp = num_gp + 1
	1499	gp_outside_of_building(6) = 1
	1500	location(num_gp,1) = (i+1) * dx
	1501	location(num_gp,2) = (j+1) * dy
	1502	location(num_gp,3) = k * dz - 0.5 * dz
	1503	ei(num_gp) = e(k,j+1,i+1)
[420]	1504	dissi(num_gp) = diss(k,j+1,i+1)
[57]	1505	de_dxi(num_gp) = de_dx(k,j+1,i+1)
	1506	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	1507	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	1508	ENDIF
	1509
	1510	IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
	1511	THEN
	1512	num_gp = num_gp + 1
	1513	gp_outside_of_building(7) = 1
	1514	location(num_gp,1) = (i+1) * dx
	1515	location(num_gp,2) = j * dy
	1516	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	1517	ei(num_gp) = e(k+1,j,i+1)
[420]	1518	dissi(num_gp) = diss(k+1,j,i+1)
[57]	1519	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	1520	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	1521	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	1522	ENDIF
	1523
	1524	IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
	1525	THEN
	1526	num_gp = num_gp + 1
	1527	gp_outside_of_building(8) = 1
	1528	location(num_gp,1) = (i+1) * dx
	1529	location(num_gp,2) = (j+1) * dy
	1530	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	1531	ei(num_gp) = e(k+1,j+1,i+1)
[420]	1532	dissi(num_gp) = diss(k+1,j+1,i+1)
[57]	1533	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	1534	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	1535	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	1536	ENDIF
	1537
	1538	!
[420]	1539	!-- If all gridpoints are situated outside of a building, then the
[57]	1540	!-- ordinary interpolation scheme can be used.
	1541	IF ( num_gp == 8 ) THEN
	1542
	1543	x = particles(n)%x - i * dx
	1544	y = particles(n)%y - j * dy
	1545	aa = x2 + y2
	1546	bb = ( dx - x )2 + y2
	1547	cc = x2 + ( dy - y )2
	1548	dd = ( dx - x )2 + ( dy - y )2
	1549	gg = aa + bb + cc + dd
	1550
	1551	e_int_l = (( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
	1552	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
	1553	) / ( 3.0 * gg )
	1554
	1555	IF ( k+1 == nzt+1 ) THEN
	1556	e_int = e_int_l
	1557	ELSE
	1558	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	1559	( gg - bb ) * e(k+1,j,i+1) + &
	1560	( gg - cc ) * e(k+1,j+1,i) + &
	1561	( gg - dd ) * e(k+1,j+1,i+1) &
	1562	) / ( 3.0 * gg )
	1563	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
	1564	( e_int_u - e_int_l )
	1565	ENDIF
	1566
	1567	!
	1568	!-- Interpolate the TKE gradient along x (adopt incides i,j,k
	1569	!-- and all position variables from above (TKE))
	1570	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	1571	( gg - bb ) * de_dx(k,j,i+1) + &
	1572	( gg - cc ) * de_dx(k,j+1,i) + &
	1573	( gg - dd ) * de_dx(k,j+1,i+1) &
	1574	) / ( 3.0 * gg )
	1575
	1576	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1577	de_dx_int = de_dx_int_l
	1578	ELSE
	1579	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	1580	( gg - bb ) * de_dx(k+1,j,i+1) + &
	1581	( gg - cc ) * de_dx(k+1,j+1,i) + &
	1582	( gg - dd ) * de_dx(k+1,j+1,i+1) &
	1583	) / ( 3.0 * gg )
	1584	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
	1585	dz * ( de_dx_int_u - de_dx_int_l )
	1586	ENDIF
	1587
	1588	!
	1589	!-- Interpolate the TKE gradient along y
	1590	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	1591	( gg - bb ) * de_dy(k,j,i+1) + &
	1592	( gg - cc ) * de_dy(k,j+1,i) + &
	1593	( gg - dd ) * de_dy(k,j+1,i+1) &
	1594	) / ( 3.0 * gg )
	1595	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1596	de_dy_int = de_dy_int_l
	1597	ELSE
	1598	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	1599	( gg - bb ) * de_dy(k+1,j,i+1) + &
	1600	( gg - cc ) * de_dy(k+1,j+1,i) + &
	1601	( gg - dd ) * de_dy(k+1,j+1,i+1) &
	1602	) / ( 3.0 * gg )
	1603	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
	1604	dz * ( de_dy_int_u - de_dy_int_l )
	1605	ENDIF
	1606
	1607	!
	1608	!-- Interpolate the TKE gradient along z
	1609	IF ( particles(n)%z < 0.5 * dz ) THEN
	1610	de_dz_int = 0.0
	1611	ELSE
	1612	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	1613	( gg - bb ) * de_dz(k,j,i+1) + &
	1614	( gg - cc ) * de_dz(k,j+1,i) + &
	1615	( gg - dd ) * de_dz(k,j+1,i+1) &
	1616	) / ( 3.0 * gg )
	1617
	1618	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	1619	de_dz_int = de_dz_int_l
	1620	ELSE
	1621	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	1622	( gg - bb ) * de_dz(k+1,j,i+1) + &
	1623	( gg - cc ) * de_dz(k+1,j+1,i) + &
	1624	( gg - dd ) * de_dz(k+1,j+1,i+1) &
	1625	) / ( 3.0 * gg )
	1626	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
	1627	dz * ( de_dz_int_u - de_dz_int_l )
	1628	ENDIF
	1629	ENDIF
	1630
	1631	!
	1632	!-- Interpolate the dissipation of TKE
	1633	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	1634	( gg - bb ) * diss(k,j,i+1) + &
	1635	( gg - cc ) * diss(k,j+1,i) + &
	1636	( gg - dd ) * diss(k,j+1,i+1) &
	1637	) / ( 3.0 * gg )
	1638
	1639	IF ( k+1 == nzt+1 ) THEN
	1640	diss_int = diss_int_l
	1641	ELSE
	1642	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	1643	( gg - bb ) * diss(k+1,j,i+1) + &
	1644	( gg - cc ) * diss(k+1,j+1,i) + &
	1645	( gg - dd ) * diss(k+1,j+1,i+1) &
	1646	) / ( 3.0 * gg )
	1647	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
	1648	( diss_int_u - diss_int_l )
	1649	ENDIF
	1650
	1651	ELSE
	1652
	1653	!
	1654	!-- If wall between gridpoint 1 and gridpoint 5, then
	1655	!-- Neumann boundary condition has to be applied
	1656	IF ( gp_outside_of_building(1) == 1 .AND. &
	1657	gp_outside_of_building(5) == 0 ) THEN
	1658	num_gp = num_gp + 1
	1659	location(num_gp,1) = i * dx + 0.5 * dx
	1660	location(num_gp,2) = j * dy
	1661	location(num_gp,3) = k * dz - 0.5 * dz
	1662	ei(num_gp) = e(k,j,i)
	1663	dissi(num_gp) = diss(k,j,i)
	1664	de_dxi(num_gp) = 0.0
	1665	de_dyi(num_gp) = de_dy(k,j,i)
[420]	1666	de_dzi(num_gp) = de_dz(k,j,i)
	1667	ENDIF
	1668
[57]	1669	IF ( gp_outside_of_building(5) == 1 .AND. &
	1670	gp_outside_of_building(1) == 0 ) THEN
	1671	num_gp = num_gp + 1
	1672	location(num_gp,1) = i * dx + 0.5 * dx
	1673	location(num_gp,2) = j * dy
	1674	location(num_gp,3) = k * dz - 0.5 * dz
	1675	ei(num_gp) = e(k,j,i+1)
	1676	dissi(num_gp) = diss(k,j,i+1)
	1677	de_dxi(num_gp) = 0.0
	1678	de_dyi(num_gp) = de_dy(k,j,i+1)
	1679	de_dzi(num_gp) = de_dz(k,j,i+1)
	1680	ENDIF
	1681
	1682	!
	1683	!-- If wall between gridpoint 5 and gridpoint 6, then
	1684	!-- then Neumann boundary condition has to be applied
	1685	IF ( gp_outside_of_building(5) == 1 .AND. &
	1686	gp_outside_of_building(6) == 0 ) THEN
	1687	num_gp = num_gp + 1
	1688	location(num_gp,1) = (i+1) * dx
	1689	location(num_gp,2) = j * dy + 0.5 * dy
	1690	location(num_gp,3) = k * dz - 0.5 * dz
	1691	ei(num_gp) = e(k,j,i+1)
	1692	dissi(num_gp) = diss(k,j,i+1)
	1693	de_dxi(num_gp) = de_dx(k,j,i+1)
	1694	de_dyi(num_gp) = 0.0
	1695	de_dzi(num_gp) = de_dz(k,j,i+1)
	1696	ENDIF
	1697
	1698	IF ( gp_outside_of_building(6) == 1 .AND. &
	1699	gp_outside_of_building(5) == 0 ) THEN
	1700	num_gp = num_gp + 1
	1701	location(num_gp,1) = (i+1) * dx
	1702	location(num_gp,2) = j * dy + 0.5 * dy
	1703	location(num_gp,3) = k * dz - 0.5 * dz
	1704	ei(num_gp) = e(k,j+1,i+1)
	1705	dissi(num_gp) = diss(k,j+1,i+1)
	1706	de_dxi(num_gp) = de_dx(k,j+1,i+1)
	1707	de_dyi(num_gp) = 0.0
	1708	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	1709	ENDIF
	1710
	1711	!
	1712	!-- If wall between gridpoint 2 and gridpoint 6, then
	1713	!-- Neumann boundary condition has to be applied
	1714	IF ( gp_outside_of_building(2) == 1 .AND. &
	1715	gp_outside_of_building(6) == 0 ) THEN
	1716	num_gp = num_gp + 1
	1717	location(num_gp,1) = i * dx + 0.5 * dx
	1718	location(num_gp,2) = (j+1) * dy
	1719	location(num_gp,3) = k * dz - 0.5 * dz
	1720	ei(num_gp) = e(k,j+1,i)
	1721	dissi(num_gp) = diss(k,j+1,i)
	1722	de_dxi(num_gp) = 0.0
	1723	de_dyi(num_gp) = de_dy(k,j+1,i)
[420]	1724	de_dzi(num_gp) = de_dz(k,j+1,i)
	1725	ENDIF
	1726
[57]	1727	IF ( gp_outside_of_building(6) == 1 .AND. &
	1728	gp_outside_of_building(2) == 0 ) THEN
	1729	num_gp = num_gp + 1
	1730	location(num_gp,1) = i * dx + 0.5 * dx
	1731	location(num_gp,2) = (j+1) * dy
	1732	location(num_gp,3) = k * dz - 0.5 * dz
	1733	ei(num_gp) = e(k,j+1,i+1)
	1734	dissi(num_gp) = diss(k,j+1,i+1)
	1735	de_dxi(num_gp) = 0.0
	1736	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	1737	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	1738	ENDIF
	1739
	1740	!
	1741	!-- If wall between gridpoint 1 and gridpoint 2, then
	1742	!-- Neumann boundary condition has to be applied
	1743	IF ( gp_outside_of_building(1) == 1 .AND. &
	1744	gp_outside_of_building(2) == 0 ) THEN
	1745	num_gp = num_gp + 1
	1746	location(num_gp,1) = i * dx
	1747	location(num_gp,2) = j * dy + 0.5 * dy
	1748	location(num_gp,3) = k * dz - 0.5 * dz
	1749	ei(num_gp) = e(k,j,i)
	1750	dissi(num_gp) = diss(k,j,i)
	1751	de_dxi(num_gp) = de_dx(k,j,i)
	1752	de_dyi(num_gp) = 0.0
	1753	de_dzi(num_gp) = de_dz(k,j,i)
	1754	ENDIF
	1755
	1756	IF ( gp_outside_of_building(2) == 1 .AND. &
	1757	gp_outside_of_building(1) == 0 ) THEN
	1758	num_gp = num_gp + 1
	1759	location(num_gp,1) = i * dx
	1760	location(num_gp,2) = j * dy + 0.5 * dy
	1761	location(num_gp,3) = k * dz - 0.5 * dz
	1762	ei(num_gp) = e(k,j+1,i)
	1763	dissi(num_gp) = diss(k,j+1,i)
	1764	de_dxi(num_gp) = de_dx(k,j+1,i)
	1765	de_dyi(num_gp) = 0.0
	1766	de_dzi(num_gp) = de_dz(k,j+1,i)
	1767	ENDIF
	1768
	1769	!
	1770	!-- If wall between gridpoint 3 and gridpoint 7, then
	1771	!-- Neumann boundary condition has to be applied
	1772	IF ( gp_outside_of_building(3) == 1 .AND. &
	1773	gp_outside_of_building(7) == 0 ) THEN
	1774	num_gp = num_gp + 1
	1775	location(num_gp,1) = i * dx + 0.5 * dx
	1776	location(num_gp,2) = j * dy
	1777	location(num_gp,3) = k * dz + 0.5 * dz
	1778	ei(num_gp) = e(k+1,j,i)
	1779	dissi(num_gp) = diss(k+1,j,i)
	1780	de_dxi(num_gp) = 0.0
	1781	de_dyi(num_gp) = de_dy(k+1,j,i)
[420]	1782	de_dzi(num_gp) = de_dz(k+1,j,i)
	1783	ENDIF
	1784
[57]	1785	IF ( gp_outside_of_building(7) == 1 .AND. &
	1786	gp_outside_of_building(3) == 0 ) THEN
	1787	num_gp = num_gp + 1
	1788	location(num_gp,1) = i * dx + 0.5 * dx
	1789	location(num_gp,2) = j * dy
	1790	location(num_gp,3) = k * dz + 0.5 * dz
	1791	ei(num_gp) = e(k+1,j,i+1)
	1792	dissi(num_gp) = diss(k+1,j,i+1)
	1793	de_dxi(num_gp) = 0.0
	1794	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	1795	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	1796	ENDIF
	1797
	1798	!
	1799	!-- If wall between gridpoint 7 and gridpoint 8, then
	1800	!-- Neumann boundary condition has to be applied
	1801	IF ( gp_outside_of_building(7) == 1 .AND. &
	1802	gp_outside_of_building(8) == 0 ) THEN
	1803	num_gp = num_gp + 1
	1804	location(num_gp,1) = (i+1) * dx
	1805	location(num_gp,2) = j * dy + 0.5 * dy
	1806	location(num_gp,3) = k * dz + 0.5 * dz
	1807	ei(num_gp) = e(k+1,j,i+1)
	1808	dissi(num_gp) = diss(k+1,j,i+1)
	1809	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	1810	de_dyi(num_gp) = 0.0
	1811	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	1812	ENDIF
	1813
	1814	IF ( gp_outside_of_building(8) == 1 .AND. &
	1815	gp_outside_of_building(7) == 0 ) THEN
	1816	num_gp = num_gp + 1
	1817	location(num_gp,1) = (i+1) * dx
	1818	location(num_gp,2) = j * dy + 0.5 * dy
	1819	location(num_gp,3) = k * dz + 0.5 * dz
	1820	ei(num_gp) = e(k+1,j+1,i+1)
	1821	dissi(num_gp) = diss(k+1,j+1,i+1)
	1822	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	1823	de_dyi(num_gp) = 0.0
	1824	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	1825	ENDIF
	1826
	1827	!
	1828	!-- If wall between gridpoint 4 and gridpoint 8, then
	1829	!-- Neumann boundary condition has to be applied
	1830	IF ( gp_outside_of_building(4) == 1 .AND. &
	1831	gp_outside_of_building(8) == 0 ) THEN
	1832	num_gp = num_gp + 1
	1833	location(num_gp,1) = i * dx + 0.5 * dx
	1834	location(num_gp,2) = (j+1) * dy
	1835	location(num_gp,3) = k * dz + 0.5 * dz
	1836	ei(num_gp) = e(k+1,j+1,i)
	1837	dissi(num_gp) = diss(k+1,j+1,i)
	1838	de_dxi(num_gp) = 0.0
	1839	de_dyi(num_gp) = de_dy(k+1,j+1,i)
[420]	1840	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	1841	ENDIF
	1842
[57]	1843	IF ( gp_outside_of_building(8) == 1 .AND. &
	1844	gp_outside_of_building(4) == 0 ) THEN
	1845	num_gp = num_gp + 1
	1846	location(num_gp,1) = i * dx + 0.5 * dx
	1847	location(num_gp,2) = (j+1) * dy
	1848	location(num_gp,3) = k * dz + 0.5 * dz
	1849	ei(num_gp) = e(k+1,j+1,i+1)
	1850	dissi(num_gp) = diss(k+1,j+1,i+1)
	1851	de_dxi(num_gp) = 0.0
	1852	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	1853	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	1854	ENDIF
	1855
	1856	!
	1857	!-- If wall between gridpoint 3 and gridpoint 4, then
	1858	!-- Neumann boundary condition has to be applied
	1859	IF ( gp_outside_of_building(3) == 1 .AND. &
	1860	gp_outside_of_building(4) == 0 ) THEN
	1861	num_gp = num_gp + 1
	1862	location(num_gp,1) = i * dx
	1863	location(num_gp,2) = j * dy + 0.5 * dy
	1864	location(num_gp,3) = k * dz + 0.5 * dz
	1865	ei(num_gp) = e(k+1,j,i)
	1866	dissi(num_gp) = diss(k+1,j,i)
	1867	de_dxi(num_gp) = de_dx(k+1,j,i)
	1868	de_dyi(num_gp) = 0.0
	1869	de_dzi(num_gp) = de_dz(k+1,j,i)
	1870	ENDIF
	1871
	1872	IF ( gp_outside_of_building(4) == 1 .AND. &
	1873	gp_outside_of_building(3) == 0 ) THEN
	1874	num_gp = num_gp + 1
	1875	location(num_gp,1) = i * dx
	1876	location(num_gp,2) = j * dy + 0.5 * dy
	1877	location(num_gp,3) = k * dz + 0.5 * dz
	1878	ei(num_gp) = e(k+1,j+1,i)
	1879	dissi(num_gp) = diss(k+1,j+1,i)
	1880	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	1881	de_dyi(num_gp) = 0.0
	1882	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	1883	ENDIF
	1884
	1885	!
	1886	!-- If wall between gridpoint 1 and gridpoint 3, then
	1887	!-- Neumann boundary condition has to be applied
	1888	!-- (only one case as only building beneath is possible)
[119]	1889	IF ( gp_outside_of_building(1) == 0 .AND. &
	1890	gp_outside_of_building(3) == 1 ) THEN
[57]	1891	num_gp = num_gp + 1
	1892	location(num_gp,1) = i * dx
	1893	location(num_gp,2) = j * dy
	1894	location(num_gp,3) = k * dz
	1895	ei(num_gp) = e(k+1,j,i)
	1896	dissi(num_gp) = diss(k+1,j,i)
	1897	de_dxi(num_gp) = de_dx(k+1,j,i)
	1898	de_dyi(num_gp) = de_dy(k+1,j,i)
	1899	de_dzi(num_gp) = 0.0
	1900	ENDIF
[799]	1901
[57]	1902	!
	1903	!-- If wall between gridpoint 5 and gridpoint 7, then
	1904	!-- Neumann boundary condition has to be applied
	1905	!-- (only one case as only building beneath is possible)
[119]	1906	IF ( gp_outside_of_building(5) == 0 .AND. &
	1907	gp_outside_of_building(7) == 1 ) THEN
[57]	1908	num_gp = num_gp + 1
	1909	location(num_gp,1) = (i+1) * dx
	1910	location(num_gp,2) = j * dy
	1911	location(num_gp,3) = k * dz
	1912	ei(num_gp) = e(k+1,j,i+1)
	1913	dissi(num_gp) = diss(k+1,j,i+1)
	1914	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	1915	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	1916	de_dzi(num_gp) = 0.0
	1917	ENDIF
	1918
	1919	!
	1920	!-- If wall between gridpoint 2 and gridpoint 4, then
	1921	!-- Neumann boundary condition has to be applied
	1922	!-- (only one case as only building beneath is possible)
[119]	1923	IF ( gp_outside_of_building(2) == 0 .AND. &
	1924	gp_outside_of_building(4) == 1 ) THEN
[57]	1925	num_gp = num_gp + 1
	1926	location(num_gp,1) = i * dx
	1927	location(num_gp,2) = (j+1) * dy
	1928	location(num_gp,3) = k * dz
	1929	ei(num_gp) = e(k+1,j+1,i)
	1930	dissi(num_gp) = diss(k+1,j+1,i)
	1931	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	1932	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	1933	de_dzi(num_gp) = 0.0
	1934	ENDIF
	1935
	1936	!
	1937	!-- If wall between gridpoint 6 and gridpoint 8, then
	1938	!-- Neumann boundary condition has to be applied
	1939	!-- (only one case as only building beneath is possible)
[119]	1940	IF ( gp_outside_of_building(6) == 0 .AND. &
	1941	gp_outside_of_building(8) == 1 ) THEN
[57]	1942	num_gp = num_gp + 1
	1943	location(num_gp,1) = (i+1) * dx
	1944	location(num_gp,2) = (j+1) * dy
	1945	location(num_gp,3) = k * dz
	1946	ei(num_gp) = e(k+1,j+1,i+1)
	1947	dissi(num_gp) = diss(k+1,j+1,i+1)
	1948	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	1949	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	1950	de_dzi(num_gp) = 0.0
	1951	ENDIF
	1952
	1953	!
	1954	!-- Carry out the interpolation
	1955	IF ( num_gp == 1 ) THEN
	1956	!
[420]	1957	!-- If only one of the gridpoints is situated outside of the
	1958	!-- building, it follows that the values at the particle
[57]	1959	!-- location are the same as the gridpoint values
[119]	1960	e_int = ei(num_gp)
	1961	diss_int = dissi(num_gp)
	1962	de_dx_int = de_dxi(num_gp)
	1963	de_dy_int = de_dyi(num_gp)
	1964	de_dz_int = de_dzi(num_gp)
[57]	1965	ELSE IF ( num_gp > 1 ) THEN
	1966
	1967	d_sum = 0.0
	1968	!
[420]	1969	!-- Evaluation of the distances between the gridpoints
	1970	!-- contributing to the interpolated values, and the particle
[57]	1971	!-- location
	1972	DO agp = 1, num_gp
	1973	d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 &
	1974	+ ( particles(n)%y-location(agp,2) )**2 &
	1975	+ ( particles(n)%z-location(agp,3) )**2
	1976	d_sum = d_sum + d_gp_pl(agp)
	1977	ENDDO
	1978
	1979	!
	1980	!-- Finally the interpolation can be carried out
	1981	e_int = 0.0
	1982	diss_int = 0.0
	1983	de_dx_int = 0.0
	1984	de_dy_int = 0.0
	1985	de_dz_int = 0.0
	1986	DO agp = 1, num_gp
	1987	e_int = e_int + ( d_sum - d_gp_pl(agp) ) * &
	1988	ei(agp) / ( (num_gp-1) * d_sum )
	1989	diss_int = diss_int + ( d_sum - d_gp_pl(agp) ) * &
	1990	dissi(agp) / ( (num_gp-1) * d_sum )
	1991	de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
	1992	de_dxi(agp) / ( (num_gp-1) * d_sum )
	1993	de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
	1994	de_dyi(agp) / ( (num_gp-1) * d_sum )
	1995	de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
	1996	de_dzi(agp) / ( (num_gp-1) * d_sum )
	1997	ENDDO
	1998
	1999	ENDIF
	2000
	2001	ENDIF
	2002
	2003	ENDIF
	2004
	2005	!
[1]	2006	!-- Vertically interpolate the horizontally averaged SGS TKE and
	2007	!-- resolved-scale velocity variances and use the interpolated values
	2008	!-- to calculate the coefficient fs, which is a measure of the ratio
	2009	!-- of the subgrid-scale turbulent kinetic energy to the total amount
	2010	!-- of turbulent kinetic energy.
	2011	IF ( k == 0 ) THEN
	2012	e_mean_int = hom(0,1,8,0)
	2013	ELSE
	2014	e_mean_int = hom(k,1,8,0) + &
	2015	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
	2016	( zu(k+1) - zu(k) ) * &
	2017	( particles(n)%z - zu(k) )
	2018	ENDIF
	2019
	2020	kw = particles(n)%z / dz
	2021
	2022	IF ( k == 0 ) THEN
	2023	aa = hom(k+1,1,30,0) * ( particles(n)%z / &
	2024	( 0.5 * ( zu(k+1) - zu(k) ) ) )
	2025	bb = hom(k+1,1,31,0) * ( particles(n)%z / &
	2026	( 0.5 * ( zu(k+1) - zu(k) ) ) )
	2027	cc = hom(kw+1,1,32,0) * ( particles(n)%z / &
	2028	( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
	2029	ELSE
	2030	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
	2031	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
	2032	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
	2033	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
	2034	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
	2035	( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
	2036	ENDIF
	2037
	2038	vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
	2039
	2040	fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
	2041	( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
	2042
	2043	!
	2044	!-- Calculate the Lagrangian timescale according to the suggestion of
	2045	!-- Weil et al. (2004).
	2046	lagr_timescale = ( 4.0 * e_int ) / &
	2047	( 3.0 * fs_int * c_0 * diss_int )
	2048
	2049	!
	2050	!-- Calculate the next particle timestep. dt_gap is the time needed to
	2051	!-- complete the current LES timestep.
	2052	dt_gap = dt_3d - particles(n)%dt_sum
	2053	dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
	2054
	2055	!
	2056	!-- The particle timestep should not be too small in order to prevent
	2057	!-- the number of particle timesteps of getting too large
	2058	IF ( dt_particle < dt_min_part .AND. dt_min_part < dt_gap ) &
	2059	THEN
	2060	dt_particle = dt_min_part
	2061	ENDIF
	2062
	2063	!
	2064	!-- Calculate the SGS velocity components
	2065	IF ( particles(n)%age == 0.0 ) THEN
	2066	!
	2067	!-- For new particles the SGS components are derived from the SGS
	2068	!-- TKE. Limit the Gaussian random number to the interval
	2069	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
	2070	!-- from becoming unrealistically large.
[824]	2071	particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * &
	2072	( random_gauss( iran_part, 5.0 ) &
[1]	2073	- 1.0 )
[824]	2074	particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * &
	2075	( random_gauss( iran_part, 5.0 ) &
[1]	2076	- 1.0 )
[824]	2077	particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * &
	2078	( random_gauss( iran_part, 5.0 ) &
[1]	2079	- 1.0 )
	2080
	2081	ELSE
[57]	2082
[1]	2083	!
[57]	2084	!-- Restriction of the size of the new timestep: compared to the
	2085	!-- previous timestep the increase must not exceed 200%
	2086
[60]	2087	dt_particle_m = particles(n)%age - particles(n)%age_m
	2088	IF ( dt_particle > 2.0 * dt_particle_m ) THEN
	2089	dt_particle = 2.0 * dt_particle_m
[57]	2090	ENDIF
	2091
	2092	!
[1]	2093	!-- For old particles the SGS components are correlated with the
	2094	!-- values from the previous timestep. Random numbers have also to
	2095	!-- be limited (see above).
[57]	2096	!-- As negative values for the subgrid TKE are not allowed, the
	2097	!-- change of the subgrid TKE with time cannot be smaller than
	2098	!-- -e_int/dt_particle. This value is used as a lower boundary
	2099	!-- value for the change of TKE
	2100
	2101	de_dt_min = - e_int / dt_particle
	2102
[60]	2103	de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m
[57]	2104
	2105	IF ( de_dt < de_dt_min ) THEN
	2106	de_dt = de_dt_min
	2107	ENDIF
	2108
[824]	2109	particles(n)%rvar1 = particles(n)%rvar1 - &
	2110	fs_int * c_0 * diss_int * particles(n)%rvar1 * &
[1]	2111	dt_particle / ( 4.0 * sgs_wfu_part * e_int ) + &
[57]	2112	( 2.0 * sgs_wfu_part * de_dt * &
[824]	2113	particles(n)%rvar1 / &
[1]	2114	( 2.0 * sgs_wfu_part * e_int ) + de_dx_int &
	2115	) * dt_particle / 2.0 + &
	2116	SQRT( fs_int * c_0 * diss_int ) * &
	2117	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	2118	SQRT( dt_particle )
	2119
[824]	2120	particles(n)%rvar2 = particles(n)%rvar2 - &
	2121	fs_int * c_0 * diss_int * particles(n)%rvar2 * &
[1]	2122	dt_particle / ( 4.0 * sgs_wfv_part * e_int ) + &
[57]	2123	( 2.0 * sgs_wfv_part * de_dt * &
[824]	2124	particles(n)%rvar2 / &
[1]	2125	( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
	2126	) * dt_particle / 2.0 + &
	2127	SQRT( fs_int * c_0 * diss_int ) * &
	2128	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	2129	SQRT( dt_particle )
	2130
[824]	2131	particles(n)%rvar3 = particles(n)%rvar3 - &
	2132	fs_int * c_0 * diss_int * particles(n)%rvar3 * &
[1]	2133	dt_particle / ( 4.0 * sgs_wfw_part * e_int ) + &
[57]	2134	( 2.0 * sgs_wfw_part * de_dt * &
[824]	2135	particles(n)%rvar3 / &
[1]	2136	( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
	2137	) * dt_particle / 2.0 + &
	2138	SQRT( fs_int * c_0 * diss_int ) * &
	2139	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	2140	SQRT( dt_particle )
	2141
	2142	ENDIF
	2143
[824]	2144	u_int = u_int + particles(n)%rvar1
	2145	v_int = v_int + particles(n)%rvar2
	2146	w_int = w_int + particles(n)%rvar3
[1]	2147
	2148	!
	2149	!-- Store the SGS TKE of the current timelevel which is needed for
	2150	!-- for calculating the SGS particle velocities at the next timestep
	2151	particles(n)%e_m = e_int
	2152
	2153	ELSE
	2154	!
	2155	!-- If no SGS velocities are used, only the particle timestep has to
	2156	!-- be set
	2157	dt_particle = dt_3d
	2158
	2159	ENDIF
	2160
	2161	!
[57]	2162	!-- Remember the old age of the particle ( needed to prevent that a
	2163	!-- particle crosses several PEs during one timestep and for the
	2164	!-- evaluation of the subgrid particle velocity fluctuations )
	2165	particles(n)%age_m = particles(n)%age
	2166
	2167
	2168	!
[1]	2169	!-- Particle advection
	2170	IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 ) THEN
	2171	!
	2172	!-- Pure passive transport (without particle inertia)
	2173	particles(n)%x = particles(n)%x + u_int * dt_particle
	2174	particles(n)%y = particles(n)%y + v_int * dt_particle
	2175	particles(n)%z = particles(n)%z + w_int * dt_particle
	2176
	2177	particles(n)%speed_x = u_int
	2178	particles(n)%speed_y = v_int
	2179	particles(n)%speed_z = w_int
	2180
	2181	ELSE
	2182	!
	2183	!-- Transport of particles with inertia
	2184	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
	2185	dt_particle
	2186	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
	2187	dt_particle
	2188	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
	2189	dt_particle
	2190
	2191	!
	2192	!-- Update of the particle velocity
	2193	dens_ratio = particle_groups(particles(n)%group)%density_ratio
	2194	IF ( cloud_droplets ) THEN
	2195	exp_arg = 4.5 * dens_ratio * molecular_viscosity / &
[420]	2196	( particles(n)%radius )*2 &
[1]	2197	( 1.0 + 0.15 * ( 2.0 * particles(n)%radius * &
	2198	SQRT( ( u_int - particles(n)%speed_x )**2 + &
	2199	( v_int - particles(n)%speed_y )**2 + &
	2200	( w_int - particles(n)%speed_z )**2 ) / &
	2201	molecular_viscosity )**0.687 &
	2202	)
	2203	exp_term = EXP( -exp_arg * dt_particle )
	2204	ELSEIF ( use_sgs_for_particles ) THEN
	2205	exp_arg = particle_groups(particles(n)%group)%exp_arg
	2206	exp_term = EXP( -exp_arg * dt_particle )
	2207	ELSE
	2208	exp_arg = particle_groups(particles(n)%group)%exp_arg
	2209	exp_term = particle_groups(particles(n)%group)%exp_term
	2210	ENDIF
	2211	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
	2212	u_int * ( 1.0 - exp_term )
[825]	2213	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
[1]	2214	v_int * ( 1.0 - exp_term )
[825]	2215	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
	2216	( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
	2217	* ( 1.0 - exp_term )
[1]	2218	ENDIF
	2219
	2220	!
	2221	!-- Increment the particle age and the total time that the particle
	2222	!-- has advanced within the particle timestep procedure
	2223	particles(n)%age = particles(n)%age + dt_particle
	2224	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
	2225
	2226	!
	2227	!-- Check whether there is still a particle that has not yet completed
	2228	!-- the total LES timestep
	2229	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 ) THEN
	2230	dt_3d_reached_l = .FALSE.
	2231	ENDIF
	2232
	2233	ENDDO ! advection loop
	2234
	2235	!
[57]	2236	!-- Particle reflection from walls
[59]	2237	CALL particle_boundary_conds
[57]	2238
	2239	!
[212]	2240	!-- User-defined actions after the calculation of the new particle position
[57]	2241	CALL user_advec_particles
	2242
	2243	!
[1]	2244	!-- Find out, if all particles on every PE have completed the LES timestep
	2245	!-- and set the switch corespondingly
	2246	#if defined( __parallel )
[622]	2247	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	2248	CALL MPI_ALLREDUCE( dt_3d_reached_l, dt_3d_reached, 1, MPI_LOGICAL, &
	2249	MPI_LAND, comm2d, ierr )
	2250	#else
	2251	dt_3d_reached = dt_3d_reached_l
[799]	2252	#endif
[1]	2253
	2254	CALL cpu_log( log_point_s(44), 'advec_part_advec', 'stop' )
	2255
	2256	!
	2257	!-- Increment time since last release
	2258	IF ( dt_3d_reached ) time_prel = time_prel + dt_3d
	2259
	2260	! WRITE ( 9, * ) '*** advec_particles: ##0.4'
[82]	2261	! CALL local_flush( 9 )
[1]	2262	! nd = 0
	2263	! DO n = 1, number_of_particles
	2264	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2265	! ENDDO
	2266	! IF ( nd /= deleted_particles ) THEN
	2267	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	2268	! CALL local_flush( 9 )
[1]	2269	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2270	! ENDIF
	2271
	2272	!
	2273	!-- If necessary, release new set of particles
	2274	IF ( time_prel >= dt_prel .AND. end_time_prel > simulated_time .AND. &
	2275	dt_3d_reached ) THEN
	2276
	2277	!
	2278	!-- Check, if particle storage must be extended
	2279	IF ( number_of_particles + number_of_initial_particles > &
	2280	maximum_number_of_particles ) THEN
[519]	2281	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	2282	message_string = 'maximum_number_of_particles ' // &
	2283	'needs to be increased ' // &
[230]	2284	'&but this is not allowed with ' // &
[519]	2285	'netcdf_data_format < 3'
[276]	2286	CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
[1]	2287	ELSE
	2288	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory dt_prel'
[82]	2289	! CALL local_flush( 9 )
[1]	2290	CALL allocate_prt_memory( number_of_initial_particles )
	2291	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory dt_prel'
[82]	2292	! CALL local_flush( 9 )
[1]	2293	ENDIF
	2294	ENDIF
	2295
	2296	!
	2297	!-- Check, if tail storage must be extended
	2298	IF ( use_particle_tails ) THEN
	2299	IF ( number_of_tails + number_of_initial_tails > &
	2300	maximum_number_of_tails ) THEN
[519]	2301	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	2302	message_string = 'maximum_number_of_tails ' // &
	2303	'needs to be increased ' // &
[230]	2304	'&but this is not allowed wi' // &
[519]	2305	'th netcdf_data_format < 3'
[276]	2306	CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
[1]	2307	ELSE
	2308	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory dt_prel'
[82]	2309	! CALL local_flush( 9 )
[1]	2310	CALL allocate_tail_memory( number_of_initial_tails )
	2311	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory dt_prel'
[82]	2312	! CALL local_flush( 9 )
[1]	2313	ENDIF
	2314	ENDIF
	2315	ENDIF
	2316
	2317	!
	2318	!-- The MOD function allows for changes in the output interval with
	2319	!-- restart runs.
	2320	time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
	2321	IF ( number_of_initial_particles /= 0 ) THEN
	2322	is = number_of_particles+1
	2323	ie = number_of_particles+number_of_initial_particles
	2324	particles(is:ie) = initial_particles(1:number_of_initial_particles)
	2325	!
	2326	!-- Add random fluctuation to particle positions. Particles should
	2327	!-- remain in the subdomain.
	2328	IF ( random_start_position ) THEN
	2329	DO n = is, ie
	2330	IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) &
	2331	THEN
	2332	particles(n)%x = particles(n)%x + &
[106]	2333	( random_function( iran_part ) - 0.5 ) *&
[1]	2334	pdx(particles(n)%group)
	2335	IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN
	2336	particles(n)%x = ( nxl - 0.4999999999 ) * dx
	2337	ELSEIF ( particles(n)%x >= ( nxr + 0.5 ) * dx ) THEN
	2338	particles(n)%x = ( nxr + 0.4999999999 ) * dx
	2339	ENDIF
	2340	ENDIF
	2341	IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) &
	2342	THEN
	2343	particles(n)%y = particles(n)%y + &
[106]	2344	( random_function( iran_part ) - 0.5 ) *&
[1]	2345	pdy(particles(n)%group)
	2346	IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN
	2347	particles(n)%y = ( nys - 0.4999999999 ) * dy
	2348	ELSEIF ( particles(n)%y >= ( nyn + 0.5 ) * dy ) THEN
	2349	particles(n)%y = ( nyn + 0.4999999999 ) * dy
	2350	ENDIF
	2351	ENDIF
	2352	IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) &
	2353	THEN
	2354	particles(n)%z = particles(n)%z + &
[106]	2355	( random_function( iran_part ) - 0.5 ) *&
[1]	2356	pdz(particles(n)%group)
	2357	ENDIF
	2358	ENDDO
	2359	ENDIF
	2360
	2361	!
	2362	!-- Set the beginning of the new particle tails and their age
	2363	IF ( use_particle_tails ) THEN
	2364	DO n = is, ie
	2365	!
	2366	!-- New particles which should have a tail, already have got a
	2367	!-- provisional tail id unequal zero (see init_particles)
	2368	IF ( particles(n)%tail_id /= 0 ) THEN
	2369	number_of_tails = number_of_tails + 1
	2370	nn = number_of_tails
	2371	particles(n)%tail_id = nn ! set the final tail id
	2372	particle_tail_coordinates(1,1,nn) = particles(n)%x
	2373	particle_tail_coordinates(1,2,nn) = particles(n)%y
	2374	particle_tail_coordinates(1,3,nn) = particles(n)%z
	2375	particle_tail_coordinates(1,4,nn) = particles(n)%color
	2376	particles(n)%tailpoints = 1
	2377	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
	2378	particle_tail_coordinates(2,1,nn) = particles(n)%x
	2379	particle_tail_coordinates(2,2,nn) = particles(n)%y
	2380	particle_tail_coordinates(2,3,nn) = particles(n)%z
	2381	particle_tail_coordinates(2,4,nn) = particles(n)%color
	2382	particle_tail_coordinates(1:2,5,nn) = 0.0
	2383	particles(n)%tailpoints = 2
	2384	ENDIF
	2385	ENDIF
	2386	ENDDO
	2387	ENDIF
	2388	! WRITE ( 9, * ) '*** advec_particles: after setting the beginning of new tails'
[82]	2389	! CALL local_flush( 9 )
[1]	2390
	2391	number_of_particles = number_of_particles + &
	2392	number_of_initial_particles
	2393	ENDIF
	2394
	2395	ENDIF
	2396
	2397	! WRITE ( 9, * ) '*** advec_particles: ##0.5'
[82]	2398	! CALL local_flush( 9 )
[1]	2399	! nd = 0
	2400	! DO n = 1, number_of_particles
	2401	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2402	! ENDDO
	2403	! IF ( nd /= deleted_particles ) THEN
	2404	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	2405	! CALL local_flush( 9 )
[1]	2406	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2407	! ENDIF
	2408
	2409	! IF ( number_of_particles /= number_of_tails ) THEN
	2410	! WRITE (9,*) '--- advec_particles: #2'
	2411	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
[82]	2412	! CALL local_flush( 9 )
[1]	2413	! ENDIF
	2414	! DO n = 1, number_of_particles
	2415	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2416	! THEN
	2417	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2418	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	2419	! CALL local_flush( 9 )
[1]	2420	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2421	! ENDIF
	2422	! ENDDO
	2423
	2424	#if defined( __parallel )
	2425	!
	2426	!-- As soon as one particle has moved beyond the boundary of the domain, it
	2427	!-- is included in the relevant transfer arrays and marked for subsequent
	2428	!-- deletion on this PE.
	2429	!-- First run for crossings in x direction. Find out first the number of
	2430	!-- particles to be transferred and allocate temporary arrays needed to store
	2431	!-- them.
	2432	!-- For a one-dimensional decomposition along y, no transfer is necessary,
	2433	!-- because the particle remains on the PE.
	2434	trlp_count = 0
	2435	trlpt_count = 0
	2436	trrp_count = 0
	2437	trrpt_count = 0
	2438	IF ( pdims(1) /= 1 ) THEN
	2439	!
	2440	!-- First calculate the storage necessary for sending and receiving the
	2441	!-- data
	2442	DO n = 1, number_of_particles
[799]	2443	i = ( particles(n)%x + 0.5 * dx ) * ddx
[1]	2444	!
	2445	!-- Above calculation does not work for indices less than zero
[799]	2446	IF ( particles(n)%x < -0.5 * dx ) i = -1
[1]	2447
	2448	IF ( i < nxl ) THEN
	2449	trlp_count = trlp_count + 1
	2450	IF ( particles(n)%tail_id /= 0 ) trlpt_count = trlpt_count + 1
	2451	ELSEIF ( i > nxr ) THEN
	2452	trrp_count = trrp_count + 1
	2453	IF ( particles(n)%tail_id /= 0 ) trrpt_count = trrpt_count + 1
	2454	ENDIF
	2455	ENDDO
	2456	IF ( trlp_count == 0 ) trlp_count = 1
	2457	IF ( trlpt_count == 0 ) trlpt_count = 1
	2458	IF ( trrp_count == 0 ) trrp_count = 1
	2459	IF ( trrpt_count == 0 ) trrpt_count = 1
	2460
	2461	ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
	2462
	2463	trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	2464	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
[57]	2465	0.0, 0, 0, 0, 0 )
[1]	2466	trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	2467	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
[57]	2468	0.0, 0, 0, 0, 0 )
[1]	2469
	2470	IF ( use_particle_tails ) THEN
	2471	ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
	2472	trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
	2473	tlength = maximum_number_of_tailpoints * 5
	2474	ENDIF
	2475
	2476	trlp_count = 0
	2477	trlpt_count = 0
	2478	trrp_count = 0
	2479	trrpt_count = 0
	2480
	2481	ENDIF
	2482
	2483	! WRITE ( 9, * ) '*** advec_particles: ##1'
[82]	2484	! CALL local_flush( 9 )
[1]	2485	! nd = 0
	2486	! DO n = 1, number_of_particles
	2487	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2488	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2489	! THEN
	2490	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2491	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	2492	! CALL local_flush( 9 )
[1]	2493	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2494	! ENDIF
	2495	! ENDDO
	2496	! IF ( nd /= deleted_particles ) THEN
	2497	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	2498	! CALL local_flush( 9 )
[1]	2499	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2500	! ENDIF
	2501
	2502	DO n = 1, number_of_particles
	2503
	2504	nn = particles(n)%tail_id
	2505
[799]	2506	i = ( particles(n)%x + 0.5 * dx ) * ddx
[1]	2507	!
	2508	!-- Above calculation does not work for indices less than zero
[799]	2509	IF ( particles(n)%x < - 0.5 * dx ) i = -1
[1]	2510
	2511	IF ( i < nxl ) THEN
	2512	IF ( i < 0 ) THEN
	2513	!
	2514	!-- Apply boundary condition along x
	2515	IF ( ibc_par_lr == 0 ) THEN
	2516	!
	2517	!-- Cyclic condition
	2518	IF ( pdims(1) == 1 ) THEN
	2519	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
	2520	particles(n)%origin_x = ( nx + 1 ) * dx + &
	2521	particles(n)%origin_x
	2522	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2523	i = particles(n)%tailpoints
	2524	particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
	2525	+ particle_tail_coordinates(1:i,1,nn)
	2526	ENDIF
	2527	ELSE
	2528	trlp_count = trlp_count + 1
	2529	trlp(trlp_count) = particles(n)
	2530	trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
	2531	trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
	2532	( nx + 1 ) * dx
	2533	particle_mask(n) = .FALSE.
	2534	deleted_particles = deleted_particles + 1
	2535
[420]	2536	IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.e-12 ) THEN
	2537	trlp(trlp_count)%x = trlp(trlp_count)%x - 1.e-10
	2538	trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x &
	2539	- 1
	2540	ENDIF
	2541
[1]	2542	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2543	trlpt_count = trlpt_count + 1
	2544	trlpt(:,:,trlpt_count) = &
	2545	particle_tail_coordinates(:,:,nn)
	2546	trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
	2547	trlpt(:,1,trlpt_count)
	2548	tail_mask(nn) = .FALSE.
	2549	deleted_tails = deleted_tails + 1
	2550	ENDIF
	2551	ENDIF
	2552
	2553	ELSEIF ( ibc_par_lr == 1 ) THEN
	2554	!
	2555	!-- Particle absorption
	2556	particle_mask(n) = .FALSE.
	2557	deleted_particles = deleted_particles + 1
	2558	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2559	tail_mask(nn) = .FALSE.
	2560	deleted_tails = deleted_tails + 1
	2561	ENDIF
	2562
	2563	ELSEIF ( ibc_par_lr == 2 ) THEN
	2564	!
	2565	!-- Particle reflection
	2566	particles(n)%x = -particles(n)%x
	2567	particles(n)%speed_x = -particles(n)%speed_x
	2568
	2569	ENDIF
	2570	ELSE
	2571	!
	2572	!-- Store particle data in the transfer array, which will be send
	2573	!-- to the neighbouring PE
	2574	trlp_count = trlp_count + 1
	2575	trlp(trlp_count) = particles(n)
	2576	particle_mask(n) = .FALSE.
	2577	deleted_particles = deleted_particles + 1
	2578
	2579	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2580	trlpt_count = trlpt_count + 1
	2581	trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
	2582	tail_mask(nn) = .FALSE.
	2583	deleted_tails = deleted_tails + 1
	2584	ENDIF
	2585	ENDIF
	2586
	2587	ELSEIF ( i > nxr ) THEN
	2588	IF ( i > nx ) THEN
	2589	!
	2590	!-- Apply boundary condition along x
	2591	IF ( ibc_par_lr == 0 ) THEN
	2592	!
	2593	!-- Cyclic condition
	2594	IF ( pdims(1) == 1 ) THEN
	2595	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
	2596	particles(n)%origin_x = particles(n)%origin_x - &
	2597	( nx + 1 ) * dx
	2598	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2599	i = particles(n)%tailpoints
	2600	particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
	2601	+ particle_tail_coordinates(1:i,1,nn)
	2602	ENDIF
	2603	ELSE
	2604	trrp_count = trrp_count + 1
	2605	trrp(trrp_count) = particles(n)
	2606	trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
	2607	trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
	2608	( nx + 1 ) * dx
	2609	particle_mask(n) = .FALSE.
	2610	deleted_particles = deleted_particles + 1
	2611
	2612	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2613	trrpt_count = trrpt_count + 1
	2614	trrpt(:,:,trrpt_count) = &
	2615	particle_tail_coordinates(:,:,nn)
	2616	trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
	2617	( nx + 1 ) * dx
	2618	tail_mask(nn) = .FALSE.
	2619	deleted_tails = deleted_tails + 1
	2620	ENDIF
	2621	ENDIF
	2622
	2623	ELSEIF ( ibc_par_lr == 1 ) THEN
	2624	!
	2625	!-- Particle absorption
	2626	particle_mask(n) = .FALSE.
	2627	deleted_particles = deleted_particles + 1
	2628	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2629	tail_mask(nn) = .FALSE.
	2630	deleted_tails = deleted_tails + 1
	2631	ENDIF
	2632
	2633	ELSEIF ( ibc_par_lr == 2 ) THEN
	2634	!
	2635	!-- Particle reflection
	2636	particles(n)%x = 2 * ( nx * dx ) - particles(n)%x
	2637	particles(n)%speed_x = -particles(n)%speed_x
	2638
	2639	ENDIF
	2640	ELSE
	2641	!
	2642	!-- Store particle data in the transfer array, which will be send
	2643	!-- to the neighbouring PE
	2644	trrp_count = trrp_count + 1
	2645	trrp(trrp_count) = particles(n)
	2646	particle_mask(n) = .FALSE.
	2647	deleted_particles = deleted_particles + 1
	2648
	2649	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2650	trrpt_count = trrpt_count + 1
	2651	trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
	2652	tail_mask(nn) = .FALSE.
	2653	deleted_tails = deleted_tails + 1
	2654	ENDIF
	2655	ENDIF
	2656
	2657	ENDIF
	2658	ENDDO
	2659
	2660	! WRITE ( 9, * ) '*** advec_particles: ##2'
[82]	2661	! CALL local_flush( 9 )
[1]	2662	! nd = 0
	2663	! DO n = 1, number_of_particles
	2664	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2665	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2666	! THEN
	2667	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2668	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	2669	! CALL local_flush( 9 )
[1]	2670	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2671	! ENDIF
	2672	! ENDDO
	2673	! IF ( nd /= deleted_particles ) THEN
	2674	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	2675	! CALL local_flush( 9 )
[1]	2676	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2677	! ENDIF
	2678
	2679	!
	2680	!-- Send left boundary, receive right boundary (but first exchange how many
	2681	!-- and check, if particle storage must be extended)
	2682	IF ( pdims(1) /= 1 ) THEN
	2683
	2684	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'start' )
	2685	CALL MPI_SENDRECV( trlp_count, 1, MPI_INTEGER, pleft, 0, &
	2686	trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
	2687	comm2d, status, ierr )
	2688
	2689	IF ( number_of_particles + trrp_count_recv > &
	2690	maximum_number_of_particles ) &
	2691	THEN
[519]	2692	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	2693	message_string = 'maximum_number_of_particles ' // &
	2694	'needs to be increased ' // &
[230]	2695	'&but this is not allowed with ' // &
[519]	2696	'netcdf-data_format < 3'
[276]	2697	CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
[1]	2698	ELSE
	2699	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trrp'
[82]	2700	! CALL local_flush( 9 )
[1]	2701	CALL allocate_prt_memory( trrp_count_recv )
	2702	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trrp'
[82]	2703	! CALL local_flush( 9 )
[1]	2704	ENDIF
	2705	ENDIF
	2706
[64]	2707	CALL MPI_SENDRECV( trlp(1)%age, trlp_count, mpi_particle_type, &
	2708	pleft, 1, particles(number_of_particles+1)%age, &
	2709	trrp_count_recv, mpi_particle_type, pright, 1, &
[1]	2710	comm2d, status, ierr )
	2711
	2712	IF ( use_particle_tails ) THEN
	2713
	2714	CALL MPI_SENDRECV( trlpt_count, 1, MPI_INTEGER, pleft, 0, &
	2715	trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
	2716	comm2d, status, ierr )
	2717
	2718	IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
	2719	THEN
[519]	2720	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	2721	message_string = 'maximum_number_of_tails ' // &
	2722	'needs to be increased ' // &
[230]	2723	'&but this is not allowed wi'// &
[519]	2724	'th netcdf_data_format < 3'
[276]	2725	CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
[1]	2726	ELSE
	2727	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trrpt'
[82]	2728	! CALL local_flush( 9 )
[1]	2729	CALL allocate_tail_memory( trrpt_count_recv )
	2730	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trrpt'
[82]	2731	! CALL local_flush( 9 )
[1]	2732	ENDIF
	2733	ENDIF
	2734
[63]	2735	CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL, &
	2736	pleft, 1, &
[1]	2737	particle_tail_coordinates(1,1,number_of_tails+1), &
	2738	trrpt_count_recv*tlength, MPI_REAL, pright, 1, &
	2739	comm2d, status, ierr )
	2740	!
	2741	!-- Update the tail ids for the transferred particles
	2742	nn = number_of_tails
	2743	DO n = number_of_particles+1, number_of_particles+trrp_count_recv
	2744	IF ( particles(n)%tail_id /= 0 ) THEN
	2745	nn = nn + 1
	2746	particles(n)%tail_id = nn
	2747	ENDIF
	2748	ENDDO
	2749
	2750	ENDIF
	2751
	2752	number_of_particles = number_of_particles + trrp_count_recv
	2753	number_of_tails = number_of_tails + trrpt_count_recv
	2754	! IF ( number_of_particles /= number_of_tails ) THEN
	2755	! WRITE (9,*) '--- advec_particles: #3'
	2756	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
[82]	2757	! CALL local_flush( 9 )
[1]	2758	! ENDIF
	2759
	2760	!
	2761	!-- Send right boundary, receive left boundary
	2762	CALL MPI_SENDRECV( trrp_count, 1, MPI_INTEGER, pright, 0, &
	2763	trlp_count_recv, 1, MPI_INTEGER, pleft, 0, &
	2764	comm2d, status, ierr )
	2765
	2766	IF ( number_of_particles + trlp_count_recv > &
	2767	maximum_number_of_particles ) &
	2768	THEN
[519]	2769	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	2770	message_string = 'maximum_number_of_particles ' // &
	2771	'needs to be increased ' // &
[230]	2772	'&but this is not allowed with '// &
[519]	2773	'netcdf_data_format < 3'
[276]	2774	CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
[1]	2775	ELSE
	2776	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trlp'
[82]	2777	! CALL local_flush( 9 )
[1]	2778	CALL allocate_prt_memory( trlp_count_recv )
	2779	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trlp'
[82]	2780	! CALL local_flush( 9 )
[1]	2781	ENDIF
	2782	ENDIF
	2783
[64]	2784	CALL MPI_SENDRECV( trrp(1)%age, trrp_count, mpi_particle_type, &
	2785	pright, 1, particles(number_of_particles+1)%age, &
	2786	trlp_count_recv, mpi_particle_type, pleft, 1, &
[1]	2787	comm2d, status, ierr )
	2788
	2789	IF ( use_particle_tails ) THEN
	2790
	2791	CALL MPI_SENDRECV( trrpt_count, 1, MPI_INTEGER, pright, 0, &
	2792	trlpt_count_recv, 1, MPI_INTEGER, pleft, 0, &
	2793	comm2d, status, ierr )
	2794
	2795	IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
	2796	THEN
[519]	2797	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	2798	message_string = 'maximum_number_of_tails ' // &
	2799	'needs to be increased ' // &
[230]	2800	'&but this is not allowed wi'// &
[519]	2801	'th netcdf_data_format < 3'
[276]	2802	CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
[1]	2803	ELSE
	2804	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trlpt'
[82]	2805	! CALL local_flush( 9 )
[1]	2806	CALL allocate_tail_memory( trlpt_count_recv )
	2807	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trlpt'
[82]	2808	! CALL local_flush( 9 )
[1]	2809	ENDIF
	2810	ENDIF
	2811
[63]	2812	CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL, &
	2813	pright, 1, &
	2814	particle_tail_coordinates(1,1,number_of_tails+1), &
[1]	2815	trlpt_count_recv*tlength, MPI_REAL, pleft, 1, &
	2816	comm2d, status, ierr )
	2817	!
	2818	!-- Update the tail ids for the transferred particles
	2819	nn = number_of_tails
	2820	DO n = number_of_particles+1, number_of_particles+trlp_count_recv
	2821	IF ( particles(n)%tail_id /= 0 ) THEN
	2822	nn = nn + 1
	2823	particles(n)%tail_id = nn
	2824	ENDIF
	2825	ENDDO
	2826
	2827	ENDIF
	2828
	2829	number_of_particles = number_of_particles + trlp_count_recv
	2830	number_of_tails = number_of_tails + trlpt_count_recv
	2831	! IF ( number_of_particles /= number_of_tails ) THEN
	2832	! WRITE (9,*) '--- advec_particles: #4'
	2833	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
[82]	2834	! CALL local_flush( 9 )
[1]	2835	! ENDIF
	2836
	2837	IF ( use_particle_tails ) THEN
	2838	DEALLOCATE( trlpt, trrpt )
	2839	ENDIF
	2840	DEALLOCATE( trlp, trrp )
	2841
	2842	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'pause' )
	2843
	2844	ENDIF
	2845
	2846	! WRITE ( 9, * ) '*** advec_particles: ##3'
[82]	2847	! CALL local_flush( 9 )
[1]	2848	! nd = 0
	2849	! DO n = 1, number_of_particles
	2850	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2851	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2852	! THEN
	2853	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2854	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	2855	! CALL local_flush( 9 )
[1]	2856	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2857	! ENDIF
	2858	! ENDDO
	2859	! IF ( nd /= deleted_particles ) THEN
	2860	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	2861	! CALL local_flush( 9 )
[1]	2862	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2863	! ENDIF
	2864
	2865	!
	2866	!-- Check whether particles have crossed the boundaries in y direction. Note
	2867	!-- that this case can also apply to particles that have just been received
	2868	!-- from the adjacent right or left PE.
	2869	!-- Find out first the number of particles to be transferred and allocate
	2870	!-- temporary arrays needed to store them.
	2871	!-- For a one-dimensional decomposition along x, no transfer is necessary,
	2872	!-- because the particle remains on the PE.
	2873	trsp_count = 0
	2874	trspt_count = 0
	2875	trnp_count = 0
	2876	trnpt_count = 0
	2877	IF ( pdims(2) /= 1 ) THEN
	2878	!
	2879	!-- First calculate the storage necessary for sending and receiving the
	2880	!-- data
	2881	DO n = 1, number_of_particles
	2882	IF ( particle_mask(n) ) THEN
[799]	2883	j = ( particles(n)%y + 0.5 * dy ) * ddy
[1]	2884	!
	2885	!-- Above calculation does not work for indices less than zero
[799]	2886	IF ( particles(n)%y < -0.5 * dy ) j = -1
[1]	2887
	2888	IF ( j < nys ) THEN
	2889	trsp_count = trsp_count + 1
	2890	IF ( particles(n)%tail_id /= 0 ) trspt_count = trspt_count+1
	2891	ELSEIF ( j > nyn ) THEN
	2892	trnp_count = trnp_count + 1
	2893	IF ( particles(n)%tail_id /= 0 ) trnpt_count = trnpt_count+1
	2894	ENDIF
	2895	ENDIF
	2896	ENDDO
	2897	IF ( trsp_count == 0 ) trsp_count = 1
	2898	IF ( trspt_count == 0 ) trspt_count = 1
	2899	IF ( trnp_count == 0 ) trnp_count = 1
	2900	IF ( trnpt_count == 0 ) trnpt_count = 1
	2901
	2902	ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
	2903
	2904	trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	2905	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
[57]	2906	0.0, 0, 0, 0, 0 )
[1]	2907	trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	2908	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
[57]	2909	0.0, 0, 0, 0, 0 )
[1]	2910
	2911	IF ( use_particle_tails ) THEN
	2912	ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
	2913	trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
	2914	tlength = maximum_number_of_tailpoints * 5
	2915	ENDIF
	2916
	2917	trsp_count = 0
	2918	trspt_count = 0
	2919	trnp_count = 0
	2920	trnpt_count = 0
	2921
	2922	ENDIF
	2923
	2924	! WRITE ( 9, * ) '*** advec_particles: ##4'
[82]	2925	! CALL local_flush( 9 )
[1]	2926	! nd = 0
	2927	! DO n = 1, number_of_particles
	2928	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	2929	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	2930	! THEN
	2931	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	2932	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	2933	! CALL local_flush( 9 )
[1]	2934	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2935	! ENDIF
	2936	! ENDDO
	2937	! IF ( nd /= deleted_particles ) THEN
	2938	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	2939	! CALL local_flush( 9 )
[1]	2940	! CALL MPI_ABORT( comm2d, 9999, ierr )
	2941	! ENDIF
	2942
	2943	DO n = 1, number_of_particles
	2944
	2945	nn = particles(n)%tail_id
	2946	!
	2947	!-- Only those particles that have not been marked as 'deleted' may be
	2948	!-- moved.
	2949	IF ( particle_mask(n) ) THEN
[799]	2950	j = ( particles(n)%y + 0.5 * dy ) * ddy
[1]	2951	!
	2952	!-- Above calculation does not work for indices less than zero
[799]	2953	IF ( particles(n)%y < -0.5 * dy ) j = -1
[1]	2954
	2955	IF ( j < nys ) THEN
	2956	IF ( j < 0 ) THEN
	2957	!
	2958	!-- Apply boundary condition along y
	2959	IF ( ibc_par_ns == 0 ) THEN
	2960	!
	2961	!-- Cyclic condition
	2962	IF ( pdims(2) == 1 ) THEN
	2963	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
	2964	particles(n)%origin_y = ( ny + 1 ) * dy + &
	2965	particles(n)%origin_y
	2966	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2967	i = particles(n)%tailpoints
	2968	particle_tail_coordinates(1:i,2,nn) = ( ny+1 ) * dy&
	2969	+ particle_tail_coordinates(1:i,2,nn)
	2970	ENDIF
	2971	ELSE
	2972	trsp_count = trsp_count + 1
	2973	trsp(trsp_count) = particles(n)
	2974	trsp(trsp_count)%y = ( ny + 1 ) * dy + &
	2975	trsp(trsp_count)%y
	2976	trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
	2977	+ ( ny + 1 ) * dy
	2978	particle_mask(n) = .FALSE.
	2979	deleted_particles = deleted_particles + 1
	2980
[420]	2981	IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.e-12 ) THEN
	2982	trsp(trsp_count)%y = trsp(trsp_count)%y - 1.e-10
	2983	trsp(trsp_count)%origin_y = &
	2984	trsp(trsp_count)%origin_y - 1
	2985	ENDIF
	2986
[1]	2987	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	2988	trspt_count = trspt_count + 1
	2989	trspt(:,:,trspt_count) = &
	2990	particle_tail_coordinates(:,:,nn)
	2991	trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
	2992	trspt(:,2,trspt_count)
	2993	tail_mask(nn) = .FALSE.
	2994	deleted_tails = deleted_tails + 1
	2995	ENDIF
	2996	ENDIF
	2997
	2998	ELSEIF ( ibc_par_ns == 1 ) THEN
	2999	!
	3000	!-- Particle absorption
	3001	particle_mask(n) = .FALSE.
	3002	deleted_particles = deleted_particles + 1
	3003	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3004	tail_mask(nn) = .FALSE.
	3005	deleted_tails = deleted_tails + 1
	3006	ENDIF
	3007
	3008	ELSEIF ( ibc_par_ns == 2 ) THEN
	3009	!
	3010	!-- Particle reflection
	3011	particles(n)%y = -particles(n)%y
	3012	particles(n)%speed_y = -particles(n)%speed_y
	3013
	3014	ENDIF
	3015	ELSE
	3016	!
	3017	!-- Store particle data in the transfer array, which will be send
	3018	!-- to the neighbouring PE
	3019	trsp_count = trsp_count + 1
	3020	trsp(trsp_count) = particles(n)
	3021	particle_mask(n) = .FALSE.
	3022	deleted_particles = deleted_particles + 1
	3023
	3024	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3025	trspt_count = trspt_count + 1
	3026	trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
	3027	tail_mask(nn) = .FALSE.
	3028	deleted_tails = deleted_tails + 1
	3029	ENDIF
	3030	ENDIF
	3031
[16]	3032	ELSEIF ( j > nyn ) THEN
[1]	3033	IF ( j > ny ) THEN
	3034	!
	3035	!-- Apply boundary condition along x
	3036	IF ( ibc_par_ns == 0 ) THEN
	3037	!
	3038	!-- Cyclic condition
	3039	IF ( pdims(2) == 1 ) THEN
	3040	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
	3041	particles(n)%origin_y = particles(n)%origin_y - &
	3042	( ny + 1 ) * dy
	3043	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3044	i = particles(n)%tailpoints
	3045	particle_tail_coordinates(1:i,2,nn) = - (ny+1) * dy&
	3046	+ particle_tail_coordinates(1:i,2,nn)
	3047	ENDIF
	3048	ELSE
	3049	trnp_count = trnp_count + 1
	3050	trnp(trnp_count) = particles(n)
	3051	trnp(trnp_count)%y = trnp(trnp_count)%y - &
	3052	( ny + 1 ) * dy
	3053	trnp(trnp_count)%origin_y = trnp(trnp_count)%origin_y &
	3054	- ( ny + 1 ) * dy
	3055	particle_mask(n) = .FALSE.
	3056	deleted_particles = deleted_particles + 1
	3057
	3058	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3059	trnpt_count = trnpt_count + 1
	3060	trnpt(:,:,trnpt_count) = &
	3061	particle_tail_coordinates(:,:,nn)
	3062	trnpt(:,2,trnpt_count) = trnpt(:,2,trnpt_count) - &
	3063	( ny + 1 ) * dy
	3064	tail_mask(nn) = .FALSE.
	3065	deleted_tails = deleted_tails + 1
	3066	ENDIF
	3067	ENDIF
	3068
	3069	ELSEIF ( ibc_par_ns == 1 ) THEN
	3070	!
	3071	!-- Particle absorption
	3072	particle_mask(n) = .FALSE.
	3073	deleted_particles = deleted_particles + 1
	3074	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3075	tail_mask(nn) = .FALSE.
	3076	deleted_tails = deleted_tails + 1
	3077	ENDIF
	3078
	3079	ELSEIF ( ibc_par_ns == 2 ) THEN
	3080	!
	3081	!-- Particle reflection
	3082	particles(n)%y = 2 * ( ny * dy ) - particles(n)%y
	3083	particles(n)%speed_y = -particles(n)%speed_y
	3084
	3085	ENDIF
	3086	ELSE
	3087	!
	3088	!-- Store particle data in the transfer array, which will be send
	3089	!-- to the neighbouring PE
	3090	trnp_count = trnp_count + 1
	3091	trnp(trnp_count) = particles(n)
	3092	particle_mask(n) = .FALSE.
	3093	deleted_particles = deleted_particles + 1
	3094
	3095	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3096	trnpt_count = trnpt_count + 1
	3097	trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
	3098	tail_mask(nn) = .FALSE.
	3099	deleted_tails = deleted_tails + 1
	3100	ENDIF
	3101	ENDIF
	3102
	3103	ENDIF
	3104	ENDIF
	3105	ENDDO
	3106
	3107	! WRITE ( 9, * ) '*** advec_particles: ##5'
[82]	3108	! CALL local_flush( 9 )
[1]	3109	! nd = 0
	3110	! DO n = 1, number_of_particles
	3111	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	3112	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	3113	! THEN
	3114	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	3115	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	3116	! CALL local_flush( 9 )
[1]	3117	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3118	! ENDIF
	3119	! ENDDO
	3120	! IF ( nd /= deleted_particles ) THEN
	3121	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	3122	! CALL local_flush( 9 )
[1]	3123	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3124	! ENDIF
	3125
	3126	!
	3127	!-- Send front boundary, receive back boundary (but first exchange how many
	3128	!-- and check, if particle storage must be extended)
	3129	IF ( pdims(2) /= 1 ) THEN
	3130
	3131	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'continue' )
	3132	CALL MPI_SENDRECV( trsp_count, 1, MPI_INTEGER, psouth, 0, &
	3133	trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
	3134	comm2d, status, ierr )
	3135
	3136	IF ( number_of_particles + trnp_count_recv > &
	3137	maximum_number_of_particles ) &
	3138	THEN
[519]	3139	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	3140	message_string = 'maximum_number_of_particles ' // &
	3141	'needs to be increased ' // &
[230]	3142	'&but this is not allowed with '// &
[519]	3143	'netcdf_data_format < 3'
[276]	3144	CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
[1]	3145	ELSE
	3146	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trnp'
[82]	3147	! CALL local_flush( 9 )
[1]	3148	CALL allocate_prt_memory( trnp_count_recv )
	3149	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trnp'
[82]	3150	! CALL local_flush( 9 )
[1]	3151	ENDIF
	3152	ENDIF
	3153
[64]	3154	CALL MPI_SENDRECV( trsp(1)%age, trsp_count, mpi_particle_type, &
	3155	psouth, 1, particles(number_of_particles+1)%age, &
	3156	trnp_count_recv, mpi_particle_type, pnorth, 1, &
[1]	3157	comm2d, status, ierr )
	3158
	3159	IF ( use_particle_tails ) THEN
	3160
[114]	3161	CALL MPI_SENDRECV( trspt_count, 1, MPI_INTEGER, psouth, 0, &
	3162	trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
[1]	3163	comm2d, status, ierr )
	3164
	3165	IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
	3166	THEN
[519]	3167	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	3168	message_string = 'maximum_number_of_tails ' // &
	3169	'needs to be increased ' // &
[230]	3170	'&but this is not allowed wi' // &
[519]	3171	'th netcdf_data_format < 3'
[276]	3172	CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
[1]	3173	ELSE
	3174	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trnpt'
[82]	3175	! CALL local_flush( 9 )
[1]	3176	CALL allocate_tail_memory( trnpt_count_recv )
	3177	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trnpt'
[82]	3178	! CALL local_flush( 9 )
[1]	3179	ENDIF
	3180	ENDIF
	3181
[63]	3182	CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL, &
	3183	psouth, 1, &
	3184	particle_tail_coordinates(1,1,number_of_tails+1), &
[1]	3185	trnpt_count_recv*tlength, MPI_REAL, pnorth, 1, &
	3186	comm2d, status, ierr )
[276]	3187
[1]	3188	!
	3189	!-- Update the tail ids for the transferred particles
	3190	nn = number_of_tails
	3191	DO n = number_of_particles+1, number_of_particles+trnp_count_recv
	3192	IF ( particles(n)%tail_id /= 0 ) THEN
	3193	nn = nn + 1
	3194	particles(n)%tail_id = nn
	3195	ENDIF
	3196	ENDDO
	3197
	3198	ENDIF
	3199
	3200	number_of_particles = number_of_particles + trnp_count_recv
	3201	number_of_tails = number_of_tails + trnpt_count_recv
	3202	! IF ( number_of_particles /= number_of_tails ) THEN
	3203	! WRITE (9,*) '--- advec_particles: #5'
	3204	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
[82]	3205	! CALL local_flush( 9 )
[1]	3206	! ENDIF
	3207
	3208	!
	3209	!-- Send back boundary, receive front boundary
	3210	CALL MPI_SENDRECV( trnp_count, 1, MPI_INTEGER, pnorth, 0, &
	3211	trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
	3212	comm2d, status, ierr )
	3213
	3214	IF ( number_of_particles + trsp_count_recv > &
	3215	maximum_number_of_particles ) &
	3216	THEN
[519]	3217	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	3218	message_string = 'maximum_number_of_particles ' // &
	3219	'needs to be increased ' // &
[230]	3220	'&but this is not allowed with ' // &
[519]	3221	'netcdf_data_format < 3'
[276]	3222	CALL message( 'advec_particles', 'PA0146', 2, 2, -1, 6, 1 )
[1]	3223	ELSE
	3224	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trsp'
[82]	3225	! CALL local_flush( 9 )
[1]	3226	CALL allocate_prt_memory( trsp_count_recv )
	3227	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trsp'
[82]	3228	! CALL local_flush( 9 )
[1]	3229	ENDIF
	3230	ENDIF
	3231
[64]	3232	CALL MPI_SENDRECV( trnp(1)%age, trnp_count, mpi_particle_type, &
	3233	pnorth, 1, particles(number_of_particles+1)%age, &
	3234	trsp_count_recv, mpi_particle_type, psouth, 1, &
[1]	3235	comm2d, status, ierr )
	3236
	3237	IF ( use_particle_tails ) THEN
	3238
[114]	3239	CALL MPI_SENDRECV( trnpt_count, 1, MPI_INTEGER, pnorth, 0, &
	3240	trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
[1]	3241	comm2d, status, ierr )
	3242
	3243	IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
	3244	THEN
[519]	3245	IF ( netcdf_output .AND. netcdf_data_format < 3 ) THEN
[274]	3246	message_string = 'maximum_number_of_tails ' // &
	3247	'needs to be increased ' // &
[230]	3248	'&but this is not allowed wi'// &
	3249	'th NetCDF output switched on'
[276]	3250	CALL message( 'advec_particles', 'PA0147', 2, 2, -1, 6, 1 )
[1]	3251	ELSE
	3252	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trspt'
[82]	3253	! CALL local_flush( 9 )
[1]	3254	CALL allocate_tail_memory( trspt_count_recv )
	3255	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trspt'
[82]	3256	! CALL local_flush( 9 )
[1]	3257	ENDIF
	3258	ENDIF
	3259
[63]	3260	CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL, &
	3261	pnorth, 1, &
	3262	particle_tail_coordinates(1,1,number_of_tails+1), &
[1]	3263	trspt_count_recv*tlength, MPI_REAL, psouth, 1, &
	3264	comm2d, status, ierr )
	3265	!
	3266	!-- Update the tail ids for the transferred particles
	3267	nn = number_of_tails
	3268	DO n = number_of_particles+1, number_of_particles+trsp_count_recv
	3269	IF ( particles(n)%tail_id /= 0 ) THEN
	3270	nn = nn + 1
	3271	particles(n)%tail_id = nn
	3272	ENDIF
	3273	ENDDO
	3274
	3275	ENDIF
	3276
	3277	number_of_particles = number_of_particles + trsp_count_recv
	3278	number_of_tails = number_of_tails + trspt_count_recv
	3279	! IF ( number_of_particles /= number_of_tails ) THEN
	3280	! WRITE (9,*) '--- advec_particles: #6'
	3281	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
[82]	3282	! CALL local_flush( 9 )
[1]	3283	! ENDIF
	3284
	3285	IF ( use_particle_tails ) THEN
	3286	DEALLOCATE( trspt, trnpt )
	3287	ENDIF
	3288	DEALLOCATE( trsp, trnp )
	3289
	3290	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'stop' )
	3291
	3292	ENDIF
	3293
	3294	! WRITE ( 9, * ) '*** advec_particles: ##6'
[82]	3295	! CALL local_flush( 9 )
[1]	3296	! nd = 0
	3297	! DO n = 1, number_of_particles
	3298	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	3299	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	3300	! THEN
	3301	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	3302	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	3303	! CALL local_flush( 9 )
[1]	3304	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3305	! ENDIF
	3306	! ENDDO
	3307	! IF ( nd /= deleted_particles ) THEN
	3308	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	3309	! CALL local_flush( 9 )
[1]	3310	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3311	! ENDIF
	3312
	3313	#else
	3314
	3315	!
	3316	!-- Apply boundary conditions
	3317	DO n = 1, number_of_particles
	3318
	3319	nn = particles(n)%tail_id
	3320
	3321	IF ( particles(n)%x < -0.5 * dx ) THEN
	3322
	3323	IF ( ibc_par_lr == 0 ) THEN
	3324	!
	3325	!-- Cyclic boundary. Relevant coordinate has to be changed.
	3326	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
	3327	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3328	i = particles(n)%tailpoints
	3329	particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
	3330	particle_tail_coordinates(1:i,1,nn)
	3331	ENDIF
	3332	ELSEIF ( ibc_par_lr == 1 ) THEN
	3333	!
	3334	!-- Particle absorption
	3335	particle_mask(n) = .FALSE.
	3336	deleted_particles = deleted_particles + 1
	3337	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3338	tail_mask(nn) = .FALSE.
	3339	deleted_tails = deleted_tails + 1
	3340	ENDIF
	3341	ELSEIF ( ibc_par_lr == 2 ) THEN
	3342	!
	3343	!-- Particle reflection
	3344	particles(n)%x = -dx - particles(n)%x
	3345	particles(n)%speed_x = -particles(n)%speed_x
	3346	ENDIF
	3347
	3348	ELSEIF ( particles(n)%x >= ( nx + 0.5 ) * dx ) THEN
	3349
	3350	IF ( ibc_par_lr == 0 ) THEN
	3351	!
	3352	!-- Cyclic boundary. Relevant coordinate has to be changed.
	3353	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
	3354	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3355	i = particles(n)%tailpoints
	3356	particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
	3357	particle_tail_coordinates(1:i,1,nn)
	3358	ENDIF
	3359	ELSEIF ( ibc_par_lr == 1 ) THEN
	3360	!
	3361	!-- Particle absorption
	3362	particle_mask(n) = .FALSE.
	3363	deleted_particles = deleted_particles + 1
	3364	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3365	tail_mask(nn) = .FALSE.
	3366	deleted_tails = deleted_tails + 1
	3367	ENDIF
	3368	ELSEIF ( ibc_par_lr == 2 ) THEN
	3369	!
	3370	!-- Particle reflection
	3371	particles(n)%x = ( nx + 1 ) * dx - particles(n)%x
	3372	particles(n)%speed_x = -particles(n)%speed_x
	3373	ENDIF
	3374
	3375	ENDIF
	3376
	3377	IF ( particles(n)%y < -0.5 * dy ) THEN
	3378
	3379	IF ( ibc_par_ns == 0 ) THEN
	3380	!
	3381	!-- Cyclic boundary. Relevant coordinate has to be changed.
	3382	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
	3383	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3384	i = particles(n)%tailpoints
	3385	particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
	3386	particle_tail_coordinates(1:i,2,nn)
	3387	ENDIF
	3388	ELSEIF ( ibc_par_ns == 1 ) THEN
	3389	!
	3390	!-- Particle absorption
	3391	particle_mask(n) = .FALSE.
	3392	deleted_particles = deleted_particles + 1
	3393	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3394	tail_mask(nn) = .FALSE.
	3395	deleted_tails = deleted_tails + 1
	3396	ENDIF
	3397	ELSEIF ( ibc_par_ns == 2 ) THEN
	3398	!
	3399	!-- Particle reflection
	3400	particles(n)%y = -dy - particles(n)%y
	3401	particles(n)%speed_y = -particles(n)%speed_y
	3402	ENDIF
	3403
	3404	ELSEIF ( particles(n)%y >= ( ny + 0.5 ) * dy ) THEN
	3405
	3406	IF ( ibc_par_ns == 0 ) THEN
	3407	!
	3408	!-- Cyclic boundary. Relevant coordinate has to be changed.
	3409	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
	3410	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3411	i = particles(n)%tailpoints
	3412	particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
	3413	particle_tail_coordinates(1:i,2,nn)
	3414	ENDIF
	3415	ELSEIF ( ibc_par_ns == 1 ) THEN
	3416	!
	3417	!-- Particle absorption
	3418	particle_mask(n) = .FALSE.
	3419	deleted_particles = deleted_particles + 1
	3420	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3421	tail_mask(nn) = .FALSE.
	3422	deleted_tails = deleted_tails + 1
	3423	ENDIF
	3424	ELSEIF ( ibc_par_ns == 2 ) THEN
	3425	!
	3426	!-- Particle reflection
	3427	particles(n)%y = ( ny + 1 ) * dy - particles(n)%y
	3428	particles(n)%speed_y = -particles(n)%speed_y
	3429	ENDIF
	3430
	3431	ENDIF
	3432	ENDDO
	3433
	3434	#endif
	3435
	3436	!
	3437	!-- Apply boundary conditions to those particles that have crossed the top or
	3438	!-- bottom boundary and delete those particles, which are older than allowed
	3439	DO n = 1, number_of_particles
	3440
	3441	nn = particles(n)%tail_id
	3442
[57]	3443	!
	3444	!-- Stop if particles have moved further than the length of one
[262]	3445	!-- PE subdomain (newly released particles have age = age_m!)
	3446	IF ( particles(n)%age /= particles(n)%age_m ) THEN
	3447	IF ( ABS(particles(n)%speed_x) > &
	3448	((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. &
	3449	ABS(particles(n)%speed_y) > &
	3450	((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN
[57]	3451
[262]	3452	WRITE( message_string, * ) 'particle too fast. n = ', n
[276]	3453	CALL message( 'advec_particles', 'PA0148', 2, 2, -1, 6, 1 )
[262]	3454	ENDIF
[57]	3455	ENDIF
	3456
[1]	3457	IF ( particles(n)%age > particle_maximum_age .AND. &
	3458	particle_mask(n) ) &
	3459	THEN
	3460	particle_mask(n) = .FALSE.
	3461	deleted_particles = deleted_particles + 1
	3462	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3463	tail_mask(nn) = .FALSE.
	3464	deleted_tails = deleted_tails + 1
	3465	ENDIF
	3466	ENDIF
	3467
	3468	IF ( particles(n)%z >= zu(nz) .AND. particle_mask(n) ) THEN
	3469	IF ( ibc_par_t == 1 ) THEN
	3470	!
	3471	!-- Particle absorption
	3472	particle_mask(n) = .FALSE.
	3473	deleted_particles = deleted_particles + 1
	3474	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3475	tail_mask(nn) = .FALSE.
	3476	deleted_tails = deleted_tails + 1
	3477	ENDIF
	3478	ELSEIF ( ibc_par_t == 2 ) THEN
	3479	!
	3480	!-- Particle reflection
	3481	particles(n)%z = 2.0 * zu(nz) - particles(n)%z
	3482	particles(n)%speed_z = -particles(n)%speed_z
	3483	IF ( use_sgs_for_particles .AND. &
[824]	3484	particles(n)%rvar3 > 0.0 ) THEN
	3485	particles(n)%rvar3 = -particles(n)%rvar3
[1]	3486	ENDIF
	3487	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3488	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
	3489	particle_tail_coordinates(1,3,nn)
	3490	ENDIF
	3491	ENDIF
	3492	ENDIF
[150]	3493	IF ( particles(n)%z < zw(0) .AND. particle_mask(n) ) THEN
[1]	3494	IF ( ibc_par_b == 1 ) THEN
	3495	!
	3496	!-- Particle absorption
	3497	particle_mask(n) = .FALSE.
	3498	deleted_particles = deleted_particles + 1
	3499	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3500	tail_mask(nn) = .FALSE.
	3501	deleted_tails = deleted_tails + 1
	3502	ENDIF
	3503	ELSEIF ( ibc_par_b == 2 ) THEN
	3504	!
	3505	!-- Particle reflection
[150]	3506	particles(n)%z = 2.0 * zw(0) - particles(n)%z
[1]	3507	particles(n)%speed_z = -particles(n)%speed_z
	3508	IF ( use_sgs_for_particles .AND. &
[824]	3509	particles(n)%rvar3 < 0.0 ) THEN
	3510	particles(n)%rvar3 = -particles(n)%rvar3
[1]	3511	ENDIF
	3512	IF ( use_particle_tails .AND. nn /= 0 ) THEN
	3513	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
	3514	particle_tail_coordinates(1,3,nn)
	3515	ENDIF
	3516	IF ( use_particle_tails .AND. nn /= 0 ) THEN
[150]	3517	particle_tail_coordinates(1,3,nn) = 2.0 * zw(0) - &
	3518	particle_tail_coordinates(1,3,nn)
[1]	3519	ENDIF
	3520	ENDIF
	3521	ENDIF
	3522	ENDDO
	3523
	3524	! WRITE ( 9, * ) '*** advec_particles: ##7'
[82]	3525	! CALL local_flush( 9 )
[1]	3526	! nd = 0
	3527	! DO n = 1, number_of_particles
	3528	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
	3529	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	3530	! THEN
	3531	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	3532	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	3533	! CALL local_flush( 9 )
[1]	3534	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3535	! ENDIF
	3536	! ENDDO
	3537	! IF ( nd /= deleted_particles ) THEN
	3538	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
[82]	3539	! CALL local_flush( 9 )
[1]	3540	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3541	! ENDIF
	3542
	3543	!
	3544	!-- Pack particles (eliminate those marked for deletion),
	3545	!-- determine new number of particles
	3546	IF ( number_of_particles > 0 .AND. deleted_particles > 0 ) THEN
	3547	nn = 0
	3548	nd = 0
	3549	DO n = 1, number_of_particles
	3550	IF ( particle_mask(n) ) THEN
	3551	nn = nn + 1
	3552	particles(nn) = particles(n)
	3553	ELSE
	3554	nd = nd + 1
	3555	ENDIF
	3556	ENDDO
	3557	! IF ( nd /= deleted_particles ) THEN
	3558	! WRITE (9,) '** advec_part nd=',nd,' deleted_particles=',deleted_particles
[82]	3559	! CALL local_flush( 9 )
[1]	3560	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3561	! ENDIF
	3562
	3563	number_of_particles = number_of_particles - deleted_particles
	3564	!
	3565	!-- Pack the tails, store the new tail ids and re-assign it to the
	3566	!-- respective
	3567	!-- particles
	3568	IF ( use_particle_tails ) THEN
	3569	nn = 0
	3570	nd = 0
	3571	DO n = 1, number_of_tails
	3572	IF ( tail_mask(n) ) THEN
	3573	nn = nn + 1
	3574	particle_tail_coordinates(:,:,nn) = &
	3575	particle_tail_coordinates(:,:,n)
	3576	new_tail_id(n) = nn
	3577	ELSE
	3578	nd = nd + 1
	3579	! WRITE (9,*) '+++ n=',n,' (of ',number_of_tails,' #oftails)'
	3580	! WRITE (9,*) ' id=',new_tail_id(n)
[82]	3581	! CALL local_flush( 9 )
[1]	3582	ENDIF
	3583	ENDDO
	3584	ENDIF
	3585
	3586	! IF ( nd /= deleted_tails .AND. use_particle_tails ) THEN
	3587	! WRITE (9,) '** advec_part nd=',nd,' deleted_tails=',deleted_tails
[82]	3588	! CALL local_flush( 9 )
[1]	3589	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3590	! ENDIF
	3591
	3592	number_of_tails = number_of_tails - deleted_tails
	3593
	3594	! nn = 0
	3595	DO n = 1, number_of_particles
	3596	IF ( particles(n)%tail_id /= 0 ) THEN
	3597	! nn = nn + 1
	3598	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id > number_of_tails ) THEN
	3599	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	3600	! WRITE (9,*) ' tail_id=',particles(n)%tail_id
	3601	! WRITE (9,*) ' new_tail_id=', new_tail_id(particles(n)%tail_id), &
	3602	! ' of (',number_of_tails,')'
[82]	3603	! CALL local_flush( 9 )
[1]	3604	! ENDIF
	3605	particles(n)%tail_id = new_tail_id(particles(n)%tail_id)
	3606	ENDIF
	3607	ENDDO
	3608
	3609	! IF ( nn /= number_of_tails .AND. use_particle_tails ) THEN
	3610	! WRITE (9,) '** advec_part #of_tails=',number_of_tails,' nn=',nn
[82]	3611	! CALL local_flush( 9 )
[1]	3612	! DO n = 1, number_of_particles
	3613	! WRITE (9,*) 'prt# ',n,' tail_id=',particles(n)%tail_id, &
	3614	! ' x=',particles(n)%x, ' y=',particles(n)%y, &
	3615	! ' z=',particles(n)%z
	3616	! ENDDO
	3617	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3618	! ENDIF
	3619
	3620	ENDIF
	3621
	3622	! IF ( number_of_particles /= number_of_tails ) THEN
	3623	! WRITE (9,*) '--- advec_particles: #7'
	3624	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
[82]	3625	! CALL local_flush( 9 )
[1]	3626	! ENDIF
	3627	! WRITE ( 9, * ) '*** advec_particles: ##8'
[82]	3628	! CALL local_flush( 9 )
[1]	3629	! DO n = 1, number_of_particles
	3630	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	3631	! THEN
	3632	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	3633	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	3634	! CALL local_flush( 9 )
[1]	3635	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3636	! ENDIF
	3637	! ENDDO
	3638
	3639	! WRITE ( 9, * ) '*** advec_particles: ##9'
[82]	3640	! CALL local_flush( 9 )
[1]	3641	! DO n = 1, number_of_particles
	3642	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	3643	! THEN
	3644	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	3645	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	3646	! CALL local_flush( 9 )
[1]	3647	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3648	! ENDIF
	3649	! ENDDO
	3650
	3651	!
	3652	!-- Accumulate the number of particles transferred between the subdomains
	3653	#if defined( __parallel )
	3654	trlp_count_sum = trlp_count_sum + trlp_count
	3655	trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
	3656	trrp_count_sum = trrp_count_sum + trrp_count
	3657	trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
	3658	trsp_count_sum = trsp_count_sum + trsp_count
	3659	trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
	3660	trnp_count_sum = trnp_count_sum + trnp_count
	3661	trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
	3662	#endif
	3663
	3664	IF ( dt_3d_reached ) EXIT
	3665
[420]	3666	!
	3667	!-- Initialize variables for the next (sub-) timestep, i.e. for marking those
	3668	!-- particles to be deleted after the timestep
	3669	particle_mask = .TRUE.
	3670	deleted_particles = 0
	3671	trlp_count_recv = 0
	3672	trnp_count_recv = 0
	3673	trrp_count_recv = 0
	3674	trsp_count_recv = 0
	3675	trlpt_count_recv = 0
	3676	trnpt_count_recv = 0
	3677	trrpt_count_recv = 0
	3678	trspt_count_recv = 0
	3679	IF ( use_particle_tails ) THEN
	3680	tail_mask = .TRUE.
	3681	ENDIF
	3682	deleted_tails = 0
	3683
[1]	3684	ENDDO ! timestep loop
	3685
	3686	!
[116]	3687	!-- Sort particles in the sequence the gridboxes are stored in the memory
	3688	time_sort_particles = time_sort_particles + dt_3d
	3689	IF ( time_sort_particles >= dt_sort_particles ) THEN
	3690	CALL sort_particles
	3691	time_sort_particles = MOD( time_sort_particles, &
	3692	MAX( dt_sort_particles, dt_3d ) )
	3693	ENDIF
	3694
[1]	3695	IF ( cloud_droplets ) THEN
	3696
	3697	CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'start' )
	3698
	3699	ql = 0.0; ql_v = 0.0; ql_vp = 0.0
	3700
	3701	!
[420]	3702	!-- Calculate the liquid water content
[1]	3703	DO i = nxl, nxr
	3704	DO j = nys, nyn
	3705	DO k = nzb, nzt+1
	3706
	3707	!
[420]	3708	!-- Calculate the total volume in the boxes (ql_v, weighting factor
	3709	!-- has to beincluded)
[1]	3710	psi = prt_start_index(k,j,i)
	3711	DO n = psi, psi+prt_count(k,j,i)-1
[420]	3712	ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
[1]	3713	particles(n)%radius**3
	3714	ENDDO
	3715
	3716	!
[420]	3717	!-- Calculate the liquid water content
	3718	IF ( ql_v(k,j,i) /= 0.0 ) THEN
	3719	ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * &
	3720	ql_v(k,j,i) / &
[1]	3721	( rho_surface * dx * dy * dz )
[420]	3722
	3723	IF ( ql(k,j,i) < 0.0 ) THEN
	3724	WRITE( message_string, * ) 'LWC out of range: ' , &
	3725	ql(k,j,i)
	3726	CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 )
	3727	ENDIF
	3728
[1]	3729	ELSE
	3730	ql(k,j,i) = 0.0
	3731	ENDIF
	3732
	3733	ENDDO
	3734	ENDDO
	3735	ENDDO
	3736
	3737	CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'stop' )
	3738
	3739	ENDIF
	3740
	3741	!
[264]	3742	!-- Set particle attributes
	3743	CALL set_particle_attributes
	3744
	3745	!
[1]	3746	!-- Set particle attributes defined by the user
	3747	CALL user_particle_attributes
	3748	! WRITE ( 9, * ) '*** advec_particles: ##10'
[82]	3749	! CALL local_flush( 9 )
[1]	3750	! DO n = 1, number_of_particles
	3751	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
	3752	! THEN
	3753	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
	3754	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
[82]	3755	! CALL local_flush( 9 )
[1]	3756	! CALL MPI_ABORT( comm2d, 9999, ierr )
	3757	! ENDIF
	3758	! ENDDO
	3759
	3760	!
	3761	!-- If necessary, add the actual particle positions to the particle tails
	3762	IF ( use_particle_tails ) THEN
	3763
	3764	distance = 0.0
	3765	DO n = 1, number_of_particles
	3766
	3767	nn = particles(n)%tail_id
	3768
	3769	IF ( nn /= 0 ) THEN
	3770	!
	3771	!-- Calculate the distance between the actual particle position and the
	3772	!-- next tailpoint
	3773	! WRITE ( 9, * ) '*** advec_particles: ##10.1 nn=',nn
[82]	3774	! CALL local_flush( 9 )
[1]	3775	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
	3776	distance = ( particle_tail_coordinates(1,1,nn) - &
	3777	particle_tail_coordinates(2,1,nn) )**2 + &
	3778	( particle_tail_coordinates(1,2,nn) - &
	3779	particle_tail_coordinates(2,2,nn) )**2 + &
	3780	( particle_tail_coordinates(1,3,nn) - &
	3781	particle_tail_coordinates(2,3,nn) )**2
	3782	ENDIF
	3783	! WRITE ( 9, * ) '*** advec_particles: ##10.2'
[82]	3784	! CALL local_flush( 9 )
[1]	3785	!
	3786	!-- First, increase the index of all existings tailpoints by one
	3787	IF ( distance >= minimum_tailpoint_distance ) THEN
	3788	DO i = particles(n)%tailpoints, 1, -1
	3789	particle_tail_coordinates(i+1,:,nn) = &
	3790	particle_tail_coordinates(i,:,nn)
	3791	ENDDO
	3792	!
	3793	!-- Increase the counter which contains the number of tailpoints.
	3794	!-- This must always be smaller than the given maximum number of
	3795	!-- tailpoints because otherwise the index bounds of
	3796	!-- particle_tail_coordinates would be exceeded
	3797	IF ( particles(n)%tailpoints < maximum_number_of_tailpoints-1 )&
	3798	THEN
	3799	particles(n)%tailpoints = particles(n)%tailpoints + 1
	3800	ENDIF
	3801	ENDIF
	3802	! WRITE ( 9, * ) '*** advec_particles: ##10.3'
[82]	3803	! CALL local_flush( 9 )
[1]	3804	!
	3805	!-- In any case, store the new point at the beginning of the tail
	3806	particle_tail_coordinates(1,1,nn) = particles(n)%x
	3807	particle_tail_coordinates(1,2,nn) = particles(n)%y
	3808	particle_tail_coordinates(1,3,nn) = particles(n)%z
	3809	particle_tail_coordinates(1,4,nn) = particles(n)%color
	3810	! WRITE ( 9, * ) '*** advec_particles: ##10.4'
[82]	3811	! CALL local_flush( 9 )
[1]	3812	!
	3813	!-- Increase the age of the tailpoints
	3814	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
	3815	particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) = &
	3816	particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) + &
	3817	dt_3d
	3818	!
	3819	!-- Delete the last tailpoint, if it has exceeded its maximum age
	3820	IF ( particle_tail_coordinates(particles(n)%tailpoints,5,nn) > &
	3821	maximum_tailpoint_age ) THEN
	3822	particles(n)%tailpoints = particles(n)%tailpoints - 1
	3823	ENDIF
	3824	ENDIF
	3825	! WRITE ( 9, * ) '*** advec_particles: ##10.5'
[82]	3826	! CALL local_flush( 9 )
[1]	3827
	3828	ENDIF
	3829
	3830	ENDDO
	3831
	3832	ENDIF
	3833	! WRITE ( 9, * ) '*** advec_particles: ##11'
[82]	3834	! CALL local_flush( 9 )
[1]	3835
	3836	!
	3837	!-- Write particle statistics on file
	3838	IF ( write_particle_statistics ) THEN
	3839	CALL check_open( 80 )
	3840	#if defined( __parallel )
	3841	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
	3842	number_of_particles, pleft, trlp_count_sum, &
	3843	trlp_count_recv_sum, pright, trrp_count_sum, &
	3844	trrp_count_recv_sum, psouth, trsp_count_sum, &
	3845	trsp_count_recv_sum, pnorth, trnp_count_sum, &
	3846	trnp_count_recv_sum, maximum_number_of_particles
	3847	CALL close_file( 80 )
	3848	#else
	3849	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
	3850	number_of_particles, maximum_number_of_particles
	3851	#endif
	3852	ENDIF
	3853
	3854	CALL cpu_log( log_point(25), 'advec_particles', 'stop' )
	3855
	3856	!
	3857	!-- Formats
	3858	8000 FORMAT (I6,1X,F7.2,4X,I6,5X,4(I3,1X,I4,'/',I4,2X),6X,I6)
	3859
	3860	END SUBROUTINE advec_particles
	3861
	3862
	3863	SUBROUTINE allocate_prt_memory( number_of_new_particles )
	3864
	3865	!------------------------------------------------------------------------------!
	3866	! Description:
	3867	! ------------
	3868	! Extend particle memory
	3869	!------------------------------------------------------------------------------!
	3870
	3871	USE particle_attributes
	3872
	3873	IMPLICIT NONE
	3874
	3875	INTEGER :: new_maximum_number, number_of_new_particles
	3876
	3877	LOGICAL, DIMENSION(:), ALLOCATABLE :: tmp_particle_mask
	3878
	3879	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles
	3880
	3881
	3882	new_maximum_number = maximum_number_of_particles + &
	3883	MAX( 5*number_of_new_particles, number_of_initial_particles )
	3884
	3885	IF ( write_particle_statistics ) THEN
	3886	CALL check_open( 80 )
	3887	WRITE ( 80, '(''*** Request: '', I7, '' new_maximum_number(prt)'')' ) &
	3888	new_maximum_number
	3889	CALL close_file( 80 )
	3890	ENDIF
	3891
	3892	ALLOCATE( tmp_particles(maximum_number_of_particles), &
	3893	tmp_particle_mask(maximum_number_of_particles) )
	3894	tmp_particles = particles
	3895	tmp_particle_mask = particle_mask
	3896
	3897	DEALLOCATE( particles, particle_mask )
	3898	ALLOCATE( particles(new_maximum_number), &
	3899	particle_mask(new_maximum_number) )
	3900	maximum_number_of_particles = new_maximum_number
	3901
	3902	particles(1:number_of_particles) = tmp_particles(1:number_of_particles)
	3903	particle_mask(1:number_of_particles) = &
	3904	tmp_particle_mask(1:number_of_particles)
	3905	particle_mask(number_of_particles+1:maximum_number_of_particles) = .TRUE.
	3906	DEALLOCATE( tmp_particles, tmp_particle_mask )
	3907
	3908	END SUBROUTINE allocate_prt_memory
	3909
	3910
	3911	SUBROUTINE allocate_tail_memory( number_of_new_tails )
	3912
	3913	!------------------------------------------------------------------------------!
	3914	! Description:
	3915	! ------------
	3916	! Extend tail memory
	3917	!------------------------------------------------------------------------------!
	3918
	3919	USE particle_attributes
	3920
	3921	IMPLICIT NONE
	3922
	3923	INTEGER :: new_maximum_number, number_of_new_tails
	3924
	3925	LOGICAL, DIMENSION(maximum_number_of_tails) :: tmp_tail_mask
	3926
	3927	REAL, DIMENSION(maximum_number_of_tailpoints,5,maximum_number_of_tails) :: &
	3928	tmp_tail
	3929
	3930
	3931	new_maximum_number = maximum_number_of_tails + &
	3932	MAX( 5*number_of_new_tails, number_of_initial_tails )
	3933
	3934	IF ( write_particle_statistics ) THEN
	3935	CALL check_open( 80 )
	3936	WRITE ( 80, '(''*** Request: '', I5, '' new_maximum_number(tails)'')' ) &
	3937	new_maximum_number
	3938	CALL close_file( 80 )
	3939	ENDIF
	3940	WRITE (9,) '** Request: ',new_maximum_number,' new_maximum_number(tails)'
[82]	3941	! CALL local_flush( 9 )
[1]	3942
	3943	tmp_tail(:,:,1:number_of_tails) = &
	3944	particle_tail_coordinates(:,:,1:number_of_tails)
	3945	tmp_tail_mask(1:number_of_tails) = tail_mask(1:number_of_tails)
	3946
	3947	DEALLOCATE( new_tail_id, particle_tail_coordinates, tail_mask )
	3948	ALLOCATE( new_tail_id(new_maximum_number), &
	3949	particle_tail_coordinates(maximum_number_of_tailpoints,5, &
	3950	new_maximum_number), &
	3951	tail_mask(new_maximum_number) )
	3952	maximum_number_of_tails = new_maximum_number
	3953
	3954	particle_tail_coordinates = 0.0
	3955	particle_tail_coordinates(:,:,1:number_of_tails) = &
	3956	tmp_tail(:,:,1:number_of_tails)
	3957	tail_mask(1:number_of_tails) = tmp_tail_mask(1:number_of_tails)
	3958	tail_mask(number_of_tails+1:maximum_number_of_tails) = .TRUE.
	3959
	3960	END SUBROUTINE allocate_tail_memory
	3961
	3962
	3963	SUBROUTINE output_particles_netcdf
	3964	#if defined( __netcdf )
	3965
	3966	USE control_parameters
	3967	USE netcdf_control
	3968	USE particle_attributes
	3969
	3970	IMPLICIT NONE
	3971
	3972
	3973	CALL check_open( 108 )
	3974
	3975	!
	3976	!-- Update the NetCDF time axis
	3977	prt_time_count = prt_time_count + 1
	3978
	3979	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, (/ simulated_time /), &
	3980	start = (/ prt_time_count /), count = (/ 1 /) )
[263]	3981	CALL handle_netcdf_error( 'output_particles_netcdf', 1 )
[1]	3982
	3983	!
	3984	!-- Output the real number of particles used
	3985	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, &
	3986	(/ number_of_particles /), &
	3987	start = (/ prt_time_count /), count = (/ 1 /) )
[263]	3988	CALL handle_netcdf_error( 'output_particles_netcdf', 2 )
[1]	3989
	3990	!
	3991	!-- Output all particle attributes
	3992	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age, &
	3993	start = (/ 1, prt_time_count /), &
	3994	count = (/ maximum_number_of_particles /) )
[263]	3995	CALL handle_netcdf_error( 'output_particles_netcdf', 3 )
[1]	3996
	3997	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize, &
	3998	start = (/ 1, prt_time_count /), &
	3999	count = (/ maximum_number_of_particles /) )
[263]	4000	CALL handle_netcdf_error( 'output_particles_netcdf', 4 )
[1]	4001
	4002	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x, &
	4003	start = (/ 1, prt_time_count /), &
	4004	count = (/ maximum_number_of_particles /) )
[263]	4005	CALL handle_netcdf_error( 'output_particles_netcdf', 5 )
[1]	4006
	4007	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y, &
	4008	start = (/ 1, prt_time_count /), &
	4009	count = (/ maximum_number_of_particles /) )
[263]	4010	CALL handle_netcdf_error( 'output_particles_netcdf', 6 )
[1]	4011
	4012	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z, &
	4013	start = (/ 1, prt_time_count /), &
	4014	count = (/ maximum_number_of_particles /) )
[263]	4015	CALL handle_netcdf_error( 'output_particles_netcdf', 7 )
[1]	4016
	4017	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius, &
	4018	start = (/ 1, prt_time_count /), &
	4019	count = (/ maximum_number_of_particles /) )
[263]	4020	CALL handle_netcdf_error( 'output_particles_netcdf', 8 )
[1]	4021
	4022	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x, &
	4023	start = (/ 1, prt_time_count /), &
	4024	count = (/ maximum_number_of_particles /) )
[263]	4025	CALL handle_netcdf_error( 'output_particles_netcdf', 9 )
[1]	4026
	4027	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y, &
	4028	start = (/ 1, prt_time_count /), &
	4029	count = (/ maximum_number_of_particles /) )
[263]	4030	CALL handle_netcdf_error( 'output_particles_netcdf', 10 )
[1]	4031
	4032	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z, &
	4033	start = (/ 1, prt_time_count /), &
	4034	count = (/ maximum_number_of_particles /) )
[263]	4035	CALL handle_netcdf_error( 'output_particles_netcdf', 11 )
[1]	4036
	4037	nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10),particles%weight_factor,&
	4038	start = (/ 1, prt_time_count /), &
	4039	count = (/ maximum_number_of_particles /) )
[263]	4040	CALL handle_netcdf_error( 'output_particles_netcdf', 12 )
[1]	4041
[274]	4042	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x, &
[1]	4043	start = (/ 1, prt_time_count /), &
	4044	count = (/ maximum_number_of_particles /) )
[263]	4045	CALL handle_netcdf_error( 'output_particles_netcdf', 13 )
[1]	4046
	4047	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y, &
	4048	start = (/ 1, prt_time_count /), &
	4049	count = (/ maximum_number_of_particles /) )
[263]	4050	CALL handle_netcdf_error( 'output_particles_netcdf', 14 )
[1]	4051
	4052	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z, &
	4053	start = (/ 1, prt_time_count /), &
	4054	count = (/ maximum_number_of_particles /) )
[263]	4055	CALL handle_netcdf_error( 'output_particles_netcdf', 15 )
[1]	4056
	4057	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%color, &
	4058	start = (/ 1, prt_time_count /), &
	4059	count = (/ maximum_number_of_particles /) )
[263]	4060	CALL handle_netcdf_error( 'output_particles_netcdf', 16 )
[1]	4061
	4062	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group, &
	4063	start = (/ 1, prt_time_count /), &
	4064	count = (/ maximum_number_of_particles /) )
[263]	4065	CALL handle_netcdf_error( 'output_particles_netcdf', 17 )
[1]	4066
	4067	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16), particles%tailpoints, &
	4068	start = (/ 1, prt_time_count /), &
	4069	count = (/ maximum_number_of_particles /) )
[263]	4070	CALL handle_netcdf_error( 'output_particles_netcdf', 18 )
[1]	4071
	4072	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id, &
	4073	start = (/ 1, prt_time_count /), &
	4074	count = (/ maximum_number_of_particles /) )
[263]	4075	CALL handle_netcdf_error( 'output_particles_netcdf', 19 )
[1]	4076
	4077	#endif
	4078	END SUBROUTINE output_particles_netcdf
	4079
	4080
	4081	SUBROUTINE write_particles
	4082
	4083	!------------------------------------------------------------------------------!
	4084	! Description:
	4085	! ------------
	4086	! Write particle data on restart file
	4087	!------------------------------------------------------------------------------!
	4088
	4089	USE control_parameters
	4090	USE particle_attributes
	4091	USE pegrid
	4092
	4093	IMPLICIT NONE
	4094
	4095	CHARACTER (LEN=10) :: particle_binary_version
[759]	4096	INTEGER :: i
	4097
[1]	4098	!
	4099	!-- First open the output unit.
	4100	IF ( myid_char == '' ) THEN
	4101	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, &
	4102	FORM='UNFORMATTED')
	4103	ELSE
	4104	IF ( myid == 0 ) CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' )
	4105	#if defined( __parallel )
	4106	!
	4107	!-- Set a barrier in order to allow that thereafter all other processors
	4108	!-- in the directory created by PE0 can open their file
	4109	CALL MPI_BARRIER( comm2d, ierr )
	4110	#endif
	4111	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, &
	4112	FORM='UNFORMATTED' )
	4113	ENDIF
	4114
[759]	4115	DO i = 0, io_blocks-1
	4116
	4117	IF ( i == io_group ) THEN
	4118
[1]	4119	!
[759]	4120	!-- Write the version number of the binary format.
	4121	!-- Attention: After changes to the following output commands the version
	4122	!-- --------- number of the variable particle_binary_version must be
	4123	!-- changed! Also, the version number and the list of arrays
	4124	!-- to be read in init_particles must be adjusted accordingly.
	4125	particle_binary_version = '3.0'
	4126	WRITE ( 90 ) particle_binary_version
[1]	4127
	4128	!
[759]	4129	!-- Write some particle parameters, the size of the particle arrays as
	4130	!-- well as other dvrp-plot variables.
	4131	WRITE ( 90 ) bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, &
	4132	maximum_number_of_particles, &
	4133	maximum_number_of_tailpoints, maximum_number_of_tails, &
	4134	number_of_initial_particles, number_of_particles, &
	4135	number_of_particle_groups, number_of_tails, &
	4136	particle_groups, time_prel, time_write_particle_data, &
	4137	uniform_particles
[1]	4138
[759]	4139	IF ( number_of_initial_particles /= 0 ) WRITE ( 90 ) initial_particles
[1]	4140
[759]	4141	WRITE ( 90 ) prt_count, prt_start_index
	4142	WRITE ( 90 ) particles
[1]	4143
[759]	4144	IF ( use_particle_tails ) THEN
	4145	WRITE ( 90 ) particle_tail_coordinates
	4146	ENDIF
[1]	4147
[759]	4148	CLOSE ( 90 )
[1]	4149
[759]	4150	ENDIF
	4151
	4152	#if defined( __parallel )
	4153	CALL MPI_BARRIER( comm2d, ierr )
	4154	#endif
	4155
	4156	ENDDO
	4157
[1]	4158	END SUBROUTINE write_particles
	4159
	4160
[799]	4161
[1]	4162	SUBROUTINE collision_efficiency( mean_r, r, e)
	4163	!------------------------------------------------------------------------------!
	4164	! Description:
	4165	! ------------
	4166	! Interpolate collision efficiency from table
	4167	!------------------------------------------------------------------------------!
	4168
	4169	IMPLICIT NONE
	4170
	4171	INTEGER :: i, j, k
	4172
	4173	LOGICAL, SAVE :: first = .TRUE.
	4174
	4175	REAL :: aa, bb, cc, dd, dx, dy, e, gg, mean_r, mean_rm, r, rm, &
	4176	x, y
	4177
	4178	REAL, DIMENSION(1:9), SAVE :: collected_r = 0.0
	4179	REAL, DIMENSION(1:19), SAVE :: collector_r = 0.0
	4180	REAL, DIMENSION(1:9,1:19), SAVE :: ef = 0.0
	4181
	4182	mean_rm = mean_r * 1.0E06
	4183	rm = r * 1.0E06
	4184
	4185	IF ( first ) THEN
	4186	collected_r = (/ 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0 /)
	4187	collector_r = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0, 150.0,&
	4188	200.0, 300.0, 400.0, 500.0, 600.0, 1000.0, 1400.0, &
	4189	1800.0, 2400.0, 3000.0 /)
	4190	ef(:,1) = (/0.017, 0.027, 0.037, 0.052, 0.052, 0.052, 0.052, 0.0, 0.0 /)
	4191	ef(:,2) = (/0.001, 0.016, 0.027, 0.060, 0.12, 0.17, 0.17, 0.17, 0.0 /)
	4192	ef(:,3) = (/0.001, 0.001, 0.02, 0.13, 0.28, 0.37, 0.54, 0.55, 0.47/)
	4193	ef(:,4) = (/0.001, 0.001, 0.02, 0.23, 0.4, 0.55, 0.7, 0.75, 0.75/)
	4194	ef(:,5) = (/0.01, 0.01, 0.03, 0.3, 0.4, 0.58, 0.73, 0.75, 0.79/)
	4195	ef(:,6) = (/0.01, 0.01, 0.13, 0.38, 0.57, 0.68, 0.80, 0.86, 0.91/)
	4196	ef(:,7) = (/0.01, 0.085, 0.23, 0.52, 0.68, 0.76, 0.86, 0.92, 0.95/)
	4197	ef(:,8) = (/0.01, 0.14, 0.32, 0.60, 0.73, 0.81, 0.90, 0.94, 0.96/)
	4198	ef(:,9) = (/0.025, 0.25, 0.43, 0.66, 0.78, 0.83, 0.92, 0.95, 0.96/)
	4199	ef(:,10)= (/0.039, 0.3, 0.46, 0.69, 0.81, 0.87, 0.93, 0.95, 0.96/)
	4200	ef(:,11)= (/0.095, 0.33, 0.51, 0.72, 0.82, 0.87, 0.93, 0.96, 0.97/)
	4201	ef(:,12)= (/0.098, 0.36, 0.51, 0.73, 0.83, 0.88, 0.93, 0.96, 0.97/)
	4202	ef(:,13)= (/0.1, 0.36, 0.52, 0.74, 0.83, 0.88, 0.93, 0.96, 0.97/)
	4203	ef(:,14)= (/0.17, 0.4, 0.54, 0.72, 0.83, 0.88, 0.94, 0.98, 1.0 /)
	4204	ef(:,15)= (/0.15, 0.37, 0.52, 0.74, 0.82, 0.88, 0.94, 0.98, 1.0 /)
	4205	ef(:,16)= (/0.11, 0.34, 0.49, 0.71, 0.83, 0.88, 0.94, 0.95, 1.0 /)
	4206	ef(:,17)= (/0.08, 0.29, 0.45, 0.68, 0.8, 0.86, 0.96, 0.94, 1.0 /)
	4207	ef(:,18)= (/0.04, 0.22, 0.39, 0.62, 0.75, 0.83, 0.92, 0.96, 1.0 /)
	4208	ef(:,19)= (/0.02, 0.16, 0.33, 0.55, 0.71, 0.81, 0.90, 0.94, 1.0 /)
	4209	ENDIF
	4210
	4211	DO k = 1, 8
	4212	IF ( collected_r(k) <= mean_rm ) i = k
	4213	ENDDO
	4214
	4215	DO k = 1, 18
	4216	IF ( collector_r(k) <= rm ) j = k
	4217	ENDDO
	4218
	4219	IF ( rm < 10.0 ) THEN
	4220	e = 0.0
	4221	ELSEIF ( mean_rm < 2.0 ) THEN
	4222	e = 0.001
	4223	ELSEIF ( mean_rm >= 25.0 ) THEN
[420]	4224	IF( j <= 2 ) e = 0.0
	4225	IF( j == 3 ) e = 0.47
[1]	4226	IF( j == 4 ) e = 0.8
	4227	IF( j == 5 ) e = 0.9
	4228	IF( j >=6 ) e = 1.0
	4229	ELSEIF ( rm >= 3000.0 ) THEN
[420]	4230	IF( i == 1 ) e = 0.02
	4231	IF( i == 2 ) e = 0.16
	4232	IF( i == 3 ) e = 0.33
	4233	IF( i == 4 ) e = 0.55
	4234	IF( i == 5 ) e = 0.71
	4235	IF( i == 6 ) e = 0.81
	4236	IF( i == 7 ) e = 0.90
	4237	IF( i >= 8 ) e = 0.94
[1]	4238	ELSE
	4239	x = mean_rm - collected_r(i)
[420]	4240	y = rm - collector_r(j)
	4241	dx = collected_r(i+1) - collected_r(i)
	4242	dy = collector_r(j+1) - collector_r(j)
[1]	4243	aa = x2 + y2
	4244	bb = ( dx - x )2 + y2
	4245	cc = x2 + ( dy - y )2
	4246	dd = ( dx - x )2 + ( dy - y )2
	4247	gg = aa + bb + cc + dd
	4248
	4249	e = ( (gg-aa)ef(i,j) + (gg-bb)ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
	4250	(gg-dd)ef(i+1,j+1) ) / (3.0gg)
	4251	ENDIF
	4252
[420]	4253	END SUBROUTINE collision_efficiency
[1]	4254
	4255
	4256
	4257	SUBROUTINE sort_particles
	4258
	4259	!------------------------------------------------------------------------------!
	4260	! Description:
	4261	! ------------
	4262	! Sort particles in the sequence the grid boxes are stored in memory
	4263	!------------------------------------------------------------------------------!
	4264
	4265	USE arrays_3d
	4266	USE control_parameters
	4267	USE cpulog
	4268	USE grid_variables
	4269	USE indices
	4270	USE interfaces
	4271	USE particle_attributes
	4272
	4273	IMPLICIT NONE
	4274
	4275	INTEGER :: i, ilow, j, k, n
	4276
[792]	4277	TYPE(particle_type), DIMENSION(:), POINTER :: particles_temp
[1]	4278
	4279
	4280	CALL cpu_log( log_point_s(47), 'sort_particles', 'start' )
	4281
	4282	!
[792]	4283	!-- Initialize counters and set pointer of the temporary array into which
	4284	!-- particles are sorted to free memory
	4285	prt_count = 0
	4286	sort_count = sort_count +1
[1]	4287
[792]	4288	SELECT CASE ( MOD( sort_count, 2 ) )
	4289
	4290	CASE ( 0 )
	4291
	4292	particles_temp => part_1
[824]	4293	! part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	4294	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	4295	! 0.0, 0, 0, 0, 0 )
[792]	4296
	4297	CASE ( 1 )
	4298
	4299	particles_temp => part_2
[824]	4300	! part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	4301	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	4302	! 0.0, 0, 0, 0, 0 )
[792]	4303
	4304	END SELECT
	4305
[1]	4306	!
	4307	!-- Count the particles per gridbox
	4308	DO n = 1, number_of_particles
	4309
	4310	i = ( particles(n)%x + 0.5 * dx ) * ddx
	4311	j = ( particles(n)%y + 0.5 * dy ) * ddy
[212]	4312	k = particles(n)%z / dz + 1 + offset_ocean_nzt
[150]	4313	! only exact if equidistant
[1]	4314
	4315	prt_count(k,j,i) = prt_count(k,j,i) + 1
	4316
	4317	IF ( i < nxl .OR. i > nxr .OR. j < nys .OR. j > nyn .OR. k < nzb+1 .OR. &
	4318	k > nzt ) THEN
[230]	4319	WRITE( message_string, * ) ' particle out of range: i=', i, ' j=', &
[274]	4320	j, ' k=', k, &
[420]	4321	' nxl=', nxl, ' nxr=', nxr, &
[274]	4322	' nys=', nys, ' nyn=', nyn, &
[230]	4323	' nzb=', nzb, ' nzt=', nzt
	4324	CALL message( 'sort_particles', 'PA0149', 1, 2, 0, 6, 0 )
[1]	4325	ENDIF
	4326
	4327	ENDDO
	4328
	4329	!
	4330	!-- Calculate the lower indices of those ranges of the particles-array
	4331	!-- containing particles which belong to the same gridpox i,j,k
	4332	ilow = 1
	4333	DO i = nxl, nxr
	4334	DO j = nys, nyn
	4335	DO k = nzb+1, nzt
	4336	prt_start_index(k,j,i) = ilow
	4337	ilow = ilow + prt_count(k,j,i)
	4338	ENDDO
	4339	ENDDO
	4340	ENDDO
	4341
	4342	!
	4343	!-- Sorting the particles
	4344	DO n = 1, number_of_particles
	4345
	4346	i = ( particles(n)%x + 0.5 * dx ) * ddx
	4347	j = ( particles(n)%y + 0.5 * dy ) * ddy
[212]	4348	k = particles(n)%z / dz + 1 + offset_ocean_nzt
[150]	4349	! only exact if equidistant
[1]	4350
	4351	particles_temp(prt_start_index(k,j,i)) = particles(n)
	4352
	4353	prt_start_index(k,j,i) = prt_start_index(k,j,i) + 1
	4354
	4355	ENDDO
	4356
[792]	4357	!
	4358	!-- Redirect the particles pointer to the sorted array
	4359	SELECT CASE ( MOD( sort_count, 2 ) )
[1]	4360
[792]	4361	CASE ( 0 )
	4362
	4363	particles => part_1
	4364
	4365	CASE ( 1 )
	4366
	4367	particles => part_2
	4368
	4369	END SELECT
	4370
[1]	4371	!
	4372	!-- Reset the index array to the actual start position
	4373	DO i = nxl, nxr
	4374	DO j = nys, nyn
	4375	DO k = nzb+1, nzt
	4376	prt_start_index(k,j,i) = prt_start_index(k,j,i) - prt_count(k,j,i)
	4377	ENDDO
	4378	ENDDO
	4379	ENDDO
	4380
	4381	CALL cpu_log( log_point_s(47), 'sort_particles', 'stop' )
	4382
	4383	END SUBROUTINE sort_particles

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

Context Navigation

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

Download in other formats: