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

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

Modularization of all lagrangian particle model code components

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