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

Last change on this file since 4018 was 4018, checked in by schwenkel, 2 years ago

bugfix for last commit

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