Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |