Home

Context Navigation

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

Last change on this file since 3163 was 3065, checked in by Giersch, 6 years ago
New vertical stretching procedure has been introduced
Property svn:keywords set to `Id`
File size: 55.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |