Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |