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

Last change on this file since 4495 was 4495, checked in by raasch, 14 months ago

restart data handling with MPI-IO added, first part

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