Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |