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

Last change on this file since 4589 was 4589, checked in by suehring, 15 months ago

remove unused variable

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