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

Last change on this file since 4275 was 4275, checked in by schwenkel, 20 months ago

move call of lpm at the end of intermediate timeloop and improve simple corrector particle interpolation method

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