Home

Context Navigation

lagrangian_particle_model_mod.f90 @ 4143

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