Home

Context Navigation

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

Last change on this file since 1369 was 1369, checked in by raasch, 11 years ago
routine description added, usage of module interfaces removed
Property svn:keywords set to `Id`
File size: 61.1 KB

Rev	Line
[1359]	1	SUBROUTINE lpm_advec (ip,jp,kp)
[849]	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[849]	20	! Current revisions:
	21	! ------------------
[1369]	22	! usage of module interfaces removed
[1360]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: lpm_advec.f90 1369 2014-04-24 05:57:38Z raasch $
	27	!
[1360]	28	! 1359 2014-04-11 17:15:14Z hoffmann
	29	! New particle structure integrated.
	30	! Kind definition added to all floating point numbers.
	31	!
[1323]	32	! 1322 2014-03-20 16:38:49Z raasch
	33	! REAL constants defined as wp_kind
	34	!
[1321]	35	! 1320 2014-03-20 08:40:49Z raasch
[1320]	36	! ONLY-attribute added to USE-statements,
	37	! kind-parameters added to all INTEGER and REAL declaration statements,
	38	! kinds are defined in new module kinds,
	39	! revision history before 2012 removed,
	40	! comment fields (!:) to be used for variable explanations added to
	41	! all variable declaration statements
[849]	42	!
[1315]	43	! 1314 2014-03-14 18:25:17Z suehring
	44	! Vertical logarithmic interpolation of horizontal particle speed for particles
	45	! between roughness height and first vertical grid level.
	46	!
[1037]	47	! 1036 2012-10-22 13:43:42Z raasch
	48	! code put under GPL (PALM 3.9)
	49	!
[850]	50	! 849 2012-03-15 10:35:09Z raasch
	51	! initial revision (former part of advec_particles)
[849]	52	!
[850]	53	!
[849]	54	! Description:
	55	! ------------
	56	! Calculation of new particle positions due to advection using a simple Euler
	57	! scheme. Particles may feel inertia effects. SGS transport can be included
	58	! using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
	59	!------------------------------------------------------------------------------!
	60
[1320]	61	USE arrays_3d, &
[1359]	62	ONLY: de_dx, de_dy, de_dz, diss, e, u, us, usws, v, vsws, w, z0, zu, &
	63	zw
[849]	64
[1359]	65	USE cpulog
	66
	67	USE pegrid
	68
[1320]	69	USE control_parameters, &
	70	ONLY: atmos_ocean_sign, cloud_droplets, dt_3d, dt_3d_reached_l, dz, &
	71	g, kappa, molecular_viscosity, prandtl_layer, topography, &
[1359]	72	u_gtrans, v_gtrans, simulated_time
[849]	73
[1320]	74	USE grid_variables, &
	75	ONLY: ddx, dx, ddy, dy
	76
	77	USE indices, &
	78	ONLY: nzb, nzb_s_inner, nzt
	79
	80	USE kinds
	81
	82	USE particle_attributes, &
[1359]	83	ONLY: block_offset, c_0, density_ratio, dt_min_part, grid_particles, &
	84	iran_part, log_z_z0, number_of_particles, number_of_sublayers, &
	85	particles, particle_groups, offset_ocean_nzt, &
	86	offset_ocean_nzt_m1, sgs_wfu_part, sgs_wfv_part, sgs_wfw_part, &
	87	use_sgs_for_particles, vertical_particle_advection, z0_av_global
[1320]	88
	89	USE statistics, &
	90	ONLY: hom
[849]	91
[1320]	92	IMPLICIT NONE
[849]	93
[1320]	94	INTEGER(iwp) :: agp !:
	95	INTEGER(iwp) :: gp_outside_of_building(1:8) !:
	96	INTEGER(iwp) :: i !:
[1359]	97	INTEGER(iwp) :: ip !:
[1320]	98	INTEGER(iwp) :: j !:
[1359]	99	INTEGER(iwp) :: jp !:
[1320]	100	INTEGER(iwp) :: k !:
[1359]	101	INTEGER(iwp) :: kp !:
[1320]	102	INTEGER(iwp) :: kw !:
	103	INTEGER(iwp) :: n !:
[1359]	104	INTEGER(iwp) :: nb !:
[1320]	105	INTEGER(iwp) :: num_gp !:
[849]	106
[1359]	107	INTEGER(iwp), DIMENSION(0:7) :: start_index !:
	108	INTEGER(iwp), DIMENSION(0:7) :: end_index !:
	109
[1320]	110	REAL(wp) :: aa !:
	111	REAL(wp) :: bb !:
	112	REAL(wp) :: cc !:
	113	REAL(wp) :: d_sum !:
	114	REAL(wp) :: d_z_p_z0 !:
	115	REAL(wp) :: dd !:
	116	REAL(wp) :: de_dx_int_l !:
	117	REAL(wp) :: de_dx_int_u !:
	118	REAL(wp) :: de_dy_int_l !:
	119	REAL(wp) :: de_dy_int_u !:
	120	REAL(wp) :: de_dt !:
	121	REAL(wp) :: de_dt_min !:
	122	REAL(wp) :: de_dz_int_l !:
	123	REAL(wp) :: de_dz_int_u !:
	124	REAL(wp) :: diss_int_l !:
	125	REAL(wp) :: diss_int_u !:
	126	REAL(wp) :: dt_gap !:
	127	REAL(wp) :: dt_particle_m !:
	128	REAL(wp) :: e_int_l !:
	129	REAL(wp) :: e_int_u !:
	130	REAL(wp) :: e_mean_int !:
	131	REAL(wp) :: exp_arg !:
	132	REAL(wp) :: exp_term !:
	133	REAL(wp) :: gg !:
	134	REAL(wp) :: height_int !:
	135	REAL(wp) :: height_p !:
	136	REAL(wp) :: lagr_timescale !:
	137	REAL(wp) :: location(1:30,1:3) !:
	138	REAL(wp) :: random_gauss !:
	139	REAL(wp) :: u_int_l !:
	140	REAL(wp) :: u_int_u !:
	141	REAL(wp) :: us_int !:
	142	REAL(wp) :: v_int_l !:
	143	REAL(wp) :: v_int_u !:
	144	REAL(wp) :: vv_int !:
	145	REAL(wp) :: w_int_l !:
	146	REAL(wp) :: w_int_u !:
	147	REAL(wp) :: x !:
	148	REAL(wp) :: y !:
	149	REAL(wp) :: z_p !:
[849]	150
[1320]	151	REAL(wp), DIMENSION(1:30) :: d_gp_pl !:
	152	REAL(wp), DIMENSION(1:30) :: de_dxi !:
	153	REAL(wp), DIMENSION(1:30) :: de_dyi !:
	154	REAL(wp), DIMENSION(1:30) :: de_dzi !:
	155	REAL(wp), DIMENSION(1:30) :: dissi !:
	156	REAL(wp), DIMENSION(1:30) :: ei !:
[849]	157
[1359]	158	REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !:
	159	REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !:
	160	REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !:
	161	REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !:
	162	REAL(wp), DIMENSION(number_of_particles) :: diss_int !:
	163	REAL(wp), DIMENSION(number_of_particles) :: dt_particle !:
	164	REAL(wp), DIMENSION(number_of_particles) :: e_int !:
	165	REAL(wp), DIMENSION(number_of_particles) :: fs_int !:
	166	REAL(wp), DIMENSION(number_of_particles) :: log_z_z0_int !:
	167	REAL(wp), DIMENSION(number_of_particles) :: u_int !:
	168	REAL(wp), DIMENSION(number_of_particles) :: v_int !:
	169	REAL(wp), DIMENSION(number_of_particles) :: w_int !:
	170	REAL(wp), DIMENSION(number_of_particles) :: xv !:
	171	REAL(wp), DIMENSION(number_of_particles) :: yv !:
	172	REAL(wp), DIMENSION(number_of_particles) :: zv !:
	173
	174	REAL(wp), DIMENSION(number_of_particles, 3) :: rg !:
	175
	176	CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
	177
[1314]	178	!
	179	!-- Determine height of Prandtl layer and distance between Prandtl-layer
	180	!-- height and horizontal mean roughness height, which are required for
	181	!-- vertical logarithmic interpolation of horizontal particle speeds
	182	!-- (for particles below first vertical grid level).
	183	z_p = zu(nzb+1) - zw(nzb)
[1359]	184	d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
[849]	185
[1359]	186	start_index = grid_particles(kp,jp,ip)%start_index
	187	end_index = grid_particles(kp,jp,ip)%end_index
[849]	188
[1359]	189	xv = particles(1:number_of_particles)%x
	190	yv = particles(1:number_of_particles)%y
	191	zv = particles(1:number_of_particles)%z
[849]	192
[1359]	193	DO nb = 0, 7
[1314]	194
[1359]	195	i = ip
	196	j = jp + block_offset(nb)%j_off
	197	k = kp + block_offset(nb)%k_off
	198
[849]	199	!
[1359]	200	!-- Interpolate u velocity-component
	201	DO n = start_index(nb), end_index(nb)
[1314]	202	!
[1359]	203	!-- Interpolation of the u velocity component onto particle position.
	204	!-- Particles are interpolation bi-linearly in the horizontal and a
	205	!-- linearly in the vertical. An exception is made for particles below
	206	!-- the first vertical grid level in case of a prandtl layer. In this
	207	!-- case the horizontal particle velocity components are determined using
	208	!-- Monin-Obukhov relations (if branch).
	209	!-- First, check if particle is located below first vertical grid level
	210	!-- (Prandtl-layer height)
	211	IF ( prandtl_layer .AND. particles(n)%z < z_p ) THEN
[1314]	212	!
[1359]	213	!-- Resolved-scale horizontal particle velocity is zero below z0.
	214	IF ( particles(n)%z < z0_av_global ) THEN
	215	u_int(n) = 0.0_wp
	216	ELSE
[1314]	217	!
[1359]	218	!-- Determine the sublayer. Further used as index.
	219	height_p = ( particles(n)%z - z0_av_global ) &
	220	* REAL( number_of_sublayers, KIND=wp ) &
	221	* d_z_p_z0
[1314]	222	!
[1359]	223	!-- Calculate LOG(z/z0) for exact particle height. Therefore,
	224	!-- interpolate linearly between precalculated logarithm.
	225	log_z_z0_int(n) = log_z_z0(INT(height_p)) &
	226	+ ( height_p - INT(height_p) ) &
	227	* ( log_z_z0(INT(height_p)+1) &
	228	- log_z_z0(INT(height_p)) &
	229	)
[1314]	230	!
[1359]	231	!-- Neutral solution is applied for all situations, e.g. also for
	232	!-- unstable and stable situations. Even though this is not exact
	233	!-- this saves a lot of CPU time since several calls of intrinsic
	234	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
	235	!-- as sensitivity studies revealed no significant effect of
	236	!-- using the neutral solution also for un/stable situations.
	237	!-- Calculated left and bottom index on u grid.
	238	us_int = 0.5_wp * ( us(j,i) + us(j,i-1) )
[1314]	239
[1359]	240	u_int = -usws(j,i) / ( us_int * kappa + 1E-10_wp ) &
	241	* log_z_z0_int(n)
[1314]	242
[1359]	243	ENDIF
	244	!
	245	!-- Particle above the first grid level. Bi-linear interpolation in the
	246	!-- horizontal and linear interpolation in the vertical direction.
[1314]	247	ELSE
	248
[1359]	249	x = xv(n) + ( 0.5_wp - i ) * dx
	250	y = yv(n) - j * dy
	251	aa = x2 + y2
	252	bb = ( dx - x )2 + y2
	253	cc = x2 + ( dy - y )2
	254	dd = ( dx - x )2 + ( dy - y )2
	255	gg = aa + bb + cc + dd
[1314]	256
[1359]	257	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
	258	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * &
	259	u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
[1314]	260
[1359]	261	IF ( k == nzt ) THEN
	262	u_int(n) = u_int_l
	263	ELSE
	264	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
	265	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * &
	266	u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
	267	u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz * &
	268	( u_int_u - u_int_l )
	269	ENDIF
[1314]	270	ENDIF
	271
[1359]	272	ENDDO
[849]	273
[1359]	274	i = ip + block_offset(nb)%i_off
	275	j = jp
	276	k = kp + block_offset(nb)%k_off
[849]	277	!
[1359]	278	!-- Same procedure for interpolation of the v velocity-component
	279	DO n = start_index(nb), end_index(nb)
	280	IF ( prandtl_layer .AND. particles(n)%z < z_p ) THEN
[849]	281
[1359]	282	IF ( particles(n)%z < z0_av_global ) THEN
[1314]	283	!
[1359]	284	!-- Resolved-scale horizontal particle velocity is zero below z0.
	285	v_int(n) = 0.0_wp
	286	ELSE
	287	!
	288	!-- Neutral solution is applied for all situations, e.g. also for
	289	!-- unstable and stable situations. Even though this is not exact
	290	!-- this saves a lot of CPU time since several calls of intrinsic
	291	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
	292	!-- as sensitivity studies revealed no significant effect of
	293	!-- using the neutral solution also for un/stable situations.
	294	!-- Calculated left and bottom index on v grid.
	295	us_int = 0.5_wp * ( us(j,i) + us(j-1,i) )
[1314]	296
[1359]	297	v_int = -vsws(j,i) / ( us_int * kappa + 1E-10_wp ) &
	298	* log_z_z0_int(n)
	299	ENDIF
	300	ELSE
	301	x = xv(n) - i * dx
	302	y = yv(n) + ( 0.5_wp - j ) * dy
	303	aa = x2 + y2
	304	bb = ( dx - x )2 + y2
	305	cc = x2 + ( dy - y )2
	306	dd = ( dx - x )2 + ( dy - y )2
	307	gg = aa + bb + cc + dd
[1314]	308
[1359]	309	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
	310	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
	311	) / ( 3.0_wp * gg ) - v_gtrans
[1314]	312
[1359]	313	IF ( k == nzt ) THEN
	314	v_int(n) = v_int_l
	315	ELSE
	316	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
	317	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
	318	) / ( 3.0_wp * gg ) - v_gtrans
	319	v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
	320	( v_int_u - v_int_l )
	321	ENDIF
[1314]	322	ENDIF
	323
[1359]	324	ENDDO
[1314]	325
[1359]	326	i = ip + block_offset(nb)%i_off
	327	j = jp + block_offset(nb)%j_off
	328	k = kp-1
[849]	329	!
[1314]	330	!-- Same procedure for interpolation of the w velocity-component
[1359]	331	DO n = start_index(nb), end_index(nb)
[849]	332
[1359]	333	IF ( vertical_particle_advection(particles(n)%group) ) THEN
[849]	334
[1359]	335	x = xv(n) - i * dx
	336	y = yv(n) - j * dy
[849]	337	aa = x2 + y2
	338	bb = ( dx - x )2 + y2
	339	cc = x2 + ( dy - y )2
	340	dd = ( dx - x )2 + ( dy - y )2
	341	gg = aa + bb + cc + dd
	342
[1359]	343	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
	344	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
	345	) / ( 3.0_wp * gg )
[849]	346
[1359]	347	IF ( k == nzt ) THEN
	348	w_int(n) = w_int_l
