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

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

Last commit fixed

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