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

Last change on this file since 4430 was 4430, checked in by suehring, 3 years ago

Lagrangian particle model: Bugfix in logarithmic interpolation of near-ground particle speed (density was not considered); Revise CFL-check when SGS particle speeds are considered; .palm.iofiles: missing output file connection for child particle statistics

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