[849]	349	ELSE
[1359]	350	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
	351	( gg-bb ) * w(k+1,j,i+1) + &
	352	( gg-cc ) * w(k+1,j+1,i) + &
	353	( gg-dd ) * w(k+1,j+1,i+1) &
	354	) / ( 3.0_wp * gg )
	355	w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
	356	( w_int_u - w_int_l )
[849]	357	ENDIF
	358
[1359]	359	ELSE
[849]	360
[1359]	361	w_int(n) = 0.0_wp
[849]	362
[1359]	363	ENDIF
	364
	365	ENDDO
	366
	367	ENDDO
	368
	369	!-- Interpolate and calculate quantities needed for calculating the SGS
	370	!-- velocities
	371	IF ( use_sgs_for_particles ) THEN
	372
	373	IF ( topography == 'flat' ) THEN
	374
	375	DO nb = 0,7
	376
	377	i = ip + block_offset(nb)%i_off
	378	j = jp + block_offset(nb)%j_off
	379	k = kp + block_offset(nb)%k_off
	380
	381	DO n = start_index(nb), end_index(nb)
[849]	382	!
[1359]	383	!-- Interpolate TKE
	384	x = xv(n) - i * dx
	385	y = yv(n) - j * dy
	386	aa = x2 + y2
	387	bb = ( dx - x )2 + y2
	388	cc = x2 + ( dy - y )2
	389	dd = ( dx - x )2 + ( dy - y )2
	390	gg = aa + bb + cc + dd
[849]	391
[1359]	392	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
	393	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
	394	) / ( 3.0_wp * gg )
	395
	396	IF ( k+1 == nzt+1 ) THEN
	397	e_int(n) = e_int_l
	398	ELSE
	399	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	400	( gg - bb ) * e(k+1,j,i+1) + &
	401	( gg - cc ) * e(k+1,j+1,i) + &
	402	( gg - dd ) * e(k+1,j+1,i+1) &
	403	) / ( 3.0_wp * gg )
	404	e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
	405	( e_int_u - e_int_l )
	406	ENDIF
[849]	407	!
[1359]	408	!-- Needed to avoid NaN particle velocities
	409	IF ( e_int(n) == 0.0_wp ) THEN
	410	e_int(n) = 1.0E-20_wp
	411	ENDIF
	412	!
	413	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
	414	!-- all position variables from above (TKE))
	415	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	416	( gg - bb ) * de_dx(k,j,i+1) + &
	417	( gg - cc ) * de_dx(k,j+1,i) + &
	418	( gg - dd ) * de_dx(k,j+1,i+1) &
	419	) / ( 3.0_wp * gg )
[849]	420
	421	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
[1359]	422	de_dx_int(n) = de_dx_int_l
[849]	423	ELSE
[1359]	424	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	425	( gg - bb ) * de_dx(k+1,j,i+1) + &
	426	( gg - cc ) * de_dx(k+1,j+1,i) + &
	427	( gg - dd ) * de_dx(k+1,j+1,i+1) &
	428	) / ( 3.0_wp * gg )
	429	de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
	430	( de_dx_int_u - de_dx_int_l )
[849]	431	ENDIF
[1359]	432	!
	433	!-- Interpolate the TKE gradient along y
	434	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	435	( gg - bb ) * de_dy(k,j,i+1) + &
	436	( gg - cc ) * de_dy(k,j+1,i) + &
	437	( gg - dd ) * de_dy(k,j+1,i+1) &
	438	) / ( 3.0_wp * gg )
	439	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	440	de_dy_int(n) = de_dy_int_l
	441	ELSE
	442	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	443	( gg - bb ) * de_dy(k+1,j,i+1) + &
	444	( gg - cc ) * de_dy(k+1,j+1,i) + &
	445	( gg - dd ) * de_dy(k+1,j+1,i+1) &
	446	) / ( 3.0_wp * gg )
	447	de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
	448	( de_dy_int_u - de_dy_int_l )
	449	ENDIF
[849]	450
	451	!
[1359]	452	!-- Interpolate the TKE gradient along z
	453	IF ( zv(n) < 0.5_wp * dz ) THEN
	454	de_dz_int(n) = 0.0_wp
	455	ELSE
	456	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	457	( gg - bb ) * de_dz(k,j,i+1) + &
	458	( gg - cc ) * de_dz(k,j+1,i) + &
	459	( gg - dd ) * de_dz(k,j+1,i+1) &
	460	) / ( 3.0_wp * gg )
[849]	461
[1359]	462	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	463	de_dz_int(n) = de_dz_int_l
	464	ELSE
	465	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	466	( gg - bb ) * de_dz(k+1,j,i+1) + &
	467	( gg - cc ) * de_dz(k+1,j+1,i) + &
	468	( gg - dd ) * de_dz(k+1,j+1,i+1) &
	469	) / ( 3.0_wp * gg )
	470	de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
	471	( de_dz_int_u - de_dz_int_l )
	472	ENDIF
	473	ENDIF
[849]	474
[1359]	475	!
	476	!-- Interpolate the dissipation of TKE
	477	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	478	( gg - bb ) * diss(k,j,i+1) + &
	479	( gg - cc ) * diss(k,j+1,i) + &
	480	( gg - dd ) * diss(k,j+1,i+1) &
	481	) / ( 3.0_wp * gg )
[849]	482
[1359]	483	IF ( k == nzt ) THEN
	484	diss_int(n) = diss_int_l
	485	ELSE
	486	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	487	( gg - bb ) * diss(k+1,j,i+1) + &
	488	( gg - cc ) * diss(k+1,j+1,i) + &
	489	( gg - dd ) * diss(k+1,j+1,i+1) &
	490	) / ( 3.0_wp * gg )
	491	diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
	492	( diss_int_u - diss_int_l )
	493	ENDIF
	494
	495	ENDDO
	496	ENDDO
	497
	498	ELSE ! non-flat topography, e.g., buildings
	499
	500	DO n = 1, number_of_particles
	501
	502	i = particles(n)%x * ddx
	503	j = particles(n)%y * ddy
	504	k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz &
	505	+ offset_ocean_nzt ! only exact if eq.dist
[849]	506	!
	507	!-- In case that there are buildings it has to be determined
	508	!-- how many of the gridpoints defining the particle box are
	509	!-- situated within a building
	510	!-- gp_outside_of_building(1): i,j,k
	511	!-- gp_outside_of_building(2): i,j+1,k
	512	!-- gp_outside_of_building(3): i,j,k+1
	513	!-- gp_outside_of_building(4): i,j+1,k+1
	514	!-- gp_outside_of_building(5): i+1,j,k
	515	!-- gp_outside_of_building(6): i+1,j+1,k
	516	!-- gp_outside_of_building(7): i+1,j,k+1
	517	!-- gp_outside_of_building(8): i+1,j+1,k+1
	518
	519	gp_outside_of_building = 0
[1359]	520	location = 0.0_wp
[849]	521	num_gp = 0
	522
	523	IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	524	num_gp = num_gp + 1
	525	gp_outside_of_building(1) = 1
	526	location(num_gp,1) = i * dx
	527	location(num_gp,2) = j * dy
[1359]	528	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	529	ei(num_gp) = e(k,j,i)
	530	dissi(num_gp) = diss(k,j,i)
	531	de_dxi(num_gp) = de_dx(k,j,i)
	532	de_dyi(num_gp) = de_dy(k,j,i)
	533	de_dzi(num_gp) = de_dz(k,j,i)
	534	ENDIF
	535
	536	IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
	537	THEN
	538	num_gp = num_gp + 1
	539	gp_outside_of_building(2) = 1
	540	location(num_gp,1) = i * dx
	541	location(num_gp,2) = (j+1) * dy
[1359]	542	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	543	ei(num_gp) = e(k,j+1,i)
	544	dissi(num_gp) = diss(k,j+1,i)
	545	de_dxi(num_gp) = de_dx(k,j+1,i)
	546	de_dyi(num_gp) = de_dy(k,j+1,i)
	547	de_dzi(num_gp) = de_dz(k,j+1,i)
	548	ENDIF
	549
	550	IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	551	num_gp = num_gp + 1
	552	gp_outside_of_building(3) = 1
	553	location(num_gp,1) = i * dx
	554	location(num_gp,2) = j * dy
[1359]	555	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	556	ei(num_gp) = e(k+1,j,i)
	557	dissi(num_gp) = diss(k+1,j,i)
	558	de_dxi(num_gp) = de_dx(k+1,j,i)
	559	de_dyi(num_gp) = de_dy(k+1,j,i)
	560	de_dzi(num_gp) = de_dz(k+1,j,i)
	561	ENDIF
	562
	563	IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
	564	THEN
	565	num_gp = num_gp + 1
	566	gp_outside_of_building(4) = 1
	567	location(num_gp,1) = i * dx
	568	location(num_gp,2) = (j+1) * dy
[1359]	569	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	570	ei(num_gp) = e(k+1,j+1,i)
	571	dissi(num_gp) = diss(k+1,j+1,i)
	572	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	573	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	574	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	575	ENDIF
	576
	577	IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
	578	THEN
	579	num_gp = num_gp + 1
	580	gp_outside_of_building(5) = 1
	581	location(num_gp,1) = (i+1) * dx
	582	location(num_gp,2) = j * dy
[1359]	583	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	584	ei(num_gp) = e(k,j,i+1)
	585	dissi(num_gp) = diss(k,j,i+1)
	586	de_dxi(num_gp) = de_dx(k,j,i+1)
	587	de_dyi(num_gp) = de_dy(k,j,i+1)
	588	de_dzi(num_gp) = de_dz(k,j,i+1)
	589	ENDIF
	590
[1359]	591	IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
[849]	592	THEN
	593	num_gp = num_gp + 1
	594	gp_outside_of_building(6) = 1
	595	location(num_gp,1) = (i+1) * dx
	596	location(num_gp,2) = (j+1) * dy
[1359]	597	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	598	ei(num_gp) = e(k,j+1,i+1)
	599	dissi(num_gp) = diss(k,j+1,i+1)
	600	de_dxi(num_gp) = de_dx(k,j+1,i+1)
	601	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	602	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	603	ENDIF
	604
	605	IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
	606	THEN
	607	num_gp = num_gp + 1
	608	gp_outside_of_building(7) = 1
	609	location(num_gp,1) = (i+1) * dx
	610	location(num_gp,2) = j * dy
[1359]	611	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	612	ei(num_gp) = e(k+1,j,i+1)
	613	dissi(num_gp) = diss(k+1,j,i+1)
	614	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	615	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	616	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	617	ENDIF
	618
[1359]	619	IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
[849]	620	THEN
	621	num_gp = num_gp + 1
	622	gp_outside_of_building(8) = 1
	623	location(num_gp,1) = (i+1) * dx
	624	location(num_gp,2) = (j+1) * dy
[1359]	625	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	626	ei(num_gp) = e(k+1,j+1,i+1)
	627	dissi(num_gp) = diss(k+1,j+1,i+1)
	628	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	629	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	630	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	631	ENDIF
	632
	633	!
	634	!-- If all gridpoints are situated outside of a building, then the
	635	!-- ordinary interpolation scheme can be used.
	636	IF ( num_gp == 8 ) THEN
	637
	638	x = particles(n)%x - i * dx
	639	y = particles(n)%y - j * dy
	640	aa = x2 + y2
	641	bb = ( dx - x )2 + y2
	642	cc = x2 + ( dy - y )2
	643	dd = ( dx - x )2 + ( dy - y )2
	644	gg = aa + bb + cc + dd
	645
[1359]	646	e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) &
	647	+ ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
	648	) / ( 3.0_wp * gg )
[849]	649
[1359]	650	IF ( k == nzt ) THEN
	651	e_int(n) = e_int_l
[849]	652	ELSE
	653	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	654	( gg - bb ) * e(k+1,j,i+1) + &
	655	( gg - cc ) * e(k+1,j+1,i) + &
	656	( gg - dd ) * e(k+1,j+1,i+1) &
[1359]	657	) / ( 3.0_wp * gg )
	658	e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
[849]	659	( e_int_u - e_int_l )
	660	ENDIF
	661	!
[1359]	662	!-- Needed to avoid NaN particle velocities
	663	IF ( e_int(n) == 0.0_wp ) THEN
	664	e_int(n) = 1.0E-20_wp
	665	ENDIF
	666	!
[849]	667	!-- Interpolate the TKE gradient along x (adopt incides i,j,k
	668	!-- and all position variables from above (TKE))
	669	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	670	( gg - bb ) * de_dx(k,j,i+1) + &
	671	( gg - cc ) * de_dx(k,j+1,i) + &
	672	( gg - dd ) * de_dx(k,j+1,i+1) &
[1359]	673	) / ( 3.0_wp * gg )
[849]	674
[1359]	675	IF ( ( k == nzt ) .OR. ( k == nzb ) ) THEN
	676	de_dx_int(n) = de_dx_int_l
[849]	677	ELSE
	678	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	679	( gg - bb ) * de_dx(k+1,j,i+1) + &
	680	( gg - cc ) * de_dx(k+1,j+1,i) + &
	681	( gg - dd ) * de_dx(k+1,j+1,i+1) &
[1359]	682	) / ( 3.0_wp * gg )
	683	de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
[849]	684	dz * ( de_dx_int_u - de_dx_int_l )
	685	ENDIF
	686
	687	!
	688	!-- Interpolate the TKE gradient along y
	689	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	690	( gg - bb ) * de_dy(k,j,i+1) + &
	691	( gg - cc ) * de_dy(k,j+1,i) + &
	692	( gg - dd ) * de_dy(k,j+1,i+1) &
[1359]	693	) / ( 3.0_wp * gg )
[849]	694	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
[1359]	695	de_dy_int(n) = de_dy_int_l
[849]	696	ELSE
	697	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	698	( gg - bb ) * de_dy(k+1,j,i+1) + &
	699	( gg - cc ) * de_dy(k+1,j+1,i) + &
	700	( gg - dd ) * de_dy(k+1,j+1,i+1) &
[1359]	701	) / ( 3.0_wp * gg )
	702	de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
[849]	703	dz * ( de_dy_int_u - de_dy_int_l )
	704	ENDIF
	705
	706	!
	707	!-- Interpolate the TKE gradient along z
[1359]	708	IF ( zv(n) < 0.5_wp * dz ) THEN
	709	de_dz_int(n) = 0.0_wp
[849]	710	ELSE
	711	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	712	( gg - bb ) * de_dz(k,j,i+1) + &
	713	( gg - cc ) * de_dz(k,j+1,i) + &
	714	( gg - dd ) * de_dz(k,j+1,i+1) &
[1359]	715	) / ( 3.0_wp * gg )
[849]	716
	717	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
[1359]	718	de_dz_int(n) = de_dz_int_l
[849]	719	ELSE
	720	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	721	( gg - bb ) * de_dz(k+1,j,i+1) + &
	722	( gg - cc ) * de_dz(k+1,j+1,i) + &
	723	( gg - dd ) * de_dz(k+1,j+1,i+1) &
[1359]	724	) / ( 3.0_wp * gg )
	725	de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
[849]	726	dz * ( de_dz_int_u - de_dz_int_l )
	727	ENDIF
	728	ENDIF
	729
	730	!
	731	!-- Interpolate the dissipation of TKE
	732	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	733	( gg - bb ) * diss(k,j,i+1) + &
	734	( gg - cc ) * diss(k,j+1,i) + &
	735	( gg - dd ) * diss(k,j+1,i+1) &
[1359]	736	) / ( 3.0_wp * gg )
[849]	737
[1359]	738	IF ( k == nzt ) THEN
	739	diss_int(n) = diss_int_l
[849]	740	ELSE
	741	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	742	( gg - bb ) * diss(k+1,j,i+1) + &
	743	( gg - cc ) * diss(k+1,j+1,i) + &
	744	( gg - dd ) * diss(k+1,j+1,i+1) &
[1359]	745	) / ( 3.0_wp * gg )
	746	diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
[849]	747	( diss_int_u - diss_int_l )
	748	ENDIF
	749
	750	ELSE
	751
	752	!
	753	!-- If wall between gridpoint 1 and gridpoint 5, then
	754	!-- Neumann boundary condition has to be applied
	755	IF ( gp_outside_of_building(1) == 1 .AND. &
	756	gp_outside_of_building(5) == 0 ) THEN
	757	num_gp = num_gp + 1
[1359]	758	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	759	location(num_gp,2) = j * dy
[1359]	760	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	761	ei(num_gp) = e(k,j,i)
	762	dissi(num_gp) = diss(k,j,i)
[1359]	763	de_dxi(num_gp) = 0.0_wp
[849]	764	de_dyi(num_gp) = de_dy(k,j,i)
	765	de_dzi(num_gp) = de_dz(k,j,i)
	766	ENDIF
	767
	768	IF ( gp_outside_of_building(5) == 1 .AND. &
	769	gp_outside_of_building(1) == 0 ) THEN
	770	num_gp = num_gp + 1
[1359]	771	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	772	location(num_gp,2) = j * dy
[1359]	773	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	774	ei(num_gp) = e(k,j,i+1)
	775	dissi(num_gp) = diss(k,j,i+1)
[1359]	776	de_dxi(num_gp) = 0.0_wp
[849]	777	de_dyi(num_gp) = de_dy(k,j,i+1)
	778	de_dzi(num_gp) = de_dz(k,j,i+1)
	779	ENDIF
	780
	781	!
	782	!-- If wall between gridpoint 5 and gridpoint 6, then
	783	!-- then Neumann boundary condition has to be applied
	784	IF ( gp_outside_of_building(5) == 1 .AND. &
	785	gp_outside_of_building(6) == 0 ) THEN
	786	num_gp = num_gp + 1
	787	location(num_gp,1) = (i+1) * dx
[1359]	788	location(num_gp,2) = j * dy + 0.5_wp * dy
	789	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	790	ei(num_gp) = e(k,j,i+1)
	791	dissi(num_gp) = diss(k,j,i+1)
	792	de_dxi(num_gp) = de_dx(k,j,i+1)
[1359]	793	de_dyi(num_gp) = 0.0_wp
[849]	794	de_dzi(num_gp) = de_dz(k,j,i+1)
	795	ENDIF
	796
	797	IF ( gp_outside_of_building(6) == 1 .AND. &
	798	gp_outside_of_building(5) == 0 ) THEN
	799	num_gp = num_gp + 1
	800	location(num_gp,1) = (i+1) * dx
[1359]	801	location(num_gp,2) = j * dy + 0.5_wp * dy
	802	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	803	ei(num_gp) = e(k,j+1,i+1)
	804	dissi(num_gp) = diss(k,j+1,i+1)
	805	de_dxi(num_gp) = de_dx(k,j+1,i+1)
[1359]	806	de_dyi(num_gp) = 0.0_wp
[849]	807	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	808	ENDIF
	809
	810	!
	811	!-- If wall between gridpoint 2 and gridpoint 6, then
	812	!-- Neumann boundary condition has to be applied
	813	IF ( gp_outside_of_building(2) == 1 .AND. &
	814	gp_outside_of_building(6) == 0 ) THEN
	815	num_gp = num_gp + 1
[1359]	816	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	817	location(num_gp,2) = (j+1) * dy
[1359]	818	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	819	ei(num_gp) = e(k,j+1,i)
	820	dissi(num_gp) = diss(k,j+1,i)
[1359]	821	de_dxi(num_gp) = 0.0_wp
[849]	822	de_dyi(num_gp) = de_dy(k,j+1,i)
	823	de_dzi(num_gp) = de_dz(k,j+1,i)
	824	ENDIF
	825
	826	IF ( gp_outside_of_building(6) == 1 .AND. &
	827	gp_outside_of_building(2) == 0 ) THEN
	828	num_gp = num_gp + 1
[1359]	829	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	830	location(num_gp,2) = (j+1) * dy
[1359]	831	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	832	ei(num_gp) = e(k,j+1,i+1)
	833	dissi(num_gp) = diss(k,j+1,i+1)
[1359]	834	de_dxi(num_gp) = 0.0_wp
[849]	835	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	836	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	837	ENDIF
	838
	839	!
	840	!-- If wall between gridpoint 1 and gridpoint 2, then
	841	!-- Neumann boundary condition has to be applied
	842	IF ( gp_outside_of_building(1) == 1 .AND. &
	843	gp_outside_of_building(2) == 0 ) THEN
	844	num_gp = num_gp + 1
	845	location(num_gp,1) = i * dx
[1359]	846	location(num_gp,2) = j * dy + 0.5_wp * dy
	847	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	848	ei(num_gp) = e(k,j,i)
	849	dissi(num_gp) = diss(k,j,i)
	850	de_dxi(num_gp) = de_dx(k,j,i)
[1359]	851	de_dyi(num_gp) = 0.0_wp
[849]	852	de_dzi(num_gp) = de_dz(k,j,i)
	853	ENDIF
	854
	855	IF ( gp_outside_of_building(2) == 1 .AND. &
	856	gp_outside_of_building(1) == 0 ) THEN
	857	num_gp = num_gp + 1
	858	location(num_gp,1) = i * dx
[1359]	859	location(num_gp,2) = j * dy + 0.5_wp * dy
	860	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	861	ei(num_gp) = e(k,j+1,i)
	862	dissi(num_gp) = diss(k,j+1,i)
	863	de_dxi(num_gp) = de_dx(k,j+1,i)
[1359]	864	de_dyi(num_gp) = 0.0_wp
[849]	865	de_dzi(num_gp) = de_dz(k,j+1,i)
	866	ENDIF
	867
	868	!
	869	!-- If wall between gridpoint 3 and gridpoint 7, then
	870	!-- Neumann boundary condition has to be applied
	871	IF ( gp_outside_of_building(3) == 1 .AND. &
	872	gp_outside_of_building(7) == 0 ) THEN
	873	num_gp = num_gp + 1
[1359]	874	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	875	location(num_gp,2) = j * dy
[1359]	876	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	877	ei(num_gp) = e(k+1,j,i)
	878	dissi(num_gp) = diss(k+1,j,i)
[1359]	879	de_dxi(num_gp) = 0.0_wp
[849]	880	de_dyi(num_gp) = de_dy(k+1,j,i)
	881	de_dzi(num_gp) = de_dz(k+1,j,i)
	882	ENDIF
	883
	884	IF ( gp_outside_of_building(7) == 1 .AND. &
	885	gp_outside_of_building(3) == 0 ) THEN
	886	num_gp = num_gp + 1
[1359]	887	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	888	location(num_gp,2) = j * dy
[1359]	889	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	890	ei(num_gp) = e(k+1,j,i+1)
	891	dissi(num_gp) = diss(k+1,j,i+1)
[1359]	892	de_dxi(num_gp) = 0.0_wp
[849]	893	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	894	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	895	ENDIF
	896
	897	!
	898	!-- If wall between gridpoint 7 and gridpoint 8, then
	899	!-- Neumann boundary condition has to be applied
	900	IF ( gp_outside_of_building(7) == 1 .AND. &
	901	gp_outside_of_building(8) == 0 ) THEN
	902	num_gp = num_gp + 1
	903	location(num_gp,1) = (i+1) * dx
[1359]	904	location(num_gp,2) = j * dy + 0.5_wp * dy
	905	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	906	ei(num_gp) = e(k+1,j,i+1)
	907	dissi(num_gp) = diss(k+1,j,i+1)
	908	de_dxi(num_gp) = de_dx(k+1,j,i+1)
[1359]	909	de_dyi(num_gp) = 0.0_wp
[849]	910	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	911	ENDIF
	912
	913	IF ( gp_outside_of_building(8) == 1 .AND. &
	914	gp_outside_of_building(7) == 0 ) THEN
	915	num_gp = num_gp + 1
	916	location(num_gp,1) = (i+1) * dx
[1359]	917	location(num_gp,2) = j * dy + 0.5_wp * dy
	918	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	919	ei(num_gp) = e(k+1,j+1,i+1)
	920	dissi(num_gp) = diss(k+1,j+1,i+1)
	921	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
[1359]	922	de_dyi(num_gp) = 0.0_wp
[849]	923	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	924	ENDIF
	925
	926	!
	927	!-- If wall between gridpoint 4 and gridpoint 8, then
	928	!-- Neumann boundary condition has to be applied
	929	IF ( gp_outside_of_building(4) == 1 .AND. &
	930	gp_outside_of_building(8) == 0 ) THEN
	931	num_gp = num_gp + 1
[1359]	932	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	933	location(num_gp,2) = (j+1) * dy
[1359]	934	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	935	ei(num_gp) = e(k+1,j+1,i)
	936	dissi(num_gp) = diss(k+1,j+1,i)
[1359]	937	de_dxi(num_gp) = 0.0_wp
[849]	938	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	939	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	940	ENDIF
	941
	942	IF ( gp_outside_of_building(8) == 1 .AND. &
	943	gp_outside_of_building(4) == 0 ) THEN
	944	num_gp = num_gp + 1
[1359]	945	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	946	location(num_gp,2) = (j+1) * dy
[1359]	947	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	948	ei(num_gp) = e(k+1,j+1,i+1)
	949	dissi(num_gp) = diss(k+1,j+1,i+1)
[1359]	950	de_dxi(num_gp) = 0.0_wp
[849]	951	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	952	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	953	ENDIF
	954
	955	!
	956	!-- If wall between gridpoint 3 and gridpoint 4, then
	957	!-- Neumann boundary condition has to be applied
	958	IF ( gp_outside_of_building(3) == 1 .AND. &
	959	gp_outside_of_building(4) == 0 ) THEN
	960	num_gp = num_gp + 1
	961	location(num_gp,1) = i * dx
[1359]	962	location(num_gp,2) = j * dy + 0.5_wp * dy
	963	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	964	ei(num_gp) = e(k+1,j,i)
	965	dissi(num_gp) = diss(k+1,j,i)
	966	de_dxi(num_gp) = de_dx(k+1,j,i)
[1359]	967	de_dyi(num_gp) = 0.0_wp
[849]	968	de_dzi(num_gp) = de_dz(k+1,j,i)
	969	ENDIF
	970
	971	IF ( gp_outside_of_building(4) == 1 .AND. &
	972	gp_outside_of_building(3) == 0 ) THEN
	973	num_gp = num_gp + 1
	974	location(num_gp,1) = i * dx
[1359]	975	location(num_gp,2) = j * dy + 0.5_wp * dy
	976	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	977	ei(num_gp) = e(k+1,j+1,i)
	978	dissi(num_gp) = diss(k+1,j+1,i)
	979	de_dxi(num_gp) = de_dx(k+1,j+1,i)
[1359]	980	de_dyi(num_gp) = 0.0_wp
[849]	981	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	982	ENDIF
	983
	984	!
	985	!-- If wall between gridpoint 1 and gridpoint 3, then
	986	!-- Neumann boundary condition has to be applied
	987	!-- (only one case as only building beneath is possible)
	988	IF ( gp_outside_of_building(1) == 0 .AND. &
	989	gp_outside_of_building(3) == 1 ) THEN
	990	num_gp = num_gp + 1
	991	location(num_gp,1) = i * dx
	992	location(num_gp,2) = j * dy
	993	location(num_gp,3) = k * dz
	994	ei(num_gp) = e(k+1,j,i)
	995	dissi(num_gp) = diss(k+1,j,i)
	996	de_dxi(num_gp) = de_dx(k+1,j,i)
	997	de_dyi(num_gp) = de_dy(k+1,j,i)
[1359]	998	de_dzi(num_gp) = 0.0_wp
[849]	999	ENDIF
	1000
	1001	!
	1002	!-- If wall between gridpoint 5 and gridpoint 7, then
	1003	!-- Neumann boundary condition has to be applied
	1004	!-- (only one case as only building beneath is possible)
	1005	IF ( gp_outside_of_building(5) == 0 .AND. &
	1006	gp_outside_of_building(7) == 1 ) THEN
	1007	num_gp = num_gp + 1
	1008	location(num_gp,1) = (i+1) * dx
	1009	location(num_gp,2) = j * dy
	1010	location(num_gp,3) = k * dz
	1011	ei(num_gp) = e(k+1,j,i+1)
	1012	dissi(num_gp) = diss(k+1,j,i+1)
	1013	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	1014	de_dyi(num_gp) = de_dy(k+1,j,i+1)
[1359]	1015	de_dzi(num_gp) = 0.0_wp
[849]	1016	ENDIF
	1017
	1018	!
	1019	!-- If wall between gridpoint 2 and gridpoint 4, then
	1020	!-- Neumann boundary condition has to be applied
	1021	!-- (only one case as only building beneath is possible)
	1022	IF ( gp_outside_of_building(2) == 0 .AND. &
	1023	gp_outside_of_building(4) == 1 ) THEN
	1024	num_gp = num_gp + 1
	1025	location(num_gp,1) = i * dx
	1026	location(num_gp,2) = (j+1) * dy
	1027	location(num_gp,3) = k * dz
	1028	ei(num_gp) = e(k+1,j+1,i)
	1029	dissi(num_gp) = diss(k+1,j+1,i)
	1030	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	1031	de_dyi(num_gp) = de_dy(k+1,j+1,i)
[1359]	1032	de_dzi(num_gp) = 0.0_wp
[849]	1033	ENDIF
	1034
	1035	!
	1036	!-- If wall between gridpoint 6 and gridpoint 8, then
	1037	!-- Neumann boundary condition has to be applied
	1038	!-- (only one case as only building beneath is possible)
	1039	IF ( gp_outside_of_building(6) == 0 .AND. &
	1040	gp_outside_of_building(8) == 1 ) THEN
	1041	num_gp = num_gp + 1
	1042	location(num_gp,1) = (i+1) * dx
	1043	location(num_gp,2) = (j+1) * dy
	1044	location(num_gp,3) = k * dz
	1045	ei(num_gp) = e(k+1,j+1,i+1)
	1046	dissi(num_gp) = diss(k+1,j+1,i+1)
	1047	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	1048	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
[1359]	1049	de_dzi(num_gp) = 0.0_wp
[849]	1050	ENDIF
	1051
	1052	!
	1053	!-- Carry out the interpolation
	1054	IF ( num_gp == 1 ) THEN
	1055	!
	1056	!-- If only one of the gridpoints is situated outside of the
	1057	!-- building, it follows that the values at the particle
	1058	!-- location are the same as the gridpoint values
[1359]	1059	e_int(n) = ei(num_gp)
	1060	diss_int(n) = dissi(num_gp)
	1061	de_dx_int(n) = de_dxi(num_gp)
	1062	de_dy_int(n) = de_dyi(num_gp)
	1063	de_dz_int(n) = de_dzi(num_gp)
[849]	1064	ELSE IF ( num_gp > 1 ) THEN
	1065
[1359]	1066	d_sum = 0.0_wp
[849]	1067	!
	1068	!-- Evaluation of the distances between the gridpoints
	1069	!-- contributing to the interpolated values, and the particle
	1070	!-- location
	1071	DO agp = 1, num_gp
	1072	d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 &
	1073	+ ( particles(n)%y-location(agp,2) )**2 &
[1359]	1074	+ ( zv(n)-location(agp,3) )**2
[849]	1075	d_sum = d_sum + d_gp_pl(agp)
	1076	ENDDO
	1077
	1078	!
	1079	!-- Finally the interpolation can be carried out
[1359]	1080	e_int(n) = 0.0_wp
	1081	diss_int(n) = 0.0_wp
	1082	de_dx_int(n) = 0.0_wp
	1083	de_dy_int(n) = 0.0_wp
	1084	de_dz_int(n) = 0.0_wp
[849]	1085	DO agp = 1, num_gp
[1359]	1086	e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1087	ei(agp) / ( (num_gp-1) * d_sum )
[1359]	1088	diss_int(n) = diss_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1089	dissi(agp) / ( (num_gp-1) * d_sum )
[1359]	1090	de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1091	de_dxi(agp) / ( (num_gp-1) * d_sum )
[1359]	1092	de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1093	de_dyi(agp) / ( (num_gp-1) * d_sum )
[1359]	1094	de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1095	de_dzi(agp) / ( (num_gp-1) * d_sum )
	1096	ENDDO
	1097
	1098	ENDIF
	1099
	1100	ENDIF
[1359]	1101	ENDDO
	1102	ENDIF
[849]	1103
[1359]	1104	DO nb = 0,7
	1105	i = ip + block_offset(nb)%i_off
	1106	j = jp + block_offset(nb)%j_off
	1107	k = kp + block_offset(nb)%k_off
[849]	1108
[1359]	1109	DO n = start_index(nb), end_index(nb)
[849]	1110	!
[1359]	1111	!-- Vertical interpolation of the horizontally averaged SGS TKE and
	1112	!-- resolved-scale velocity variances and use the interpolated values
	1113	!-- to calculate the coefficient fs, which is a measure of the ratio
	1114	!-- of the subgrid-scale turbulent kinetic energy to the total amount
	1115	!-- of turbulent kinetic energy.
	1116	IF ( k == 0 ) THEN
	1117	e_mean_int = hom(0,1,8,0)
	1118	ELSE
	1119	e_mean_int = hom(k,1,8,0) + &
	1120	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
	1121	( zu(k+1) - zu(k) ) * &
	1122	( zv(n) - zu(k) )
	1123	ENDIF
[849]	1124
[1359]	1125	! kw = particles(n)%z / dz
	1126	kw = kp-1 ! ok for ocean??? ( + offset_ocean_nzt_m1 ???)
[849]	1127
[1359]	1128	IF ( k == 0 ) THEN
	1129	aa = hom(k+1,1,30,0) * ( zv(n) / &
	1130	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
	1131	bb = hom(k+1,1,31,0) * ( zv(n) / &
	1132	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
	1133	cc = hom(kw+1,1,32,0) * ( zv(n) / &
	1134	( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
	1135	ELSE
	1136	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
	1137	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
	1138	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
	1139	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
	1140	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) * &
	1141	( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
	1142	ENDIF
[849]	1143
[1359]	1144	vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
	1145	!
	1146	!-- Needed to avoid NaN particle velocities. The value of 1.0 is just
	1147	!-- an educated guess for the given case.
	1148	IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp ) THEN
	1149	fs_int(n) = 1.0_wp
	1150	ELSE
	1151	fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int / &
	1152	( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
	1153	ENDIF
[849]	1154
[1359]	1155	ENDDO
	1156	ENDDO
[849]	1157
[1359]	1158	DO n = 1, number_of_particles
	1159
	1160	rg(n,1) = random_gauss( iran_part, 5.0_wp )
	1161	rg(n,2) = random_gauss( iran_part, 5.0_wp )
	1162	rg(n,3) = random_gauss( iran_part, 5.0_wp )
	1163
	1164	ENDDO
	1165
	1166	DO n = 1, number_of_particles
[849]	1167	!
	1168	!-- Calculate the Lagrangian timescale according to Weil et al. (2004).
[1359]	1169	lagr_timescale = ( 4.0_wp * e_int(n) ) / &
	1170	( 3.0_wp * fs_int(n) * c_0 * diss_int(n) )
[849]	1171
	1172	!
	1173	!-- Calculate the next particle timestep. dt_gap is the time needed to
	1174	!-- complete the current LES timestep.
	1175	dt_gap = dt_3d - particles(n)%dt_sum
[1359]	1176	dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
[849]	1177
	1178	!
	1179	!-- The particle timestep should not be too small in order to prevent
	1180	!-- the number of particle timesteps of getting too large
[1359]	1181	IF ( dt_particle(n) < dt_min_part .AND. dt_min_part < dt_gap ) THEN
	1182	dt_particle(n) = dt_min_part
[849]	1183	ENDIF
	1184
	1185	!
	1186	!-- Calculate the SGS velocity components
[1359]	1187	IF ( particles(n)%age == 0.0_wp ) THEN
[849]	1188	!
	1189	!-- For new particles the SGS components are derived from the SGS
	1190	!-- TKE. Limit the Gaussian random number to the interval
	1191	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
	1192	!-- from becoming unrealistically large.
[1359]	1193	particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wfu_part * e_int(n) ) * &
	1194	( rg(n,1) - 1.0_wp )
	1195	particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wfv_part * e_int(n) ) * &
	1196	( rg(n,2) - 1.0_wp )
	1197	particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wfw_part * e_int(n) ) * &
	1198	( rg(n,3) - 1.0_wp )
[849]	1199
	1200	ELSE
	1201	!
	1202	!-- Restriction of the size of the new timestep: compared to the
	1203	!-- previous timestep the increase must not exceed 200%
	1204
	1205	dt_particle_m = particles(n)%age - particles(n)%age_m
[1359]	1206	IF ( dt_particle(n) > 2.0_wp * dt_particle_m ) THEN
	1207	dt_particle(n) = 2.0_wp * dt_particle_m
[849]	1208	ENDIF
	1209
	1210	!
	1211	!-- For old particles the SGS components are correlated with the
	1212	!-- values from the previous timestep. Random numbers have also to
	1213	!-- be limited (see above).
	1214	!-- As negative values for the subgrid TKE are not allowed, the
	1215	!-- change of the subgrid TKE with time cannot be smaller than
[1359]	1216	!-- -e_int(n)/dt_particle. This value is used as a lower boundary
[849]	1217	!-- value for the change of TKE
	1218
[1359]	1219	de_dt_min = - e_int(n) / dt_particle(n)
[849]	1220
[1359]	1221	de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
[849]	1222
	1223	IF ( de_dt < de_dt_min ) THEN
	1224	de_dt = de_dt_min
	1225	ENDIF
	1226
[1359]	1227	particles(n)%rvar1 = particles(n)%rvar1 - fs_int(n) * c_0 * &
	1228	diss_int(n) * particles(n)%rvar1 * dt_particle(n) / &
	1229	( 4.0_wp * sgs_wfu_part * e_int(n) ) + &
	1230	( 2.0_wp * sgs_wfu_part * de_dt * &
	1231	particles(n)%rvar1 / &
	1232	( 2.0_wp * sgs_wfu_part * e_int(n) ) + &
	1233	de_dx_int(n) &
	1234	) * dt_particle(n) / 2.0_wp + &
	1235	SQRT( fs_int(n) * c_0 * diss_int(n) ) * &
	1236	( rg(n,1) - 1.0_wp ) * &
	1237	SQRT( dt_particle(n) )
[849]	1238
[1359]	1239	particles(n)%rvar2 = particles(n)%rvar2 - fs_int(n) * c_0 * &
	1240	diss_int(n) * particles(n)%rvar2 * dt_particle(n) / &
	1241	( 4.0_wp * sgs_wfv_part * e_int(n) ) + &
	1242	( 2.0_wp * sgs_wfv_part * de_dt * &
	1243	particles(n)%rvar2 / &
	1244	( 2.0_wp * sgs_wfv_part * e_int(n) ) + &
	1245	de_dy_int(n) &
	1246	) * dt_particle(n) / 2.0_wp + &
	1247	SQRT( fs_int(n) * c_0 * diss_int(n) ) * &
	1248	( rg(n,2) - 1.0_wp ) * &
	1249	SQRT( dt_particle(n) )
[849]	1250
[1359]	1251	particles(n)%rvar3 = particles(n)%rvar3 - fs_int(n) * c_0 * &
	1252	diss_int(n) * particles(n)%rvar3 * dt_particle(n) / &
	1253	( 4.0_wp * sgs_wfw_part * e_int(n) ) + &
	1254	( 2.0_wp * sgs_wfw_part * de_dt * &
	1255	particles(n)%rvar3 / &
	1256	( 2.0_wp * sgs_wfw_part * e_int(n) ) + &
	1257	de_dz_int(n) &
	1258	) * dt_particle(n) / 2.0_wp + &
	1259	SQRT( fs_int(n) * c_0 * diss_int(n) ) * &
	1260	( rg(n,3) - 1.0_wp ) * &
	1261	SQRT( dt_particle(n) )
[849]	1262
	1263	ENDIF
[1359]	1264	u_int(n) = u_int(n) + particles(n)%rvar1
	1265	v_int(n) = v_int(n) + particles(n)%rvar2
	1266	w_int(n) = w_int(n) + particles(n)%rvar3
[849]	1267
	1268	!
	1269	!-- Store the SGS TKE of the current timelevel which is needed for
	1270	!-- for calculating the SGS particle velocities at the next timestep
[1359]	1271	particles(n)%e_m = e_int(n)
	1272	ENDDO
[849]	1273
[1359]	1274	ELSE
[849]	1275	!
[1359]	1276	!-- If no SGS velocities are used, only the particle timestep has to
	1277	!-- be set
	1278	dt_particle = dt_3d
[849]	1279
[1359]	1280	ENDIF
[849]	1281	!
[1359]	1282	!-- Store the old age of the particle ( needed to prevent that a
	1283	!-- particle crosses several PEs during one timestep, and for the
	1284	!-- evaluation of the subgrid particle velocity fluctuations )
	1285	particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
[849]	1286
[1359]	1287	dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
[849]	1288
[1359]	1289	IF ( ANY( dens_ratio == 0.0_wp ) ) THEN
	1290	DO n = 1, number_of_particles
	1291
[849]	1292	!
[1359]	1293	!-- Particle advection
	1294	IF ( dens_ratio(n) == 0.0_wp ) THEN
[849]	1295	!
[1359]	1296	!-- Pure passive transport (without particle inertia)
	1297	particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
	1298	particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
	1299	particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
[849]	1300
[1359]	1301	particles(n)%speed_x = u_int(n)
	1302	particles(n)%speed_y = v_int(n)
	1303	particles(n)%speed_z = w_int(n)
[849]	1304
[1359]	1305	ELSE
[849]	1306	!
[1359]	1307	!-- Transport of particles with inertia
	1308	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
	1309	dt_particle(n)
	1310	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
	1311	dt_particle(n)
	1312	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
	1313	dt_particle(n)
[849]	1314
	1315	!
[1359]	1316	!-- Update of the particle velocity
	1317	IF ( cloud_droplets ) THEN
	1318	exp_arg = 4.5_wp * dens_ratio(n) * molecular_viscosity / &
	1319	( particles(n)%radius )*2 &
	1320	( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius &
	1321	* SQRT( ( u_int(n) - particles(n)%speed_x )**2 + &
	1322	( v_int(n) - particles(n)%speed_y )**2 + &
	1323	( w_int(n) - particles(n)%speed_z )**2 ) &
	1324	/ molecular_viscosity )**0.687_wp &
	1325	)
	1326
	1327	exp_term = EXP( -exp_arg * dt_particle(n) )
	1328	ELSEIF ( use_sgs_for_particles ) THEN
	1329	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1330	exp_term = EXP( -exp_arg * dt_particle(n) )
	1331	ELSE
	1332	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1333	exp_term = particle_groups(particles(n)%group)%exp_term
	1334	ENDIF
	1335	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
	1336	u_int(n) * ( 1.0_wp - exp_term )
	1337	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
	1338	v_int(n) * ( 1.0_wp - exp_term )
	1339	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
	1340	( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
	1341	g / exp_arg ) * ( 1.0_wp - exp_term )
	1342	ENDIF
	1343
	1344	ENDDO
	1345
	1346	ELSE
	1347
	1348	DO n = 1, number_of_particles
	1349
	1350	!-- Transport of particles with inertia
	1351	particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
	1352	particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
	1353	particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
	1354	!
[849]	1355	!-- Update of the particle velocity
	1356	IF ( cloud_droplets ) THEN
[1359]	1357
	1358	exp_arg = 4.5_wp * dens_ratio(n) * molecular_viscosity / &
	1359	( particles(n)%radius )*2 &
	1360	( 1.0_wp + 0.15_wp * ( 2.0_wp * particles(n)%radius * &
	1361	SQRT( ( u_int(n) - particles(n)%speed_x )**2 + &
	1362	( v_int(n) - particles(n)%speed_y )**2 + &
	1363	( w_int(n) - particles(n)%speed_z )**2 ) / &
	1364	molecular_viscosity )**0.687_wp &
[849]	1365	)
[1359]	1366
	1367	exp_term = EXP( -exp_arg * dt_particle(n) )
[849]	1368	ELSEIF ( use_sgs_for_particles ) THEN
	1369	exp_arg = particle_groups(particles(n)%group)%exp_arg
[1359]	1370	exp_term = EXP( -exp_arg * dt_particle(n) )
[849]	1371	ELSE
	1372	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1373	exp_term = particle_groups(particles(n)%group)%exp_term
	1374	ENDIF
	1375	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
[1359]	1376	u_int(n) * ( 1.0_wp - exp_term )
[849]	1377	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
[1359]	1378	v_int(n) * ( 1.0_wp - exp_term )
[849]	1379	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
[1359]	1380	( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
	1381	exp_arg ) * ( 1.0_wp - exp_term )
	1382	ENDDO
[849]	1383
[1359]	1384	ENDIF
	1385
	1386	DO n = 1, number_of_particles
[849]	1387	!
	1388	!-- Increment the particle age and the total time that the particle
	1389	!-- has advanced within the particle timestep procedure
[1359]	1390	particles(n)%age = particles(n)%age + dt_particle(n)
	1391	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
[849]	1392
	1393	!
	1394	!-- Check whether there is still a particle that has not yet completed
	1395	!-- the total LES timestep
[1359]	1396	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) THEN
[849]	1397	dt_3d_reached_l = .FALSE.
	1398	ENDIF
	1399
	1400	ENDDO
	1401
[1359]	1402	CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
[849]	1403
	1404	END SUBROUTINE lpm_advec

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |