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

Last change on this file since 4444 was 4444, checked in by raasch, 3 years ago

bugfix: cpp-directives for serial mode added

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