Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |