Home

Context Navigation

lagrangian_particle_model_mod.f90 @ 4598

Last change on this file since 4598 was 4589, checked in by suehring, 10 months ago
remove unused variable
Property svn:keywords set to `Id`
File size: 362.6 KB

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