Home

Context Navigation

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

Last change on this file since 236 was 230, checked in by heinze, 16 years ago
Output of messages replaced by message handling routine
Property svn:keywords set to `Id`
File size: 165.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |