Home

Context Navigation

source: palm/trunk/SOURCE/lpm_advec.f90 @ 3038

Last change on this file since 3038 was 2969, checked in by thiele, 7 years ago
Bugfix in Interpolation indices
Property svn:keywords set to `Id`
File size: 55.2 KB

Line
1	!> @file lpm_advec.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-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: lpm_advec.f90 2969 2018-04-13 11:55:09Z gronemeier $
27	! Bugfix in Interpolation indices.
28	!
29	! 2886 2018-03-14 11:51:53Z thiele
30	! Bugfix in passive particle SGS Model:
31	! Sometimes the added SGS velocities would lead to a violation of the CFL
32	! criterion for single particles. For this a check was added after the
33	! calculation of SGS velocities.
34	!
35	! 2718 2018-01-02 08:49:38Z maronga
36	! Corrected "Former revisions" section
37	!
38	! 2701 2017-12-15 15:40:50Z suehring
39	! Changes from last commit documented
40	!
41	! 2698 2017-12-14 18:46:24Z suehring
42	! Particle interpolations at walls in case of SGS velocities revised and not
43	! required parts are removed. (responsible Philipp Thiele)
44	! Bugfix in get_topography_top_index
45	!
46	! 2696 2017-12-14 17:12:51Z kanani
47	! Change in file header (GPL part)
48	!
49	! 2630 2017-11-20 12:58:20Z schwenkel
50	! Removed indices ilog and jlog which are no longer needed since particle box
51	! locations are identical to scalar boxes and topography.
52	!
53	! 2628 2017-11-20 12:40:38Z raasch
54	! bugfix in logarithmic interpolation of v-component (usws was used by mistake)
55	!
56	! 2606 2017-11-10 10:36:31Z schwenkel
57	! Changed particle box locations: center of particle box now coincides
58	! with scalar grid point of same index.
59	! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
60	! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
61	! lpm_sort -> lpm_sort_timeloop_done
62	!
63	! 2417 2017-09-06 15:22:27Z suehring
64	! Particle loops adapted for sub-box structure, i.e. for each sub-box the
65	! particle loop runs from start_index up to end_index instead from 1 to
66	! number_of_particles. This way, it is possible to skip unnecessary
67	! computations for particles that already completed the LES timestep.
68	!
69	! 2318 2017-07-20 17:27:44Z suehring
70	! Get topography top index via Function call
71	!
72	! 2317 2017-07-20 17:27:19Z suehring
73	!
74	! 2232 2017-05-30 17:47:52Z suehring
75	! Adjustments to new topography and surface concept
76	!
77	! 2100 2017-01-05 16:40:16Z suehring
78	! Prevent extremely large SGS-velocities in regions where TKE is zero, e.g.
79	! at the begin of simulations and/or in non-turbulent regions.
80	!
81	! 2000 2016-08-20 18:09:15Z knoop
82	! Forced header and separation lines into 80 columns
83	!
84	! 1936 2016-06-13 13:37:44Z suehring
85	! Formatting adjustments
86	!
87	! 1929 2016-06-09 16:25:25Z suehring
88	! Put stochastic equation in an extra subroutine.
89	! Set flag for stochastic equation to communicate whether a particle is near
90	! topography. This case, memory and drift term are disabled in the Weil equation.
91	!
92	! Enable vertical logarithmic interpolation also above topography. This case,
93	! set a lower limit for the friction velocity, as it can become very small
94	! in narrow street canyons, leading to too large particle speeds.
95	!
96	! 1888 2016-04-21 12:20:49Z suehring
97	! Bugfix concerning logarithmic interpolation of particle speed
98	!
99	! 1822 2016-04-07 07:49:42Z hoffmann
100	! Random velocity fluctuations for particles added. Terminal fall velocity
101	! for droplets is calculated from a parameterization (which is better than
102	! the previous, physically correct calculation, which demands a very short
103	! time step that is not used in the model).
104	!
105	! Unused variables deleted.
106	!
107	! 1691 2015-10-26 16:17:44Z maronga
108	! Renamed prandtl_layer to constant_flux_layer.
109	!
110	! 1685 2015-10-08 07:32:13Z raasch
111	! TKE check for negative values (so far, only zero value was checked)
112	! offset_ocean_nzt_m1 removed
113	!
114	! 1682 2015-10-07 23:56:08Z knoop
115	! Code annotations made doxygen readable
116	!
117	! 1583 2015-04-15 12:16:27Z suehring
118	! Bugfix: particle advection within Prandtl-layer in case of Galilei
119	! transformation.
120	!
121	! 1369 2014-04-24 05:57:38Z raasch
122	! usage of module interfaces removed
123	!
124	! 1359 2014-04-11 17:15:14Z hoffmann
125	! New particle structure integrated.
126	! Kind definition added to all floating point numbers.
127	!
128	! 1322 2014-03-20 16:38:49Z raasch
129	! REAL constants defined as wp_kind
130	!
131	! 1320 2014-03-20 08:40:49Z raasch
132	! ONLY-attribute added to USE-statements,
133	! kind-parameters added to all INTEGER and REAL declaration statements,
134	! kinds are defined in new module kinds,
135	! revision history before 2012 removed,
136	! comment fields (!:) to be used for variable explanations added to
137	! all variable declaration statements
138	!
139	! 1314 2014-03-14 18:25:17Z suehring
140	! Vertical logarithmic interpolation of horizontal particle speed for particles
141	! between roughness height and first vertical grid level.
142	!
143	! 1036 2012-10-22 13:43:42Z raasch
144	! code put under GPL (PALM 3.9)
145	!
146	! 849 2012-03-15 10:35:09Z raasch
147	! initial revision (former part of advec_particles)
148	!
149	!
150	! Description:
151	! ------------
152	!> Calculation of new particle positions due to advection using a simple Euler
153	!> scheme. Particles may feel inertia effects. SGS transport can be included
154	!> using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
155	!------------------------------------------------------------------------------!
156	SUBROUTINE lpm_advec (ip,jp,kp)
157
158
159	USE arrays_3d, &
160	ONLY: de_dx, de_dy, de_dz, diss, e, km, u, v, w, zu, zw
161
162	USE cpulog
163
164	USE pegrid
165
166	USE control_parameters, &
167	ONLY: atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d, &
168	dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans
169
170	USE grid_variables, &
171	ONLY: ddx, dx, ddy, dy
172
173	USE indices, &
174	ONLY: nzb, nzt, wall_flags_0
175
176	USE kinds
177
178	USE particle_attributes, &
179	ONLY: block_offset, c_0, dt_min_part, grid_particles, &
180	iran_part, log_z_z0, number_of_particles, number_of_sublayers, &
181	particles, particle_groups, offset_ocean_nzt, sgs_wf_part, &
182	use_sgs_for_particles, vertical_particle_advection, z0_av_global
183
184	USE statistics, &
185	ONLY: hom
186
187	USE surface_mod, &
188	ONLY: get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
189
190	IMPLICIT NONE
191
192	LOGICAL :: subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
193
194	INTEGER(iwp) :: agp !< loop variable
195	INTEGER(iwp) :: gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography
196	INTEGER(iwp) :: i !< index variable along x
197	INTEGER(iwp) :: ip !< index variable along x
198	INTEGER(iwp) :: j !< index variable along y
199	INTEGER(iwp) :: jp !< index variable along y
200	INTEGER(iwp) :: k !< index variable along z
201	INTEGER(iwp) :: k_wall !< vertical index of topography top
202	INTEGER(iwp) :: kp !< index variable along z
203	INTEGER(iwp) :: kw !< index variable along z
204	INTEGER(iwp) :: n !< loop variable over all particles in a grid box
205	INTEGER(iwp) :: nb !< block number particles are sorted in
206	INTEGER(iwp) :: num_gp !< number of adjacent grid points inside topography
207	INTEGER(iwp) :: surf_start !< Index on surface data-type for current grid box
208
209	INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current block
210	INTEGER(iwp), DIMENSION(0:7) :: end_index !< start particle index for current block
211
212	REAL(wp) :: aa !< dummy argument for horizontal particle interpolation
213	REAL(wp) :: bb !< dummy argument for horizontal particle interpolation
214	REAL(wp) :: cc !< dummy argument for horizontal particle interpolation
215	REAL(wp) :: d_sum !< dummy argument for horizontal particle interpolation in case of topography
216	REAL(wp) :: d_z_p_z0 !< inverse of interpolation length for logarithmic interpolation
217	REAL(wp) :: dd !< dummy argument for horizontal particle interpolation
218	REAL(wp) :: de_dx_int_l !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
219	REAL(wp) :: de_dx_int_u !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
220	REAL(wp) :: de_dy_int_l !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
221	REAL(wp) :: de_dy_int_u !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
222	REAL(wp) :: de_dt !< temporal derivative of TKE experienced by the particle
223	REAL(wp) :: de_dt_min !< lower level for temporal TKE derivative
224	REAL(wp) :: de_dz_int_l !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
225	REAL(wp) :: de_dz_int_u !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
226	REAL(wp) :: diameter !< diamter of droplet
227	REAL(wp) :: diss_int_l !< x/y-interpolated dissipation at particle position at lower vertical level
228	REAL(wp) :: diss_int_u !< x/y-interpolated dissipation at particle position at upper vertical level
229	REAL(wp) :: dt_particle_m !< previous particle time step
230	REAL(wp) :: dz_temp !<
231	REAL(wp) :: e_int_l !< x/y-interpolated TKE at particle position at lower vertical level
232	REAL(wp) :: e_int_u !< x/y-interpolated TKE at particle position at upper vertical level
233	REAL(wp) :: e_mean_int !< horizontal mean TKE at particle height
234	REAL(wp) :: exp_arg !<
235	REAL(wp) :: exp_term !<
236	REAL(wp) :: gg !< dummy argument for horizontal particle interpolation
237	REAL(wp) :: height_p !< dummy argument for logarithmic interpolation
238	REAL(wp) :: location(1:30,1:3) !< wall locations
239	REAL(wp) :: log_z_z0_int !< logarithmus used for surface_layer interpolation
240	REAL(wp) :: random_gauss !<
241	REAL(wp) :: RL !< Lagrangian autocorrelation coefficient
242	REAL(wp) :: rg1 !< Gaussian distributed random number
243	REAL(wp) :: rg2 !< Gaussian distributed random number
244	REAL(wp) :: rg3 !< Gaussian distributed random number
245	REAL(wp) :: sigma !< velocity standard deviation
246	REAL(wp) :: u_int_l !< x/y-interpolated u-component at particle position at lower vertical level
247	REAL(wp) :: u_int_u !< x/y-interpolated u-component at particle position at upper vertical level
248	REAL(wp) :: us_int !< friction velocity at particle grid box
249	REAL(wp) :: usws_int !< surface momentum flux (u component) at particle grid box
250	REAL(wp) :: v_int_l !< x/y-interpolated v-component at particle position at lower vertical level
251	REAL(wp) :: v_int_u !< x/y-interpolated v-component at particle position at upper vertical level
252	REAL(wp) :: vsws_int !< surface momentum flux (u component) at particle grid box
253	REAL(wp) :: vv_int !<
254	REAL(wp) :: w_int_l !< x/y-interpolated w-component at particle position at lower vertical level
255	REAL(wp) :: w_int_u !< x/y-interpolated w-component at particle position at upper vertical level
256	REAL(wp) :: w_s !< terminal velocity of droplets
257	REAL(wp) :: x !< dummy argument for horizontal particle interpolation
258	REAL(wp) :: y !< dummy argument for horizontal particle interpolation
259	REAL(wp) :: z_p !< surface layer height (0.5 dz)
260
261	REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter for fall velocity
262	REAL(wp), PARAMETER :: b_rog = 10.43_wp !< parameter for fall velocity
263	REAL(wp), PARAMETER :: c_rog = 0.6_wp !< parameter for fall velocity
264	REAL(wp), PARAMETER :: k_cap_rog = 4.0_wp !< parameter for fall velocity
265	REAL(wp), PARAMETER :: k_low_rog = 12.0_wp !< parameter for fall velocity
266	REAL(wp), PARAMETER :: d0_rog = 0.745_wp !< separation diameter
267
268	REAL(wp), DIMENSION(1:30) :: d_gp_pl !< dummy argument for particle interpolation scheme in case of topography
269	REAL(wp), DIMENSION(1:30) :: de_dxi !< horizontal TKE gradient along x at adjacent wall
270	REAL(wp), DIMENSION(1:30) :: de_dyi !< horizontal TKE gradient along y at adjacent wall
271	REAL(wp), DIMENSION(1:30) :: de_dzi !< horizontal TKE gradient along z at adjacent wall
272	REAL(wp), DIMENSION(1:30) :: dissi !< dissipation at adjacent wall
273	REAL(wp), DIMENSION(1:30) :: ei !< TKE at adjacent wall
274
275	REAL(wp), DIMENSION(number_of_particles) :: term_1_2 !< flag to communicate whether a particle is near topography or not
276	REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !<
277	REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< horizontal TKE gradient along x at particle position
278	REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< horizontal TKE gradient along y at particle position
279	REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< horizontal TKE gradient along z at particle position
280	REAL(wp), DIMENSION(number_of_particles) :: diss_int !< dissipation at particle position
281	REAL(wp), DIMENSION(number_of_particles) :: dt_gap !< remaining time until particle time integration reaches LES time
282	REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< particle time step
283	REAL(wp), DIMENSION(number_of_particles) :: e_int !< TKE at particle position
284	REAL(wp), DIMENSION(number_of_particles) :: fs_int !< weighting factor for subgrid-scale particle speed
285	REAL(wp), DIMENSION(number_of_particles) :: lagr_timescale !< Lagrangian timescale
286	REAL(wp), DIMENSION(number_of_particles) :: rvar1_temp !<
287	REAL(wp), DIMENSION(number_of_particles) :: rvar2_temp !<
288	REAL(wp), DIMENSION(number_of_particles) :: rvar3_temp !<
289	REAL(wp), DIMENSION(number_of_particles) :: u_int !< u-component of particle speed
290	REAL(wp), DIMENSION(number_of_particles) :: v_int !< v-component of particle speed
291	REAL(wp), DIMENSION(number_of_particles) :: w_int !< w-component of particle speed
292	REAL(wp), DIMENSION(number_of_particles) :: xv !< x-position
293	REAL(wp), DIMENSION(number_of_particles) :: yv !< y-position
294	REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position
295
296	REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers
297
298	CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
299
300	!
301	!-- Determine height of Prandtl layer and distance between Prandtl-layer
302	!-- height and horizontal mean roughness height, which are required for
303	!-- vertical logarithmic interpolation of horizontal particle speeds
304	!-- (for particles below first vertical grid level).
305	z_p = zu(nzb+1) - zw(nzb)
306	d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
307
308	start_index = grid_particles(kp,jp,ip)%start_index
309	end_index = grid_particles(kp,jp,ip)%end_index
310
311	xv = particles(1:number_of_particles)%x
312	yv = particles(1:number_of_particles)%y
313	zv = particles(1:number_of_particles)%z
314
315	DO nb = 0, 7
316	!
317	!-- Interpolate u velocity-component
318	i = ip
319	j = jp + block_offset(nb)%j_off
320	k = kp + block_offset(nb)%k_off
321
322	DO n = start_index(nb), end_index(nb)
323	!
324	!-- Interpolation of the u velocity component onto particle position.
325	!-- Particles are interpolation bi-linearly in the horizontal and a
326	!-- linearly in the vertical. An exception is made for particles below
327	!-- the first vertical grid level in case of a prandtl layer. In this
328	!-- case the horizontal particle velocity components are determined using
329	!-- Monin-Obukhov relations (if branch).
330	!-- First, check if particle is located below first vertical grid level
331	!-- above topography (Prandtl-layer height)
332	!-- Determine vertical index of topography top
333	k_wall = get_topography_top_index_ji( jp, ip, 's' )
334
335	IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN
336	!
337	!-- Resolved-scale horizontal particle velocity is zero below z0.
338	IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN
339	u_int(n) = 0.0_wp
340	ELSE
341	!
342	!-- Determine the sublayer. Further used as index.
343	height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
344	* REAL( number_of_sublayers, KIND=wp ) &
345	* d_z_p_z0
346	!
347	!-- Calculate LOG(z/z0) for exact particle height. Therefore,
348	!-- interpolate linearly between precalculated logarithm.
349	log_z_z0_int = log_z_z0(INT(height_p)) &
350	+ ( height_p - INT(height_p) ) &
351	* ( log_z_z0(INT(height_p)+1) &
352	- log_z_z0(INT(height_p)) &
353	)
354	!
355	!-- Get friction velocity and momentum flux from new surface data
356	!-- types.
357	IF ( surf_def_h(0)%start_index(jp,ip) <= &
358	surf_def_h(0)%end_index(jp,ip) ) THEN
359	surf_start = surf_def_h(0)%start_index(jp,ip)
360	!-- Limit friction velocity. In narrow canyons or holes the
361	!-- friction velocity can become very small, resulting in a too
362	!-- large particle speed.
363	us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
364	usws_int = surf_def_h(0)%usws(surf_start)
365	ELSEIF ( surf_lsm_h%start_index(jp,ip) <= &
366	surf_lsm_h%end_index(jp,ip) ) THEN
367	surf_start = surf_lsm_h%start_index(jp,ip)
368	us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
369	usws_int = surf_lsm_h%usws(surf_start)
370	ELSEIF ( surf_usm_h%start_index(jp,ip) <= &
371	surf_usm_h%end_index(jp,ip) ) THEN
372	surf_start = surf_usm_h%start_index(jp,ip)
373	us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp )
374	usws_int = surf_usm_h%usws(surf_start)
375	ENDIF
376
377	!
378	!-- Neutral solution is applied for all situations, e.g. also for
379	!-- unstable and stable situations. Even though this is not exact
380	!-- this saves a lot of CPU time since several calls of intrinsic
381	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
382	!-- as sensitivity studies revealed no significant effect of
383	!-- using the neutral solution also for un/stable situations.
384	u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp ) &
385	* log_z_z0_int - u_gtrans
386
387	ENDIF
388	!
389	!-- Particle above the first grid level. Bi-linear interpolation in the
390	!-- horizontal and linear interpolation in the vertical direction.
391	ELSE
392
393	x = xv(n) - i * dx
394	y = yv(n) + ( 0.5_wp - j ) * dy
395	aa = x2 + y2
396	bb = ( dx - x )2 + y2
397	cc = x2 + ( dy - y )2
398	dd = ( dx - x )2 + ( dy - y )2
399	gg = aa + bb + cc + dd
400
401	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
402	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * &
403	u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
404
405	IF ( k == nzt ) THEN
406	u_int(n) = u_int_l
407	ELSE
408	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
409	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * &
410	u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
411	u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz * &
412	( u_int_u - u_int_l )
413	ENDIF
414
415	ENDIF
416
417	ENDDO
418	!
419	!-- Same procedure for interpolation of the v velocity-component
420	i = ip + block_offset(nb)%i_off
421	j = jp
422	k = kp + block_offset(nb)%k_off
423
424	DO n = start_index(nb), end_index(nb)
425
426	!
427	!-- Determine vertical index of topography top
428	k_wall = get_topography_top_index_ji( jp,ip, 's' )
429
430	IF ( constant_flux_layer .AND. zv(n) - zw(k_wall) < z_p ) THEN
431	IF ( zv(n) - zw(k_wall) < z0_av_global ) THEN
432	!
433	!-- Resolved-scale horizontal particle velocity is zero below z0.
434	v_int(n) = 0.0_wp
435	ELSE
436	!
437	!-- Determine the sublayer. Further used as index. Please note,
438	!-- logarithmus can not be reused from above, as in in case of
439	!-- topography particle on u-grid can be above surface-layer height,
440	!-- whereas it can be below on v-grid.
441	height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
442	* REAL( number_of_sublayers, KIND=wp ) &
443	* d_z_p_z0
444	!
445	!-- Calculate LOG(z/z0) for exact particle height. Therefore,
446	!-- interpolate linearly between precalculated logarithm.
447	log_z_z0_int = log_z_z0(INT(height_p)) &
448	+ ( height_p - INT(height_p) ) &
449	* ( log_z_z0(INT(height_p)+1) &
450	- log_z_z0(INT(height_p)) &
451	)
452	!
453	!-- Get friction velocity and momentum flux from new surface data
454	!-- types.
455	IF ( surf_def_h(0)%start_index(jp,ip) <= &
456	surf_def_h(0)%end_index(jp,ip) ) THEN
457	surf_start = surf_def_h(0)%start_index(jp,ip)
458	!-- Limit friction velocity. In narrow canyons or holes the
459	!-- friction velocity can become very small, resulting in a too
460	!-- large particle speed.
461	us_int = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
462	vsws_int = surf_def_h(0)%vsws(surf_start)
463	ELSEIF ( surf_lsm_h%start_index(jp,ip) <= &
464	surf_lsm_h%end_index(jp,ip) ) THEN
465	surf_start = surf_lsm_h%start_index(jp,ip)
466	us_int = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
467	vsws_int = surf_lsm_h%vsws(surf_start)
468	ELSEIF ( surf_usm_h%start_index(jp,ip) <= &
469	surf_usm_h%end_index(jp,ip) ) THEN
470	surf_start = surf_usm_h%start_index(jp,ip)
471	us_int = MAX( surf_usm_h%us(surf_start), 0.01_wp )
472	vsws_int = surf_usm_h%vsws(surf_start)
473	ENDIF
474	!
475	!-- Neutral solution is applied for all situations, e.g. also for
476	!-- unstable and stable situations. Even though this is not exact
477	!-- this saves a lot of CPU time since several calls of intrinsic
478	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
479	!-- as sensitivity studies revealed no significant effect of
480	!-- using the neutral solution also for un/stable situations.
481	v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp ) &
482	* log_z_z0_int - v_gtrans
483
484	ENDIF
485
486	ELSE
487	x = xv(n) + ( 0.5_wp - i ) * dx
488	y = yv(n) - j * dy
489	aa = x2 + y2
490	bb = ( dx - x )2 + y2
491	cc = x2 + ( dy - y )2
492	dd = ( dx - x )2 + ( dy - y )2
493	gg = aa + bb + cc + dd
494
495	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
496	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
497	) / ( 3.0_wp * gg ) - v_gtrans
498
499	IF ( k == nzt ) THEN
500	v_int(n) = v_int_l
501	ELSE
502	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
503	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
504	) / ( 3.0_wp * gg ) - v_gtrans
505	v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
506	( v_int_u - v_int_l )
507	ENDIF
508
509	ENDIF
510
511	ENDDO
512	!
513	!-- Same procedure for interpolation of the w velocity-component
514	i = ip + block_offset(nb)%i_off
515	j = jp + block_offset(nb)%j_off
516	k = kp - 1
517
518	DO n = start_index(nb), end_index(nb)
519
520	IF ( vertical_particle_advection(particles(n)%group) ) THEN
521
522	x = xv(n) + ( 0.5_wp - i ) * dx
523	y = yv(n) + ( 0.5_wp - j ) * dy
524	aa = x2 + y2
525	bb = ( dx - x )2 + y2
526	cc = x2 + ( dy - y )2
527	dd = ( dx - x )2 + ( dy - y )2
528	gg = aa + bb + cc + dd
529
530	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
531	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
532	) / ( 3.0_wp * gg )
533
534	IF ( k == nzt ) THEN
535	w_int(n) = w_int_l
536	ELSE
537	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
538	( gg-bb ) * w(k+1,j,i+1) + &
539	( gg-cc ) * w(k+1,j+1,i) + &
540	( gg-dd ) * w(k+1,j+1,i+1) &
541	) / ( 3.0_wp * gg )
542	w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
543	( w_int_u - w_int_l )
544	ENDIF
545
546	ELSE
547
548	w_int(n) = 0.0_wp
549
550	ENDIF
551
552	ENDDO
553
554	ENDDO
555
556	!-- Interpolate and calculate quantities needed for calculating the SGS
557	!-- velocities
558	IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN
559
560	DO nb = 0,7
561
562	subbox_at_wall = .FALSE.
563	!
564	!-- In case of topography check if subbox is adjacent to a wall
565	IF ( .NOT. topography == 'flat' ) THEN
566	i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) )
567	j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) )
568	k = kp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 0 ) )
569	IF ( .NOT. BTEST(wall_flags_0(k, jp, ip), 0) .OR. &
570	.NOT. BTEST(wall_flags_0(kp, j, ip), 0) .OR. &
571	.NOT. BTEST(wall_flags_0(kp, jp, i ), 0) ) &
572	THEN
573	subbox_at_wall = .TRUE.
574	ENDIF
575	ENDIF
576	IF ( subbox_at_wall ) THEN
577	e_int(start_index(nb):end_index(nb)) = e(kp,jp,ip)
578	diss_int(start_index(nb):end_index(nb)) = diss(kp,jp,ip)
579	de_dx_int(start_index(nb):end_index(nb)) = de_dx(kp,jp,ip)
580	de_dy_int(start_index(nb):end_index(nb)) = de_dy(kp,jp,ip)
581	de_dz_int(start_index(nb):end_index(nb)) = de_dz(kp,jp,ip)
582	!
583	!-- Set flag for stochastic equation.
584	term_1_2(start_index(nb):end_index(nb)) = 0.0_wp
585	ELSE
586	i = ip + block_offset(nb)%i_off
587	j = jp + block_offset(nb)%j_off
588	k = kp + block_offset(nb)%k_off
589
590	DO n = start_index(nb), end_index(nb)
591	!
592	!-- Interpolate TKE
593	x = xv(n) + ( 0.5_wp - i ) * dx
594	y = yv(n) + ( 0.5_wp - j ) * dy
595	aa = x2 + y2
596	bb = ( dx - x )2 + y2
597	cc = x2 + ( dy - y )2
598	dd = ( dx - x )2 + ( dy - y )2
599	gg = aa + bb + cc + dd
600
601	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
602	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
603	) / ( 3.0_wp * gg )
604
605	IF ( k+1 == nzt+1 ) THEN
606	e_int(n) = e_int_l
607	ELSE
608	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
609	( gg - bb ) * e(k+1,j,i+1) + &
610	( gg - cc ) * e(k+1,j+1,i) + &
611	( gg - dd ) * e(k+1,j+1,i+1) &
612	) / ( 3.0_wp * gg )
613	e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
614	( e_int_u - e_int_l )
615	ENDIF
616	!
617	!-- Needed to avoid NaN particle velocities (this might not be
618	!-- required any more)
619	IF ( e_int(n) <= 0.0_wp ) THEN
620	e_int(n) = 1.0E-20_wp
621	ENDIF
622	!
623	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
624	!-- all position variables from above (TKE))
625	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
626	( gg - bb ) * de_dx(k,j,i+1) + &
627	( gg - cc ) * de_dx(k,j+1,i) + &
628	( gg - dd ) * de_dx(k,j+1,i+1) &
629	) / ( 3.0_wp * gg )
630
631	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
632	de_dx_int(n) = de_dx_int_l
633	ELSE
634	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
635	( gg - bb ) * de_dx(k+1,j,i+1) + &
636	( gg - cc ) * de_dx(k+1,j+1,i) + &
637	( gg - dd ) * de_dx(k+1,j+1,i+1) &
638	) / ( 3.0_wp * gg )
639	de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
640	( de_dx_int_u - de_dx_int_l )
641	ENDIF
642	!
643	!-- Interpolate the TKE gradient along y
644	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
645	( gg - bb ) * de_dy(k,j,i+1) + &
646	( gg - cc ) * de_dy(k,j+1,i) + &
647	( gg - dd ) * de_dy(k,j+1,i+1) &
648	) / ( 3.0_wp * gg )
649	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
650	de_dy_int(n) = de_dy_int_l
651	ELSE
652	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
653	( gg - bb ) * de_dy(k+1,j,i+1) + &
654	( gg - cc ) * de_dy(k+1,j+1,i) + &
655	( gg - dd ) * de_dy(k+1,j+1,i+1) &
656	) / ( 3.0_wp * gg )
657	de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
658	( de_dy_int_u - de_dy_int_l )
659	ENDIF
660
661	!
662	!-- Interpolate the TKE gradient along z
663	IF ( zv(n) < 0.5_wp * dz ) THEN
664	de_dz_int(n) = 0.0_wp
665	ELSE
666	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
667	( gg - bb ) * de_dz(k,j,i+1) + &
668	( gg - cc ) * de_dz(k,j+1,i) + &
669	( gg - dd ) * de_dz(k,j+1,i+1) &
670	) / ( 3.0_wp * gg )
671
672	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
673	de_dz_int(n) = de_dz_int_l
674	ELSE
675	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
676	( gg - bb ) * de_dz(k+1,j,i+1) + &
677	( gg - cc ) * de_dz(k+1,j+1,i) + &
678	( gg - dd ) * de_dz(k+1,j+1,i+1) &
679	) / ( 3.0_wp * gg )
680	de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
681	( de_dz_int_u - de_dz_int_l )
682	ENDIF
683	ENDIF
684
685	!
686	!-- Interpolate the dissipation of TKE
687	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
688	( gg - bb ) * diss(k,j,i+1) + &
689	( gg - cc ) * diss(k,j+1,i) + &
690	( gg - dd ) * diss(k,j+1,i+1) &
691	) / ( 3.0_wp * gg )
692
693	IF ( k == nzt ) THEN
694	diss_int(n) = diss_int_l
695	ELSE
696	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
697	( gg - bb ) * diss(k+1,j,i+1) + &
698	( gg - cc ) * diss(k+1,j+1,i) + &
699	( gg - dd ) * diss(k+1,j+1,i+1) &
700	) / ( 3.0_wp * gg )
701	diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
702	( diss_int_u - diss_int_l )
703	ENDIF
704
705	!
706	!-- Set flag for stochastic equation.
707	term_1_2(n) = 1.0_wp
708	ENDDO
709	ENDIF
710	ENDDO
711
712	DO nb = 0,7
713	i = ip + block_offset(nb)%i_off
714	j = jp + block_offset(nb)%j_off
715	k = kp + block_offset(nb)%k_off
716
717	DO n = start_index(nb), end_index(nb)
718	!
719	!-- Vertical interpolation of the horizontally averaged SGS TKE and
720	!-- resolved-scale velocity variances and use the interpolated values
721	!-- to calculate the coefficient fs, which is a measure of the ratio
722	!-- of the subgrid-scale turbulent kinetic energy to the total amount
723	!-- of turbulent kinetic energy.
724	IF ( k == 0 ) THEN
725	e_mean_int = hom(0,1,8,0)
726	ELSE
727	e_mean_int = hom(k,1,8,0) + &
728	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
729	( zu(k+1) - zu(k) ) * &
730	( zv(n) - zu(k) )
731	ENDIF
732
733	kw = kp - 1
734
735	IF ( k == 0 ) THEN
736	aa = hom(k+1,1,30,0) * ( zv(n) / &
737	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
738	bb = hom(k+1,1,31,0) * ( zv(n) / &
739	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
740	cc = hom(kw+1,1,32,0) * ( zv(n) / &
741	( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
742	ELSE
743	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
744	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
745	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
746	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
747	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) * &
748	( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
749	ENDIF
750
751	vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
752	!
753	!-- Needed to avoid NaN particle velocities. The value of 1.0 is just
754	!-- an educated guess for the given case.
755	IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp ) THEN
756	fs_int(n) = 1.0_wp
757	ELSE
758	fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int / &
759	( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
760	ENDIF
761
762	ENDDO
763	ENDDO
764
765	DO nb = 0, 7
766	DO n = start_index(nb), end_index(nb)
767	rg(n,1) = random_gauss( iran_part, 5.0_wp )
768	rg(n,2) = random_gauss( iran_part, 5.0_wp )
769	rg(n,3) = random_gauss( iran_part, 5.0_wp )
770	ENDDO
771	ENDDO
772
773	DO nb = 0, 7
774	DO n = start_index(nb), end_index(nb)
775
776	!
777	!-- Calculate the Lagrangian timescale according to Weil et al. (2004).
778	lagr_timescale(n) = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
779	( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
780
781	!
782	!-- Calculate the next particle timestep. dt_gap is the time needed to
783	!-- complete the current LES timestep.
784	dt_gap(n) = dt_3d - particles(n)%dt_sum
785	dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale(n), dt_gap(n) )
786	particles(n)%aux1 = lagr_timescale(n)
787	particles(n)%aux2 = dt_gap(n)
788	!
789	!-- The particle timestep should not be too small in order to prevent
790	!-- the number of particle timesteps of getting too large
791	IF ( dt_particle(n) < dt_min_part .AND. dt_min_part < dt_gap(n) ) THEN
792	dt_particle(n) = dt_min_part
793	ENDIF
794	rvar1_temp(n) = particles(n)%rvar1
795	rvar2_temp(n) = particles(n)%rvar2
796	rvar3_temp(n) = particles(n)%rvar3
797	!
798	!-- Calculate the SGS velocity components
799	IF ( particles(n)%age == 0.0_wp ) THEN
800	!
801	!-- For new particles the SGS components are derived from the SGS
802	!-- TKE. Limit the Gaussian random number to the interval
803	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
804	!-- from becoming unrealistically large.
805	rvar1_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) &
806	+ 1E-20_wp ) * ( rg(n,1) - 1.0_wp )
807	rvar2_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) &
808	+ 1E-20_wp ) * ( rg(n,2) - 1.0_wp )
809	rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n) &
810	+ 1E-20_wp ) * ( rg(n,3) - 1.0_wp )
811
812	ELSE
813	!
814	!-- Restriction of the size of the new timestep: compared to the
815	!-- previous timestep the increase must not exceed 200%. First,
816	!-- check if age > age_m, in order to prevent that particles get zero
817	!-- timestep.
818	dt_particle_m = MERGE( dt_particle(n), &
819	particles(n)%age - particles(n)%age_m, &
820	particles(n)%age - particles(n)%age_m < &
821	1E-8_wp )
822	IF ( dt_particle(n) > 2.0_wp * dt_particle_m ) THEN
823	dt_particle(n) = 2.0_wp * dt_particle_m
824	ENDIF
825
826	!-- For old particles the SGS components are correlated with the
827	!-- values from the previous timestep. Random numbers have also to
828	!-- be limited (see above).
829	!-- As negative values for the subgrid TKE are not allowed, the
830	!-- change of the subgrid TKE with time cannot be smaller than
831	!-- -e_int(n)/dt_particle. This value is used as a lower boundary
832	!-- value for the change of TKE
833	de_dt_min = - e_int(n) / dt_particle(n)
834
835	de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
836
837	IF ( de_dt < de_dt_min ) THEN
838	de_dt = de_dt_min
839	ENDIF
840
841	CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),&
842	de_dx_int(n), de_dt, diss_int(n), &
843	dt_particle(n), rg(n,1), term_1_2(n) )
844
845	CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),&
846	de_dy_int(n), de_dt, diss_int(n), &
847	dt_particle(n), rg(n,2), term_1_2(n) )
848
849	CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),&
850	de_dz_int(n), de_dt, diss_int(n), &
851	dt_particle(n), rg(n,3), term_1_2(n) )
852
853	ENDIF
854
855	ENDDO
856	ENDDO
857	!
858	!-- Check if the added SGS velocities result in a violation of the CFL-
859	!-- criterion. If yes choose a smaller timestep based on the new velocities
860	!-- and calculate SGS velocities again
861	dz_temp = zw(kp)-zw(kp-1)
862
863	DO nb = 0, 7
864	DO n = start_index(nb), end_index(nb)
865	IF ( .NOT. particles(n)%age == 0.0_wp .AND. &
866	(ABS( u_int(n) + rvar1_temp(n) ) > (dx/dt_particle(n)) .OR. &
867	ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n)) .OR. &
868	ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN
869
870	dt_particle(n) = 0.9_wp * MIN( &
871	( dx / ABS( u_int(n) + rvar1_temp(n) ) ), &
872	( dy / ABS( v_int(n) + rvar2_temp(n) ) ), &
873	( dz_temp / ABS( w_int(n) + rvar3_temp(n) ) ) )
874
875	!
876	!-- Reset temporary SGS velocites to "current" ones
877	rvar1_temp(n) = particles(n)%rvar1
878	rvar2_temp(n) = particles(n)%rvar2
879	rvar3_temp(n) = particles(n)%rvar3
880
881	de_dt_min = - e_int(n) / dt_particle(n)
882
883	de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
884
885	IF ( de_dt < de_dt_min ) THEN
886	de_dt = de_dt_min
887	ENDIF
888
889	CALL weil_stochastic_eq(rvar1_temp(n), fs_int(n), e_int(n),&
890	de_dx_int(n), de_dt, diss_int(n), &
891	dt_particle(n), rg(n,1), term_1_2(n) )
892
893	CALL weil_stochastic_eq(rvar2_temp(n), fs_int(n), e_int(n),&
894	de_dy_int(n), de_dt, diss_int(n), &
895	dt_particle(n), rg(n,2), term_1_2(n) )
896
897	CALL weil_stochastic_eq(rvar3_temp(n), fs_int(n), e_int(n),&
898	de_dz_int(n), de_dt, diss_int(n), &
899	dt_particle(n), rg(n,3), term_1_2(n) )
900	ENDIF
901
902	!
903	!-- Update particle velocites
904	particles(n)%rvar1 = rvar1_temp(n)
905	particles(n)%rvar2 = rvar2_temp(n)
906	particles(n)%rvar3 = rvar3_temp(n)
907	u_int(n) = u_int(n) + particles(n)%rvar1
908	v_int(n) = v_int(n) + particles(n)%rvar2
909	w_int(n) = w_int(n) + particles(n)%rvar3
910	!
911	!-- Store the SGS TKE of the current timelevel which is needed for
912	!-- for calculating the SGS particle velocities at the next timestep
913	particles(n)%e_m = e_int(n)
914	ENDDO
915	ENDDO
916
917	ELSE
918	!
919	!-- If no SGS velocities are used, only the particle timestep has to
920	!-- be set
921	dt_particle = dt_3d
922
923	ENDIF
924
925	dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
926
927	IF ( ANY( dens_ratio == 0.0_wp ) ) THEN
928	DO nb = 0, 7
929	DO n = start_index(nb), end_index(nb)
930
931	!
932	!-- Particle advection
933	IF ( dens_ratio(n) == 0.0_wp ) THEN
934	!
935	!-- Pure passive transport (without particle inertia)
936	particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
937	particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
938	particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
939
940	particles(n)%speed_x = u_int(n)
941	particles(n)%speed_y = v_int(n)
942	particles(n)%speed_z = w_int(n)
943
944	ELSE
945	!
946	!-- Transport of particles with inertia
947	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
948	dt_particle(n)
949	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
950	dt_particle(n)
951	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
952	dt_particle(n)
953
954	!
955	!-- Update of the particle velocity
956	IF ( cloud_droplets ) THEN
957	!
958	!-- Terminal velocity is computed for vertical direction (Rogers et
959	!-- al., 1993, J. Appl. Meteorol.)
960	diameter = particles(n)%radius * 2000.0_wp !diameter in mm
961	IF ( diameter <= d0_rog ) THEN
962	w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
963	ELSE
964	w_s = a_rog - b_rog * EXP( -c_rog * diameter )
965	ENDIF
966
967	!
968	!-- If selected, add random velocities following Soelch and Kaercher
969	!-- (2010, Q. J. R. Meteorol. Soc.)
970	IF ( use_sgs_for_particles ) THEN
971	lagr_timescale(n) = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
972	RL = EXP( -1.0_wp * dt_3d / lagr_timescale(n) )
973	sigma = SQRT( e(kp,jp,ip) )
974
975	rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
976	rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
977	rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
978
979	particles(n)%rvar1 = RL * particles(n)%rvar1 + &
980	SQRT( 1.0_wp - RL*2 ) sigma * rg1
981	particles(n)%rvar2 = RL * particles(n)%rvar2 + &
982	SQRT( 1.0_wp - RL*2 ) sigma * rg2
983	particles(n)%rvar3 = RL * particles(n)%rvar3 + &
984	SQRT( 1.0_wp - RL*2 ) sigma * rg3
985
986	particles(n)%speed_x = u_int(n) + particles(n)%rvar1
987	particles(n)%speed_y = v_int(n) + particles(n)%rvar2
988	particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
989	ELSE
990	particles(n)%speed_x = u_int(n)
991	particles(n)%speed_y = v_int(n)
992	particles(n)%speed_z = w_int(n) - w_s
993	ENDIF
994
995	ELSE
996
997	IF ( use_sgs_for_particles ) THEN
998	exp_arg = particle_groups(particles(n)%group)%exp_arg
999	exp_term = EXP( -exp_arg * dt_particle(n) )
1000	ELSE
1001	exp_arg = particle_groups(particles(n)%group)%exp_arg
1002	exp_term = particle_groups(particles(n)%group)%exp_term
1003	ENDIF
1004	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
1005	u_int(n) * ( 1.0_wp - exp_term )
1006	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
1007	v_int(n) * ( 1.0_wp - exp_term )
1008	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
1009	( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
1010	g / exp_arg ) * ( 1.0_wp - exp_term )
1011	ENDIF
1012
1013	ENDIF
1014	ENDDO
1015	ENDDO
1016
1017	ELSE
1018
1019	DO nb = 0, 7
1020	DO n = start_index(nb), end_index(nb)
1021	!
1022	!-- Transport of particles with inertia
1023	particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
1024	particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
1025	particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
1026	!
1027	!-- Update of the particle velocity
1028	IF ( cloud_droplets ) THEN
1029	!
1030	!-- Terminal velocity is computed for vertical direction (Rogers et al.,
1031	!-- 1993, J. Appl. Meteorol.)
1032	diameter = particles(n)%radius * 2000.0_wp !diameter in mm
1033	IF ( diameter <= d0_rog ) THEN
1034	w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
1035	ELSE
1036	w_s = a_rog - b_rog * EXP( -c_rog * diameter )
1037	ENDIF
1038
1039	!
1040	!-- If selected, add random velocities following Soelch and Kaercher
1041	!-- (2010, Q. J. R. Meteorol. Soc.)
1042	IF ( use_sgs_for_particles ) THEN
1043	lagr_timescale(n) = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
1044	RL = EXP( -1.0_wp * dt_3d / lagr_timescale(n) )
1045	sigma = SQRT( e(kp,jp,ip) )
1046
1047	rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1048	rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1049	rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
1050
1051	particles(n)%rvar1 = RL * particles(n)%rvar1 + &
1052	SQRT( 1.0_wp - RL*2 ) sigma * rg1
1053	particles(n)%rvar2 = RL * particles(n)%rvar2 + &
1054	SQRT( 1.0_wp - RL*2 ) sigma * rg2
1055	particles(n)%rvar3 = RL * particles(n)%rvar3 + &
1056	SQRT( 1.0_wp - RL*2 ) sigma * rg3
1057
1058	particles(n)%speed_x = u_int(n) + particles(n)%rvar1
1059	particles(n)%speed_y = v_int(n) + particles(n)%rvar2
1060	particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
1061	ELSE
1062	particles(n)%speed_x = u_int(n)
1063	particles(n)%speed_y = v_int(n)
1064	particles(n)%speed_z = w_int(n) - w_s
1065	ENDIF
1066
1067	ELSE
1068
1069	IF ( use_sgs_for_particles ) THEN
1070	exp_arg = particle_groups(particles(n)%group)%exp_arg
1071	exp_term = EXP( -exp_arg * dt_particle(n) )
1072	ELSE
1073	exp_arg = particle_groups(particles(n)%group)%exp_arg
1074	exp_term = particle_groups(particles(n)%group)%exp_term
1075	ENDIF
1076	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
1077	u_int(n) * ( 1.0_wp - exp_term )
1078	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
1079	v_int(n) * ( 1.0_wp - exp_term )
1080	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
1081	( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
1082	exp_arg ) * ( 1.0_wp - exp_term )
1083	ENDIF
1084	ENDDO
1085	ENDDO
1086
1087	ENDIF
1088
1089	!
1090	!-- Store the old age of the particle ( needed to prevent that a
1091	!-- particle crosses several PEs during one timestep, and for the
1092	!-- evaluation of the subgrid particle velocity fluctuations )
1093	particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
1094
1095	DO nb = 0, 7
1096	DO n = start_index(nb), end_index(nb)
1097	!
1098	!-- Increment the particle age and the total time that the particle
1099	!-- has advanced within the particle timestep procedure
1100	particles(n)%age = particles(n)%age + dt_particle(n)
1101	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
1102
1103	!
1104	!-- Check whether there is still a particle that has not yet completed
1105	!-- the total LES timestep
1106	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) THEN
1107	dt_3d_reached_l = .FALSE.
1108	ENDIF
1109
1110	ENDDO
1111	ENDDO
1112
1113	CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
1114
1115
1116	END SUBROUTINE lpm_advec
1117
1118	! Description:
1119	! ------------
1120	!> Calculation of subgrid-scale particle speed using the stochastic model
1121	!> of Weil et al. (2004, JAS, 61, 2877-2887).
1122	!------------------------------------------------------------------------------!
1123	SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n, &
1124	dt_n, rg_n, fac )
1125
1126	USE kinds
1127
1128	USE particle_attributes, &
1129	ONLY: c_0, sgs_wf_part
1130
1131	IMPLICIT NONE
1132
1133	REAL(wp) :: a1 !< dummy argument
1134	REAL(wp) :: dedt_n !< time derivative of TKE at particle position
1135	REAL(wp) :: dedxi_n !< horizontal derivative of TKE at particle position
1136	REAL(wp) :: diss_n !< dissipation at particle position
1137	REAL(wp) :: dt_n !< particle timestep
1138	REAL(wp) :: e_n !< TKE at particle position
1139	REAL(wp) :: fac !< flag to identify adjacent topography
1140	REAL(wp) :: fs_n !< weighting factor to prevent that subgrid-scale particle speed becomes too large
1141	REAL(wp) :: sgs_w !< constant (1/3)
1142	REAL(wp) :: rg_n !< random number
1143	REAL(wp) :: term1 !< memory term
1144	REAL(wp) :: term2 !< drift correction term
1145	REAL(wp) :: term3 !< random term
1146	REAL(wp) :: v_sgs !< subgrid-scale velocity component
1147
1148	!-- At first, limit TKE to a small non-zero number, in order to prevent
1149	!-- the occurrence of extremely large SGS-velocities in case TKE is zero,
1150	!-- (could occur at the simulation begin).
1151	e_n = MAX( e_n, 1E-20_wp )
1152	!
1153	!-- Please note, terms 1 and 2 (drift and memory term, respectively) are
1154	!-- multiplied by a flag to switch of both terms near topography.
1155	!-- This is necessary, as both terms may cause a subgrid-scale velocity build up
1156	!-- if particles are trapped in regions with very small TKE, e.g. in narrow street
1157	!-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are
1158	!-- disabled if one of the adjacent grid points belongs to topography.
1159	!-- Moreover, in this case, the previous subgrid-scale component is also set
1160	!-- to zero.
1161
1162	a1 = fs_n * c_0 * diss_n
1163	!
1164	!-- Memory term
1165	term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp ) &
1166	* fac
1167	!
1168	!-- Drift correction term
1169	term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n &
1170	* fac
1171	!
1172	!-- Random term
1173	term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n )
1174	!
1175	!-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous
1176	!-- subgrid-scale velocity component is set to zero, in order to prevent a
1177	!-- velocity build-up.
1178	!-- This case, set also previous subgrid-scale component to zero.
1179	v_sgs = v_sgs * fac + term1 + term2 + term3
1180
1181	END SUBROUTINE weil_stochastic_eq

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |