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

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

added restart with MPI-IO for reading local arrays

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