Home

Context Navigation

lagrangian_particle_model_mod.f90 @ 4444

Last change on this file since 4444 was 4444, checked in by raasch, 5 years ago
bugfix: cpp-directives for serial mode added
Property svn:keywords set to `Id`
File size: 356.7 KB

Line
1	!> @file lagrangian_particle_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2020 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: lagrangian_particle_model_mod.f90 4444 2020-03-05 15:59:50Z raasch $
27	! bugfix: cpp-directives for serial mode added
28	!
29	! 4430 2020-02-27 18:02:20Z suehring
30	! - Bugfix in logarithmic interpolation of near-ground particle speed (density
31	! was not considered).
32	! - Revise CFL-check when SGS particle speeds are considered.
33	! - In nested case with SGS particle speeds in the child domain, do not give
34	! warning that particles are on domain boundaries. At the end of the particle
35	! time integration these will be transferred to the parent domain anyhow.
36	!
37	! 4360 2020-01-07 11:25:50Z suehring
38	! Introduction of wall_flags_total_0, which currently sets bits based on static
39	! topography information used in wall_flags_static_0
40	!
41	! 4336 2019-12-13 10:12:05Z raasch
42	! bugfix: wrong header output of particle group features (density ratio) in case
43	! of restarts corrected
44	!
45	! 4329 2019-12-10 15:46:36Z motisi
46	! Renamed wall_flags_0 to wall_flags_static_0
47	!
48	! 4282 2019-10-29 16:18:46Z schwenkel
49	! Bugfix of particle timeseries in case of more than one particle group
50	!
51	! 4277 2019-10-28 16:53:23Z schwenkel
52	! Bugfix: Added first_call_lpm in use statement
53	!
54	! 4276 2019-10-28 16:03:29Z schwenkel
55	! Modularize lpm: Move conditions in time intergration to module
56	!
57	! 4275 2019-10-28 15:34:55Z schwenkel
58	! Change call of simple predictor corrector method, i.e. two divergence free
59	! velocitiy fields are now used.
60	!
61	! 4232 2019-09-20 09:34:22Z knoop
62	! Removed INCLUDE "mpif.h", as it is not needed because of USE pegrid
63	!
64	! 4195 2019-08-28 13:44:27Z schwenkel
65	! Bugfix for simple_corrector interpolation method in case of ocean runs and
66	! output particle advection interpolation method into header
67	!
68	! 4182 2019-08-22 15:20:23Z scharf
69	! Corrected "Former revisions" section
70	!
71	! 4168 2019-08-16 13:50:17Z suehring
72	! Replace function get_topography_top_index by topo_top_ind
73	!
74	! 4145 2019-08-06 09:55:22Z schwenkel
75	! Some reformatting
76	!
77	! 4144 2019-08-06 09:11:47Z raasch
78	! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
79	!
80	! 4143 2019-08-05 15:14:53Z schwenkel
81	! Rename variable and change select case to if statement
82	!
83	! 4122 2019-07-26 13:11:56Z schwenkel
84	! Implement reset method as bottom boundary condition
85	!
86	! 4121 2019-07-26 10:01:22Z schwenkel
87	! Implementation of an simple method for interpolating the velocities to
88	! particle position
89	!
90	! 4114 2019-07-23 14:09:27Z schwenkel
91	! Bugfix: Added working precision for if statement
92	!
93	! 4054 2019-06-27 07:42:18Z raasch
94	! bugfix for calculating the minimum particle time step
95	!
96	! 4044 2019-06-19 12:28:27Z schwenkel
97	! Bugfix in case of grid strecting: corrected calculation of k-Index
98	!
99	! 4043 2019-06-18 16:59:00Z schwenkel
100	! Remove min_nr_particle, Add lpm_droplet_interactions_ptq into module
101	!
102	! 4028 2019-06-13 12:21:37Z schwenkel
103	! Further modularization of particle code components
104	!
105	! 4020 2019-06-06 14:57:48Z schwenkel
106	! Removing submodules
107	!
108	! 4018 2019-06-06 13:41:50Z eckhard
109	! Bugfix for former revision
110	!
111	! 4017 2019-06-06 12:16:46Z schwenkel
112	! Modularization of all lagrangian particle model code components
113	!
114	! 3655 2019-01-07 16:51:22Z knoop
115	! bugfix to guarantee correct particle releases in case that the release
116	! interval is smaller than the model timestep
117	!
118	! Revision 1.1 1999/11/25 16:16:06 raasch
119	! Initial revision
120	!
121	!
122	! Description:
123	! ------------
124	!> The embedded LPM allows for studying transport and dispersion processes within
125	!> turbulent flows. This model including passive particles that do not show any
126	!> feedback on the turbulent flow. Further also particles with inertia and
127	!> cloud droplets ca be simulated explicitly.
128	!>
129	!> @todo test lcm
130	!> implement simple interpolation method for subgrid scale velocites
131	!> @note <Enter notes on the module>
132	!> @bug <Enter bug on the module>
133	!------------------------------------------------------------------------------!
134	MODULE lagrangian_particle_model_mod
135
136	USE, INTRINSIC :: ISO_C_BINDING
137
138	USE arrays_3d, &
139	ONLY: de_dx, de_dy, de_dz, &
140	d_exner, &
141	drho_air_zw, &
142	dzw, zu, zw, ql_c, ql_v, ql_vp, hyp, &
143	pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p
144
145	USE averaging, &
146	ONLY: ql_c_av, pr_av, pc_av, ql_vp_av, ql_v_av
147
148	USE basic_constants_and_equations_mod, &
149	ONLY: molecular_weight_of_solute, molecular_weight_of_water, magnus, &
150	pi, rd_d_rv, rho_l, r_v, rho_s, vanthoff, l_v, kappa, g, lv_d_cp
151
152	USE control_parameters, &
153	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
154	child_domain, &
155	cloud_droplets, constant_flux_layer, current_timestep_number, &
156	dt_3d, dt_3d_reached, first_call_lpm, humidity, &
157	dt_3d_reached_l, dt_dopts, dz, initializing_actions, &
158	intermediate_timestep_count, intermediate_timestep_count_max, &
159	message_string, molecular_viscosity, ocean_mode, &
160	particle_maximum_age, iran, &
161	simulated_time, topography, dopts_time_count, &
162	time_since_reference_point, rho_surface, u_gtrans, v_gtrans, &
163	dz_stretch_level, dz_stretch_level_start
164
165	USE cpulog, &
166	ONLY: cpu_log, log_point, log_point_s
167
168	USE indices, &
169	ONLY: nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb, &
170	nzb_max, nzt,nbgp, ngp_2dh_outer, &
171	topo_top_ind, &
172	wall_flags_total_0
173
174	USE kinds
175
176	USE pegrid
177
178	USE particle_attributes
179
180	#if defined( __parallel )
181	USE pmc_particle_interface, &
182	ONLY: pmcp_c_get_particle_from_parent, pmcp_p_fill_particle_win, &
183	pmcp_c_send_particle_to_parent, pmcp_p_empty_particle_win, &
184	pmcp_p_delete_particles_in_fine_grid_area, pmcp_g_init, &
185	pmcp_g_print_number_of_particles
186	#endif
187
188	USE pmc_interface, &
189	ONLY: nested_run
190
191	USE grid_variables, &
192	ONLY: ddx, dx, ddy, dy
193
194	USE netcdf_interface, &
195	ONLY: netcdf_data_format, netcdf_deflate, dopts_num, id_set_pts, &
196	id_var_dopts, id_var_time_pts, nc_stat, &
197	netcdf_handle_error
198
199	USE random_function_mod, &
200	ONLY: random_function
201
202	USE statistics, &
203	ONLY: hom
204
205	USE surface_mod, &
206	ONLY: bc_h, &
207	surf_def_h, &
208	surf_lsm_h, &
209	surf_usm_h
210
211	#if defined( __parallel ) && !defined( __mpifh )
212	USE MPI
213	#endif
214
215	#if defined( __netcdf )
216	USE NETCDF
217	#endif
218
219	IMPLICIT NONE
220
221	CHARACTER(LEN=15) :: aero_species = 'nacl' !< aerosol species
222	CHARACTER(LEN=15) :: aero_type = 'maritime' !< aerosol type
223	CHARACTER(LEN=15) :: bc_par_lr = 'cyclic' !< left/right boundary condition
224	CHARACTER(LEN=15) :: bc_par_ns = 'cyclic' !< north/south boundary condition
225	CHARACTER(LEN=15) :: bc_par_b = 'reflect' !< bottom boundary condition
226	CHARACTER(LEN=15) :: bc_par_t = 'absorb' !< top boundary condition
227	CHARACTER(LEN=15) :: collision_kernel = 'none' !< collision kernel
228
229	CHARACTER(LEN=5) :: splitting_function = 'gamma' !< function for calculation critical weighting factor
230	CHARACTER(LEN=5) :: splitting_mode = 'const' !< splitting mode
231
232	CHARACTER(LEN=25) :: particle_advection_interpolation = 'trilinear' !< interpolation method for calculatin the particle
233
234	INTEGER(iwp) :: deleted_particles = 0 !< number of deleted particles per time step
235	INTEGER(iwp) :: i_splitting_mode !< dummy for splitting mode
236	INTEGER(iwp) :: iran_part = -1234567 !< number for random generator
237	INTEGER(iwp) :: max_number_particles_per_gridbox = 100 !< namelist parameter (see documentation)
238	INTEGER(iwp) :: isf !< dummy for splitting function
239	INTEGER(iwp) :: number_particles_per_gridbox = -1 !< namelist parameter (see documentation)
240	INTEGER(iwp) :: number_of_sublayers = 20 !< number of sublayers for particle velocities betwenn surface and first grid level
241	INTEGER(iwp) :: offset_ocean_nzt = 0 !< in case of oceans runs, the vertical index calculations need an offset
242	INTEGER(iwp) :: offset_ocean_nzt_m1 = 0 !< in case of oceans runs, the vertical index calculations need an offset
243	INTEGER(iwp) :: particles_per_point = 1 !< namelist parameter (see documentation)
244	INTEGER(iwp) :: radius_classes = 20 !< namelist parameter (see documentation)
245	INTEGER(iwp) :: splitting_factor = 2 !< namelist parameter (see documentation)
246	INTEGER(iwp) :: splitting_factor_max = 5 !< namelist parameter (see documentation)
247	INTEGER(iwp) :: step_dealloc = 100 !< namelist parameter (see documentation)
248	INTEGER(iwp) :: total_number_of_particles !< total number of particles in the whole model domain
249	INTEGER(iwp) :: trlp_count_sum !< parameter for particle exchange of PEs
250	INTEGER(iwp) :: trlp_count_recv_sum !< parameter for particle exchange of PEs
251	INTEGER(iwp) :: trrp_count_sum !< parameter for particle exchange of PEs
252	INTEGER(iwp) :: trrp_count_recv_sum !< parameter for particle exchange of PEs
253	INTEGER(iwp) :: trsp_count_sum !< parameter for particle exchange of PEs
254	INTEGER(iwp) :: trsp_count_recv_sum !< parameter for particle exchange of PEs
255	INTEGER(iwp) :: trnp_count_sum !< parameter for particle exchange of PEs
256	INTEGER(iwp) :: trnp_count_recv_sum !< parameter for particle exchange of PEs
257
258	LOGICAL :: lagrangian_particle_model = .FALSE. !< namelist parameter (see documentation)
259	LOGICAL :: curvature_solution_effects = .FALSE. !< namelist parameter (see documentation)
260	LOGICAL :: deallocate_memory = .TRUE. !< namelist parameter (see documentation)
261	LOGICAL :: hall_kernel = .FALSE. !< flag for collision kernel
262	LOGICAL :: merging = .FALSE. !< namelist parameter (see documentation)
263	LOGICAL :: random_start_position = .FALSE. !< namelist parameter (see documentation)
264	LOGICAL :: read_particles_from_restartfile = .TRUE. !< namelist parameter (see documentation)
265	LOGICAL :: seed_follows_topography = .FALSE. !< namelist parameter (see documentation)
266	LOGICAL :: splitting = .FALSE. !< namelist parameter (see documentation)
267	LOGICAL :: use_kernel_tables = .FALSE. !< parameter, which turns on the use of precalculated collision kernels
268	LOGICAL :: write_particle_statistics = .FALSE. !< namelist parameter (see documentation)
269	LOGICAL :: interpolation_simple_predictor = .FALSE. !< flag for simple particle advection interpolation with predictor step
270	LOGICAL :: interpolation_simple_corrector = .FALSE. !< flag for simple particle advection interpolation with corrector step
271	LOGICAL :: interpolation_trilinear = .FALSE. !< flag for trilinear particle advection interpolation
272
273	LOGICAL, DIMENSION(max_number_of_particle_groups) :: vertical_particle_advection = .TRUE. !< Switch for vertical particle transport
274
275	REAL(wp) :: aero_weight = 1.0_wp !< namelist parameter (see documentation)
276	REAL(wp) :: dt_min_part = 0.0002_wp !< minimum particle time step when SGS velocities are used (s)
277	REAL(wp) :: dt_prel = 9999999.9_wp !< namelist parameter (see documentation)
278	REAL(wp) :: dt_write_particle_data = 9999999.9_wp !< namelist parameter (see documentation)
279	REAL(wp) :: end_time_prel = 9999999.9_wp !< namelist parameter (see documentation)
280	REAL(wp) :: initial_weighting_factor = 1.0_wp !< namelist parameter (see documentation)
281	REAL(wp) :: last_particle_release_time = 0.0_wp !< last time of particle release
282	REAL(wp) :: log_sigma(3) = 1.0_wp !< namelist parameter (see documentation)
283	REAL(wp) :: na(3) = 0.0_wp !< namelist parameter (see documentation)
284	REAL(wp) :: number_concentration = -1.0_wp !< namelist parameter (see documentation)
285	REAL(wp) :: radius_merge = 1.0E-7_wp !< namelist parameter (see documentation)
286	REAL(wp) :: radius_split = 40.0E-6_wp !< namelist parameter (see documentation)
287	REAL(wp) :: rm(3) = 1.0E-6_wp !< namelist parameter (see documentation)
288	REAL(wp) :: sgs_wf_part !< parameter for sgs
289	REAL(wp) :: time_write_particle_data = 0.0_wp !< write particle data at current time on file
290	REAL(wp) :: weight_factor_merge = -1.0_wp !< namelist parameter (see documentation)
291	REAL(wp) :: weight_factor_split = -1.0_wp !< namelist parameter (see documentation)
292	REAL(wp) :: z0_av_global !< horizontal mean value of z0
293
294	REAL(wp) :: rclass_lbound !<
295	REAL(wp) :: rclass_ubound !<
296
297	REAL(wp), PARAMETER :: c_0 = 3.0_wp !< parameter for lagrangian timescale
298
299	REAL(wp), DIMENSION(max_number_of_particle_groups) :: density_ratio = 9999999.9_wp !< namelist parameter (see documentation)
300	REAL(wp), DIMENSION(max_number_of_particle_groups) :: pdx = 9999999.9_wp !< namelist parameter (see documentation)
301	REAL(wp), DIMENSION(max_number_of_particle_groups) :: pdy = 9999999.9_wp !< namelist parameter (see documentation)
302	REAL(wp), DIMENSION(max_number_of_particle_groups) :: pdz = 9999999.9_wp !< namelist parameter (see documentation)
303	REAL(wp), DIMENSION(max_number_of_particle_groups) :: psb = 9999999.9_wp !< namelist parameter (see documentation)
304	REAL(wp), DIMENSION(max_number_of_particle_groups) :: psl = 9999999.9_wp !< namelist parameter (see documentation)
305	REAL(wp), DIMENSION(max_number_of_particle_groups) :: psn = 9999999.9_wp !< namelist parameter (see documentation)
306	REAL(wp), DIMENSION(max_number_of_particle_groups) :: psr = 9999999.9_wp !< namelist parameter (see documentation)
307	REAL(wp), DIMENSION(max_number_of_particle_groups) :: pss = 9999999.9_wp !< namelist parameter (see documentation)
308	REAL(wp), DIMENSION(max_number_of_particle_groups) :: pst = 9999999.9_wp !< namelist parameter (see documentation).
309	REAL(wp), DIMENSION(max_number_of_particle_groups) :: radius = 9999999.9_wp !< namelist parameter (see documentation)
310
311	REAL(wp), DIMENSION(:), ALLOCATABLE :: log_z_z0 !< Precalculate LOG(z/z0)
312
313	INTEGER(iwp), PARAMETER :: NR_2_direction_move = 10000 !<
314
315	#if defined( __parallel )
316	INTEGER(iwp) :: nr_move_north !<
317	INTEGER(iwp) :: nr_move_south !<
318
319	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: move_also_north
320	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: move_also_south
321	#endif
322
323	REAL(wp) :: epsilon_collision !<
324	REAL(wp) :: urms !<
325
326	REAL(wp), DIMENSION(:), ALLOCATABLE :: epsclass !< dissipation rate class
327	REAL(wp), DIMENSION(:), ALLOCATABLE :: radclass !< radius class
328	REAL(wp), DIMENSION(:), ALLOCATABLE :: winf !<
329
330	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ec !<
331	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ecf !<
332	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gck !<
333	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hkernel !<
334	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hwratio !<
335
336	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ckernel !<
337	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_t !< u value of old timelevel t
338	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: v_t !< v value of old timelevel t
339	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: w_t !< w value of old timelevel t
340
341
342	INTEGER(iwp), PARAMETER :: PHASE_INIT = 1 !<
343	INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2 !<
344
345	SAVE
346
347	PRIVATE
348
349	PUBLIC lpm_parin, &
350	lpm_header, &
351	lpm_init_arrays,&
352	lpm_init, &
353	lpm_actions, &
354	lpm_data_output_ptseries, &
355	lpm_interaction_droplets_ptq, &
356	lpm_rrd_local_particles, &
357	lpm_wrd_local, &
358	lpm_rrd_global, &
359	lpm_wrd_global, &
360	lpm_rrd_local, &
361	lpm_check_parameters
362
363	PUBLIC lagrangian_particle_model
364
365	INTERFACE lpm_check_parameters
366	MODULE PROCEDURE lpm_check_parameters
367	END INTERFACE lpm_check_parameters
368
369	INTERFACE lpm_parin
370	MODULE PROCEDURE lpm_parin
371	END INTERFACE lpm_parin
372
373	INTERFACE lpm_header
374	MODULE PROCEDURE lpm_header
375	END INTERFACE lpm_header
376
377	INTERFACE lpm_init_arrays
378	MODULE PROCEDURE lpm_init_arrays
379	END INTERFACE lpm_init_arrays
380
381	INTERFACE lpm_init
382	MODULE PROCEDURE lpm_init
383	END INTERFACE lpm_init
384
385	INTERFACE lpm_actions
386	MODULE PROCEDURE lpm_actions
387	END INTERFACE lpm_actions
388
389	INTERFACE lpm_data_output_ptseries
390	MODULE PROCEDURE lpm_data_output_ptseries
391	END INTERFACE
392
393	INTERFACE lpm_rrd_local_particles
394	MODULE PROCEDURE lpm_rrd_local_particles
395	END INTERFACE lpm_rrd_local_particles
396
397	INTERFACE lpm_rrd_global
398	MODULE PROCEDURE lpm_rrd_global
399	END INTERFACE lpm_rrd_global
400
401	INTERFACE lpm_rrd_local
402	MODULE PROCEDURE lpm_rrd_local
403	END INTERFACE lpm_rrd_local
404
405	INTERFACE lpm_wrd_local
406	MODULE PROCEDURE lpm_wrd_local
407	END INTERFACE lpm_wrd_local
408
409	INTERFACE lpm_wrd_global
410	MODULE PROCEDURE lpm_wrd_global
411	END INTERFACE lpm_wrd_global
412
413	INTERFACE lpm_advec
414	MODULE PROCEDURE lpm_advec
415	END INTERFACE lpm_advec
416
417	INTERFACE lpm_calc_liquid_water_content
418	MODULE PROCEDURE lpm_calc_liquid_water_content
419	END INTERFACE
420
421	INTERFACE lpm_interaction_droplets_ptq
422	MODULE PROCEDURE lpm_interaction_droplets_ptq
423	MODULE PROCEDURE lpm_interaction_droplets_ptq_ij
424	END INTERFACE lpm_interaction_droplets_ptq
425
426	INTERFACE lpm_boundary_conds
427	MODULE PROCEDURE lpm_boundary_conds
428	END INTERFACE lpm_boundary_conds
429
430	INTERFACE lpm_droplet_condensation
431	MODULE PROCEDURE lpm_droplet_condensation
432	END INTERFACE
433
434	INTERFACE lpm_droplet_collision
435	MODULE PROCEDURE lpm_droplet_collision
436	END INTERFACE lpm_droplet_collision
437
438	INTERFACE lpm_init_kernels
439	MODULE PROCEDURE lpm_init_kernels
440	END INTERFACE lpm_init_kernels
441
442	INTERFACE lpm_splitting
443	MODULE PROCEDURE lpm_splitting
444	END INTERFACE lpm_splitting
445
446	INTERFACE lpm_merging
447	MODULE PROCEDURE lpm_merging
448	END INTERFACE lpm_merging
449
450	INTERFACE lpm_exchange_horiz
451	MODULE PROCEDURE lpm_exchange_horiz
452	END INTERFACE lpm_exchange_horiz
453
454	INTERFACE lpm_move_particle
455	MODULE PROCEDURE lpm_move_particle
456	END INTERFACE lpm_move_particle
457
458	INTERFACE realloc_particles_array
459	MODULE PROCEDURE realloc_particles_array
460	END INTERFACE realloc_particles_array
461
462	INTERFACE dealloc_particles_array
463	MODULE PROCEDURE dealloc_particles_array
464	END INTERFACE dealloc_particles_array
465
466	INTERFACE lpm_sort_and_delete
467	MODULE PROCEDURE lpm_sort_and_delete
468	END INTERFACE lpm_sort_and_delete
469
470	INTERFACE lpm_sort_timeloop_done
471	MODULE PROCEDURE lpm_sort_timeloop_done
472	END INTERFACE lpm_sort_timeloop_done
473
474	INTERFACE lpm_pack
475	MODULE PROCEDURE lpm_pack
476	END INTERFACE lpm_pack
477
478	CONTAINS
479
480
481	!------------------------------------------------------------------------------!
482	! Description:
483	! ------------
484	!> Parin for &particle_parameters for the Lagrangian particle model
485	!------------------------------------------------------------------------------!
486	SUBROUTINE lpm_parin
487
488	CHARACTER (LEN=80) :: line !<
489
490	NAMELIST /particles_par/ &
491	aero_species, &
492	aero_type, &
493	aero_weight, &
494	alloc_factor, &
495	bc_par_b, &
496	bc_par_lr, &
497	bc_par_ns, &
498	bc_par_t, &
499	collision_kernel, &
500	curvature_solution_effects, &
501	deallocate_memory, &
502	density_ratio, &
503	dissipation_classes, &
504	dt_dopts, &
505	dt_min_part, &
506	dt_prel, &
507	dt_write_particle_data, &
508	end_time_prel, &
509	initial_weighting_factor, &
510	log_sigma, &
511	max_number_particles_per_gridbox, &
512	merging, &
513	na, &
514	number_concentration, &
515	number_of_particle_groups, &
516	number_particles_per_gridbox, &
517	particles_per_point, &
518	particle_advection_start, &
519	particle_advection_interpolation, &
520	particle_maximum_age, &
521	pdx, &
522	pdy, &
523	pdz, &
524	psb, &
525	psl, &
526	psn, &
527	psr, &
528	pss, &
529	pst, &
530	radius, &
531	radius_classes, &
532	radius_merge, &
533	radius_split, &
534	random_start_position, &
535	read_particles_from_restartfile, &
536	rm, &
537	seed_follows_topography, &
538	splitting, &
539	splitting_factor, &
540	splitting_factor_max, &
541	splitting_function, &
542	splitting_mode, &
543	step_dealloc, &
544	use_sgs_for_particles, &
545	vertical_particle_advection, &
546	weight_factor_merge, &
547	weight_factor_split, &
548	write_particle_statistics
549
550	NAMELIST /particle_parameters/ &
551	aero_species, &
552	aero_type, &
553	aero_weight, &
554	alloc_factor, &
555	bc_par_b, &
556	bc_par_lr, &
557	bc_par_ns, &
558	bc_par_t, &
559	collision_kernel, &
560	curvature_solution_effects, &
561	deallocate_memory, &
562	density_ratio, &
563	dissipation_classes, &
564	dt_dopts, &
565	dt_min_part, &
566	dt_prel, &
567	dt_write_particle_data, &
568	end_time_prel, &
569	initial_weighting_factor, &
570	log_sigma, &
571	max_number_particles_per_gridbox, &
572	merging, &
573	na, &
574	number_concentration, &
575	number_of_particle_groups, &
576	number_particles_per_gridbox, &
577	particles_per_point, &
578	particle_advection_start, &
579	particle_advection_interpolation, &
580	particle_maximum_age, &
581	pdx, &
582	pdy, &
583	pdz, &
584	psb, &
585	psl, &
586	psn, &
587	psr, &
588	pss, &
589	pst, &
590	radius, &
591	radius_classes, &
592	radius_merge, &
593	radius_split, &
594	random_start_position, &
595	read_particles_from_restartfile, &
596	rm, &
597	seed_follows_topography, &
598	splitting, &
599	splitting_factor, &
600	splitting_factor_max, &
601	splitting_function, &
602	splitting_mode, &
603	step_dealloc, &
604	use_sgs_for_particles, &
605	vertical_particle_advection, &
606	weight_factor_merge, &
607	weight_factor_split, &
608	write_particle_statistics
609
610	!
611	!-- Position the namelist-file at the beginning (it was already opened in
612	!-- parin), search for the namelist-group of the package and position the
613	!-- file at this line. Do the same for each optionally used package.
614	line = ' '
615
616	!
617	!-- Try to find particles package
618	REWIND ( 11 )
619	line = ' '
620	DO WHILE ( INDEX( line, '&particle_parameters' ) == 0 )
621	READ ( 11, '(A)', END=12 ) line
622	ENDDO
623	BACKSPACE ( 11 )
624	!
625	!-- Read user-defined namelist
626	READ ( 11, particle_parameters, ERR = 10 )
627	!
628	!-- Set flag that indicates that particles are switched on
629	particle_advection = .TRUE.
630
631	GOTO 14
632
633	10 BACKSPACE( 11 )
634	READ( 11 , '(A)') line
635	CALL parin_fail_message( 'particle_parameters', line )
636	!
637	!-- Try to find particles package (old namelist)
638	12 REWIND ( 11 )
639	line = ' '
640	DO WHILE ( INDEX( line, '&particles_par' ) == 0 )
641	READ ( 11, '(A)', END=14 ) line
642	ENDDO
643	BACKSPACE ( 11 )
644	!
645	!-- Read user-defined namelist
646	READ ( 11, particles_par, ERR = 13, END = 14 )
647
648	message_string = 'namelist particles_par is deprecated and will be ' // &
649	'removed in near future. Please use namelist ' // &
650	'particle_parameters instead'
651	CALL message( 'package_parin', 'PA0487', 0, 1, 0, 6, 0 )
652
653	!
654	!-- Set flag that indicates that particles are switched on
655	particle_advection = .TRUE.
656
657	GOTO 14
658
659	13 BACKSPACE( 11 )
660	READ( 11 , '(A)') line
661	CALL parin_fail_message( 'particles_par', line )
662
663	14 CONTINUE
664
665	END SUBROUTINE lpm_parin
666
667	!------------------------------------------------------------------------------!
668	! Description:
669	! ------------
670	!> Writes used particle attributes in header file.
671	!------------------------------------------------------------------------------!
672	SUBROUTINE lpm_header ( io )
673
674	CHARACTER (LEN=40) :: output_format !< netcdf format
675
676	INTEGER(iwp) :: i !<
677	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
678
679
680	IF ( humidity .AND. cloud_droplets ) THEN
681	WRITE ( io, 433 )
682	IF ( curvature_solution_effects ) WRITE ( io, 434 )
683	IF ( collision_kernel /= 'none' ) THEN
684	WRITE ( io, 435 ) TRIM( collision_kernel )
685	IF ( collision_kernel(6:9) == 'fast' ) THEN
686	WRITE ( io, 436 ) radius_classes, dissipation_classes
687	ENDIF
688	ELSE
689	WRITE ( io, 437 )
690	ENDIF
691	ENDIF
692
693	IF ( particle_advection ) THEN
694	!
695	!-- Particle attributes
696	WRITE ( io, 480 ) particle_advection_start, TRIM(particle_advection_interpolation), &
697	dt_prel, bc_par_lr, &
698	bc_par_ns, bc_par_b, bc_par_t, particle_maximum_age, &
699	end_time_prel
700	IF ( use_sgs_for_particles ) WRITE ( io, 488 ) dt_min_part
701	IF ( random_start_position ) WRITE ( io, 481 )
702	IF ( seed_follows_topography ) WRITE ( io, 496 )
703	IF ( particles_per_point > 1 ) WRITE ( io, 489 ) particles_per_point
704	WRITE ( io, 495 ) total_number_of_particles
705	IF ( dt_write_particle_data /= 9999999.9_wp ) THEN
706	WRITE ( io, 485 ) dt_write_particle_data
707	IF ( netcdf_data_format > 1 ) THEN
708	output_format = 'netcdf (64 bit offset) and binary'
709	ELSE
710	output_format = 'netcdf and binary'
711	ENDIF
712	IF ( netcdf_deflate == 0 ) THEN
713	WRITE ( io, 344 ) output_format
714	ELSE
715	WRITE ( io, 354 ) TRIM( output_format ), netcdf_deflate
716	ENDIF
717	ENDIF
718	IF ( dt_dopts /= 9999999.9_wp ) WRITE ( io, 494 ) dt_dopts
719	IF ( write_particle_statistics ) WRITE ( io, 486 )
720
721	WRITE ( io, 487 ) number_of_particle_groups
722
723	DO i = 1, number_of_particle_groups
724	WRITE ( io, 490 ) i, radius(i)
725	IF ( density_ratio(i) /= 0.0_wp ) THEN
726	WRITE ( io, 491 ) density_ratio(i)
727	ELSE
728	WRITE ( io, 492 )
729	ENDIF
730	WRITE ( io, 493 ) psl(i), psr(i), pss(i), psn(i), psb(i), pst(i), &
731	pdx(i), pdy(i), pdz(i)
732	IF ( .NOT. vertical_particle_advection(i) ) WRITE ( io, 482 )
733	ENDDO
734
735	ENDIF
736
737	344 FORMAT (' Output format: ',A/)
738	354 FORMAT (' Output format: ',A, ' compressed with level: ',I1/)
739
740	433 FORMAT (' Cloud droplets treated explicitly using the Lagrangian part', &
741	'icle model')
742	434 FORMAT (' Curvature and solution effecs are considered for growth of', &
743	' droplets < 1.0E-6 m')
744	435 FORMAT (' Droplet collision is handled by ',A,'-kernel')
745	436 FORMAT (' Fast kernel with fixed radius- and dissipation classes ', &
746	'are used'/ &
747	' number of radius classes: ',I3,' interval ', &
748	'[1.0E-6,2.0E-4] m'/ &
749	' number of dissipation classes: ',I2,' interval ', &
750	'[0,1000] cm2/s3')
751	437 FORMAT (' Droplet collision is switched off')
752
753	480 FORMAT (' Particles:'/ &
754	' ---------'// &
755	' Particle advection is active (switched on at t = ', F7.1, &
756	' s)'/ &
757	' Interpolation of particle velocities is done by using ', A, &
758	' method'/ &
759	' Start of new particle generations every ',F6.1,' s'/ &
760	' Boundary conditions: left/right: ', A, ' north/south: ', A/&
761	' bottom: ', A, ' top: ', A/&
762	' Maximum particle age: ',F9.1,' s'/ &
763	' Advection stopped at t = ',F9.1,' s'/)
764	481 FORMAT (' Particles have random start positions'/)
765	482 FORMAT (' Particles are advected only horizontally'/)
766	485 FORMAT (' Particle data are written on file every ', F9.1, ' s')
767	486 FORMAT (' Particle statistics are written on file'/)
768	487 FORMAT (' Number of particle groups: ',I2/)
769	488 FORMAT (' SGS velocity components are used for particle advection'/ &
770	' minimum timestep for advection:', F8.5/)
771	489 FORMAT (' Number of particles simultaneously released at each ', &
772	'point: ', I5/)
773	490 FORMAT (' Particle group ',I2,':'/ &
774	' Particle radius: ',E10.3, 'm')
775	491 FORMAT (' Particle inertia is activated'/ &
776	' density_ratio (rho_fluid/rho_particle) =',F6.3/)
777	492 FORMAT (' Particles are advected only passively (no inertia)'/)
778	493 FORMAT (' Boundaries of particle source: x:',F8.1,' - ',F8.1,' m'/&
779	' y:',F8.1,' - ',F8.1,' m'/&
780	' z:',F8.1,' - ',F8.1,' m'/&
781	' Particle distances: dx = ',F8.1,' m dy = ',F8.1, &
782	' m dz = ',F8.1,' m'/)
783	494 FORMAT (' Output of particle time series in NetCDF format every ', &
784	F8.2,' s'/)
785	495 FORMAT (' Number of particles in total domain: ',I10/)
786	496 FORMAT (' Initial vertical particle positions are interpreted ', &
787	'as relative to the given topography')
788
789	END SUBROUTINE lpm_header
790
791	!------------------------------------------------------------------------------!
792	! Description:
793	! ------------
794	!> Writes used particle attributes in header file.
795	!------------------------------------------------------------------------------!
796	SUBROUTINE lpm_check_parameters
797
798	!
799	!-- Collision kernels:
800	SELECT CASE ( TRIM( collision_kernel ) )
801
802	CASE ( 'hall', 'hall_fast' )
803	hall_kernel = .TRUE.
804
805	CASE ( 'wang', 'wang_fast' )
806	wang_kernel = .TRUE.
807
808	CASE ( 'none' )
809
810
811	CASE DEFAULT
812	message_string = 'unknown collision kernel: collision_kernel = "' // &
813	TRIM( collision_kernel ) // '"'
814	CALL message( 'lpm_check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
815
816	END SELECT
817	IF ( collision_kernel(6:9) == 'fast' ) use_kernel_tables = .TRUE.
818
819	!
820	!-- Subgrid scale velocites with the simple interpolation method for resolved
821	!-- velocites is not implemented for passive particles. However, for cloud
822	!-- it can be combined as the sgs-velocites for active particles are
823	!-- calculated differently, i.e. no subboxes are needed.
824	IF ( .NOT. TRIM( particle_advection_interpolation ) == 'trilinear' .AND. &
825	use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN
826	message_string = 'subrgrid scale velocities in combination with ' // &
827	'simple interpolation method is not ' // &
828	'implemented'
829	CALL message( 'lpm_check_parameters', 'PA0659', 1, 2, 0, 6, 0 )
830	ENDIF
831
832	IF ( nested_run .AND. cloud_droplets ) THEN
833	message_string = 'nested runs in combination with cloud droplets ' // &
834	'is not implemented'
835	CALL message( 'lpm_check_parameters', 'PA0687', 1, 2, 0, 6, 0 )
836	ENDIF
837
838
839	END SUBROUTINE lpm_check_parameters
840
841	!------------------------------------------------------------------------------!
842	! Description:
843	! ------------
844	!> Initialize arrays for lpm
845	!------------------------------------------------------------------------------!
846	SUBROUTINE lpm_init_arrays
847
848	IF ( cloud_droplets ) THEN
849	!
850	!-- Liquid water content, change in liquid water content
851	ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
852	ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
853	!-- Real volume of particles (with weighting), volume of particles
854	ALLOCATE ( ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
855	ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
856	ENDIF
857
858
859	ALLOCATE( u_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
860	v_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
861	w_t(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
862	!
863	!-- Initialize values with current time step
864	u_t = u
865	v_t = v
866	w_t = w
867	!
868	!-- Initial assignment of the pointers
869	IF ( cloud_droplets ) THEN
870	ql => ql_1
871	ql_c => ql_2
872	ENDIF
873
874	END SUBROUTINE lpm_init_arrays
875
876	!------------------------------------------------------------------------------!
877	! Description:
878	! ------------
879	!> Initialize Lagrangian particle model
880	!------------------------------------------------------------------------------!
881	SUBROUTINE lpm_init
882
883	INTEGER(iwp) :: i !<
884	INTEGER(iwp) :: j !<
885	INTEGER(iwp) :: k !<
886
887	REAL(wp) :: div !<
888	REAL(wp) :: height_int !<
889	REAL(wp) :: height_p !<
890	REAL(wp) :: z_p !<
891	REAL(wp) :: z0_av_local !<
892
893	!
894	!-- In case of oceans runs, the vertical index calculations need an offset,
895	!-- because otherwise the k indices will become negative
896	IF ( ocean_mode ) THEN
897	offset_ocean_nzt = nzt
898	offset_ocean_nzt_m1 = nzt - 1
899	ENDIF
900
901	!
902	!-- Define block offsets for dividing a gridcell in 8 sub cells
903	!-- See documentation for List of subgrid boxes
904	!-- See pack_and_sort in lpm_pack_arrays.f90 for assignment of the subgrid boxes
905	block_offset(0) = block_offset_def ( 0, 0, 0)
906	block_offset(1) = block_offset_def ( 0, 0,-1)
907	block_offset(2) = block_offset_def ( 0,-1, 0)
908	block_offset(3) = block_offset_def ( 0,-1,-1)
909	block_offset(4) = block_offset_def (-1, 0, 0)
910	block_offset(5) = block_offset_def (-1, 0,-1)
911	block_offset(6) = block_offset_def (-1,-1, 0)
912	block_offset(7) = block_offset_def (-1,-1,-1)
913	!
914	!-- Check the number of particle groups.
915	IF ( number_of_particle_groups > max_number_of_particle_groups ) THEN
916	WRITE( message_string, * ) 'max_number_of_particle_groups =', &
917	max_number_of_particle_groups , &
918	'&number_of_particle_groups reset to ', &
919	max_number_of_particle_groups
920	CALL message( 'lpm_init', 'PA0213', 0, 1, 0, 6, 0 )
921	number_of_particle_groups = max_number_of_particle_groups
922	ENDIF
923	!
924	!-- Check if downward-facing walls exist. This case, reflection boundary
925	!-- conditions (as well as subgrid-scale velocities) may do not work
926	!-- propably (not realized so far).
927	IF ( surf_def_h(1)%ns >= 1 ) THEN
928	WRITE( message_string, * ) 'Overhanging topography do not work '// &
929	'with particles'
930	CALL message( 'lpm_init', 'PA0212', 0, 1, 0, 6, 0 )
931
932	ENDIF
933
934	!
935	!-- Set default start positions, if necessary
936	IF ( psl(1) == 9999999.9_wp ) psl(1) = 0.0_wp
937	IF ( psr(1) == 9999999.9_wp ) psr(1) = ( nx +1 ) * dx
938	IF ( pss(1) == 9999999.9_wp ) pss(1) = 0.0_wp
939	IF ( psn(1) == 9999999.9_wp ) psn(1) = ( ny +1 ) * dy
940	IF ( psb(1) == 9999999.9_wp ) psb(1) = zu(nz/2)
941	IF ( pst(1) == 9999999.9_wp ) pst(1) = psb(1)
942
943	IF ( pdx(1) == 9999999.9_wp .OR. pdx(1) == 0.0_wp ) pdx(1) = dx
944	IF ( pdy(1) == 9999999.9_wp .OR. pdy(1) == 0.0_wp ) pdy(1) = dy
945	IF ( pdz(1) == 9999999.9_wp .OR. pdz(1) == 0.0_wp ) pdz(1) = zu(2) - zu(1)
946
947	!
948	!-- If number_particles_per_gridbox is set, the parametres pdx, pdy and pdz are
949	!-- calculated diagnostically. Therfore an isotropic distribution is prescribed.
950	IF ( number_particles_per_gridbox /= -1 .AND. &
951	number_particles_per_gridbox >= 1 ) THEN
952	pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) / &
953	REAL(number_particles_per_gridbox))**0.3333333_wp
954	!
955	!-- Ensure a smooth value (two significant digits) of distance between
956	!-- particles (pdx, pdy, pdz).
957	div = 1000.0_wp
958	DO WHILE ( pdx(1) < div )
959	div = div / 10.0_wp
960	ENDDO
961	pdx(1) = NINT( pdx(1) * 100.0_wp / div ) * div / 100.0_wp
962	pdy(1) = pdx(1)
963	pdz(1) = pdx(1)
964
965	ENDIF
966
967	DO j = 2, number_of_particle_groups
968	IF ( psl(j) == 9999999.9_wp ) psl(j) = psl(j-1)
969	IF ( psr(j) == 9999999.9_wp ) psr(j) = psr(j-1)
970	IF ( pss(j) == 9999999.9_wp ) pss(j) = pss(j-1)
971	IF ( psn(j) == 9999999.9_wp ) psn(j) = psn(j-1)
972	IF ( psb(j) == 9999999.9_wp ) psb(j) = psb(j-1)
973	IF ( pst(j) == 9999999.9_wp ) pst(j) = pst(j-1)
974	IF ( pdx(j) == 9999999.9_wp .OR. pdx(j) == 0.0_wp ) pdx(j) = pdx(j-1)
975	IF ( pdy(j) == 9999999.9_wp .OR. pdy(j) == 0.0_wp ) pdy(j) = pdy(j-1)
976	IF ( pdz(j) == 9999999.9_wp .OR. pdz(j) == 0.0_wp ) pdz(j) = pdz(j-1)
977	ENDDO
978
979	!
980	!-- Allocate arrays required for calculating particle SGS velocities.
981	!-- Initialize prefactor required for stoachastic Weil equation.
982	IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN
983	ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
984	de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
985	de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
986
987	de_dx = 0.0_wp
988	de_dy = 0.0_wp
989	de_dz = 0.0_wp
990
991	sgs_wf_part = 1.0_wp / 3.0_wp
992	ENDIF
993
994	!
995	!-- Allocate array required for logarithmic vertical interpolation of
996	!-- horizontal particle velocities between the surface and the first vertical
997	!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
998	!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
999	!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
1000	!-- To obtain exact height levels of particles, linear interpolation is applied
1001	!-- (see lpm_advec.f90).
1002	IF ( constant_flux_layer ) THEN
1003
1004	ALLOCATE ( log_z_z0(0:number_of_sublayers) )
1005	z_p = zu(nzb+1) - zw(nzb)
1006
1007	!
1008	!-- Calculate horizontal mean value of z0 used for logartihmic
1009	!-- interpolation. Note: this is not exact for heterogeneous z0.
1010	!-- However, sensitivity studies showed that the effect is
1011	!-- negligible.
1012	z0_av_local = SUM( surf_def_h(0)%z0 ) + SUM( surf_lsm_h%z0 ) + &
1013	SUM( surf_usm_h%z0 )
1014	z0_av_global = 0.0_wp
1015
1016	#if defined( __parallel )
1017	CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, &
1018	comm2d, ierr )
1019	#else
1020	z0_av_global = z0_av_local
1021	#endif
1022
1023	z0_av_global = z0_av_global / ( ( ny + 1 ) * ( nx + 1 ) )
1024	!
1025	!-- Horizontal wind speed is zero below and at z0
1026	log_z_z0(0) = 0.0_wp
1027	!
1028	!-- Calculate vertical depth of the sublayers
1029	height_int = ( z_p - z0_av_global ) / REAL( number_of_sublayers, KIND=wp )
1030	!
1031	!-- Precalculate LOG(z/z0)
1032	height_p = z0_av_global
1033	DO k = 1, number_of_sublayers
1034
1035	height_p = height_p + height_int
1036	log_z_z0(k) = LOG( height_p / z0_av_global )
1037
1038	ENDDO
1039
1040	ENDIF
1041
1042	!
1043	!-- Check which particle interpolation method should be used
1044	IF ( TRIM( particle_advection_interpolation ) == 'trilinear' ) THEN
1045	interpolation_simple_corrector = .FALSE.
1046	interpolation_simple_predictor = .FALSE.
1047	interpolation_trilinear = .TRUE.
1048	ELSEIF ( TRIM( particle_advection_interpolation ) == 'simple_corrector' ) THEN
1049	interpolation_simple_corrector = .TRUE.
1050	interpolation_simple_predictor = .FALSE.
1051	interpolation_trilinear = .FALSE.
1052	ELSEIF ( TRIM( particle_advection_interpolation ) == 'simple_predictor' ) THEN
1053	interpolation_simple_corrector = .FALSE.
1054	interpolation_simple_predictor = .TRUE.
1055	interpolation_trilinear = .FALSE.
1056	ENDIF
1057
1058	!
1059	!-- Check boundary condition and set internal variables
1060	SELECT CASE ( bc_par_b )
1061
1062	CASE ( 'absorb' )
1063	ibc_par_b = 1
1064
1065	CASE ( 'reflect' )
1066	ibc_par_b = 2
1067
1068	CASE ( 'reset' )
1069	ibc_par_b = 3
1070
1071	CASE DEFAULT
1072	WRITE( message_string, * ) 'unknown boundary condition ', &
1073	'bc_par_b = "', TRIM( bc_par_b ), '"'
1074	CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
1075
1076	END SELECT
1077	SELECT CASE ( bc_par_t )
1078
1079	CASE ( 'absorb' )
1080	ibc_par_t = 1
1081
1082	CASE ( 'reflect' )
1083	ibc_par_t = 2
1084
1085	CASE ( 'nested' )
1086	ibc_par_t = 3
1087
1088	CASE DEFAULT
1089	WRITE( message_string, * ) 'unknown boundary condition ', &
1090	'bc_par_t = "', TRIM( bc_par_t ), '"'
1091	CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
1092
1093	END SELECT
1094	SELECT CASE ( bc_par_lr )
1095
1096	CASE ( 'cyclic' )
1097	ibc_par_lr = 0
1098
1099	CASE ( 'absorb' )
1100	ibc_par_lr = 1
1101
1102	CASE ( 'reflect' )
1103	ibc_par_lr = 2
1104
1105	CASE ( 'nested' )
1106	ibc_par_lr = 3
1107
1108	CASE DEFAULT
1109	WRITE( message_string, * ) 'unknown boundary condition ', &
1110	'bc_par_lr = "', TRIM( bc_par_lr ), '"'
1111	CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
1112
1113	END SELECT
1114	SELECT CASE ( bc_par_ns )
1115
1116	CASE ( 'cyclic' )
1117	ibc_par_ns = 0
1118
1119	CASE ( 'absorb' )
1120	ibc_par_ns = 1
1121
1122	CASE ( 'reflect' )
1123	ibc_par_ns = 2
1124
1125	CASE ( 'nested' )
1126	ibc_par_ns = 3
1127
1128	CASE DEFAULT
1129	WRITE( message_string, * ) 'unknown boundary condition ', &
1130	'bc_par_ns = "', TRIM( bc_par_ns ), '"'
1131	CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
1132
1133	END SELECT
1134	SELECT CASE ( splitting_mode )
1135
1136	CASE ( 'const' )
1137	i_splitting_mode = 1
1138
1139	CASE ( 'cl_av' )
1140	i_splitting_mode = 2
1141
1142	CASE ( 'gb_av' )
1143	i_splitting_mode = 3
1144
1145	CASE DEFAULT
1146	WRITE( message_string, * ) 'unknown splitting_mode = "', &
1147	TRIM( splitting_mode ), '"'
1148	CALL message( 'lpm_init', 'PA0146', 1, 2, 0, 6, 0 )
1149
1150	END SELECT
1151	SELECT CASE ( splitting_function )
1152
1153	CASE ( 'gamma' )
1154	isf = 1
1155
1156	CASE ( 'log' )
1157	isf = 2
1158
1159	CASE ( 'exp' )
1160	isf = 3
1161
1162	CASE DEFAULT
1163	WRITE( message_string, * ) 'unknown splitting function = "', &
1164	TRIM( splitting_function ), '"'
1165	CALL message( 'lpm_init', 'PA0147', 1, 2, 0, 6, 0 )
1166
1167	END SELECT
1168	!
1169	!-- Initialize collision kernels
1170	IF ( collision_kernel /= 'none' ) CALL lpm_init_kernels
1171	!
1172	!-- For the first model run of a possible job chain initialize the
1173	!-- particles, otherwise read the particle data from restart file.
1174	IF ( TRIM( initializing_actions ) == 'read_restart_data' &
1175	.AND. read_particles_from_restartfile ) THEN
1176	CALL lpm_rrd_local_particles
1177	ELSE
1178	!
1179	!-- Allocate particle arrays and set attributes of the initial set of
1180	!-- particles, which can be also periodically released at later times.
1181	ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
1182	grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
1183
1184	number_of_particles = 0
1185	prt_count = 0
1186	!
1187	!-- initialize counter for particle IDs
1188	grid_particles%id_counter = 1
1189	!
1190	!-- Initialize all particles with dummy values (otherwise errors may
1191	!-- occur within restart runs). The reason for this is still not clear
1192	!-- and may be presumably caused by errors in the respective user-interface.
1193	zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
1194	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
1195	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
1196	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
1197	0, 0, 0_idp, .FALSE., -1 )
1198
1199	particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
1200	!
1201	!-- Set values for the density ratio and radius for all particle
1202	!-- groups, if necessary
1203	IF ( density_ratio(1) == 9999999.9_wp ) density_ratio(1) = 0.0_wp
1204	IF ( radius(1) == 9999999.9_wp ) radius(1) = 0.0_wp
1205	DO i = 2, number_of_particle_groups
1206	IF ( density_ratio(i) == 9999999.9_wp ) THEN
1207	density_ratio(i) = density_ratio(i-1)
1208	ENDIF
1209	IF ( radius(i) == 9999999.9_wp ) radius(i) = radius(i-1)
1210	ENDDO
1211
1212	DO i = 1, number_of_particle_groups
1213	IF ( density_ratio(i) /= 0.0_wp .AND. radius(i) == 0 ) THEN
1214	WRITE( message_string, * ) 'particle group #', i, ' has a', &
1215	'density ratio /= 0 but radius = 0'
1216	CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 )
1217	ENDIF
1218	particle_groups(i)%density_ratio = density_ratio(i)
1219	particle_groups(i)%radius = radius(i)
1220	ENDDO
1221	!
1222	!-- Set a seed value for the random number generator to be exclusively
1223	!-- used for the particle code. The generated random numbers should be
1224	!-- different on the different PEs.
1225	iran_part = iran_part + myid
1226	!
1227	!-- Create the particle set, and set the initial particles
1228	CALL lpm_create_particle( phase_init )
1229	last_particle_release_time = particle_advection_start
1230	!
1231	!-- User modification of initial particles
1232	CALL user_lpm_init
1233	!
1234	!-- Open file for statistical informations about particle conditions
1235	IF ( write_particle_statistics ) THEN
1236	CALL check_open( 80 )
1237	WRITE ( 80, 8000 ) current_timestep_number, simulated_time, &
1238	number_of_particles
1239	CALL close_file( 80 )
1240	ENDIF
1241
1242	ENDIF
1243
1244	#if defined( __parallel )
1245	IF ( nested_run ) CALL pmcp_g_init
1246	#endif
1247
1248	!
1249	!-- To avoid programm abort, assign particles array to the local version of
1250	!-- first grid cell
1251	number_of_particles = prt_count(nzb+1,nys,nxl)
1252	particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles)
1253	!
1254	!-- Formats
1255	8000 FORMAT (I6,1X,F7.2,4X,I10,71X,I10)
1256
1257	END SUBROUTINE lpm_init
1258
1259	!------------------------------------------------------------------------------!
1260	! Description:
1261	! ------------
1262	!> Create Lagrangian particles
1263	!------------------------------------------------------------------------------!
1264	SUBROUTINE lpm_create_particle (phase)
1265
1266	INTEGER(iwp) :: alloc_size !< relative increase of allocated memory for particles
1267	INTEGER(iwp) :: i !< loop variable ( particle groups )
1268	INTEGER(iwp) :: ip !< index variable along x
1269	INTEGER(iwp) :: j !< loop variable ( particles per point )
1270	INTEGER(iwp) :: jp !< index variable along y
1271	INTEGER(iwp) :: k !< index variable along z
1272	INTEGER(iwp) :: k_surf !< index of surface grid point
1273	INTEGER(iwp) :: kp !< index variable along z
1274	INTEGER(iwp) :: loop_stride !< loop variable for initialization
1275	INTEGER(iwp) :: n !< loop variable ( number of particles )
1276	INTEGER(iwp) :: new_size !< new size of allocated memory for particles
1277
1278	INTEGER(iwp), INTENT(IN) :: phase !< mode of inititialization
1279
1280	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_count !< start address of new particle
1281	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: local_start !< start address of new particle
1282
1283	LOGICAL :: first_stride !< flag for initialization
1284
1285	REAL(wp) :: pos_x !< increment for particle position in x
1286	REAL(wp) :: pos_y !< increment for particle position in y
1287	REAL(wp) :: pos_z !< increment for particle position in z
1288	REAL(wp) :: rand_contr !< dummy argument for random position
1289
1290	TYPE(particle_type),TARGET :: tmp_particle !< temporary particle used for initialization
1291
1292
1293	!
1294	!-- Calculate particle positions and store particle attributes, if
1295	!-- particle is situated on this PE
1296	DO loop_stride = 1, 2
1297	first_stride = (loop_stride == 1)
1298	IF ( first_stride ) THEN
1299	local_count = 0 ! count number of particles
1300	ELSE
1301	local_count = prt_count ! Start address of new particles
1302	ENDIF
1303
1304	!
1305	!-- Calculate initial_weighting_factor diagnostically
1306	IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN
1307	initial_weighting_factor = number_concentration * &
1308	pdx(1) * pdy(1) * pdz(1)
1309	END IF
1310
1311	n = 0
1312	DO i = 1, number_of_particle_groups
1313	pos_z = psb(i)
1314	DO WHILE ( pos_z <= pst(i) )
1315	IF ( pos_z >= zw(0) .AND. pos_z < zw(nzt) ) THEN
1316	pos_y = pss(i)
1317	DO WHILE ( pos_y <= psn(i) )
1318	IF ( pos_y >= nys * dy .AND. &
1319	pos_y < ( nyn + 1 ) * dy ) THEN
1320	pos_x = psl(i)
1321	xloop: DO WHILE ( pos_x <= psr(i) )
1322	IF ( pos_x >= nxl * dx .AND. &
1323	pos_x < ( nxr + 1) * dx ) THEN
1324	DO j = 1, particles_per_point
1325	n = n + 1
1326	tmp_particle%x = pos_x
1327	tmp_particle%y = pos_y
1328	tmp_particle%z = pos_z
1329	tmp_particle%age = 0.0_wp
1330	tmp_particle%age_m = 0.0_wp
1331	tmp_particle%dt_sum = 0.0_wp
1332	tmp_particle%e_m = 0.0_wp
1333	tmp_particle%rvar1 = 0.0_wp
1334	tmp_particle%rvar2 = 0.0_wp
1335	tmp_particle%rvar3 = 0.0_wp
1336	tmp_particle%speed_x = 0.0_wp
1337	tmp_particle%speed_y = 0.0_wp
1338	tmp_particle%speed_z = 0.0_wp
1339	tmp_particle%origin_x = pos_x
1340	tmp_particle%origin_y = pos_y
1341	tmp_particle%origin_z = pos_z
1342	IF ( curvature_solution_effects ) THEN
1343	tmp_particle%aux1 = 0.0_wp ! dry aerosol radius
1344	tmp_particle%aux2 = dt_3d ! last Rosenbrock timestep
1345	ELSE
1346	tmp_particle%aux1 = 0.0_wp ! free to use
1347	tmp_particle%aux2 = 0.0_wp ! free to use
1348	ENDIF
1349	tmp_particle%radius = particle_groups(i)%radius
1350	tmp_particle%weight_factor = initial_weighting_factor
1351	tmp_particle%class = 1
1352	tmp_particle%group = i
1353	tmp_particle%id = 0_idp
1354	tmp_particle%particle_mask = .TRUE.
1355	tmp_particle%block_nr = -1
1356	!
1357	!-- Determine the grid indices of the particle position
1358	ip = INT( tmp_particle%x * ddx )
1359	jp = INT( tmp_particle%y * ddy )
1360	!
1361	!-- In case of stretching the actual k index is found iteratively
1362	IF ( dz_stretch_level /= -9999999.9_wp .OR. &
1363	dz_stretch_level_start(1) /= -9999999.9_wp ) THEN
1364	kp = MINLOC( ABS( tmp_particle%z - zu ), DIM = 1 ) - 1
1365	ELSE
1366	kp = INT( tmp_particle%z / dz(1) + 1 + offset_ocean_nzt )
1367	ENDIF
1368	!
1369	!-- Determine surface level. Therefore, check for
1370	!-- upward-facing wall on w-grid.
1371	k_surf = topo_top_ind(jp,ip,3)
1372	IF ( seed_follows_topography ) THEN
1373	!
1374	!-- Particle height is given relative to topography
1375	kp = kp + k_surf
1376	tmp_particle%z = tmp_particle%z + zw(k_surf)
1377	!-- Skip particle release if particle position is
1378	!-- above model top, or within topography in case
1379	!-- of overhanging structures.
1380	IF ( kp > nzt .OR. &
1381	.NOT. BTEST( wall_flags_total_0(kp,jp,ip), 0 ) ) THEN
1382	pos_x = pos_x + pdx(i)
1383	CYCLE xloop
1384	ENDIF
1385	!
1386	!-- Skip particle release if particle position is
1387	!-- below surface, or within topography in case
1388	!-- of overhanging structures.
1389	ELSEIF ( .NOT. seed_follows_topography .AND. &
1390	tmp_particle%z <= zw(k_surf) .OR. &
1391	.NOT. BTEST( wall_flags_total_0(kp,jp,ip), 0 ) )&
1392	THEN
1393	pos_x = pos_x + pdx(i)
1394	CYCLE xloop
1395	ENDIF
1396
1397	local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1
1398
1399	IF ( .NOT. first_stride ) THEN
1400	IF ( ip < nxl .OR. jp < nys .OR. kp < nzb+1 ) THEN
1401	write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1
1402	ENDIF
1403	IF ( ip > nxr .OR. jp > nyn .OR. kp > nzt ) THEN
1404	write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt
1405	ENDIF
1406	grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
1407	ENDIF
1408	ENDDO
1409	ENDIF
1410	pos_x = pos_x + pdx(i)
1411	ENDDO xloop
1412	ENDIF
1413	pos_y = pos_y + pdy(i)
1414	ENDDO
1415	ENDIF
1416
1417	pos_z = pos_z + pdz(i)
1418	ENDDO
1419	ENDDO
1420
1421	IF ( first_stride ) THEN
1422	DO ip = nxl, nxr
1423	DO jp = nys, nyn
1424	DO kp = nzb+1, nzt
1425	IF ( phase == PHASE_INIT ) THEN
1426	IF ( local_count(kp,jp,ip) > 0 ) THEN
1427	alloc_size = MAX( INT( local_count(kp,jp,ip) * &
1428	( 1.0_wp + alloc_factor / 100.0_wp ) ), &
1429	1 )
1430	ELSE
1431	alloc_size = 1
1432	ENDIF
1433	ALLOCATE(grid_particles(kp,jp,ip)%particles(1:alloc_size))
1434	DO n = 1, alloc_size
1435	grid_particles(kp,jp,ip)%particles(n) = zero_particle
1436	ENDDO
1437	ELSEIF ( phase == PHASE_RELEASE ) THEN
1438	IF ( local_count(kp,jp,ip) > 0 ) THEN
1439	new_size = local_count(kp,jp,ip) + prt_count(kp,jp,ip)
1440	alloc_size = MAX( INT( new_size * ( 1.0_wp + &
1441	alloc_factor / 100.0_wp ) ), 1 )
1442	IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) ) THEN
1443	CALL realloc_particles_array( ip, jp, kp, alloc_size )
1444	ENDIF
1445	ENDIF
1446	ENDIF
1447	ENDDO
1448	ENDDO
1449	ENDDO
1450	ENDIF
1451
1452	ENDDO
1453
1454	local_start = prt_count+1
1455	prt_count = local_count
1456	!
1457	!-- Calculate particle IDs
1458	DO ip = nxl, nxr
1459	DO jp = nys, nyn
1460	DO kp = nzb+1, nzt
1461	number_of_particles = prt_count(kp,jp,ip)
1462	IF ( number_of_particles <= 0 ) CYCLE
1463	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
1464
1465	DO n = local_start(kp,jp,ip), number_of_particles !only new particles
1466
1467	particles(n)%id = 10000_idp*3 grid_particles(kp,jp,ip)%id_counter + &
1468	10000_idp*2 kp + 10000_idp * jp + ip
1469	!
1470	!-- Count the number of particles that have been released before
1471	grid_particles(kp,jp,ip)%id_counter = &
1472	grid_particles(kp,jp,ip)%id_counter + 1
1473
1474	ENDDO
1475
1476	ENDDO
1477	ENDDO
1478	ENDDO
1479	!
1480	!-- Initialize aerosol background spectrum
1481	IF ( curvature_solution_effects ) THEN
1482	CALL lpm_init_aerosols( local_start )
1483	ENDIF
1484	!
1485	!-- Add random fluctuation to particle positions.
1486	IF ( random_start_position ) THEN
1487	DO ip = nxl, nxr
1488	DO jp = nys, nyn
1489	DO kp = nzb+1, nzt
1490	number_of_particles = prt_count(kp,jp,ip)
1491	IF ( number_of_particles <= 0 ) CYCLE
1492	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
1493	!
1494	!-- Move only new particles. Moreover, limit random fluctuation
1495	!-- in order to prevent that particles move more than one grid box,
1496	!-- which would lead to problems concerning particle exchange
1497	!-- between processors in case pdx/pdy are larger than dx/dy,
1498	!-- respectively.
1499	DO n = local_start(kp,jp,ip), number_of_particles
1500	IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN
1501	rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
1502	pdx(particles(n)%group)
1503	particles(n)%x = particles(n)%x + &
1504	MERGE( rand_contr, SIGN( dx, rand_contr ), &
1505	ABS( rand_contr ) < dx &
1506	)
1507	ENDIF
1508	IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN
1509	rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
1510	pdy(particles(n)%group)
1511	particles(n)%y = particles(n)%y + &
1512	MERGE( rand_contr, SIGN( dy, rand_contr ), &
1513	ABS( rand_contr ) < dy &
1514	)
1515	ENDIF
1516	IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN
1517	rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
1518	pdz(particles(n)%group)
1519	particles(n)%z = particles(n)%z + &
1520	MERGE( rand_contr, SIGN( dzw(kp), rand_contr ), &
1521	ABS( rand_contr ) < dzw(kp) &
1522	)
1523	ENDIF
1524	ENDDO
1525	!
1526	!-- Identify particles located outside the model domain and reflect
1527	!-- or absorb them if necessary.
1528	CALL lpm_boundary_conds( 'bottom/top', i, j, k )
1529	!
1530	!-- Furthermore, remove particles located in topography. Note, as
1531	!-- the particle speed is still zero at this point, wall
1532	!-- reflection boundary conditions will not work in this case.
1533	particles => &
1534	grid_particles(kp,jp,ip)%particles(1:number_of_particles)
1535	DO n = local_start(kp,jp,ip), number_of_particles
1536	i = particles(n)%x * ddx
1537	j = particles(n)%y * ddy
1538	k = particles(n)%z / dz(1) + 1 + offset_ocean_nzt
1539	DO WHILE( zw(k) < particles(n)%z )
1540	k = k + 1
1541	ENDDO
1542	DO WHILE( zw(k-1) > particles(n)%z )
1543	k = k - 1
1544	ENDDO
1545	!
1546	!-- Check if particle is within topography
1547	IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 0 ) ) THEN
1548	particles(n)%particle_mask = .FALSE.
1549	deleted_particles = deleted_particles + 1
1550	ENDIF
1551
1552	ENDDO
1553	ENDDO
1554	ENDDO
1555	ENDDO
1556	!
1557	!-- Exchange particles between grid cells and processors
1558	CALL lpm_move_particle
1559	CALL lpm_exchange_horiz
1560
1561	ENDIF
1562	!
1563	!-- In case of random_start_position, delete particles identified by
1564	!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks,
1565	!-- which is needed for a fast interpolation of the LES fields on the particle
1566	!-- position.
1567	CALL lpm_sort_and_delete
1568	!
1569	!-- Determine the current number of particles
1570	DO ip = nxl, nxr
1571	DO jp = nys, nyn
1572	DO kp = nzb+1, nzt
1573	number_of_particles = number_of_particles &
1574	+ prt_count(kp,jp,ip)
1575	ENDDO
1576	ENDDO
1577	ENDDO
1578	!
1579	!-- Calculate the number of particles of the total domain
1580	#if defined( __parallel )
1581	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1582	CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
1583	MPI_INTEGER, MPI_SUM, comm2d, ierr )
1584	#else
1585	total_number_of_particles = number_of_particles
1586	#endif
1587
1588	RETURN
1589
1590	END SUBROUTINE lpm_create_particle
1591
1592
1593	!------------------------------------------------------------------------------!
1594	! Description:
1595	! ------------
1596	!> This routine initialize the particles as aerosols with physio-chemical
1597	!> properties.
1598	!------------------------------------------------------------------------------!
1599	SUBROUTINE lpm_init_aerosols(local_start)
1600
1601	REAL(wp) :: afactor !< curvature effects
1602	REAL(wp) :: bfactor !< solute effects
1603	REAL(wp) :: dlogr !< logarithmic width of radius bin
1604	REAL(wp) :: e_a !< vapor pressure
1605	REAL(wp) :: e_s !< saturation vapor pressure
1606	REAL(wp) :: rmin = 0.005e-6_wp !< minimum aerosol radius
1607	REAL(wp) :: rmax = 10.0e-6_wp !< maximum aerosol radius
1608	REAL(wp) :: r_mid !< mean radius of bin
1609	REAL(wp) :: r_l !< left radius of bin
1610	REAL(wp) :: r_r !< right radius of bin
1611	REAL(wp) :: sigma !< surface tension
1612	REAL(wp) :: t_int !< temperature
1613
1614	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: local_start !<
1615
1616	INTEGER(iwp) :: n !<
1617	INTEGER(iwp) :: ip !<
1618	INTEGER(iwp) :: jp !<
1619	INTEGER(iwp) :: kp !<
1620
1621	!
1622	!-- Set constants for different aerosol species
1623	IF ( TRIM( aero_species ) == 'nacl' ) THEN
1624	molecular_weight_of_solute = 0.05844_wp
1625	rho_s = 2165.0_wp
1626	vanthoff = 2.0_wp
1627	ELSEIF ( TRIM( aero_species ) == 'c3h4o4' ) THEN
1628	molecular_weight_of_solute = 0.10406_wp
1629	rho_s = 1600.0_wp
1630	vanthoff = 1.37_wp
1631	ELSEIF ( TRIM( aero_species ) == 'nh4o3' ) THEN
1632	molecular_weight_of_solute = 0.08004_wp
1633	rho_s = 1720.0_wp
1634	vanthoff = 2.31_wp
1635	ELSE
1636	WRITE( message_string, * ) 'unknown aerosol species ', &
1637	'aero_species = "', TRIM( aero_species ), '"'
1638	CALL message( 'lpm_init', 'PA0470', 1, 2, 0, 6, 0 )
1639	ENDIF
1640	!
1641	!-- The following typical aerosol spectra are taken from Jaenicke (1993):
1642	!-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
1643	IF ( TRIM( aero_type ) == 'polar' ) THEN
1644	na = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6_wp
1645	rm = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6_wp
1646	log_sigma = (/ 0.245, 0.300, 0.291 /)
1647	ELSEIF ( TRIM( aero_type ) == 'background' ) THEN
1648	na = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6_wp
1649	rm = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6_wp
1650	log_sigma = (/ 0.645, 0.253, 0.425 /)
1651	ELSEIF ( TRIM( aero_type ) == 'maritime' ) THEN
1652	na = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6_wp
1653	rm = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6_wp
1654	log_sigma = (/ 0.657, 0.210, 0.396 /)
1655	ELSEIF ( TRIM( aero_type ) == 'continental' ) THEN
1656	na = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6_wp
1657	rm = (/ 0.01, 0.058, 0.9 /) * 1.0E-6_wp
1658	log_sigma = (/ 0.161, 0.217, 0.380 /)
1659	ELSEIF ( TRIM( aero_type ) == 'desert' ) THEN
1660	na = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6_wp
1661	rm = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6_wp
1662	log_sigma = (/ 0.247, 0.770, 0.438 /)
1663	ELSEIF ( TRIM( aero_type ) == 'rural' ) THEN
1664	na = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6_wp
1665	rm = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6_wp
1666	log_sigma = (/ 0.225, 0.557, 0.266 /)
1667	ELSEIF ( TRIM( aero_type ) == 'urban' ) THEN
1668	na = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6_wp
1669	rm = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6_wp
1670	log_sigma = (/ 0.245, 0.666, 0.337 /)
1671	ELSEIF ( TRIM( aero_type ) == 'user' ) THEN
1672	CONTINUE
1673	ELSE
1674	WRITE( message_string, * ) 'unknown aerosol type ', &
1675	'aero_type = "', TRIM( aero_type ), '"'
1676	CALL message( 'lpm_init', 'PA0459', 1, 2, 0, 6, 0 )
1677	ENDIF
1678
1679	DO ip = nxl, nxr
1680	DO jp = nys, nyn
1681	DO kp = nzb+1, nzt
1682
1683	number_of_particles = prt_count(kp,jp,ip)
1684	IF ( number_of_particles <= 0 ) CYCLE
1685	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
1686
1687	dlogr = ( LOG10(rmax) - LOG10(rmin) ) / ( number_of_particles - local_start(kp,jp,ip) + 1 )
1688	!
1689	!-- Initialize the aerosols with a predefined spectral distribution
1690	!-- of the dry radius (logarithmically increasing bins) and a varying
1691	!-- weighting factor
1692	DO n = local_start(kp,jp,ip), number_of_particles !only new particles
1693
1694	r_l = 10.0*( LOG10( rmin ) + (n-1) dlogr )
1695	r_r = 10.0*( LOG10( rmin ) + n dlogr )
1696	r_mid = SQRT( r_l * r_r )
1697
1698	particles(n)%aux1 = r_mid
1699	particles(n)%weight_factor = &
1700	( na(1) / ( SQRT( 2.0_wp * pi ) * log_sigma(1) ) * &
1701	EXP( - LOG10( r_mid / rm(1) )*2 / ( 2.0_wp log_sigma(1)**2 ) ) + &
1702	na(2) / ( SQRT( 2.0_wp * pi ) * log_sigma(2) ) * &
1703	EXP( - LOG10( r_mid / rm(2) )*2 / ( 2.0_wp log_sigma(2)**2 ) ) + &
1704	na(3) / ( SQRT( 2.0_wp * pi ) * log_sigma(3) ) * &
1705	EXP( - LOG10( r_mid / rm(3) )*2 / ( 2.0_wp log_sigma(3)**2 ) ) &
1706	) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dzw(kp) )
1707
1708	!
1709	!-- Multiply weight_factor with the namelist parameter aero_weight
1710	!-- to increase or decrease the number of simulated aerosols
1711	particles(n)%weight_factor = particles(n)%weight_factor * aero_weight
1712
1713	IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) &
1714	> random_function( iran_part ) ) THEN
1715	particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0_wp
1716	ELSE
1717	particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp)
1718	ENDIF
1719	!
1720	!-- Unnecessary particles will be deleted
1721	IF ( particles(n)%weight_factor <= 0.0_wp ) particles(n)%particle_mask = .FALSE.
1722
1723	ENDDO
1724	!
1725	!-- Set particle radius to equilibrium radius based on the environmental
1726	!-- supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids
1727	!-- the sometimes lengthy growth toward their equilibrium radius within
1728	!-- the simulation.
1729	t_int = pt(kp,jp,ip) * exner(kp)
1730
1731	e_s = magnus( t_int )
1732	e_a = q(kp,jp,ip) * hyp(kp) / ( q(kp,jp,ip) + rd_d_rv )
1733
1734	sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
1735	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
1736
1737	bfactor = vanthoff * molecular_weight_of_water * &
1738	rho_s / ( molecular_weight_of_solute * rho_l )
1739	!
1740	!-- The formula is only valid for subsaturated environments. For
1741	!-- supersaturations higher than -5 %, the supersaturation is set to -5%.
1742	IF ( e_a / e_s >= 0.95_wp ) e_a = 0.95_wp * e_s
1743
1744	DO n = local_start(kp,jp,ip), number_of_particles !only new particles
1745	!
1746	!-- For details on this equation, see Eq. (14) of Khvorostyanov and
1747	!-- Curry (2007, JGR)
1748	particles(n)%radius = bfactor*0.3333333_wp &
1749	particles(n)%aux1 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
1750	( 1.0_wp + ( afactor / ( 3.0_wp * bfactor*0.3333333_wp &
1751	particles(n)%aux1 ) ) / &
1752	( 1.0_wp - e_a / e_s )**0.6666666_wp &
1753	)
1754
1755	ENDDO
1756
1757	ENDDO
1758	ENDDO
1759	ENDDO
1760
1761	END SUBROUTINE lpm_init_aerosols
1762
1763
1764	!------------------------------------------------------------------------------!
1765	! Description:
1766	! ------------
1767	!> Calculates quantities required for considering the SGS velocity fluctuations
1768	!> in the particle transport by a stochastic approach. The respective
1769	!> quantities are: SGS-TKE gradients and horizontally averaged profiles of the
1770	!> SGS TKE and the resolved-scale velocity variances.
1771	!------------------------------------------------------------------------------!
1772	SUBROUTINE lpm_init_sgs_tke
1773
1774	USE statistics, &
1775	ONLY: flow_statistics_called, hom, sums, sums_l
1776
1777	INTEGER(iwp) :: i !< index variable along x
1778	INTEGER(iwp) :: j !< index variable along y
1779	INTEGER(iwp) :: k !< index variable along z
1780	INTEGER(iwp) :: m !< running index for the surface elements
1781
1782	REAL(wp) :: flag1 !< flag to mask topography
1783
1784	!
1785	!-- TKE gradient along x and y
1786	DO i = nxl, nxr
1787	DO j = nys, nyn
1788	DO k = nzb, nzt+1
1789
1790	IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 ) .AND. &
1791	BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
1792	BTEST( wall_flags_total_0(k,j,i+1), 0 ) ) &
1793	THEN
1794	de_dx(k,j,i) = 2.0_wp * sgs_wf_part * &
1795	( e(k,j,i+1) - e(k,j,i) ) * ddx
1796	ELSEIF ( BTEST( wall_flags_total_0(k,j,i-1), 0 ) .AND. &
1797	BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
1798	.NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) ) &
1799	THEN
1800	de_dx(k,j,i) = 2.0_wp * sgs_wf_part * &
1801	( e(k,j,i) - e(k,j,i-1) ) * ddx
1802	ELSEIF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 22 ) .AND. &
1803	.NOT. BTEST( wall_flags_total_0(k,j,i+1), 22 ) ) &
1804	THEN
1805	de_dx(k,j,i) = 0.0_wp
1806	ELSEIF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 22 ) .AND. &
1807	.NOT. BTEST( wall_flags_total_0(k,j,i), 22 ) ) &
1808	THEN
1809	de_dx(k,j,i) = 0.0_wp
1810	ELSE
1811	de_dx(k,j,i) = sgs_wf_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
1812	ENDIF
1813
1814	IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 ) .AND. &
1815	BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
1816	BTEST( wall_flags_total_0(k,j+1,i), 0 ) ) &
1817	THEN
1818	de_dy(k,j,i) = 2.0_wp * sgs_wf_part * &
1819	( e(k,j+1,i) - e(k,j,i) ) * ddy
1820	ELSEIF ( BTEST( wall_flags_total_0(k,j-1,i), 0 ) .AND. &
1821	BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. &
1822	.NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) ) &
1823	THEN
1824	de_dy(k,j,i) = 2.0_wp * sgs_wf_part * &
1825	( e(k,j,i) - e(k,j-1,i) ) * ddy
1826	ELSEIF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 22 ) .AND. &
1827	.NOT. BTEST( wall_flags_total_0(k,j+1,i), 22 ) ) &
1828	THEN
1829	de_dy(k,j,i) = 0.0_wp
1830	ELSEIF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 22 ) .AND. &
1831	.NOT. BTEST( wall_flags_total_0(k,j,i), 22 ) ) &
1832	THEN
1833	de_dy(k,j,i) = 0.0_wp
1834	ELSE
1835	de_dy(k,j,i) = sgs_wf_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
1836	ENDIF
1837
1838	ENDDO
1839	ENDDO
1840	ENDDO
1841
1842	!
1843	!-- TKE gradient along z at topograhy and including bottom and top boundary conditions
1844	DO i = nxl, nxr
1845	DO j = nys, nyn
1846	DO k = nzb+1, nzt-1
1847	!
1848	!-- Flag to mask topography
1849	flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
1850
1851	de_dz(k,j,i) = 2.0_wp * sgs_wf_part * &
1852	( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1) - zu(k-1) ) &
1853	* flag1
1854	ENDDO
1855	!
1856	!-- upward-facing surfaces
1857	DO m = bc_h(0)%start_index(j,i), bc_h(0)%end_index(j,i)
1858	k = bc_h(0)%k(m)
1859	de_dz(k,j,i) = 2.0_wp * sgs_wf_part * &
1860	( e(k+1,j,i) - e(k,j,i) ) / ( zu(k+1) - zu(k) )
1861	ENDDO
1862	!
1863	!-- downward-facing surfaces
1864	DO m = bc_h(1)%start_index(j,i), bc_h(1)%end_index(j,i)
1865	k = bc_h(1)%k(m)
1866	de_dz(k,j,i) = 2.0_wp * sgs_wf_part * &
1867	( e(k,j,i) - e(k-1,j,i) ) / ( zu(k) - zu(k-1) )
1868	ENDDO
1869
1870	de_dz(nzb,j,i) = 0.0_wp
1871	de_dz(nzt,j,i) = 0.0_wp
1872	de_dz(nzt+1,j,i) = 0.0_wp
1873	ENDDO
1874	ENDDO
1875	!
1876	!-- Ghost point exchange
1877	CALL exchange_horiz( de_dx, nbgp )
1878	CALL exchange_horiz( de_dy, nbgp )
1879	CALL exchange_horiz( de_dz, nbgp )
1880	CALL exchange_horiz( diss, nbgp )
1881	!
1882	!-- Set boundary conditions at non-periodic boundaries. Note, at non-period
1883	!-- boundaries zero-gradient boundary conditions are set for the subgrid TKE.
1884	!-- Thus, TKE gradients normal to the respective lateral boundaries are zero,
1885	!-- while tangetial TKE gradients then must be the same as within the prognostic
1886	!-- domain.
1887	IF ( bc_dirichlet_l ) THEN
1888	de_dx(:,:,-1) = 0.0_wp
1889	de_dy(:,:,-1) = de_dy(:,:,0)
1890	de_dz(:,:,-1) = de_dz(:,:,0)
1891	ENDIF
1892	IF ( bc_dirichlet_r ) THEN
1893	de_dx(:,:,nxr+1) = 0.0_wp
1894	de_dy(:,:,nxr+1) = de_dy(:,:,nxr)
1895	de_dz(:,:,nxr+1) = de_dz(:,:,nxr)
1896	ENDIF
1897	IF ( bc_dirichlet_n ) THEN
1898	de_dx(:,nyn+1,:) = de_dx(:,nyn,:)
1899	de_dy(:,nyn+1,:) = 0.0_wp
1900	de_dz(:,nyn+1,:) = de_dz(:,nyn,:)
1901	ENDIF
1902	IF ( bc_dirichlet_s ) THEN
1903	de_dx(:,nys-1,:) = de_dx(:,nys,:)
1904	de_dy(:,nys-1,:) = 0.0_wp
1905	de_dz(:,nys-1,:) = de_dz(:,nys,:)
1906	ENDIF
1907	!
1908	!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
1909	!-- velocity variances (they may have been already calculated in routine
1910	!-- flow_statistics).
1911	IF ( .NOT. flow_statistics_called ) THEN
1912
1913	!
1914	!-- First calculate horizontally averaged profiles of the horizontal
1915	!-- velocities.
1916	sums_l(:,1,0) = 0.0_wp
1917	sums_l(:,2,0) = 0.0_wp
1918
1919	DO i = nxl, nxr
1920	DO j = nys, nyn
1921	DO k = nzb, nzt+1
1922	!
1923	!-- Flag indicating vicinity of wall
1924	flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 24 ) )
1925
1926	sums_l(k,1,0) = sums_l(k,1,0) + u(k,j,i) * flag1
1927	sums_l(k,2,0) = sums_l(k,2,0) + v(k,j,i) * flag1
1928	ENDDO
1929	ENDDO
1930	ENDDO
1931
1932	#if defined( __parallel )
1933	!
1934	!-- Compute total sum from local sums
1935	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1936	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
1937	MPI_REAL, MPI_SUM, comm2d, ierr )
1938	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1939	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
1940	MPI_REAL, MPI_SUM, comm2d, ierr )
1941	#else
1942	sums(:,1) = sums_l(:,1,0)
1943	sums(:,2) = sums_l(:,2,0)
1944	#endif
1945
1946	!
1947	!-- Final values are obtained by division by the total number of grid
1948	!-- points used for the summation.
1949	hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0) ! u
1950	hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0) ! v
1951
1952	!
1953	!-- Now calculate the profiles of SGS TKE and the resolved-scale
1954	!-- velocity variances
1955	sums_l(:,8,0) = 0.0_wp
1956	sums_l(:,30,0) = 0.0_wp
1957	sums_l(:,31,0) = 0.0_wp
1958	sums_l(:,32,0) = 0.0_wp
1959	DO i = nxl, nxr
1960	DO j = nys, nyn
1961	DO k = nzb, nzt+1
1962	!
1963	!-- Flag indicating vicinity of wall
1964	flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 24 ) )
1965
1966	sums_l(k,8,0) = sums_l(k,8,0) + e(k,j,i) * flag1
1967	sums_l(k,30,0) = sums_l(k,30,0) + ( u(k,j,i) - hom(k,1,1,0) )*2 flag1
1968	sums_l(k,31,0) = sums_l(k,31,0) + ( v(k,j,i) - hom(k,1,2,0) )*2 flag1
1969	sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)*2 flag1
1970	ENDDO
1971	ENDDO
1972	ENDDO
1973
1974	#if defined( __parallel )
1975	!
1976	!-- Compute total sum from local sums
1977	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1978	CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
1979	MPI_REAL, MPI_SUM, comm2d, ierr )
1980	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1981	CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
1982	MPI_REAL, MPI_SUM, comm2d, ierr )
1983	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1984	CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
1985	MPI_REAL, MPI_SUM, comm2d, ierr )
1986	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1987	CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
1988	MPI_REAL, MPI_SUM, comm2d, ierr )
1989
1990	#else
1991	sums(:,8) = sums_l(:,8,0)
1992	sums(:,30) = sums_l(:,30,0)
1993	sums(:,31) = sums_l(:,31,0)
1994	sums(:,32) = sums_l(:,32,0)
1995	#endif
1996
1997	!
1998	!-- Final values are obtained by division by the total number of grid
1999	!-- points used for the summation.
2000	hom(:,1,8,0) = sums(:,8) / ngp_2dh_outer(:,0) ! e
2001	hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0) ! u*2
2002	hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0) ! v*2
2003	hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0) ! w*2
2004
2005	ENDIF
2006
2007	END SUBROUTINE lpm_init_sgs_tke
2008
2009
2010	!------------------------------------------------------------------------------!
2011	! Description:
2012	! ------------
2013	!> Sobroutine control lpm actions, i.e. all actions during one time step.
2014	!------------------------------------------------------------------------------!
2015	SUBROUTINE lpm_actions( location )
2016
2017	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
2018
2019	INTEGER(iwp) :: i !<
2020	INTEGER(iwp) :: ie !<
2021	INTEGER(iwp) :: is !<
2022	INTEGER(iwp) :: j !<
2023	INTEGER(iwp) :: je !<
2024	INTEGER(iwp) :: js !<
2025	INTEGER(iwp), SAVE :: lpm_count = 0 !<
2026	INTEGER(iwp) :: k !<
2027	INTEGER(iwp) :: ke !<
2028	INTEGER(iwp) :: ks !<
2029	INTEGER(iwp) :: m !<
2030	INTEGER(iwp), SAVE :: steps = 0 !<
2031
2032	LOGICAL :: first_loop_stride !<
2033
2034
2035	SELECT CASE ( location )
2036
2037	CASE ( 'after_pressure_solver' )
2038	!
2039	!-- The particle model is executed if particle advection start is reached and only at the end
2040	!-- of the intermediate time step loop.
2041	IF ( time_since_reference_point >= particle_advection_start &
2042	.AND. intermediate_timestep_count == intermediate_timestep_count_max ) &
2043	THEN
2044	CALL cpu_log( log_point(25), 'lpm', 'start' )
2045	!
2046	!-- Write particle data at current time on file.
2047	!-- This has to be done here, before particles are further processed,
2048	!-- because they may be deleted within this timestep (in case that
2049	!-- dt_write_particle_data = dt_prel = particle_maximum_age).
2050	time_write_particle_data = time_write_particle_data + dt_3d
2051	IF ( time_write_particle_data >= dt_write_particle_data ) THEN
2052
2053	CALL lpm_data_output_particles
2054	!
2055	!-- The MOD function allows for changes in the output interval with restart
2056	!-- runs.
2057	time_write_particle_data = MOD( time_write_particle_data, &
2058	MAX( dt_write_particle_data, dt_3d ) )
2059	ENDIF
2060
2061	!
2062	!-- Initialize arrays for marking those particles to be deleted after the
2063	!-- (sub-) timestep
2064	deleted_particles = 0
2065
2066	!
2067	!-- Initialize variables used for accumulating the number of particles
2068	!-- xchanged between the subdomains during all sub-timesteps (if sgs
2069	!-- velocities are included). These data are output further below on the
2070	!-- particle statistics file.
2071	trlp_count_sum = 0
2072	trlp_count_recv_sum = 0
2073	trrp_count_sum = 0
2074	trrp_count_recv_sum = 0
2075	trsp_count_sum = 0
2076	trsp_count_recv_sum = 0
2077	trnp_count_sum = 0
2078	trnp_count_recv_sum = 0
2079	!
2080	!-- Calculate exponential term used in case of particle inertia for each
2081	!-- of the particle groups
2082	DO m = 1, number_of_particle_groups
2083	IF ( particle_groups(m)%density_ratio /= 0.0_wp ) THEN
2084	particle_groups(m)%exp_arg = &
2085	4.5_wp * particle_groups(m)%density_ratio * &
2086	molecular_viscosity / ( particle_groups(m)%radius )**2
2087
2088	particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
2089	dt_3d )
2090	ENDIF
2091	ENDDO
2092	!
2093	!-- If necessary, release new set of particles
2094	IF ( ( simulated_time - last_particle_release_time ) >= dt_prel .AND. &
2095	end_time_prel > simulated_time ) THEN
2096	DO WHILE ( ( simulated_time - last_particle_release_time ) >= dt_prel )
2097	CALL lpm_create_particle( PHASE_RELEASE )
2098	last_particle_release_time = last_particle_release_time + dt_prel
2099	ENDDO
2100	ENDIF
2101	!
2102	!-- Reset summation arrays
2103	IF ( cloud_droplets ) THEN
2104	ql_c = 0.0_wp
2105	ql_v = 0.0_wp
2106	ql_vp = 0.0_wp
2107	ENDIF
2108
2109	first_loop_stride = .TRUE.
2110	grid_particles(:,:,:)%time_loop_done = .TRUE.
2111	!
2112	!-- Timestep loop for particle advection.
2113	!-- This loop has to be repeated until the advection time of every particle
2114	!-- (within the total domain!) has reached the LES timestep (dt_3d).
2115	!-- In case of including the SGS velocities, the particle timestep may be
2116	!-- smaller than the LES timestep (because of the Lagrangian timescale
2117	!-- restriction) and particles may require to undergo several particle
2118	!-- timesteps, before the LES timestep is reached. Because the number of these
2119	!-- particle timesteps to be carried out is unknown at first, these steps are
2120	!-- carried out in the following infinite loop with exit condition.
2121	DO
2122	CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' )
2123	CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
2124
2125	!
2126	!-- If particle advection includes SGS velocity components, calculate the
2127	!-- required SGS quantities (i.e. gradients of the TKE, as well as
2128	!-- horizontally averaged profiles of the SGS TKE and the resolved-scale
2129	!-- velocity variances)
2130	IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN
2131	CALL lpm_init_sgs_tke
2132	ENDIF
2133	!
2134	!-- In case SGS-particle speed is considered, particles may carry out
2135	!-- several particle timesteps. In order to prevent unnecessary
2136	!-- treatment of particles that already reached the final time level,
2137	!-- particles are sorted into contiguous blocks of finished and
2138	!-- not-finished particles, in addition to their already sorting
2139	!-- according to their sub-boxes.
2140	IF ( .NOT. first_loop_stride .AND. use_sgs_for_particles ) &
2141	CALL lpm_sort_timeloop_done
2142	DO i = nxl, nxr
2143	DO j = nys, nyn
2144	DO k = nzb+1, nzt
2145
2146	number_of_particles = prt_count(k,j,i)
2147	!
2148	!-- If grid cell gets empty, flag must be true
2149	IF ( number_of_particles <= 0 ) THEN
2150	grid_particles(k,j,i)%time_loop_done = .TRUE.
2151	CYCLE
2152	ENDIF
2153
2154	IF ( .NOT. first_loop_stride .AND. &
2155	grid_particles(k,j,i)%time_loop_done ) CYCLE
2156
2157	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
2158
2159	particles(1:number_of_particles)%particle_mask = .TRUE.
2160	!
2161	!-- Initialize the variable storing the total time that a particle
2162	!-- has advanced within the timestep procedure
2163	IF ( first_loop_stride ) THEN
2164	particles(1:number_of_particles)%dt_sum = 0.0_wp
2165	ENDIF
2166	!
2167	!-- Particle (droplet) growth by condensation/evaporation and
2168	!-- collision
2169	IF ( cloud_droplets .AND. first_loop_stride) THEN
2170	!
2171	!-- Droplet growth by condensation / evaporation
2172	CALL lpm_droplet_condensation(i,j,k)
2173	!
2174	!-- Particle growth by collision
2175	IF ( collision_kernel /= 'none' ) THEN
2176	CALL lpm_droplet_collision(i,j,k)
2177	ENDIF
2178
2179	ENDIF
2180	!
2181	!-- Initialize the switch used for the loop exit condition checked
2182	!-- at the end of this loop. If at least one particle has failed to
2183	!-- reach the LES timestep, this switch will be set false in
2184	!-- lpm_advec.
2185	dt_3d_reached_l = .TRUE.
2186
2187	!
2188	!-- Particle advection
2189	CALL lpm_advec( i, j, k )
2190	!
2191	!-- Particle reflection from walls. Only applied if the particles
2192	!-- are in the vertical range of the topography. (Here, some
2193	!-- optimization is still possible.)
2194	IF ( topography /= 'flat' .AND. k < nzb_max + 2 ) THEN
2195	CALL lpm_boundary_conds( 'walls', i, j, k )
2196	ENDIF
2197	!
2198	!-- User-defined actions after the calculation of the new particle
2199	!-- position
2200	CALL user_lpm_advec( i, j, k )
2201	!
2202	!-- Apply boundary conditions to those particles that have crossed
2203	!-- the top or bottom boundary and delete those particles, which are
2204	!-- older than allowed
2205	CALL lpm_boundary_conds( 'bottom/top', i, j, k )
2206	!
2207	!--- If not all particles of the actual grid cell have reached the
2208	!-- LES timestep, this cell has to do another loop iteration. Due to
2209	!-- the fact that particles can move into neighboring grid cells,
2210	!-- these neighbor cells also have to perform another loop iteration.
2211	!-- Please note, this realization does not work properly if
2212	!-- particles move into another subdomain.
2213	IF ( .NOT. dt_3d_reached_l ) THEN
2214	ks = MAX(nzb+1,k-1)
2215	ke = MIN(nzt,k+1)
2216	js = MAX(nys,j-1)
2217	je = MIN(nyn,j+1)
2218	is = MAX(nxl,i-1)
2219	ie = MIN(nxr,i+1)
2220	grid_particles(ks:ke,js:je,is:ie)%time_loop_done = .FALSE.
2221	ELSE
2222	grid_particles(k,j,i)%time_loop_done = .TRUE.
2223	ENDIF
2224
2225	ENDDO
2226	ENDDO
2227	ENDDO
2228	steps = steps + 1
2229	dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done)
2230	!
2231	!-- Find out, if all particles on every PE have completed the LES timestep
2232	!-- and set the switch corespondingly
2233	#if defined( __parallel )
2234	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2235	CALL MPI_ALLREDUCE( dt_3d_reached_l, dt_3d_reached, 1, MPI_LOGICAL, &
2236	MPI_LAND, comm2d, ierr )
2237	#else
2238	dt_3d_reached = dt_3d_reached_l
2239	#endif
2240	CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' )
2241
2242	!
2243	!-- Apply splitting and merging algorithm
2244	IF ( cloud_droplets ) THEN
2245	IF ( splitting ) THEN
2246	CALL lpm_splitting
2247	ENDIF
2248	IF ( merging ) THEN
2249	CALL lpm_merging
2250	ENDIF
2251	ENDIF
2252	!
2253	!-- Move Particles local to PE to a different grid cell
2254	CALL lpm_move_particle
2255	!
2256	!-- Horizontal boundary conditions including exchange between subdmains
2257	CALL lpm_exchange_horiz
2258
2259	!
2260	!-- IF .FALSE., lpm_sort_and_delete is done inside pcmp
2261	IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run ) THEN
2262	!
2263	!-- Pack particles (eliminate those marked for deletion),
2264	!-- determine new number of particles
2265	CALL lpm_sort_and_delete
2266
2267	!-- Initialize variables for the next (sub-) timestep, i.e., for marking
2268	!-- those particles to be deleted after the timestep
2269	deleted_particles = 0
2270	ENDIF
2271
2272	IF ( dt_3d_reached ) EXIT
2273
2274	first_loop_stride = .FALSE.
2275	ENDDO ! timestep loop
2276
2277	#if defined( __parallel )
2278	!
2279	!-- in case of nested runs do the transfer of particles after every full model time step
2280	IF ( nested_run ) THEN
2281	CALL particles_from_parent_to_child
2282	CALL particles_from_child_to_parent
2283	CALL pmcp_p_delete_particles_in_fine_grid_area
2284
2285	CALL lpm_sort_and_delete
2286
2287	deleted_particles = 0
2288	ENDIF
2289	#endif
2290
2291	!
2292	!-- Calculate the new liquid water content for each grid box
2293	IF ( cloud_droplets ) CALL lpm_calc_liquid_water_content
2294
2295	!
2296	!-- At the end all arrays are exchanged
2297	IF ( cloud_droplets ) THEN
2298	CALL exchange_horiz( ql, nbgp )
2299	CALL exchange_horiz( ql_c, nbgp )
2300	CALL exchange_horiz( ql_v, nbgp )
2301	CALL exchange_horiz( ql_vp, nbgp )
2302	ENDIF
2303
2304	!
2305	!-- Deallocate unused memory
2306	IF ( deallocate_memory .AND. lpm_count == step_dealloc ) THEN
2307	CALL dealloc_particles_array
2308	lpm_count = 0
2309	ELSEIF ( deallocate_memory ) THEN
2310	lpm_count = lpm_count + 1
2311	ENDIF
2312
2313	!
2314	!-- Write particle statistics (in particular the number of particles
2315	!-- exchanged between the subdomains) on file
2316	IF ( write_particle_statistics ) CALL lpm_write_exchange_statistics
2317	!
2318	!-- Execute Interactions of condnesation and evaporation to humidity and
2319	!-- temperature field
2320	IF ( cloud_droplets ) THEN
2321	CALL lpm_interaction_droplets_ptq
2322	CALL exchange_horiz( pt, nbgp )
2323	CALL exchange_horiz( q, nbgp )
2324	ENDIF
2325
2326	CALL cpu_log( log_point(25), 'lpm', 'stop' )
2327
2328	! !
2329	! !-- Output of particle time series
2330	! IF ( particle_advection ) THEN
2331	! IF ( time_dopts >= dt_dopts .OR. &
2332	! ( time_since_reference_point >= particle_advection_start .AND. &
2333	! first_call_lpm ) ) THEN
2334	! CALL lpm_data_output_ptseries
2335	! time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
2336	! ENDIF
2337	! ENDIF
2338
2339	!
2340	!-- Set this switch to .false. @todo: maybe find better solution.
2341	first_call_lpm = .FALSE.
2342	ENDIF! ENDIF statement of lpm_actions('after_pressure_solver')
2343
2344	CASE ( 'after_integration' )
2345	!
2346	!-- Call at the end of timestep routine to save particle velocities fields
2347	!-- for the next timestep
2348	CALL lpm_swap_timelevel_for_particle_advection
2349
2350	CASE DEFAULT
2351	CONTINUE
2352
2353	END SELECT
2354
2355	END SUBROUTINE lpm_actions
2356
2357
2358	#if defined( __parallel )
2359	!------------------------------------------------------------------------------!
2360	! Description:
2361	! ------------
2362	!
2363	!------------------------------------------------------------------------------!
2364	SUBROUTINE particles_from_parent_to_child
2365
2366	CALL pmcp_c_get_particle_from_parent ! Child actions
2367	CALL pmcp_p_fill_particle_win ! Parent actions
2368
2369	RETURN
2370
2371	END SUBROUTINE particles_from_parent_to_child
2372
2373
2374	!------------------------------------------------------------------------------!
2375	! Description:
2376	! ------------
2377	!
2378	!------------------------------------------------------------------------------!
2379	SUBROUTINE particles_from_child_to_parent
2380
2381	CALL pmcp_c_send_particle_to_parent ! Child actions
2382	CALL pmcp_p_empty_particle_win ! Parent actions
2383
2384	RETURN
2385
2386	END SUBROUTINE particles_from_child_to_parent
2387	#endif
2388
2389	!------------------------------------------------------------------------------!
2390	! Description:
2391	! ------------
2392	!> This routine write exchange statistics of the lpm in a ascii file.
2393	!------------------------------------------------------------------------------!
2394	SUBROUTINE lpm_write_exchange_statistics
2395
2396	INTEGER(iwp) :: ip !<
2397	INTEGER(iwp) :: jp !<
2398	INTEGER(iwp) :: kp !<
2399	INTEGER(iwp) :: tot_number_of_particles !<
2400
2401	!
2402	!-- Determine the current number of particles
2403	number_of_particles = 0
2404	DO ip = nxl, nxr
2405	DO jp = nys, nyn
2406	DO kp = nzb+1, nzt
2407	number_of_particles = number_of_particles &
2408	+ prt_count(kp,jp,ip)
2409	ENDDO
2410	ENDDO
2411	ENDDO
2412
2413	CALL check_open( 80 )
2414	#if defined( __parallel )
2415	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
2416	number_of_particles, pleft, trlp_count_sum, &
2417	trlp_count_recv_sum, pright, trrp_count_sum, &
2418	trrp_count_recv_sum, psouth, trsp_count_sum, &
2419	trsp_count_recv_sum, pnorth, trnp_count_sum, &
2420	trnp_count_recv_sum
2421	#else
2422	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
2423	number_of_particles
2424	#endif
2425	CALL close_file( 80 )
2426
2427	IF ( number_of_particles > 0 ) THEN
2428	WRITE(9,*) 'number_of_particles ', number_of_particles, &
2429	current_timestep_number + 1, simulated_time + dt_3d
2430	ENDIF
2431
2432	#if defined( __parallel )
2433	CALL MPI_ALLREDUCE( number_of_particles, tot_number_of_particles, 1, &
2434	MPI_INTEGER, MPI_SUM, comm2d, ierr )
2435	#else
2436	tot_number_of_particles = number_of_particles
2437	#endif
2438
2439	#if defined( __parallel )
2440	IF ( nested_run ) THEN
2441	CALL pmcp_g_print_number_of_particles( simulated_time+dt_3d, &
2442	tot_number_of_particles)
2443	ENDIF
2444	#endif
2445
2446	!
2447	!-- Formats
2448	8000 FORMAT (I6,1X,F7.2,4X,I10,5X,4(I3,1X,I4,'/',I4,2X),6X,I10)
2449
2450
2451	END SUBROUTINE lpm_write_exchange_statistics
2452
2453
2454	!------------------------------------------------------------------------------!
2455	! Description:
2456	! ------------
2457	!> Write particle data in FORTRAN binary and/or netCDF format
2458	!------------------------------------------------------------------------------!
2459	SUBROUTINE lpm_data_output_particles
2460
2461	INTEGER(iwp) :: ip !<
2462	INTEGER(iwp) :: jp !<
2463	INTEGER(iwp) :: kp !<
2464
2465	CALL cpu_log( log_point_s(40), 'lpm_data_output', 'start' )
2466
2467	!
2468	!-- Attention: change version number for unit 85 (in routine check_open)
2469	!-- whenever the output format for this unit is changed!
2470	CALL check_open( 85 )
2471
2472	WRITE ( 85 ) simulated_time
2473	WRITE ( 85 ) prt_count
2474
2475	DO ip = nxl, nxr
2476	DO jp = nys, nyn
2477	DO kp = nzb+1, nzt
2478	number_of_particles = prt_count(kp,jp,ip)
2479	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
2480	IF ( number_of_particles <= 0 ) CYCLE
2481	WRITE ( 85 ) particles
2482	ENDDO
2483	ENDDO
2484	ENDDO
2485
2486	CALL close_file( 85 )
2487
2488
2489	#if defined( __netcdf )
2490	! !
2491	! !-- Output in netCDF format
2492	! CALL check_open( 108 )
2493	!
2494	! !
2495	! !-- Update the NetCDF time axis
2496	! prt_time_count = prt_time_count + 1
2497	!
2498	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, &
2499	! (/ simulated_time /), &
2500	! start = (/ prt_time_count /), count = (/ 1 /) )
2501	! CALL netcdf_handle_error( 'lpm_data_output_particles', 1 )
2502	!
2503	! !
2504	! !-- Output the real number of particles used
2505	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, &
2506	! (/ number_of_particles /), &
2507	! start = (/ prt_time_count /), count = (/ 1 /) )
2508	! CALL netcdf_handle_error( 'lpm_data_output_particles', 2 )
2509	!
2510	! !
2511	! !-- Output all particle attributes
2512	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age, &
2513	! start = (/ 1, prt_time_count /), &
2514	! count = (/ maximum_number_of_particles /) )
2515	! CALL netcdf_handle_error( 'lpm_data_output_particles', 3 )
2516	!
2517	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%user, &
2518	! start = (/ 1, prt_time_count /), &
2519	! count = (/ maximum_number_of_particles /) )
2520	! CALL netcdf_handle_error( 'lpm_data_output_particles', 4 )
2521	!
2522	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x, &
2523	! start = (/ 1, prt_time_count /), &
2524	! count = (/ maximum_number_of_particles /) )
2525	! CALL netcdf_handle_error( 'lpm_data_output_particles', 5 )
2526	!
2527	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y, &
2528	! start = (/ 1, prt_time_count /), &
2529	! count = (/ maximum_number_of_particles /) )
2530	! CALL netcdf_handle_error( 'lpm_data_output_particles', 6 )
2531	!
2532	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z, &
2533	! start = (/ 1, prt_time_count /), &
2534	! count = (/ maximum_number_of_particles /) )
2535	! CALL netcdf_handle_error( 'lpm_data_output_particles', 7 )
2536	!
2537	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius, &
2538	! start = (/ 1, prt_time_count /), &
2539	! count = (/ maximum_number_of_particles /) )
2540	! CALL netcdf_handle_error( 'lpm_data_output_particles', 8 )
2541	!
2542	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x, &
2543	! start = (/ 1, prt_time_count /), &
2544	! count = (/ maximum_number_of_particles /) )
2545	! CALL netcdf_handle_error( 'lpm_data_output_particles', 9 )
2546	!
2547	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y, &
2548	! start = (/ 1, prt_time_count /), &
2549	! count = (/ maximum_number_of_particles /) )
2550	! CALL netcdf_handle_error( 'lpm_data_output_particles', 10 )
2551	!
2552	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z, &
2553	! start = (/ 1, prt_time_count /), &
2554	! count = (/ maximum_number_of_particles /) )
2555	! CALL netcdf_handle_error( 'lpm_data_output_particles', 11 )
2556	!
2557	! nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10), &
2558	! particles%weight_factor, &
2559	! start = (/ 1, prt_time_count /), &
2560	! count = (/ maximum_number_of_particles /) )
2561	! CALL netcdf_handle_error( 'lpm_data_output_particles', 12 )
2562	!
2563	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x, &
2564	! start = (/ 1, prt_time_count /), &
2565	! count = (/ maximum_number_of_particles /) )
2566	! CALL netcdf_handle_error( 'lpm_data_output_particles', 13 )
2567	!
2568	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y, &
2569	! start = (/ 1, prt_time_count /), &
2570	! count = (/ maximum_number_of_particles /) )
2571	! CALL netcdf_handle_error( 'lpm_data_output_particles', 14 )
2572	!
2573	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z, &
2574	! start = (/ 1, prt_time_count /), &
2575	! count = (/ maximum_number_of_particles /) )
2576	! CALL netcdf_handle_error( 'lpm_data_output_particles', 15 )
2577	!
2578	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%class, &
2579	! start = (/ 1, prt_time_count /), &
2580	! count = (/ maximum_number_of_particles /) )
2581	! CALL netcdf_handle_error( 'lpm_data_output_particles', 16 )
2582	!
2583	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group, &
2584	! start = (/ 1, prt_time_count /), &
2585	! count = (/ maximum_number_of_particles /) )
2586	! CALL netcdf_handle_error( 'lpm_data_output_particles', 17 )
2587	!
2588	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16), &
2589	! particles%id2, &
2590	! start = (/ 1, prt_time_count /), &
2591	! count = (/ maximum_number_of_particles /) )
2592	! CALL netcdf_handle_error( 'lpm_data_output_particles', 18 )
2593	!
2594	! nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%id1, &
2595	! start = (/ 1, prt_time_count /), &
2596	! count = (/ maximum_number_of_particles /) )
2597	! CALL netcdf_handle_error( 'lpm_data_output_particles', 19 )
2598	!
2599	#endif
2600
2601	CALL cpu_log( log_point_s(40), 'lpm_data_output', 'stop' )
2602
2603	END SUBROUTINE lpm_data_output_particles
2604
2605	!------------------------------------------------------------------------------!
2606	! Description:
2607	! ------------
2608	!> This routine calculates and provide particle timeseries output.
2609	!------------------------------------------------------------------------------!
2610	SUBROUTINE lpm_data_output_ptseries
2611
2612	INTEGER(iwp) :: i !<
2613	INTEGER(iwp) :: inum !<
2614	INTEGER(iwp) :: j !<
2615	INTEGER(iwp) :: jg !<
2616	INTEGER(iwp) :: k !<
2617	INTEGER(iwp) :: n !<
2618
2619	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pts_value !<
2620	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pts_value_l !<
2621
2622
2623	CALL cpu_log( log_point(36), 'data_output_ptseries', 'start' )
2624
2625	IF ( myid == 0 ) THEN
2626	!
2627	!-- Open file for time series output in NetCDF format
2628	dopts_time_count = dopts_time_count + 1
2629	CALL check_open( 109 )
2630	#if defined( __netcdf )
2631	!
2632	!-- Update the particle time series time axis
2633	nc_stat = NF90_PUT_VAR( id_set_pts, id_var_time_pts, &
2634	(/ time_since_reference_point /), &
2635	start = (/ dopts_time_count /), count = (/ 1 /) )
2636	CALL netcdf_handle_error( 'data_output_ptseries', 391 )
2637	#endif
2638
2639	ENDIF
2640
2641	ALLOCATE( pts_value(0:number_of_particle_groups,dopts_num), &
2642	pts_value_l(0:number_of_particle_groups,dopts_num) )
2643
2644	pts_value_l = 0.0_wp
2645	pts_value_l(:,16) = 9999999.9_wp ! for calculation of minimum radius
2646
2647	!
2648	!-- Calculate or collect the particle time series quantities for all particles
2649	!-- and seperately for each particle group (if there is more than one group)
2650	DO i = nxl, nxr
2651	DO j = nys, nyn
2652	DO k = nzb, nzt
2653	number_of_particles = prt_count(k,j,i)
2654	IF (number_of_particles <= 0) CYCLE
2655	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
2656	DO n = 1, number_of_particles
2657
2658	IF ( particles(n)%particle_mask ) THEN ! Restrict analysis to active particles
2659
2660	pts_value_l(0,1) = pts_value_l(0,1) + 1.0_wp ! total # of particles
2661	pts_value_l(0,2) = pts_value_l(0,2) + &
2662	( particles(n)%x - particles(n)%origin_x ) ! mean x
2663	pts_value_l(0,3) = pts_value_l(0,3) + &
2664	( particles(n)%y - particles(n)%origin_y ) ! mean y
2665	pts_value_l(0,4) = pts_value_l(0,4) + &
2666	( particles(n)%z - particles(n)%origin_z ) ! mean z
2667	pts_value_l(0,5) = pts_value_l(0,5) + particles(n)%z ! mean z (absolute)
2668	pts_value_l(0,6) = pts_value_l(0,6) + particles(n)%speed_x ! mean u
2669	pts_value_l(0,7) = pts_value_l(0,7) + particles(n)%speed_y ! mean v
2670	pts_value_l(0,8) = pts_value_l(0,8) + particles(n)%speed_z ! mean w
2671	pts_value_l(0,9) = pts_value_l(0,9) + particles(n)%rvar1 ! mean sgsu
2672	pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
2673	pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
2674	IF ( particles(n)%speed_z > 0.0_wp ) THEN
2675	pts_value_l(0,12) = pts_value_l(0,12) + 1.0_wp ! # of upward moving prts
2676	pts_value_l(0,13) = pts_value_l(0,13) + &
2677	particles(n)%speed_z ! mean w upw.
2678	ELSE
2679	pts_value_l(0,14) = pts_value_l(0,14) + &
2680	particles(n)%speed_z ! mean w down
2681	ENDIF
2682	pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius ! mean rad
2683	pts_value_l(0,16) = MIN( pts_value_l(0,16), particles(n)%radius ) ! minrad
2684	pts_value_l(0,17) = MAX( pts_value_l(0,17), particles(n)%radius ) ! maxrad
2685	pts_value_l(0,18) = pts_value_l(0,18) + 1.0_wp
2686	pts_value_l(0,19) = pts_value_l(0,18) + 1.0_wp
2687	!
2688	!-- Repeat the same for the respective particle group
2689	IF ( number_of_particle_groups > 1 ) THEN
2690	jg = particles(n)%group
2691
2692	pts_value_l(jg,1) = pts_value_l(jg,1) + 1.0_wp
2693	pts_value_l(jg,2) = pts_value_l(jg,2) + &
2694	( particles(n)%x - particles(n)%origin_x )
2695	pts_value_l(jg,3) = pts_value_l(jg,3) + &
2696	( particles(n)%y - particles(n)%origin_y )
2697	pts_value_l(jg,4) = pts_value_l(jg,4) + &
2698	( particles(n)%z - particles(n)%origin_z )
2699	pts_value_l(jg,5) = pts_value_l(jg,5) + particles(n)%z
2700	pts_value_l(jg,6) = pts_value_l(jg,6) + particles(n)%speed_x
2701	pts_value_l(jg,7) = pts_value_l(jg,7) + particles(n)%speed_y
2702	pts_value_l(jg,8) = pts_value_l(jg,8) + particles(n)%speed_z
2703	pts_value_l(jg,9) = pts_value_l(jg,9) + particles(n)%rvar1
2704	pts_value_l(jg,10) = pts_value_l(jg,10) + particles(n)%rvar2
2705	pts_value_l(jg,11) = pts_value_l(jg,11) + particles(n)%rvar3
2706	IF ( particles(n)%speed_z > 0.0_wp ) THEN
2707	pts_value_l(jg,12) = pts_value_l(jg,12) + 1.0_wp
2708	pts_value_l(jg,13) = pts_value_l(jg,13) + particles(n)%speed_z
2709	ELSE
2710	pts_value_l(jg,14) = pts_value_l(jg,14) + particles(n)%speed_z
2711	ENDIF
2712	pts_value_l(jg,15) = pts_value_l(jg,15) + particles(n)%radius
2713	pts_value_l(jg,16) = MIN( pts_value_l(jg,16), particles(n)%radius )
2714	pts_value_l(jg,17) = MAX( pts_value_l(jg,17), particles(n)%radius )
2715	pts_value_l(jg,18) = pts_value_l(jg,18) + 1.0_wp
2716	pts_value_l(jg,19) = pts_value_l(jg,19) + 1.0_wp
2717	ENDIF
2718
2719	ENDIF
2720
2721	ENDDO
2722
2723	ENDDO
2724	ENDDO
2725	ENDDO
2726
2727
2728	#if defined( __parallel )
2729	!
2730	!-- Sum values of the subdomains
2731	inum = number_of_particle_groups + 1
2732
2733	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2734	CALL MPI_ALLREDUCE( pts_value_l(0,1), pts_value(0,1), 15*inum, MPI_REAL, &
2735	MPI_SUM, comm2d, ierr )
2736	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2737	CALL MPI_ALLREDUCE( pts_value_l(0,16), pts_value(0,16), inum, MPI_REAL, &
2738	MPI_MIN, comm2d, ierr )
2739	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2740	CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, &
2741	MPI_MAX, comm2d, ierr )
2742	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2743	CALL MPI_ALLREDUCE( pts_value_l(0,18), pts_value(0,18), inum, MPI_REAL, &
2744	MPI_MAX, comm2d, ierr )
2745	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2746	CALL MPI_ALLREDUCE( pts_value_l(0,19), pts_value(0,19), inum, MPI_REAL, &
2747	MPI_MIN, comm2d, ierr )
2748	#else
2749	pts_value(:,1:19) = pts_value_l(:,1:19)
2750	#endif
2751
2752	!
2753	!-- Normalize the above calculated quantities (except min/max values) with the
2754	!-- total number of particles
2755	IF ( number_of_particle_groups > 1 ) THEN
2756	inum = number_of_particle_groups
2757	ELSE
2758	inum = 0
2759	ENDIF
2760
2761	DO j = 0, inum
2762
2763	IF ( pts_value(j,1) > 0.0_wp ) THEN
2764
2765	pts_value(j,2:15) = pts_value(j,2:15) / pts_value(j,1)
2766	IF ( pts_value(j,12) > 0.0_wp .AND. pts_value(j,12) < 1.0_wp ) THEN
2767	pts_value(j,13) = pts_value(j,13) / pts_value(j,12)
2768	pts_value(j,14) = pts_value(j,14) / ( 1.0_wp - pts_value(j,12) )
2769	ELSEIF ( pts_value(j,12) == 0.0_wp ) THEN
2770	pts_value(j,13) = -1.0_wp
2771	ELSE
2772	pts_value(j,14) = -1.0_wp
2773	ENDIF
2774
2775	ENDIF
2776
2777	ENDDO
2778
2779	!
2780	!-- Calculate higher order moments of particle time series quantities,
2781	!-- seperately for each particle group (if there is more than one group)
2782	DO i = nxl, nxr
2783	DO j = nys, nyn
2784	DO k = nzb, nzt
2785	number_of_particles = prt_count(k,j,i)
2786	IF (number_of_particles <= 0) CYCLE
2787	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
2788	DO n = 1, number_of_particles
2789
2790	pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%x - &
2791	particles(n)%origin_x - pts_value(0,2) )*2 ! x2
2792	pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%y - &
2793	particles(n)%origin_y - pts_value(0,3) )*2 ! y2
2794	pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%z - &
2795	particles(n)%origin_z - pts_value(0,4) )*2 ! z2
2796	pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_x - &
2797	pts_value(0,6) )*2 ! u2
2798	pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%speed_y - &
2799	pts_value(0,7) )*2 ! v2
2800	pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
2801	pts_value(0,8) )*2 ! w2
2802	pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
2803	pts_value(0,9) )**2 ! u"2
2804	pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
2805	pts_value(0,10) )**2 ! v"2
2806	pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
2807	pts_value(0,11) )**2 ! w"2
2808	!
2809	!-- Repeat the same for the respective particle group
2810	IF ( number_of_particle_groups > 1 ) THEN
2811	jg = particles(n)%group
2812
2813	pts_value_l(jg,20) = pts_value_l(jg,20) + ( particles(n)%x - &
2814	particles(n)%origin_x - pts_value(jg,2) )**2
2815	pts_value_l(jg,21) = pts_value_l(jg,21) + ( particles(n)%y - &
2816	particles(n)%origin_y - pts_value(jg,3) )**2
2817	pts_value_l(jg,22) = pts_value_l(jg,22) + ( particles(n)%z - &
2818	particles(n)%origin_z - pts_value(jg,4) )**2
2819	pts_value_l(jg,23) = pts_value_l(jg,23) + ( particles(n)%speed_x - &
2820	pts_value(jg,6) )**2
2821	pts_value_l(jg,24) = pts_value_l(jg,24) + ( particles(n)%speed_y - &
2822	pts_value(jg,7) )**2
2823	pts_value_l(jg,25) = pts_value_l(jg,25) + ( particles(n)%speed_z - &
2824	pts_value(jg,8) )**2
2825	pts_value_l(jg,26) = pts_value_l(jg,26) + ( particles(n)%rvar1 - &
2826	pts_value(jg,9) )**2
2827	pts_value_l(jg,27) = pts_value_l(jg,27) + ( particles(n)%rvar2 - &
2828	pts_value(jg,10) )**2
2829	pts_value_l(jg,28) = pts_value_l(jg,28) + ( particles(n)%rvar3 - &
2830	pts_value(jg,11) )**2
2831	ENDIF
2832
2833	ENDDO
2834	ENDDO
2835	ENDDO
2836	ENDDO
2837
2838	pts_value_l(0,29) = ( number_of_particles - pts_value(0,1) / numprocs )**2
2839	! variance of particle numbers
2840	IF ( number_of_particle_groups > 1 ) THEN
2841	DO j = 1, number_of_particle_groups
2842	pts_value_l(j,29) = ( pts_value_l(j,1) - &
2843	pts_value(j,1) / numprocs )**2
2844	ENDDO
2845	ENDIF
2846
2847	#if defined( __parallel )
2848	!
2849	!-- Sum values of the subdomains
2850	inum = number_of_particle_groups + 1
2851
2852	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2853	CALL MPI_ALLREDUCE( pts_value_l(0,20), pts_value(0,20), inum*10, MPI_REAL, &
2854	MPI_SUM, comm2d, ierr )
2855	#else
2856	pts_value(:,20:29) = pts_value_l(:,20:29)
2857	#endif
2858
2859	!
2860	!-- Normalize the above calculated quantities with the total number of
2861	!-- particles
2862	IF ( number_of_particle_groups > 1 ) THEN
2863	inum = number_of_particle_groups
2864	ELSE
2865	inum = 0
2866	ENDIF
2867
2868	DO j = 0, inum
2869
2870	IF ( pts_value(j,1) > 0.0_wp ) THEN
2871	pts_value(j,20:28) = pts_value(j,20:28) / pts_value(j,1)
2872	ENDIF
2873	pts_value(j,29) = pts_value(j,29) / numprocs
2874
2875	ENDDO
2876
2877	#if defined( __netcdf )
2878	!
2879	!-- Output particle time series quantities in NetCDF format
2880	IF ( myid == 0 ) THEN
2881	DO j = 0, inum
2882	DO i = 1, dopts_num
2883	nc_stat = NF90_PUT_VAR( id_set_pts, id_var_dopts(i,j), &
2884	(/ pts_value(j,i) /), &
2885	start = (/ dopts_time_count /), &
2886	count = (/ 1 /) )
2887	CALL netcdf_handle_error( 'data_output_ptseries', 392 )
2888	ENDDO
2889	ENDDO
2890	ENDIF
2891	#endif
2892
2893	DEALLOCATE( pts_value, pts_value_l )
2894
2895	CALL cpu_log( log_point(36), 'data_output_ptseries', 'stop' )
2896
2897	END SUBROUTINE lpm_data_output_ptseries
2898
2899
2900	!------------------------------------------------------------------------------!
2901	! Description:
2902	! ------------
2903	!> This routine reads the respective restart data for the lpm.
2904	!------------------------------------------------------------------------------!
2905	SUBROUTINE lpm_rrd_local_particles
2906
2907	CHARACTER (LEN=10) :: particle_binary_version !<
2908	CHARACTER (LEN=10) :: version_on_file !<
2909
2910	INTEGER(iwp) :: alloc_size !<
2911	INTEGER(iwp) :: ip !<
2912	INTEGER(iwp) :: jp !<
2913	INTEGER(iwp) :: kp !<
2914
2915	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles !<
2916
2917	!
2918	!-- Read particle data from previous model run.
2919	!-- First open the input unit.
2920	IF ( myid_char == '' ) THEN
2921	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN'//myid_char, &
2922	FORM='UNFORMATTED' )
2923	ELSE
2924	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN/'//myid_char, &
2925	FORM='UNFORMATTED' )
2926	ENDIF
2927
2928	!
2929	!-- First compare the version numbers
2930	READ ( 90 ) version_on_file
2931	particle_binary_version = '4.0'
2932	IF ( TRIM( version_on_file ) /= TRIM( particle_binary_version ) ) THEN
2933	message_string = 'version mismatch concerning data from prior ' // &
2934	'run &version on file = "' // &
2935	TRIM( version_on_file ) // &
2936	'&version in program = "' // &
2937	TRIM( particle_binary_version ) // '"'
2938	CALL message( 'lpm_read_restart_file', 'PA0214', 1, 2, 0, 6, 0 )
2939	ENDIF
2940
2941	!
2942	!-- If less particles are stored on the restart file than prescribed by
2943	!-- 1, the remainder is initialized by zero_particle to avoid
2944	!-- errors.
2945	zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
2946	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
2947	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
2948	0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &
2949	0, 0, 0_idp, .FALSE., -1 )
2950	!
2951	!-- Read some particle parameters and the size of the particle arrays,
2952	!-- allocate them and read their contents.
2953	READ ( 90 ) bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, &
2954	last_particle_release_time, number_of_particle_groups, &
2955	particle_groups, time_write_particle_data
2956
2957	ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
2958	grid_particles(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2959
2960	READ ( 90 ) prt_count
2961
2962	DO ip = nxl, nxr
2963	DO jp = nys, nyn
2964	DO kp = nzb+1, nzt
2965
2966	number_of_particles = prt_count(kp,jp,ip)
2967	IF ( number_of_particles > 0 ) THEN
2968	alloc_size = MAX( INT( number_of_particles * &
2969	( 1.0_wp + alloc_factor / 100.0_wp ) ), &
2970	1 )
2971	ELSE
2972	alloc_size = 1
2973	ENDIF
2974
2975	ALLOCATE( grid_particles(kp,jp,ip)%particles(1:alloc_size) )
2976
2977	IF ( number_of_particles > 0 ) THEN
2978	ALLOCATE( tmp_particles(1:number_of_particles) )
2979	READ ( 90 ) tmp_particles
2980	grid_particles(kp,jp,ip)%particles(1:number_of_particles) = tmp_particles
2981	DEALLOCATE( tmp_particles )
2982	IF ( number_of_particles < alloc_size ) THEN
2983	grid_particles(kp,jp,ip)%particles(number_of_particles+1:alloc_size) &
2984	= zero_particle
2985	ENDIF
2986	ELSE
2987	grid_particles(kp,jp,ip)%particles(1:alloc_size) = zero_particle
2988	ENDIF
2989
2990	ENDDO
2991	ENDDO
2992	ENDDO
2993
2994	CLOSE ( 90 )
2995	!
2996	!-- Must be called to sort particles into blocks, which is needed for a fast
2997	!-- interpolation of the LES fields on the particle position.
2998	CALL lpm_sort_and_delete
2999
3000
3001	END SUBROUTINE lpm_rrd_local_particles
3002
3003
3004	SUBROUTINE lpm_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
3005	nxr_on_file, nynf, nync, nyn_on_file, nysf, &
3006	nysc, nys_on_file, tmp_3d, found )
3007
3008
3009	USE control_parameters, &
3010	ONLY: length, restart_string
3011
3012	INTEGER(iwp) :: k !<
3013	INTEGER(iwp) :: nxlc !<
3014	INTEGER(iwp) :: nxlf !<
3015	INTEGER(iwp) :: nxl_on_file !<
3016	INTEGER(iwp) :: nxrc !<
3017	INTEGER(iwp) :: nxrf !<
3018	INTEGER(iwp) :: nxr_on_file !<
3019	INTEGER(iwp) :: nync !<
3020	INTEGER(iwp) :: nynf !<
3021	INTEGER(iwp) :: nyn_on_file !<
3022	INTEGER(iwp) :: nysc !<
3023	INTEGER(iwp) :: nysf !<
3024	INTEGER(iwp) :: nys_on_file !<
3025
3026	LOGICAL, INTENT(OUT) :: found
3027
3028	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
3029
3030
3031	found = .TRUE.
3032
3033	SELECT CASE ( restart_string(1:length) )
3034
3035	CASE ( 'iran' ) ! matching random numbers is still unresolved issue
3036	IF ( k == 1 ) READ ( 13 ) iran, iran_part
3037
3038	CASE ( 'pc_av' )
3039	IF ( .NOT. ALLOCATED( pc_av ) ) THEN
3040	ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3041	ENDIF
3042	IF ( k == 1 ) READ ( 13 ) tmp_3d
3043	pc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3044	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3045
3046	CASE ( 'pr_av' )
3047	IF ( .NOT. ALLOCATED( pr_av ) ) THEN
3048	ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3049	ENDIF
3050	IF ( k == 1 ) READ ( 13 ) tmp_3d
3051	pr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3052	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3053
3054	CASE ( 'ql_c_av' )
3055	IF ( .NOT. ALLOCATED( ql_c_av ) ) THEN
3056	ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3057	ENDIF
3058	IF ( k == 1 ) READ ( 13 ) tmp_3d
3059	ql_c_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3060	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3061
3062	CASE ( 'ql_v_av' )
3063	IF ( .NOT. ALLOCATED( ql_v_av ) ) THEN
3064	ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3065	ENDIF
3066	IF ( k == 1 ) READ ( 13 ) tmp_3d
3067	ql_v_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3068	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3069
3070	CASE ( 'ql_vp_av' )
3071	IF ( .NOT. ALLOCATED( ql_vp_av ) ) THEN
3072	ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3073	ENDIF
3074	IF ( k == 1 ) READ ( 13 ) tmp_3d
3075	ql_vp_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3076	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3077
3078	CASE DEFAULT
3079
3080	found = .FALSE.
3081
3082	END SELECT
3083
3084
3085	END SUBROUTINE lpm_rrd_local
3086
3087	!------------------------------------------------------------------------------!
3088	! Description:
3089	! ------------
3090	!> This routine writes the respective restart data for the lpm.
3091	!------------------------------------------------------------------------------!
3092	SUBROUTINE lpm_wrd_local
3093
3094	CHARACTER (LEN=10) :: particle_binary_version !<
3095
3096	INTEGER(iwp) :: ip !<
3097	INTEGER(iwp) :: jp !<
3098	INTEGER(iwp) :: kp !<
3099	!
3100	!-- First open the output unit.
3101	IF ( myid_char == '' ) THEN
3102	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, &
3103	FORM='UNFORMATTED')
3104	ELSE
3105	IF ( myid == 0 ) CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' )
3106	#if defined( __parallel )
3107	!
3108	!-- Set a barrier in order to allow that thereafter all other processors
3109	!-- in the directory created by PE0 can open their file
3110	CALL MPI_BARRIER( comm2d, ierr )
3111	#endif
3112	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, &
3113	FORM='UNFORMATTED' )
3114	ENDIF
3115
3116	!
3117	!-- Write the version number of the binary format.
3118	!-- Attention: After changes to the following output commands the version
3119	!-- --------- number of the variable particle_binary_version must be
3120	!-- changed! Also, the version number and the list of arrays
3121	!-- to be read in lpm_read_restart_file must be adjusted
3122	!-- accordingly.
3123	particle_binary_version = '4.0'
3124	WRITE ( 90 ) particle_binary_version
3125
3126	!
3127	!-- Write some particle parameters, the size of the particle arrays
3128	WRITE ( 90 ) bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, &
3129	last_particle_release_time, number_of_particle_groups, &
3130	particle_groups, time_write_particle_data
3131
3132	WRITE ( 90 ) prt_count
3133
3134	DO ip = nxl, nxr
3135	DO jp = nys, nyn
3136	DO kp = nzb+1, nzt
3137	number_of_particles = prt_count(kp,jp,ip)
3138	particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
3139	IF ( number_of_particles <= 0 ) CYCLE
3140	WRITE ( 90 ) particles
3141	ENDDO
3142	ENDDO
3143	ENDDO
3144
3145	CLOSE ( 90 )
3146
3147	#if defined( __parallel )
3148	CALL MPI_BARRIER( comm2d, ierr )
3149	#endif
3150
3151	CALL wrd_write_string( 'iran' )
3152	WRITE ( 14 ) iran, iran_part
3153
3154
3155	END SUBROUTINE lpm_wrd_local
3156
3157
3158	!------------------------------------------------------------------------------!
3159	! Description:
3160	! ------------
3161	!> This routine writes the respective restart data for the lpm.
3162	!------------------------------------------------------------------------------!
3163	SUBROUTINE lpm_wrd_global
3164
3165	CALL wrd_write_string( 'curvature_solution_effects' )
3166	WRITE ( 14 ) curvature_solution_effects
3167
3168	CALL wrd_write_string( 'interpolation_simple_corrector' )
3169	WRITE ( 14 ) interpolation_simple_corrector
3170
3171	CALL wrd_write_string( 'interpolation_simple_predictor' )
3172	WRITE ( 14 ) interpolation_simple_predictor
3173
3174	CALL wrd_write_string( 'interpolation_trilinear' )
3175	WRITE ( 14 ) interpolation_trilinear
3176
3177	END SUBROUTINE lpm_wrd_global
3178
3179
3180	!------------------------------------------------------------------------------!
3181	! Description:
3182	! ------------
3183	!> This routine writes the respective restart data for the lpm.
3184	!------------------------------------------------------------------------------!
3185	SUBROUTINE lpm_rrd_global( found )
3186
3187	USE control_parameters, &
3188	ONLY: length, restart_string
3189
3190	LOGICAL, INTENT(OUT) :: found
3191
3192	found = .TRUE.
3193
3194	SELECT CASE ( restart_string(1:length) )
3195
3196	CASE ( 'curvature_solution_effects' )
3197	READ ( 13 ) curvature_solution_effects
3198
3199	CASE ( 'interpolation_simple_corrector' )
3200	READ ( 13 ) interpolation_simple_corrector
3201
3202	CASE ( 'interpolation_simple_predictor' )
3203	READ ( 13 ) interpolation_simple_predictor
3204
3205	CASE ( 'interpolation_trilinear' )
3206	READ ( 13 ) interpolation_trilinear
3207
3208	! CASE ( 'global_paramter' )
3209	! READ ( 13 ) global_parameter
3210	! CASE ( 'global_array' )
3211	! IF ( .NOT. ALLOCATED( global_array ) ) ALLOCATE( global_array(1:10) )
3212	! READ ( 13 ) global_array
3213
3214	CASE DEFAULT
3215
3216	found = .FALSE.
3217
3218	END SELECT
3219
3220	END SUBROUTINE lpm_rrd_global
3221
3222
3223	!------------------------------------------------------------------------------!
3224	! Description:
3225	! ------------
3226	!> This is a submodule of the lagrangian particle model. It contains all
3227	!> dynamic processes of the lpm. This includes the advection (resolved and sub-
3228	!> grid scale) as well as the boundary conditions of particles. As a next step
3229	!> this submodule should be excluded as an own file.
3230	!------------------------------------------------------------------------------!
3231	SUBROUTINE lpm_advec (ip,jp,kp)
3232
3233	LOGICAL :: subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
3234
3235	INTEGER(iwp) :: i !< index variable along x
3236	INTEGER(iwp) :: i_next !< index variable along x
3237	INTEGER(iwp) :: ip !< index variable along x
3238	INTEGER(iwp) :: iteration_steps = 1 !< amount of iterations steps for corrector step
3239	INTEGER(iwp) :: j !< index variable along y
3240	INTEGER(iwp) :: j_next !< index variable along y
3241	INTEGER(iwp) :: jp !< index variable along y
3242	INTEGER(iwp) :: k !< index variable along z
3243	INTEGER(iwp) :: k_wall !< vertical index of topography top
3244	INTEGER(iwp) :: kp !< index variable along z
3245	INTEGER(iwp) :: k_next !< index variable along z
3246	INTEGER(iwp) :: kw !< index variable along z
3247	INTEGER(iwp) :: kkw !< index variable along z
3248	INTEGER(iwp) :: n !< loop variable over all particles in a grid box
3249	INTEGER(iwp) :: nb !< block number particles are sorted in
3250	INTEGER(iwp) :: particle_end !< end index for partilce loop
3251	INTEGER(iwp) :: particle_start !< start index for particle loop
3252	INTEGER(iwp) :: surf_start !< Index on surface data-type for current grid box
3253	INTEGER(iwp) :: subbox_end !< end index for loop over subboxes in particle advection
3254	INTEGER(iwp) :: subbox_start !< start index for loop over subboxes in particle advection
3255	INTEGER(iwp) :: nn !< loop variable over iterations steps
3256
3257	INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current block
3258	INTEGER(iwp), DIMENSION(0:7) :: end_index !< start particle index for current block
3259
3260	REAL(wp) :: aa !< dummy argument for horizontal particle interpolation
3261	REAL(wp) :: alpha !< interpolation facor for x-direction
3262
3263	REAL(wp) :: bb !< dummy argument for horizontal particle interpolation
3264	REAL(wp) :: beta !< interpolation facor for y-direction
3265	REAL(wp) :: cc !< dummy argument for horizontal particle interpolation
3266	REAL(wp) :: d_z_p_z0 !< inverse of interpolation length for logarithmic interpolation
3267	REAL(wp) :: dd !< dummy argument for horizontal particle interpolation
3268	REAL(wp) :: de_dx_int_l !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
3269	REAL(wp) :: de_dx_int_u !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
3270	REAL(wp) :: de_dy_int_l !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
3271	REAL(wp) :: de_dy_int_u !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
3272	REAL(wp) :: de_dt !< temporal derivative of TKE experienced by the particle
3273	REAL(wp) :: de_dt_min !< lower level for temporal TKE derivative
3274	REAL(wp) :: de_dz_int_l !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
3275	REAL(wp) :: de_dz_int_u !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
3276	REAL(wp) :: diameter !< diamter of droplet
3277	REAL(wp) :: diss_int_l !< x/y-interpolated dissipation at particle position at lower vertical level
3278	REAL(wp) :: diss_int_u !< x/y-interpolated dissipation at particle position at upper vertical level
3279	REAL(wp) :: dt_particle_m !< previous particle time step
3280	REAL(wp) :: dz_temp !< dummy for the vertical grid spacing
3281	REAL(wp) :: e_int_l !< x/y-interpolated TKE at particle position at lower vertical level
3282	REAL(wp) :: e_int_u !< x/y-interpolated TKE at particle position at upper vertical level
3283	REAL(wp) :: e_mean_int !< horizontal mean TKE at particle height
3284	REAL(wp) :: exp_arg !< argument in the exponent - particle radius
3285	REAL(wp) :: exp_term !< exponent term
3286	REAL(wp) :: gamma !< interpolation facor for z-direction
3287	REAL(wp) :: gg !< dummy argument for horizontal particle interpolation
3288	REAL(wp) :: height_p !< dummy argument for logarithmic interpolation
3289	REAL(wp) :: log_z_z0_int !< logarithmus used for surface_layer interpolation
3290	REAL(wp) :: random_gauss !< Gaussian-distributed random number used for SGS particle advection
3291	REAL(wp) :: RL !< Lagrangian autocorrelation coefficient
3292	REAL(wp) :: rg1 !< Gaussian distributed random number
3293	REAL(wp) :: rg2 !< Gaussian distributed random number
3294	REAL(wp) :: rg3 !< Gaussian distributed random number
3295	REAL(wp) :: sigma !< velocity standard deviation
3296	REAL(wp) :: u_int_l !< x/y-interpolated u-component at particle position at lower vertical level
3297	REAL(wp) :: u_int_u !< x/y-interpolated u-component at particle position at upper vertical level
3298	REAL(wp) :: unext !< calculated particle u-velocity of corrector step
3299	REAL(wp) :: us_int !< friction velocity at particle grid box
3300	REAL(wp) :: usws_int !< surface momentum flux (u component) at particle grid box
3301	REAL(wp) :: v_int_l !< x/y-interpolated v-component at particle position at lower vertical level
3302	REAL(wp) :: v_int_u !< x/y-interpolated v-component at particle position at upper vertical level
3303	REAL(wp) :: vsws_int !< surface momentum flux (u component) at particle grid box
3304	REAL(wp) :: vnext !< calculated particle v-velocity of corrector step
3305	REAL(wp) :: vv_int !< dummy to compute interpolated mean SGS TKE, used to scale SGS advection
3306	REAL(wp) :: w_int_l !< x/y-interpolated w-component at particle position at lower vertical level
3307	REAL(wp) :: w_int_u !< x/y-interpolated w-component at particle position at upper vertical level
3308	REAL(wp) :: wnext !< calculated particle w-velocity of corrector step
3309	REAL(wp) :: w_s !< terminal velocity of droplets
3310	REAL(wp) :: x !< dummy argument for horizontal particle interpolation
3311	REAL(wp) :: xp !< calculated particle position in x of predictor step
3312	REAL(wp) :: y !< dummy argument for horizontal particle interpolation
3313	REAL(wp) :: yp !< calculated particle position in y of predictor step
3314	REAL(wp) :: z_p !< surface layer height (0.5 dz)
3315	REAL(wp) :: zp !< calculated particle position in z of predictor step
3316
3317	REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter for fall velocity
3318	REAL(wp), PARAMETER :: b_rog = 10.43_wp !< parameter for fall velocity
3319	REAL(wp), PARAMETER :: c_rog = 0.6_wp !< parameter for fall velocity
3320	REAL(wp), PARAMETER :: k_cap_rog = 4.0_wp !< parameter for fall velocity
3321	REAL(wp), PARAMETER :: k_low_rog = 12.0_wp !< parameter for fall velocity
3322	REAL(wp), PARAMETER :: d0_rog = 0.745_wp !< separation diameter
3323
3324	REAL(wp), DIMENSION(number_of_particles) :: term_1_2 !< flag to communicate whether a particle is near topography or not
3325	REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !< ratio between the density of the fluid and the density of the particles
3326	REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< horizontal TKE gradient along x at particle position
3327	REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< horizontal TKE gradient along y at particle position
3328	REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< horizontal TKE gradient along z at particle position
3329	REAL(wp), DIMENSION(number_of_particles) :: diss_int !< dissipation at particle position
3330	REAL(wp), DIMENSION(number_of_particles) :: dt_gap !< remaining time until particle time integration reaches LES time
3331	REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< particle time step
3332	REAL(wp), DIMENSION(number_of_particles) :: e_int !< TKE at particle position
3333	REAL(wp), DIMENSION(number_of_particles) :: fs_int !< weighting factor for subgrid-scale particle speed
3334	REAL(wp), DIMENSION(number_of_particles) :: lagr_timescale !< Lagrangian timescale
3335	REAL(wp), DIMENSION(number_of_particles) :: rvar1_temp !< SGS particle velocity - u-component
3336	REAL(wp), DIMENSION(number_of_particles) :: rvar2_temp !< SGS particle velocity - v-component
3337	REAL(wp), DIMENSION(number_of_particles) :: rvar3_temp !< SGS particle velocity - w-component
3338	REAL(wp), DIMENSION(number_of_particles) :: u_int !< u-component of particle speed
3339	REAL(wp), DIMENSION(number_of_particles) :: v_int !< v-component of particle speed
3340	REAL(wp), DIMENSION(number_of_particles) :: w_int !< w-component of particle speed
3341	REAL(wp), DIMENSION(number_of_particles) :: xv !< x-position
3342	REAL(wp), DIMENSION(number_of_particles) :: yv !< y-position
3343	REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position
3344
3345	REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers
3346
3347	CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
3348	!
3349	!-- Determine height of Prandtl layer and distance between Prandtl-layer
3350	!-- height and horizontal mean roughness height, which are required for
3351	!-- vertical logarithmic interpolation of horizontal particle speeds
3352	!-- (for particles below first vertical grid level).
3353	z_p = zu(nzb+1) - zw(nzb)
3354	d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
3355
3356	xv = particles(1:number_of_particles)%x
3357	yv = particles(1:number_of_particles)%y
3358	zv = particles(1:number_of_particles)%z
3359	dt_particle = dt_3d
3360
3361	!
3362	!-- This case uses a simple interpolation method for the particle velocites,
3363	!-- and applying a predictor-corrector method. @note the current time divergence
3364	!-- free time step is denoted with u_t etc.; the velocities of the time level of
3365	!-- t+1 wit u,v, and w, as the model is called after swap timelevel
3366	!-- @attention: for the corrector step the velocities of t(n+1) are required.
3367	!-- Therefore the particle code is executed at the end of the time intermediate
3368	!-- timestep routine. This interpolation method is described in more detail
3369	!-- in Grabowski et al., 2018 (GMD).
3370	IF ( interpolation_simple_corrector ) THEN
3371	!
3372	!-- Predictor step
3373	kkw = kp - 1
3374	DO n = 1, number_of_particles
3375
3376	alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
3377	u_int(n) = u_t(kp,jp,ip) * ( 1.0_wp - alpha ) + u_t(kp,jp,ip+1) * alpha
3378
3379	beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
3380	v_int(n) = v_t(kp,jp,ip) * ( 1.0_wp - beta ) + v_t(kp,jp+1,ip) * beta
3381
3382	gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / &
3383	( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
3384	w_int(n) = w_t(kkw,jp,ip) * ( 1.0_wp - gamma ) + w_t(kkw+1,jp,ip) * gamma
3385
3386	ENDDO
3387	!
3388	!-- Corrector step
3389	DO n = 1, number_of_particles
3390
3391	IF ( .NOT. particles(n)%particle_mask ) CYCLE
3392
3393	DO nn = 1, iteration_steps
3394
3395	!
3396	!-- Guess new position
3397	xp = particles(n)%x + u_int(n) * dt_particle(n)
3398	yp = particles(n)%y + v_int(n) * dt_particle(n)
3399	zp = particles(n)%z + w_int(n) * dt_particle(n)
3400	!
3401	!-- x direction
3402	i_next = FLOOR( xp * ddx , KIND=iwp)
3403	alpha = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp )
3404	!
3405	!-- y direction
3406	j_next = FLOOR( yp * ddy )
3407	beta = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp )
3408	!
3409	!-- z_direction
3410	k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) + offset_ocean_nzt ), nzt ), 0)
3411	gamma = MAX( MIN( ( zp - zw(k_next) ) / &
3412	( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp )
3413	!
3414	!-- Calculate part of the corrector step
3415	unext = u(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) + &
3416	u(k_next+1, j_next, i_next+1) * alpha
3417
3418	vnext = v(k_next+1, j_next, i_next) * ( 1.0_wp - beta ) + &
3419	v(k_next+1, j_next+1, i_next ) * beta
3420
3421	wnext = w(k_next, j_next, i_next) * ( 1.0_wp - gamma ) + &
3422	w(k_next+1, j_next, i_next ) * gamma
3423
3424	!
3425	!-- Calculate interpolated particle velocity with predictor
3426	!-- corrector step. u_int, v_int and w_int describes the part of
3427	!-- the predictor step. unext, vnext and wnext is the part of the
3428	!-- corrector step. The resulting new position is set below. The
3429	!-- implementation is based on Grabowski et al., 2018 (GMD).
3430	u_int(n) = 0.5_wp * ( u_int(n) + unext )
3431	v_int(n) = 0.5_wp * ( v_int(n) + vnext )
3432	w_int(n) = 0.5_wp * ( w_int(n) + wnext )
3433
3434	ENDDO
3435	ENDDO
3436	!
3437	!-- This case uses a simple interpolation method for the particle velocites,
3438	!-- and applying a predictor.
3439	ELSEIF ( interpolation_simple_predictor ) THEN
3440	!
3441	!-- The particle position for the w velociy is based on the value of kp and kp-1
3442	kkw = kp - 1
3443	DO n = 1, number_of_particles
3444	IF ( .NOT. particles(n)%particle_mask ) CYCLE
3445
3446	alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
3447	u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
3448
3449	beta = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
3450	v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
3451
3452	gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) / &
3453	( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
3454	w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
3455	ENDDO
3456	!
3457	!-- The trilinear interpolation.
3458	ELSEIF ( interpolation_trilinear ) THEN
3459
3460	start_index = grid_particles(kp,jp,ip)%start_index
3461	end_index = grid_particles(kp,jp,ip)%end_index
3462
3463	DO nb = 0, 7
3464	!
3465	!-- Interpolate u velocity-component
3466	i = ip
3467	j = jp + block_offset(nb)%j_off
3468	k = kp + block_offset(nb)%k_off
3469
3470	DO n = start_index(nb), end_index(nb)
3471	!
3472	!-- Interpolation of the u velocity component onto particle position.
3473	!-- Particles are interpolation bi-linearly in the horizontal and a
3474	!-- linearly in the vertical. An exception is made for particles below
3475	!-- the first vertical grid level in case of a prandtl layer. In this
3476	!-- case the horizontal particle velocity components are determined using
3477	!-- Monin-Obukhov relations (if branch).
3478	!-- First, check if particle is located below first vertical grid level
3479	!-- above topography (Prandtl-layer height)
3480	!-- Determine vertical index of topography top
3481	k_wall = topo_top_ind(jp,ip,0)
3482
3483	IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN
3484	!
3485	!-- Resolved-scale horizontal particle velocity is zero below z0.
3486	IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN
3487	u_int(n) = 0.0_wp
3488	ELSE
3489	!
3490	!-- Determine the sublayer. Further used as index.
3491	height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
3492	* REAL( number_of_sublayers, KIND=wp ) &
3493	* d_z_p_z0
3494	!
3495	!-- Calculate LOG(z/z0) for exact particle height. Therefore,
3496	!-- interpolate linearly between precalculated logarithm.
3497	log_z_z0_int = log_z_z0(INT(height_p)) &
3498	+ ( height_p - INT(height_p) ) &
3499	* ( log_z_z0(INT(height_p)+1) &
3500	- log_z_z0(INT(height_p)) &
3501	)
3502	!
3503	!-- Get friction velocity and momentum flux from new surface data
3504	!-- types.
3505	IF ( surf_def_h(0)%start_index(jp,ip) <= &
3506	surf_def_h(0)%end_index(jp,ip) ) THEN
3507	surf_start = surf_def_h(0)%start_index(jp,ip)
3508	!-- Limit friction velocity. In narrow canyons or holes the
3509	!-- friction velocity can become very small, resulting in a too
3510	!-- large particle speed.
3511	us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
3512	usws_int = surf_def_h(0)%usws(surf_start) &
3513	* drho_air_zw(k_wall)
3514	ELSEIF ( surf_lsm_h%start_index(jp,ip) <= &
3515	surf_lsm_h%end_index(jp,ip) ) THEN
3516	surf_start = surf_lsm_h%start_index(jp,ip)
3517	us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
3518	usws_int = surf_lsm_h%usws(surf_start) &
3519	* drho_air_zw(k_wall)
3520	ELSEIF ( surf_usm_h%start_index(jp,ip) <= &
3521	surf_usm_h%end_index(jp,ip) ) THEN
3522	surf_start = surf_usm_h%start_index(jp,ip)
3523	us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp )
3524	usws_int = surf_usm_h%usws(surf_start) &
3525	* drho_air_zw(k_wall)
3526	ENDIF
3527	!
3528	!-- Neutral solution is applied for all situations, e.g. also for
3529	!-- unstable and stable situations. Even though this is not exact
3530	!-- this saves a lot of CPU time since several calls of intrinsic
3531	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
3532	!-- as sensitivity studies revealed no significant effect of
3533	!-- using the neutral solution also for un/stable situations.
3534	u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) &
3535	* log_z_z0_int - u_gtrans
3536	ENDIF
3537	!
3538	!-- Particle above the first grid level. Bi-linear interpolation in the
3539	!-- horizontal and linear interpolation in the vertical direction.
3540	ELSE
3541	x = xv(n) - i * dx
3542	y = yv(n) + ( 0.5_wp - j ) * dy
3543	aa = x2 + y2
3544	bb = ( dx - x )2 + y2
3545	cc = x2 + ( dy - y )2
3546	dd = ( dx - x )2 + ( dy - y )2
3547	gg = aa + bb + cc + dd
3548
3549	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
3550	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * &
3551	u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
3552
3553	IF ( k == nzt ) THEN
3554	u_int(n) = u_int_l
3555	ELSE
3556	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
3557	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * &
3558	u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
3559	u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * &
3560	( u_int_u - u_int_l )
3561	ENDIF
3562	ENDIF
3563	ENDDO
3564	!
3565	!-- Same procedure for interpolation of the v velocity-component
3566	i = ip + block_offset(nb)%i_off
3567	j = jp
3568	k = kp + block_offset(nb)%k_off
3569
3570	DO n = start_index(nb), end_index(nb)
3571	!
3572	!-- Determine vertical index of topography top
3573	k_wall = topo_top_ind(jp,ip,0)
3574
3575	IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN
3576	IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN
3577	!
3578	!-- Resolved-scale horizontal particle velocity is zero below z0.
3579	v_int(n) = 0.0_wp
3580	ELSE
3581	!
3582	!-- Determine the sublayer. Further used as index. Please note,
3583	!-- logarithmus can not be reused from above, as in in case of
3584	!-- topography particle on u-grid can be above surface-layer height,
3585	!-- whereas it can be below on v-grid.
3586	height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
3587	* REAL( number_of_sublayers, KIND=wp ) &
3588	* d_z_p_z0
3589	!
3590	!-- Calculate LOG(z/z0) for exact particle height. Therefore,
3591	!-- interpolate linearly between precalculated logarithm.
3592	log_z_z0_int = log_z_z0(INT(height_p)) &
3593	+ ( height_p - INT(height_p) ) &
3594	* ( log_z_z0(INT(height_p)+1) &
3595	- log_z_z0(INT(height_p)) &
3596	)
3597	!
3598	!-- Get friction velocity and momentum flux from new surface data
3599	!-- types.
3600	IF ( surf_def_h(0)%start_index(jp,ip) <= &
3601	surf_def_h(0)%end_index(jp,ip) ) THEN
3602	surf_start = surf_def_h(0)%start_index(jp,ip)
3603	!-- Limit friction velocity. In narrow canyons or holes the
3604	!-- friction velocity can become very small, resulting in a too
3605	!-- large particle speed.
3606	us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
3607	vsws_int = surf_def_h(0)%vsws(surf_start) &
3608	* drho_air_zw(k_wall)
3609	ELSEIF ( surf_lsm_h%start_index(jp,ip) <= &
3610	surf_lsm_h%end_index(jp,ip) ) THEN
3611	surf_start = surf_lsm_h%start_index(jp,ip)
3612	us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
3613	vsws_int = surf_lsm_h%vsws(surf_start) &
3614	* drho_air_zw(k_wall)
3615	ELSEIF ( surf_usm_h%start_index(jp,ip) <= &
3616	surf_usm_h%end_index(jp,ip) ) THEN
3617	surf_start = surf_usm_h%start_index(jp,ip)
3618	us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp )
3619	vsws_int = surf_usm_h%vsws(surf_start) &
3620	* drho_air_zw(k_wall)
3621	ENDIF
3622	!
3623	!-- Neutral solution is applied for all situations, e.g. also for
3624	!-- unstable and stable situations. Even though this is not exact
3625	!-- this saves a lot of CPU time since several calls of intrinsic
3626	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
3627	!-- as sensitivity studies revealed no significant effect of
3628	!-- using the neutral solution also for un/stable situations.
3629	v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) &
3630	* log_z_z0_int - v_gtrans
3631
3632	ENDIF
3633	ELSE
3634	x = xv(n) + ( 0.5_wp - i ) * dx
3635	y = yv(n) - j * dy
3636	aa = x2 + y2
3637	bb = ( dx - x )2 + y2
3638	cc = x2 + ( dy - y )2
3639	dd = ( dx - x )2 + ( dy - y )2
3640	gg = aa + bb + cc + dd
3641
3642	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
3643	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
3644	) / ( 3.0_wp * gg ) - v_gtrans
3645
3646	IF ( k == nzt ) THEN
3647	v_int(n) = v_int_l
3648	ELSE
3649	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
3650	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
3651	) / ( 3.0_wp * gg ) - v_gtrans
3652	v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * &
3653	( v_int_u - v_int_l )
3654	ENDIF
3655	ENDIF
3656	ENDDO
3657	!
3658	!-- Same procedure for interpolation of the w velocity-component
3659	i = ip + block_offset(nb)%i_off
3660	j = jp + block_offset(nb)%j_off
3661	k = kp - 1
3662
3663	DO n = start_index(nb), end_index(nb)
3664	IF ( vertical_particle_advection(particles(n)%group) ) THEN
3665	x = xv(n) + ( 0.5_wp - i ) * dx
3666	y = yv(n) + ( 0.5_wp - j ) * dy
3667	aa = x2 + y2
3668	bb = ( dx - x )2 + y2
3669	cc = x2 + ( dy - y )2
3670	dd = ( dx - x )2 + ( dy - y )2
3671	gg = aa + bb + cc + dd
3672
3673	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
3674	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
3675	) / ( 3.0_wp * gg )
3676
3677	IF ( k == nzt ) THEN
3678	w_int(n) = w_int_l
3679	ELSE
3680	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
3681	( gg-bb ) * w(k+1,j,i+1) + &
3682	( gg-cc ) * w(k+1,j+1,i) + &
3683	( gg-dd ) * w(k+1,j+1,i+1) &
3684	) / ( 3.0_wp * gg )
3685	w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) * &
3686	( w_int_u - w_int_l )
3687	ENDIF
3688	ELSE
3689	w_int(n) = 0.0_wp
3690	ENDIF
3691	ENDDO
3692	ENDDO
3693	ENDIF
3694
3695	!-- Interpolate and calculate quantities needed for calculating the SGS
3696	!-- velocities
3697	IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN
3698
3699	DO nb = 0,7
3700
3701	subbox_at_wall = .FALSE.
3702	!
3703	!-- In case of topography check if subbox is adjacent to a wall
3704	IF ( .NOT. topography == 'flat' ) THEN
3705	i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) )
3706	j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) )
3707	k = kp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 0 ) )
3708	IF ( .NOT. BTEST(wall_flags_total_0(k, jp, ip), 0) .OR. &
3709	.NOT. BTEST(wall_flags_total_0(kp, j, ip), 0) .OR. &
3710	.NOT. BTEST(wall_flags_total_0(kp, jp, i ), 0) ) &
3711	THEN
3712	subbox_at_wall = .TRUE.
3713	ENDIF
3714	ENDIF
3715	IF ( subbox_at_wall ) THEN
3716	e_int(start_index(nb):end_index(nb)) = e(kp,jp,ip)
3717	diss_int(start_index(nb):end_index(nb)) = diss(kp,jp,ip)
3718	de_dx_int(start_index(nb):end_index(nb)) = de_dx(kp,jp,ip)
3719	de_dy_int(start_index(nb):end_index(nb)) = de_dy(kp,jp,ip)
3720	de_dz_int(start_index(nb):end_index(nb)) = de_dz(kp,jp,ip)
3721	!
3722	!-- Set flag for stochastic equation.
3723	term_1_2(start_index(nb):end_index(nb)) = 0.0_wp
3724	ELSE
3725	i = ip + block_offset(nb)%i_off
3726	j = jp + block_offset(nb)%j_off
3727	k = kp + block_offset(nb)%k_off
3728
3729	DO n = start_index(nb), end_index(nb)
3730	!
3731	!-- Interpolate TKE
3732	x = xv(n) + ( 0.5_wp - i ) * dx
3733	y = yv(n) + ( 0.5_wp - j ) * dy
3734	aa = x2 + y2
3735	bb = ( dx - x )2 + y2
3736	cc = x2 + ( dy - y )2
3737	dd = ( dx - x )2 + ( dy - y )2
3738	gg = aa + bb + cc + dd
3739
3740	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
3741	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
3742	) / ( 3.0_wp * gg )
3743
3744	IF ( k+1 == nzt+1 ) THEN
3745	e_int(n) = e_int_l
3746	ELSE
3747	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
3748	( gg - bb ) * e(k+1,j,i+1) + &
3749	( gg - cc ) * e(k+1,j+1,i) + &
3750	( gg - dd ) * e(k+1,j+1,i+1) &
3751	) / ( 3.0_wp * gg )
3752	e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * &
3753	( e_int_u - e_int_l )
3754	ENDIF
3755	!
3756	!-- Needed to avoid NaN particle velocities (this might not be
3757	!-- required any more)
3758	IF ( e_int(n) <= 0.0_wp ) THEN
3759	e_int(n) = 1.0E-20_wp
3760	ENDIF
3761	!
3762	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
3763	!-- all position variables from above (TKE))
3764	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
3765	( gg - bb ) * de_dx(k,j,i+1) + &
3766	( gg - cc ) * de_dx(k,j+1,i) + &
3767	( gg - dd ) * de_dx(k,j+1,i+1) &
3768	) / ( 3.0_wp * gg )
3769
3770	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
3771	de_dx_int(n) = de_dx_int_l
3772	ELSE
3773	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
3774	( gg - bb ) * de_dx(k+1,j,i+1) + &
3775	( gg - cc ) * de_dx(k+1,j+1,i) + &
3776	( gg - dd ) * de_dx(k+1,j+1,i+1) &
3777	) / ( 3.0_wp * gg )
3778	de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * &
3779	( de_dx_int_u - de_dx_int_l )
3780	ENDIF
3781	!
3782	!-- Interpolate the TKE gradient along y
3783	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
3784	( gg - bb ) * de_dy(k,j,i+1) + &
3785	( gg - cc ) * de_dy(k,j+1,i) + &
3786	( gg - dd ) * de_dy(k,j+1,i+1) &
3787	) / ( 3.0_wp * gg )
3788	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
3789	de_dy_int(n) = de_dy_int_l
3790	ELSE
3791	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
3792	( gg - bb ) * de_dy(k+1,j,i+1) + &
3793	( gg - cc ) * de_dy(k+1,j+1,i) + &
3794	( gg - dd ) * de_dy(k+1,j+1,i+1) &
3795	) / ( 3.0_wp * gg )
3796	de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * &
3797	( de_dy_int_u - de_dy_int_l )
3798	ENDIF
3799
3800	!
3801	!-- Interpolate the TKE gradient along z
3802	IF ( zv(n) < 0.5_wp * dz(1) ) THEN
3803	de_dz_int(n) = 0.0_wp
3804	ELSE
3805	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
3806	( gg - bb ) * de_dz(k,j,i+1) + &
3807	( gg - cc ) * de_dz(k,j+1,i) + &
3808	( gg - dd ) * de_dz(k,j+1,i+1) &
3809	) / ( 3.0_wp * gg )
3810
3811	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
3812	de_dz_int(n) = de_dz_int_l
3813	ELSE
3814	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
3815	( gg - bb ) * de_dz(k+1,j,i+1) + &
3816	( gg - cc ) * de_dz(k+1,j+1,i) + &
3817	( gg - dd ) * de_dz(k+1,j+1,i+1) &
3818	) / ( 3.0_wp * gg )
3819	de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * &
3820	( de_dz_int_u - de_dz_int_l )
3821	ENDIF
3822	ENDIF
3823
3824	!
3825	!-- Interpolate the dissipation of TKE
3826	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
3827	( gg - bb ) * diss(k,j,i+1) + &
3828	( gg - cc ) * diss(k,j+1,i) + &
3829	( gg - dd ) * diss(k,j+1,i+1) &
3830	) / ( 3.0_wp * gg )
3831
3832	IF ( k == nzt ) THEN
3833	diss_int(n) = diss_int_l
3834	ELSE
3835	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
3836	( gg - bb ) * diss(k+1,j,i+1) + &
3837	( gg - cc ) * diss(k+1,j+1,i) + &
3838	( gg - dd ) * diss(k+1,j+1,i+1) &
3839	) / ( 3.0_wp * gg )
3840	diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dzw(k+1) * &
3841	( diss_int_u - diss_int_l )
3842	ENDIF
3843
3844	!
3845	!-- Set flag for stochastic equation.
3846	term_1_2(n) = 1.0_wp
3847	ENDDO
3848	ENDIF
3849	ENDDO
3850
3851	DO nb = 0,7
3852	i = ip + block_offset(nb)%i_off
3853	j = jp + block_offset(nb)%j_off
3854	k = kp + block_offset(nb)%k_off
3855
3856	DO n = start_index(nb), end_index(nb)
3857	!
3858	!-- Vertical interpolation of the horizontally averaged SGS TKE and
3859	!-- resolved-scale velocity variances and use the interpolated values
3860	!-- to calculate the coefficient fs, which is a measure of the ratio
3861	!-- of the subgrid-scale turbulent kinetic energy to the total amount
3862	!-- of turbulent kinetic energy.
3863	IF ( k == 0 ) THEN
3864	e_mean_int = hom(0,1,8,0)
3865	ELSE
3866	e_mean_int = hom(k,1,8,0) + &
3867	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
3868	( zu(k+1) - zu(k) ) * &
3869	( zv(n) - zu(k) )
3870	ENDIF
3871
3872	kw = kp - 1
3873
3874	IF ( k == 0 ) THEN
3875	aa = hom(k+1,1,30,0) * ( zv(n) / &
3876	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
3877	bb = hom(k+1,1,31,0) * ( zv(n) / &
3878	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
3879	cc = hom(kw+1,1,32,0) * ( zv(n) / &
3880	( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
3881	ELSE
3882	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
3883	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
3884	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
3885	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
3886	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) * &
3887	( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
3888	ENDIF
3889
3890	vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
3891	!
3892	!-- Needed to avoid NaN particle velocities. The value of 1.0 is just
3893	!-- an educated guess for the given case.
3894	IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp ) THEN
3895	fs_int(n) = 1.0_wp
3896	ELSE
3897	fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int / &
3898	( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
3899	ENDIF
3900
3901	ENDDO
3902	ENDDO
3903
3904	DO nb = 0, 7
3905	DO n = start_index(nb), end_index(nb)
3906	rg(n,1) = random_gauss( iran_part, 5.0_wp )
3907	rg(n,2) = random_gauss( iran_part, 5.0_wp )
3908	rg(n,3) = random_gauss( iran_part, 5.0_wp )
3909	ENDDO
3910	ENDDO
3911
3912	DO nb = 0, 7
3913	DO n = start_index(nb), end_index(nb)
3914
3915	!
3916	!-- Calculate the Lagrangian timescale according to Weil et al. (2004).
3917	lagr_timescale(n) = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
3918	( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
3919
3920	!
3921	!-- Calculate the next particle timestep. dt_gap is the time needed to
3922	!-- complete the current LES timestep.
3923	dt_gap(n) = dt_3d - particles(n)%dt_sum
3924	dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale(n), dt_gap(n) )
3925	particles(n)%aux1 = lagr_timescale(n)
3926	particles(n)%aux2 = dt_gap(n)
3927	!
3928	!-- The particle timestep should not be too small in order to prevent
3929	!-- the number of particle timesteps of getting too large
3930	IF ( dt_particle(n) < dt_min_part ) THEN
3931	IF ( dt_min_part < dt_gap(n) ) THEN
3932	dt_particle(n) = dt_min_part
3933	ELSE
3934	dt_particle(n) = dt_gap(n)
3935	ENDIF
3936	ENDIF
3937
3938	rvar1_temp(n) = particles(n)%rvar1
3939	rvar2_temp(n) = particles(n)%rvar2
3940	rvar3_temp(n) = particles(n)%rvar3
3941	!
3942	!-- Calculate the SGS velocity components
3943	IF ( particles(n)%age == 0.0_wp ) THEN
3944	!
3945	!-- For new particles the SGS components are derived from the SGS
3946	!-- TKE. Limit the Gaussian random number to the interval
3947	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
3948	!-- from becoming unrealistically large.
3949	rvar1_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) &
3950	+ 1E-20_wp ) * ( rg(n,1) - 1.0_wp )
3951	rvar2_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) &
3952	+ 1E-20_wp ) * ( rg(n,2) - 1.0_wp )
3953	rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) &
3954	+ 1E-20_wp ) * ( rg(n,3) - 1.0_wp )
3955	ELSE
3956	!
3957	!-- Restriction of the size of the new timestep: compared to the
3958	!-- previous timestep the increase must not exceed 200%. First,
3959	!-- check if age > age_m, in order to prevent that particles get zero
3960	!-- timestep.
3961	dt_particle_m = MERGE( dt_particle(n), &
3962	particles(n)%age - particles(n)%age_m, &
3963	particles(n)%age - particles(n)%age_m < &
3964	1E-8_wp )
3965	IF ( dt_particle(n) > 2.0_wp * dt_particle_m ) THEN
3966	dt_particle(n) = 2.0_wp * dt_particle_m
3967	ENDIF
3968
3969	!-- For old particles the SGS components are correlated with the
3970	!-- values from the previous timestep. Random numbers have also to
3971	!-- be limited (see above).
3972	!-- As negative values for the subgrid TKE are not allowed, the
3973	!-- change of the subgrid TKE with time cannot be smaller than
3974	!-- -e_int(n)/dt_particle. This value is used as a lower boundary
3975	!-- value for the change of TKE
3976	de_dt_min = - e_int(n) / dt_particle(n)
3977
3978	de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
3979
3980	IF ( de_dt < de_dt_min ) THEN
3981	de_dt = de_dt_min
3982	ENDIF
3983
3984	CALL weil_stochastic_eq( rvar1_temp(n), fs_int(n), e_int(n), &
3985	de_dx_int(n), de_dt, diss_int(n), &
3986	dt_particle(n), rg(n,1), term_1_2(n) )
3987
3988	CALL weil_stochastic_eq( rvar2_temp(n), fs_int(n), e_int(n), &
3989	de_dy_int(n), de_dt, diss_int(n), &
3990	dt_particle(n), rg(n,2), term_1_2(n) )
3991
3992	CALL weil_stochastic_eq( rvar3_temp(n), fs_int(n), e_int(n), &
3993	de_dz_int(n), de_dt, diss_int(n), &
3994	dt_particle(n), rg(n,3), term_1_2(n) )
3995
3996	ENDIF
3997
3998	ENDDO
3999	ENDDO
4000	!
4001	!-- Check if the added SGS velocities result in a violation of the CFL-
4002	!-- criterion. If yes, limt the SGS particle speed to match the
4003	!-- CFL criterion. Note, a re-calculation of the SGS particle speed with
4004	!-- smaller timestep does not necessarily fulfill the CFL criterion as the
4005	!-- new SGS speed can be even larger (due to the random term with scales with
4006	!-- the square-root of dt_particle, for small dt the random contribution increases).
4007	!-- Thus, we would need to re-calculate the SGS speeds as long as they would
4008	!-- fulfill the requirements, which could become computationally expensive,
4009	!-- Hence, we just limit them.
4010	dz_temp = zw(kp)-zw(kp-1)
4011
4012	DO nb = 0, 7
4013	DO n = start_index(nb), end_index(nb)
4014	IF ( ABS( u_int(n) + rvar1_temp(n) ) > ( dx / dt_particle(n) ) .OR. &
4015	ABS( v_int(n) + rvar2_temp(n) ) > ( dy / dt_particle(n) ) .OR. &
4016	ABS( w_int(n) + rvar3_temp(n) ) > ( dz_temp / dt_particle(n) ) ) THEN
4017	!
4018	!-- If total speed exceeds the allowed speed according to CFL
4019	!-- criterion, limit the SGS speed to
4020	!-- dx_i / dt_particle - u_resolved_i, considering a safty factor.
4021	rvar1_temp(n) = MERGE( rvar1_temp(n), &
4022	0.9_wp * &
4023	SIGN( dx / dt_particle(n) &
4024	- ABS( u_int(n) ), rvar1_temp(n) ), &
4025	ABS( u_int(n) + rvar1_temp(n) ) < &
4026	( dx / dt_particle(n) ) )
4027	rvar2_temp(n) = MERGE( rvar2_temp(n), &
4028	0.9_wp * &
4029	SIGN( dy / dt_particle(n) &
4030	- ABS( v_int(n) ), rvar2_temp(n) ), &
4031	ABS( v_int(n) + rvar2_temp(n) ) < &
4032	( dy / dt_particle(n) ) )
4033	rvar3_temp(n) = MERGE( rvar3_temp(n), &
4034	0.9_wp * &
4035	SIGN( zw(kp)-zw(kp-1) / dt_particle(n) &
4036	- ABS( w_int(n) ), rvar3_temp(n) ), &
4037	ABS( w_int(n) + rvar3_temp(n) ) < &
4038	( zw(kp)-zw(kp-1) / dt_particle(n) ) )
4039	ENDIF
4040	!
4041	!-- Update particle velocites
4042	particles(n)%rvar1 = rvar1_temp(n)
4043	particles(n)%rvar2 = rvar2_temp(n)
4044	particles(n)%rvar3 = rvar3_temp(n)
4045	u_int(n) = u_int(n) + particles(n)%rvar1
4046	v_int(n) = v_int(n) + particles(n)%rvar2
4047	w_int(n) = w_int(n) + particles(n)%rvar3
4048	!
4049	!-- Store the SGS TKE of the current timelevel which is needed for
4050	!-- for calculating the SGS particle velocities at the next timestep
4051	particles(n)%e_m = e_int(n)
4052	ENDDO
4053	ENDDO
4054
4055	ELSE
4056	!
4057	!-- If no SGS velocities are used, only the particle timestep has to
4058	!-- be set
4059	dt_particle = dt_3d
4060
4061	ENDIF
4062
4063	dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
4064	IF ( ANY( dens_ratio == 0.0_wp ) ) THEN
4065	!
4066	!-- Decide whether the particle loop runs over the subboxes or only over 1,
4067	!-- number_of_particles. This depends on the selected interpolation method.
4068	!-- If particle interpolation method is not trilinear, then the sorting within
4069	!-- subboxes is not required. However, therefore the index start_index(nb) and
4070	!-- end_index(nb) are not defined and the loops are still over
4071	!-- number_of_particles. @todo find a more generic way to write this loop or
4072	!-- delet