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

Last change on this file since 4616 was 4616, checked in by schwenkel, 14 months ago

Bugfix in case of strechting: k-calculation limited lower bound of 1

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