Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |