Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |