source: palm/trunk/SOURCE/lagrangian_particle_model_mod.f90 @ 4883

Last change on this file since 4883 was 4883, checked in by hellstea, 7 months ago

user switch for particle coupling added

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