Home

Context Navigation

source: palm/trunk/SOURCE/advec_ws.f90 @ 2292

Last change on this file since 2292 was 2292, checked in by schwenkel, 8 years ago
implementation of new bulk microphysics scheme
Property svn:keywords set to `Id`
File size: 278.6 KB

Line
1	!> @file advec_ws.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: advec_ws.f90 2292 2017-06-20 09:51:42Z schwenkel $
27	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
28	! includes two more prognostic equations for cloud drop concentration (nc)
29	! and cloud water content (qc).
30	!
31	! 2233 2017-05-30 18:08:54Z suehring
32	!
33	! 2232 2017-05-30 17:47:52Z suehring
34	! Rename wall_flags_0 and wall_flags_00 into advc_flags_1 and advc_flags_2,
35	! respectively.
36	! Set advc_flags_1/2 on basis of wall_flags_0/00 instead of nzb_s/u/v/w_inner.
37	! Setting advc_flags_1/2 also for downward-facing walls
38	!
39	! 2200 2017-04-11 11:37:51Z suehring
40	! monotonic_adjustment removed
41	!
42	! 2118 2017-01-17 16:38:49Z raasch
43	! OpenACC version of subroutines removed
44	!
45	! 2037 2016-10-26 11:15:40Z knoop
46	! Anelastic approximation implemented
47	!
48	! 2000 2016-08-20 18:09:15Z knoop
49	! Forced header and separation lines into 80 columns
50	!
51	! 1996 2016-08-18 11:42:29Z suehring
52	! Bugfix concerning calculation of turbulent of turbulent fluxes
53	!
54	! 1960 2016-07-12 16:34:24Z suehring
55	! Separate humidity and passive scalar
56	!
57	! 1942 2016-06-14 12:18:18Z suehring
58	! Initialization of flags for ws-scheme moved from init_grid.
59	!
60	! 1873 2016-04-18 14:50:06Z maronga
61	! Module renamed (removed _mod)
62	!
63	!
64	! 1850 2016-04-08 13:29:27Z maronga
65	! Module renamed
66	!
67	!
68	! 1822 2016-04-07 07:49:42Z hoffmann
69	! icloud_scheme removed, microphysics_seifert added
70	!
71	! 1682 2015-10-07 23:56:08Z knoop
72	! Code annotations made doxygen readable
73	!
74	! 1630 2015-08-26 16:57:23Z suehring
75	!
76	!
77	! 1629 2015-08-26 16:56:11Z suehring
78	! Bugfix concerning wall_flags at left and south PE boundaries
79	!
80	! 1581 2015-04-10 13:45:59Z suehring
81	!
82	!
83	! 1580 2015-04-10 13:43:49Z suehring
84	! Bugfix: statistical evaluation of scalar fluxes in case of monotonic limiter
85	!
86	! 1567 2015-03-10 17:57:55Z suehring
87	! Bugfixes in monotonic limiter.
88	!
89	! 2015-03-09 13:10:37Z heinze
90	! Bugfix: REAL constants provided with KIND-attribute in call of
91	! intrinsic functions like MAX and MIN
92	!
93	! 1557 2015-03-05 16:43:04Z suehring
94	! Enable monotone advection for scalars using monotonic limiter
95	!
96	! 1374 2014-04-25 12:55:07Z raasch
97	! missing variables added to ONLY list
98	!
99	! 1361 2014-04-16 15:17:48Z hoffmann
100	! accelerator and vector version for qr and nr added
101	!
102	! 1353 2014-04-08 15:21:23Z heinze
103	! REAL constants provided with KIND-attribute,
104	! module kinds added
105	! some formatting adjustments
106	!
107	! 1322 2014-03-20 16:38:49Z raasch
108	! REAL constants defined as wp-kind
109	!
110	! 1320 2014-03-20 08:40:49Z raasch
111	! ONLY-attribute added to USE-statements,
112	! kind-parameters added to all INTEGER and REAL declaration statements,
113	! kinds are defined in new module kinds,
114	! old module precision_kind is removed,
115	! revision history before 2012 removed,
116	! comment fields (!:) to be used for variable explanations added to
117	! all variable declaration statements
118	!
119	! 1257 2013-11-08 15:18:40Z raasch
120	! accelerator loop directives removed
121	!
122	! 1221 2013-09-10 08:59:13Z raasch
123	! wall_flags_00 introduced, which holds bits 32-...
124	!
125	! 1128 2013-04-12 06:19:32Z raasch
126	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
127	! j_north
128	!
129	! 1115 2013-03-26 18:16:16Z hoffmann
130	! calculation of qr and nr is restricted to precipitation
131	!
132	! 1053 2012-11-13 17:11:03Z hoffmann
133	! necessary expansions according to the two new prognostic equations (nr, qr)
134	! of the two-moment cloud physics scheme:
135	! +flux_l_, flux_s_, diss_l_, diss_s_, sums_ws*s_ws_l
136	!
137	! 1036 2012-10-22 13:43:42Z raasch
138	! code put under GPL (PALM 3.9)
139	!
140	! 1027 2012-10-15 17:18:39Z suehring
141	! Bugfix in calculation indices k_mm, k_pp in accelerator version
142	!
143	! 1019 2012-09-28 06:46:45Z raasch
144	! small change in comment lines
145	!
146	! 1015 2012-09-27 09:23:24Z raasch
147	! accelerator versions (*_acc) added
148	!
149	! 1010 2012-09-20 07:59:54Z raasch
150	! cpp switch __nopointer added for pointer free version
151	!
152	! 888 2012-04-20 15:03:46Z suehring
153	! Number of IBITS() calls with identical arguments is reduced.
154	!
155	! 862 2012-03-26 14:21:38Z suehring
156	! ws-scheme also work with topography in combination with vector version.
157	! ws-scheme also work with outflow boundaries in combination with
158	! vector version.
159	! Degradation of the applied order of scheme is now steered by multiplying with
160	! Integer advc_flags_1. 2nd order scheme, WS3 and WS5 are calculated on each
161	! grid point and mulitplied with the appropriate flag.
162	! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
163	! term derived according to the 4th and 6th order terms is applied. It turns
164	! out that diss_2nd does not provide sufficient dissipation near walls.
165	! Therefore, the function diss_2nd is removed.
166	! Near walls a divergence correction is necessary to overcome numerical
167	! instabilities due to too less divergence reduction of the velocity field.
168	! boundary_flags and logicals steering the degradation are removed.
169	! Empty SUBROUTINE local_diss removed.
170	! Further formatting adjustments.
171	!
172	! 801 2012-01-10 17:30:36Z suehring
173	! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
174	! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
175	! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
176	! dimension.
177	!
178	! Initial revision
179	!
180	! 411 2009-12-11 12:31:43 Z suehring
181	!
182	! Description:
183	! ------------
184	!> Advection scheme for scalars and momentum using the flux formulation of
185	!> Wicker and Skamarock 5th order. Additionally the module contains of a
186	!> routine using for initialisation and steering of the statical evaluation.
187	!> The computation of turbulent fluxes takes place inside the advection
188	!> routines.
189	!> Near non-cyclic boundaries the order of the applied advection scheme is
190	!> degraded.
191	!> A divergence correction is applied. It is necessary for topography, since
192	!> the divergence is not sufficiently reduced, resulting in erroneous fluxes and
193	!> partly numerical instabilities.
194	!-----------------------------------------------------------------------------!
195	MODULE advec_ws
196
197
198
199	PRIVATE
200	PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init, &
201	ws_init_flags, ws_statistics
202
203	INTERFACE ws_init
204	MODULE PROCEDURE ws_init
205	END INTERFACE ws_init
206
207	INTERFACE ws_init_flags
208	MODULE PROCEDURE ws_init_flags
209	END INTERFACE ws_init_flags
210
211	INTERFACE ws_statistics
212	MODULE PROCEDURE ws_statistics
213	END INTERFACE ws_statistics
214
215	INTERFACE advec_s_ws
216	MODULE PROCEDURE advec_s_ws
217	MODULE PROCEDURE advec_s_ws_ij
218	END INTERFACE advec_s_ws
219
220	INTERFACE advec_u_ws
221	MODULE PROCEDURE advec_u_ws
222	MODULE PROCEDURE advec_u_ws_ij
223	END INTERFACE advec_u_ws
224
225	INTERFACE advec_v_ws
226	MODULE PROCEDURE advec_v_ws
227	MODULE PROCEDURE advec_v_ws_ij
228	END INTERFACE advec_v_ws
229
230	INTERFACE advec_w_ws
231	MODULE PROCEDURE advec_w_ws
232	MODULE PROCEDURE advec_w_ws_ij
233	END INTERFACE advec_w_ws
234
235	CONTAINS
236
237
238	!------------------------------------------------------------------------------!
239	! Description:
240	! ------------
241	!> Initialization of WS-scheme
242	!------------------------------------------------------------------------------!
243	SUBROUTINE ws_init
244
245	USE arrays_3d, &
246	ONLY: diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q, &
247	diss_l_qc, diss_l_qr, diss_l_s, diss_l_sa, diss_l_u, &
248	diss_l_v, diss_l_w, flux_l_e, flux_l_nc, flux_l_nr, &
249	flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s, &
250	flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e, &
251	diss_s_nc, diss_s_nr, diss_s_pt, diss_s_q, diss_s_qc, &
252	diss_s_qr, diss_s_s, diss_s_sa, diss_s_u, diss_s_v, &
253	diss_s_w, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, &
254	flux_s_q, flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, &
255	flux_s_u, flux_s_v, flux_s_w
256
257	USE constants, &
258	ONLY: adv_mom_1, adv_mom_3, adv_mom_5, adv_sca_1, adv_sca_3, &
259	adv_sca_5
260
261	USE control_parameters, &
262	ONLY: cloud_physics, humidity, loop_optimization, &
263	passive_scalar, microphysics_morrison, microphysics_seifert, &
264	ocean, ws_scheme_mom, ws_scheme_sca
265
266	USE indices, &
267	ONLY: nyn, nys, nzb, nzt
268
269	USE kinds
270
271	USE pegrid
272
273	USE statistics, &
274	ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
275	sums_wsnrs_ws_l,sums_wspts_ws_l, sums_wsqcs_ws_l, &
276	sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l, &
277	sums_wssas_ws_l, sums_wsss_ws_l, sums_wsus_ws_l, &
278	sums_wsvs_ws_l
279
280
281	!
282	!-- Set the appropriate factors for scalar and momentum advection.
283	adv_sca_5 = 1.0_wp / 60.0_wp
284	adv_sca_3 = 1.0_wp / 12.0_wp
285	adv_sca_1 = 1.0_wp / 2.0_wp
286	adv_mom_5 = 1.0_wp / 120.0_wp
287	adv_mom_3 = 1.0_wp / 24.0_wp
288	adv_mom_1 = 1.0_wp / 4.0_wp
289	!
290	!-- Arrays needed for statical evaluation of fluxes.
291	IF ( ws_scheme_mom ) THEN
292
293	ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), &
294	sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), &
295	sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
296	sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
297	sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
298
299	sums_wsus_ws_l = 0.0_wp
300	sums_wsvs_ws_l = 0.0_wp
301	sums_us2_ws_l = 0.0_wp
302	sums_vs2_ws_l = 0.0_wp
303	sums_ws2_ws_l = 0.0_wp
304
305	ENDIF
306
307	IF ( ws_scheme_sca ) THEN
308
309	ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
310	sums_wspts_ws_l = 0.0_wp
311
312	IF ( humidity ) THEN
313	ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
314	sums_wsqs_ws_l = 0.0_wp
315	ENDIF
316
317	IF ( passive_scalar ) THEN
318	ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
319	sums_wsss_ws_l = 0.0_wp
320	ENDIF
321
322	IF ( cloud_physics .AND. microphysics_morrison ) THEN
323	ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
324	ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
325	sums_wsqcs_ws_l = 0.0_wp
326	sums_wsncs_ws_l = 0.0_wp
327	ENDIF
328
329	IF ( cloud_physics .AND. microphysics_seifert ) THEN
330	ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
331	ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
332	sums_wsqrs_ws_l = 0.0_wp
333	sums_wsnrs_ws_l = 0.0_wp
334	ENDIF
335
336	IF ( ocean ) THEN
337	ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
338	sums_wssas_ws_l = 0.0_wp
339	ENDIF
340
341	ENDIF
342
343	!
344	!-- Arrays needed for reasons of speed optimization for cache version.
345	!-- For the vector version the buffer arrays are not necessary,
346	!-- because the the fluxes can swapped directly inside the loops of the
347	!-- advection routines.
348	IF ( loop_optimization /= 'vector' ) THEN
349
350	IF ( ws_scheme_mom ) THEN
351
352	ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &
353	flux_s_v(nzb+1:nzt,0:threads_per_task-1), &
354	flux_s_w(nzb+1:nzt,0:threads_per_task-1), &
355	diss_s_u(nzb+1:nzt,0:threads_per_task-1), &
356	diss_s_v(nzb+1:nzt,0:threads_per_task-1), &
357	diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
358	ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
359	flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
360	flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
361	diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
362	diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
363	diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
364
365	ENDIF
366
367	IF ( ws_scheme_sca ) THEN
368
369	ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), &
370	flux_s_e(nzb+1:nzt,0:threads_per_task-1), &
371	diss_s_pt(nzb+1:nzt,0:threads_per_task-1), &
372	diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
373	ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
374	flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
375	diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
376	diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
377
378	IF ( humidity ) THEN
379	ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), &
380	diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
381	ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
382	diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
383	ENDIF
384
385	IF ( passive_scalar ) THEN
386	ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1), &
387	diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
388	ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
389	diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
390	ENDIF
391
392	IF ( cloud_physics .AND. microphysics_morrison ) THEN
393	ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1), &
394	diss_s_qc(nzb+1:nzt,0:threads_per_task-1), &
395	flux_s_nc(nzb+1:nzt,0:threads_per_task-1), &
396	diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
397	ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
398	diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
399	flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
400	diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
401	ENDIF
402
403	IF ( cloud_physics .AND. microphysics_seifert ) THEN
404	ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), &
405	diss_s_qr(nzb+1:nzt,0:threads_per_task-1), &
406	flux_s_nr(nzb+1:nzt,0:threads_per_task-1), &
407	diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
408	ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
409	diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
410	flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
411	diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
412	ENDIF
413
414	IF ( ocean ) THEN
415	ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), &
416	diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
417	ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
418	diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
419	ENDIF
420
421	ENDIF
422
423	ENDIF
424
425	END SUBROUTINE ws_init
426
427	!------------------------------------------------------------------------------!
428	! Description:
429	! ------------
430	!> Initialization of flags for WS-scheme used to degrade the order of the scheme
431	!> near walls.
432	!------------------------------------------------------------------------------!
433	SUBROUTINE ws_init_flags
434
435	USE control_parameters, &
436	ONLY: inflow_l, inflow_n, inflow_r, inflow_s, momentum_advec, &
437	nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, &
438	outflow_l, outflow_n, outflow_r, outflow_s, scalar_advec
439
440	USE indices, &
441	ONLY: advc_flags_1, advc_flags_2, nbgp, nxl, nxlu, nxr, nyn, nys, &
442	nysv, nzb, nzt, wall_flags_0
443
444	USE kinds
445
446	IMPLICIT NONE
447
448	INTEGER(iwp) :: i !< index variable along x
449	INTEGER(iwp) :: j !< index variable along y
450	INTEGER(iwp) :: k !< index variable along z
451	INTEGER(iwp) :: k_mm !< dummy index along z
452	INTEGER(iwp) :: k_pp !< dummy index along z
453	INTEGER(iwp) :: k_ppp !< dummy index along z
454
455	LOGICAL :: flag_set !< steering variable for advection flags
456
457
458	IF ( scalar_advec == 'ws-scheme' ) THEN
459	!
460	!-- Set flags to steer the degradation of the advection scheme in advec_ws
461	!-- near topography, inflow- and outflow boundaries as well as bottom and
462	!-- top of model domain. advc_flags_1 remains zero for all non-prognostic
463	!-- grid points.
464	DO i = nxl, nxr
465	DO j = nys, nyn
466	DO k = nzb+1, nzt
467	!
468	!-- scalar - x-direction
469	!-- WS1 (0), WS3 (1), WS5 (2)
470	IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0) &
471	.OR. .NOT. BTEST(wall_flags_0(k,j,i+2),0) &
472	.OR. .NOT. BTEST(wall_flags_0(k,j,i-1),0) ) &
473	.OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
474	.AND. i == nxl ) &
475	.OR. ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
476	.AND. i == nxr ) ) &
477	THEN
478	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 0 )
479	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0) &
480	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
481	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
482	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
483	) .OR. &
484	( .NOT. BTEST(wall_flags_0(k,j,i-2),0) &
485	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
486	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
487	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
488	) &
489	.OR. &
490	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
491	.AND. i == nxr-1 ) .OR. &
492	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
493	.AND. i == nxlu ) ) &
494	THEN
495	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 1 )
496	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0) &
497	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
498	.AND. BTEST(wall_flags_0(k,j,i+3),0) &
499	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
500	.AND. BTEST(wall_flags_0(k,j,i-2),0) ) &
501	THEN
502	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 2 )
503	ENDIF
504	!
505	!-- scalar - y-direction
506	!-- WS1 (3), WS3 (4), WS5 (5)
507	IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0) &
508	.OR. .NOT. BTEST(wall_flags_0(k,j+2,i),0) &
509	.OR. .NOT. BTEST(wall_flags_0(k,j-1,i),0)) &
510	.OR. ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
511	.AND. j == nys ) &
512	.OR. ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
513	.AND. j == nyn ) ) &
514	THEN
515	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 3 )
516	!
517	!-- WS3
518	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0) &
519	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
520	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
521	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
522	) .OR. &
523	( .NOT. BTEST(wall_flags_0(k,j-2,i),0) &
524	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
525	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
526	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
527	) &
528	.OR. &
529	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
530	.AND. j == nysv ) .OR. &
531	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
532	.AND. j == nyn-1 ) ) &
533	THEN
534	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 4 )
535	!
536	!-- WS5
537	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0) &
538	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
539	.AND. BTEST(wall_flags_0(k,j+3,i),0) &
540	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
541	.AND. BTEST(wall_flags_0(k,j-2,i),0) ) &
542	THEN
543	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 5 )
544	ENDIF
545	!
546	!-- scalar - z-direction
547	!-- WS1 (6), WS3 (7), WS5 (8)
548	IF ( k == nzb+1 ) THEN
549	k_mm = nzb
550	ELSE
551	k_mm = k - 2
552	ENDIF
553	IF ( k > nzt-1 ) THEN
554	k_pp = nzt+1
555	ELSE
556	k_pp = k + 2
557	ENDIF
558	IF ( k > nzt-2 ) THEN
559	k_ppp = nzt+1
560	ELSE
561	k_ppp = k + 3
562	ENDIF
563
564	flag_set = .FALSE.
565	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0) .AND. &
566	BTEST(wall_flags_0(k,j,i),0) .OR. &
567	.NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND. &
568	BTEST(wall_flags_0(k,j,i),0) .OR. &
569	k == nzt ) &
570	THEN
571	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 6 )
572	flag_set = .TRUE.
573	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0) .OR. &
574	.NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND. &
575	BTEST(wall_flags_0(k-1,j,i),0) .AND. &
576	BTEST(wall_flags_0(k,j,i),0) .AND. &
577	BTEST(wall_flags_0(k+1,j,i),0) .AND. &
578	BTEST(wall_flags_0(k_pp,j,i),0) .OR. &
579	k == nzt - 1 ) &
580	THEN
581	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 7 )
582	flag_set = .TRUE.
583	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0) &
584	.AND. BTEST(wall_flags_0(k-1,j,i),0) &
585	.AND. BTEST(wall_flags_0(k,j,i),0) &
586	.AND. BTEST(wall_flags_0(k+1,j,i),0) &
587	.AND. BTEST(wall_flags_0(k_pp,j,i),0) &
588	.AND. BTEST(wall_flags_0(k_ppp,j,i),0) &
589	.AND. .NOT. flag_set ) &
590	THEN
591	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 8 )
592	ENDIF
593
594	ENDDO
595	ENDDO
596	ENDDO
597	ENDIF
598
599	IF ( momentum_advec == 'ws-scheme' ) THEN
600	!
601	!-- Set advc_flags_1 to steer the degradation of the advection scheme in advec_ws
602	!-- near topography, inflow- and outflow boundaries as well as bottom and
603	!-- top of model domain. advc_flags_1 remains zero for all non-prognostic
604	!-- grid points.
605	DO i = nxl, nxr
606	DO j = nys, nyn
607	DO k = nzb+1, nzt
608	!
609	!-- At first, set flags to WS1.
610	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
611	!-- in order to handle the left/south flux.
612	!-- near vertical walls.
613	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
614	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
615	!
616	!-- u component - x-direction
617	!-- WS1 (9), WS3 (10), WS5 (11)
618	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1) .OR. &
619	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
620	.AND. i <= nxlu ) .OR. &
621	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
622	.AND. i == nxr ) ) &
623	THEN
624	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
625	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1) .AND. &
626	BTEST(wall_flags_0(k,j,i+1),1) .OR. &
627	.NOT. BTEST(wall_flags_0(k,j,i-1),1) ) &
628	.OR. &
629	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
630	.AND. i == nxr-1 ) .OR. &
631	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
632	.AND. i == nxlu+1) ) &
633	THEN
634	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 10 )
635	!
636	!-- Clear flag for WS1
637	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
638	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1) .AND. &
639	BTEST(wall_flags_0(k,j,i+2),1) .AND. &
640	BTEST(wall_flags_0(k,j,i-1),1) ) &
641	THEN
642	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 11 )
643	!
644	!-- Clear flag for WS1
645	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
646	ENDIF
647	!
648	!-- u component - y-direction
649	!-- WS1 (12), WS3 (13), WS5 (14)
650	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1) .OR. &
651	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
652	.AND. j == nys ) .OR. &
653	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
654	.AND. j == nyn ) ) &
655	THEN
656	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
657	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1) .AND. &
658	BTEST(wall_flags_0(k,j+1,i),1) .OR. &
659	.NOT. BTEST(wall_flags_0(k,j-1,i),1) ) &
660	.OR. &
661	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
662	.AND. j == nysv ) .OR. &
663	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
664	.AND. j == nyn-1 ) ) &
665	THEN
666	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 13 )
667	!
668	!-- Clear flag for WS1
669	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
670	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1) .AND. &
671	BTEST(wall_flags_0(k,j+2,i),1) .AND. &
672	BTEST(wall_flags_0(k,j-1,i),1) ) &
673	THEN
674	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 14 )
675	!
676	!-- Clear flag for WS1
677	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
678	ENDIF
679	!
680	!-- u component - z-direction
681	!-- WS1 (15), WS3 (16), WS5 (17)
682	IF ( k == nzb+1 ) THEN
683	k_mm = nzb
684	ELSE
685	k_mm = k - 2
686	ENDIF
687	IF ( k > nzt-1 ) THEN
688	k_pp = nzt+1
689	ELSE
690	k_pp = k + 2
691	ENDIF
692	IF ( k > nzt-2 ) THEN
693	k_ppp = nzt+1
694	ELSE
695	k_ppp = k + 3
696	ENDIF
697
698	flag_set = .FALSE.
699	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1) .AND. &
700	BTEST(wall_flags_0(k,j,i),1) .OR. &
701	.NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
702	BTEST(wall_flags_0(k,j,i),1) .OR. &
703	k == nzt ) &
704	THEN
705	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 15 )
706	flag_set = .TRUE.
707	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1) .OR. &
708	.NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND. &
709	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
710	BTEST(wall_flags_0(k,j,i),1) .AND. &
711	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
712	BTEST(wall_flags_0(k_pp,j,i),1) .OR. &
713	k == nzt - 1 ) &
714	THEN
715	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 16 )
716	flag_set = .TRUE.
717	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1) .AND. &
718	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
719	BTEST(wall_flags_0(k,j,i),1) .AND. &
720	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
721	BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
722	BTEST(wall_flags_0(k_ppp,j,i),1) .AND. &
723	.NOT. flag_set ) &
724	THEN
725	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 17 )
726	ENDIF
727
728	ENDDO
729	ENDDO
730	ENDDO
731
732	DO i = nxl, nxr
733	DO j = nys, nyn
734	DO k = nzb+1, nzt
735	!
736	!-- At first, set flags to WS1.
737	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
738	!-- in order to handle the left/south flux.
739	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
740	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
741	!
742	!-- v component - x-direction
743	!-- WS1 (18), WS3 (19), WS5 (20)
744	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2) .OR. &
745	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
746	.AND. i == nxl ) .OR. &
747	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
748	.AND. i == nxr ) ) &
749	THEN
750	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
751	!
752	!-- WS3
753	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2) .AND. &
754	BTEST(wall_flags_0(k,j,i+1),2) ) .OR. &
755	.NOT. BTEST(wall_flags_0(k,j,i-1),2) &
756	.OR. &
757	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
758	.AND. i == nxr-1 ) .OR. &
759	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
760	.AND. i == nxlu ) ) &
761	THEN
762	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 19 )
763	!
764	!-- Clear flag for WS1
765	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
766	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2) .AND. &
767	BTEST(wall_flags_0(k,j,i+2),2) .AND. &
768	BTEST(wall_flags_0(k,j,i-1),2) ) &
769	THEN
770	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 20 )
771	!
772	!-- Clear flag for WS1
773	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
774	ENDIF
775	!
776	!-- v component - y-direction
777	!-- WS1 (21), WS3 (22), WS5 (23)
778	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR. &
779	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
780	.AND. j <= nysv ) .OR. &
781	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
782	.AND. j == nyn ) ) &
783	THEN
784	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
785	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2) .AND. &
786	BTEST(wall_flags_0(k,j+1,i),2) .OR. &
787	.NOT. BTEST(wall_flags_0(k,j-1,i),2) ) &
788	.OR. &
789	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
790	.AND. j == nysv+1) .OR. &
791	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
792	.AND. j == nyn-1 ) ) &
793	THEN
794	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 22 )
795	!
796	!-- Clear flag for WS1
797	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
798	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2) .AND. &
799	BTEST(wall_flags_0(k,j+2,i),2) .AND. &
800	BTEST(wall_flags_0(k,j-1,i),2) ) &
801	THEN
802	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 23 )
803	!
804	!-- Clear flag for WS1
805	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
806	ENDIF
807	!
808	!-- v component - z-direction
809	!-- WS1 (24), WS3 (25), WS5 (26)
810	IF ( k == nzb+1 ) THEN
811	k_mm = nzb
812	ELSE
813	k_mm = k - 2
814	ENDIF
815	IF ( k > nzt-1 ) THEN
816	k_pp = nzt+1
817	ELSE
818	k_pp = k + 2
819	ENDIF
820	IF ( k > nzt-2 ) THEN
821	k_ppp = nzt+1
822	ELSE
823	k_ppp = k + 3
824	ENDIF
825
826	flag_set = .FALSE.
827	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2) .AND. &
828	BTEST(wall_flags_0(k,j,i),2) .OR. &
829	.NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
830	BTEST(wall_flags_0(k,j,i),2) .OR. &
831	k == nzt ) &
832	THEN
833	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 24 )
834	flag_set = .TRUE.
835	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2) .OR. &
836	.NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND. &
837	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
838	BTEST(wall_flags_0(k,j,i),2) .AND. &
839	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
840	BTEST(wall_flags_0(k_pp,j,i),2) .OR. &
841	k == nzt - 1 ) &
842	THEN
843	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 25 )
844	flag_set = .TRUE.
845	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2) .AND. &
846	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
847	BTEST(wall_flags_0(k,j,i),2) .AND. &
848	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
849	BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
850	BTEST(wall_flags_0(k_ppp,j,i),2) .AND. &
851	.NOT. flag_set ) &
852	THEN
853	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 26 )
854	ENDIF
855
856	ENDDO
857	ENDDO
858	ENDDO
859	DO i = nxl, nxr
860	DO j = nys, nyn
861	DO k = nzb+1, nzt
862	!
863	!-- At first, set flags to WS1.
864	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
865	!-- in order to handle the left/south flux.
866	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
867	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
868	!
869	!-- w component - x-direction
870	!-- WS1 (27), WS3 (28), WS5 (29)
871	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR. &
872	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
873	.AND. i == nxl ) .OR. &
874	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
875	.AND. i == nxr ) ) &
876	THEN
877	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
878	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3) .AND. &
879	BTEST(wall_flags_0(k,j,i+1),3) .OR. &
880	.NOT. BTEST(wall_flags_0(k,j,i-1),3) ) &
881	.OR. &
882	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
883	.AND. i == nxr-1 ) .OR. &
884	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
885	.AND. i == nxlu ) ) &
886	THEN
887	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 28 )
888	!
889	!-- Clear flag for WS1
890	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
891	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3) .AND. &
892	BTEST(wall_flags_0(k,j,i+2),3) .AND. &
893	BTEST(wall_flags_0(k,j,i-1),3) ) &
894	THEN
895	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i),29 )
896	!
897	!-- Clear flag for WS1
898	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
899	ENDIF
900	!
901	!-- w component - y-direction
902	!-- WS1 (30), WS3 (31), WS5 (32)
903	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR. &
904	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
905	.AND. j == nys ) .OR. &
906	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
907	.AND. j == nyn ) ) &
908	THEN
909	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
910	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3) .AND. &
911	BTEST(wall_flags_0(k,j+1,i),3) .OR. &
912	.NOT. BTEST(wall_flags_0(k,j-1,i),3) ) &
913	.OR. &
914	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
915	.AND. j == nysv ) .OR. &
916	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
917	.AND. j == nyn-1 ) ) &
918	THEN
919	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 31 )
920	!
921	!-- Clear flag for WS1
922	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
923	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3) .AND. &
924	BTEST(wall_flags_0(k,j+2,i),3) .AND. &
925	BTEST(wall_flags_0(k,j-1,i),3) ) &
926	THEN
927	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 0 )
928	!
929	!-- Clear flag for WS1
930	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
931	ENDIF
932	!
933	!-- w component - z-direction
934	!-- WS1 (33), WS3 (34), WS5 (35)
935	flag_set = .FALSE.
936	IF ( k == nzb+1 ) THEN
937	k_mm = nzb
938	ELSE
939	k_mm = k - 2
940	ENDIF
941	IF ( k > nzt-1 ) THEN
942	k_pp = nzt+1
943	ELSE
944	k_pp = k + 2
945	ENDIF
946	IF ( k > nzt-2 ) THEN
947	k_ppp = nzt+1
948	ELSE
949	k_ppp = k + 3
950	ENDIF
951
952	IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
953	.NOT. BTEST(wall_flags_0(k,j,i),3) .AND. &
954	BTEST(wall_flags_0(k+1,j,i),3) ) .OR. &
955	( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
956	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
957	( .NOT. BTEST(wall_flags_0(k+1,j,i),3) .AND. &
958	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
959	k == nzt ) &
960	THEN
961	!
962	!-- Please note, at k == nzb_w_inner(j,i) a flag is explictely
963	!-- set, although this is not a prognostic level. However,
964	!-- contrary to the advection of u,v and s this is necessary
965	!-- because flux_t(nzb_w_inner(j,i)) is used for the tendency
966	!-- at k == nzb_w_inner(j,i)+1.
967	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 1 )
968	flag_set = .TRUE.
969	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3) .OR. &
970	.NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND. &
971	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
972	BTEST(wall_flags_0(k,j,i),3) .AND. &
973	BTEST(wall_flags_0(k+1,j,i),3) .OR. &
974	k == nzt - 1 ) &
975	THEN
976	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 2 )
977	flag_set = .TRUE.
978	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3) .AND. &
979	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
980	BTEST(wall_flags_0(k,j,i),3) .AND. &
981	BTEST(wall_flags_0(k+1,j,i),3) .AND. &
982	BTEST(wall_flags_0(k_pp,j,i),3) .AND. &
983	BTEST(wall_flags_0(k_ppp,j,i),3) .AND. &
984	.NOT. flag_set ) &
985	THEN
986	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 3 )
987	ENDIF
988
989	ENDDO
990	ENDDO
991	ENDDO
992
993	ENDIF
994
995
996	!
997	!-- Exchange 3D integer wall_flags.
998	IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' &
999	) THEN
1000	!
1001	!-- Exchange ghost points for advection flags
1002	CALL exchange_horiz_int( advc_flags_1, nbgp )
1003	CALL exchange_horiz_int( advc_flags_2, nbgp )
1004	!
1005	!-- Set boundary flags at inflow and outflow boundary in case of
1006	!-- non-cyclic boundary conditions.
1007	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN
1008	advc_flags_1(:,:,nxl-1) = advc_flags_1(:,:,nxl)
1009	advc_flags_2(:,:,nxl-1) = advc_flags_2(:,:,nxl)
1010	ENDIF
1011
1012	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN
1013	advc_flags_1(:,:,nxr+1) = advc_flags_1(:,:,nxr)
1014	advc_flags_2(:,:,nxr+1) = advc_flags_2(:,:,nxr)
1015	ENDIF
1016
1017	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN
1018	advc_flags_1(:,nyn+1,:) = advc_flags_1(:,nyn,:)
1019	advc_flags_2(:,nyn+1,:) = advc_flags_2(:,nyn,:)
1020	ENDIF
1021
1022	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN
1023	advc_flags_1(:,nys-1,:) = advc_flags_1(:,nys,:)
1024	advc_flags_2(:,nys-1,:) = advc_flags_2(:,nys,:)
1025	ENDIF
1026
1027	ENDIF
1028
1029
1030	END SUBROUTINE ws_init_flags
1031
1032
1033	!------------------------------------------------------------------------------!
1034	! Description:
1035	! ------------
1036	!> Initialize variables used for storing statistic quantities (fluxes, variances)
1037	!------------------------------------------------------------------------------!
1038	SUBROUTINE ws_statistics
1039
1040	USE control_parameters, &
1041	ONLY: cloud_physics, humidity, passive_scalar, ocean, &
1042	microphysics_morrison, microphysics_seifert, ws_scheme_mom, &
1043	ws_scheme_sca
1044
1045	USE kinds
1046
1047	USE statistics, &
1048	ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
1049	sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqcs_ws_l, &
1050	sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l, &
1051	sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l
1052
1053
1054	IMPLICIT NONE
1055
1056	!
1057	!-- The arrays needed for statistical evaluation are set to to 0 at the
1058	!-- beginning of prognostic_equations.
1059	IF ( ws_scheme_mom ) THEN
1060	sums_wsus_ws_l = 0.0_wp
1061	sums_wsvs_ws_l = 0.0_wp
1062	sums_us2_ws_l = 0.0_wp
1063	sums_vs2_ws_l = 0.0_wp
1064	sums_ws2_ws_l = 0.0_wp
1065	ENDIF
1066
1067	IF ( ws_scheme_sca ) THEN
1068	sums_wspts_ws_l = 0.0_wp
1069	IF ( humidity ) sums_wsqs_ws_l = 0.0_wp
1070	IF ( passive_scalar ) sums_wsss_ws_l = 0.0_wp
1071	IF ( cloud_physics .AND. microphysics_morrison ) THEN
1072	sums_wsqcs_ws_l = 0.0_wp
1073	sums_wsncs_ws_l = 0.0_wp
1074	ENDIF
1075	IF ( cloud_physics .AND. microphysics_seifert ) THEN
1076	sums_wsqrs_ws_l = 0.0_wp
1077	sums_wsnrs_ws_l = 0.0_wp
1078	ENDIF
1079	IF ( ocean ) sums_wssas_ws_l = 0.0_wp
1080
1081	ENDIF
1082
1083	END SUBROUTINE ws_statistics
1084
1085
1086	!------------------------------------------------------------------------------!
1087	! Description:
1088	! ------------
1089	!> Scalar advection - Call for grid point i,j
1090	!------------------------------------------------------------------------------!
1091	SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local, &
1092	swap_diss_y_local, swap_flux_x_local, &
1093	swap_diss_x_local, i_omp, tn )
1094
1095	USE arrays_3d, &
1096	ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
1097
1098	USE constants, &
1099	ONLY: adv_sca_1, adv_sca_3, adv_sca_5
1100
1101	USE control_parameters, &
1102	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
1103
1104	USE grid_variables, &
1105	ONLY: ddx, ddy
1106
1107	USE indices, &
1108	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
1109	nzt, advc_flags_1
1110
1111	USE kinds
1112
1113	USE pegrid
1114
1115	USE statistics, &
1116	ONLY: hom, sums_wsncs_ws_l, sums_wsnrs_ws_l, sums_wspts_ws_l, &
1117	sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l, &
1118	sums_wssas_ws_l, sums_wsss_ws_l, weight_substep
1119
1120
1121	IMPLICIT NONE
1122
1123	CHARACTER (LEN = *), INTENT(IN) :: sk_char !<
1124
1125	INTEGER(iwp) :: i !<
1126	INTEGER(iwp) :: ibit0 !<
1127	INTEGER(iwp) :: ibit1 !<
1128	INTEGER(iwp) :: ibit2 !<
1129	INTEGER(iwp) :: ibit3 !<
1130	INTEGER(iwp) :: ibit4 !<
1131	INTEGER(iwp) :: ibit5 !<
1132	INTEGER(iwp) :: ibit6 !<
1133	INTEGER(iwp) :: ibit7 !<
1134	INTEGER(iwp) :: ibit8 !<
1135	INTEGER(iwp) :: i_omp !<
1136	INTEGER(iwp) :: j !<
1137	INTEGER(iwp) :: k !<
1138	INTEGER(iwp) :: k_mm !<
1139	INTEGER(iwp) :: k_mmm !<
1140	INTEGER(iwp) :: k_pp !<
1141	INTEGER(iwp) :: k_ppp !<
1142	INTEGER(iwp) :: tn !<
1143
1144	REAL(wp) :: diss_d !<
1145	REAL(wp) :: div !<
1146	REAL(wp) :: flux_d !<
1147	REAL(wp) :: u_comp !<
1148	REAL(wp) :: v_comp !<
1149
1150	#if defined( __nopointer )
1151	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !<
1152	#else
1153	REAL(wp), DIMENSION(:,:,:), POINTER :: sk !<
1154	#endif
1155	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
1156	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
1157	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
1158	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
1159	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
1160	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
1161
1162	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !<
1163	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !<
1164
1165	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !<
1166	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !<
1167
1168
1169	!
1170	!-- Compute southside fluxes of the respective PE bounds.
1171	IF ( j == nys ) THEN
1172	!
1173	!-- Up to the top of the highest topography.
1174	DO k = nzb+1, nzb_max
1175
1176	ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
1177	ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
1178	ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
1179
1180	v_comp = v(k,j,i) - v_gtrans
1181	swap_flux_y_local(k,tn) = v_comp * ( &
1182	( 37.0_wp * ibit5 * adv_sca_5 &
1183	+ 7.0_wp * ibit4 * adv_sca_3 &
1184	+ ibit3 * adv_sca_1 &
1185	) * &
1186	( sk(k,j,i) + sk(k,j-1,i) ) &
1187	- ( 8.0_wp * ibit5 * adv_sca_5 &
1188	+ ibit4 * adv_sca_3 &
1189	) * &
1190	( sk(k,j+1,i) + sk(k,j-2,i) ) &
1191	+ ( ibit5 * adv_sca_5 &
1192	) * &
1193	( sk(k,j+2,i) + sk(k,j-3,i) ) &
1194	)
1195
1196	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1197	( 10.0_wp * ibit5 * adv_sca_5 &
1198	+ 3.0_wp * ibit4 * adv_sca_3 &
1199	+ ibit3 * adv_sca_1 &
1200	) * &
1201	( sk(k,j,i) - sk(k,j-1,i) ) &
1202	- ( 5.0_wp * ibit5 * adv_sca_5 &
1203	+ ibit4 * adv_sca_3 &
1204	) * &
1205	( sk(k,j+1,i) - sk(k,j-2,i) ) &
1206	+ ( ibit5 * adv_sca_5 &
1207	) * &
1208	( sk(k,j+2,i) - sk(k,j-3,i) ) &
1209	)
1210
1211	ENDDO
1212	!
1213	!-- Above to the top of the highest topography. No degradation necessary.
1214	DO k = nzb_max+1, nzt
1215
1216	v_comp = v(k,j,i) - v_gtrans
1217	swap_flux_y_local(k,tn) = v_comp * ( &
1218	37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &
1219	- 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
1220	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
1221	) * adv_sca_5
1222	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1223	10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &
1224	- 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
1225	+ sk(k,j+2,i) - sk(k,j-3,i) &
1226	) * adv_sca_5
1227
1228	ENDDO
1229
1230	ENDIF
1231	!
1232	!-- Compute leftside fluxes of the respective PE bounds.
1233	IF ( i == i_omp ) THEN
1234
1235	DO k = nzb+1, nzb_max
1236
1237	ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
1238	ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
1239	ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
1240
1241	u_comp = u(k,j,i) - u_gtrans
1242	swap_flux_x_local(k,j,tn) = u_comp * ( &
1243	( 37.0_wp * ibit2 * adv_sca_5 &
1244	+ 7.0_wp * ibit1 * adv_sca_3 &
1245	+ ibit0 * adv_sca_1 &
1246	) * &
1247	( sk(k,j,i) + sk(k,j,i-1) ) &
1248	- ( 8.0_wp * ibit2 * adv_sca_5 &
1249	+ ibit1 * adv_sca_3 &
1250	) * &
1251	( sk(k,j,i+1) + sk(k,j,i-2) ) &
1252	+ ( ibit2 * adv_sca_5 &
1253	) * &
1254	( sk(k,j,i+2) + sk(k,j,i-3) ) &
1255	)
1256
1257	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1258	( 10.0_wp * ibit2 * adv_sca_5 &
1259	+ 3.0_wp * ibit1 * adv_sca_3 &
1260	+ ibit0 * adv_sca_1 &
1261	) * &
1262	( sk(k,j,i) - sk(k,j,i-1) ) &
1263	- ( 5.0_wp * ibit2 * adv_sca_5 &
1264	+ ibit1 * adv_sca_3 &
1265	) * &
1266	( sk(k,j,i+1) - sk(k,j,i-2) ) &
1267	+ ( ibit2 * adv_sca_5 &
1268	) * &
1269	( sk(k,j,i+2) - sk(k,j,i-3) ) &
1270	)
1271
1272	ENDDO
1273
1274	DO k = nzb_max+1, nzt
1275
1276	u_comp = u(k,j,i) - u_gtrans
1277	swap_flux_x_local(k,j,tn) = u_comp * ( &
1278	37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &
1279	- 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
1280	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
1281	) * adv_sca_5
1282
1283	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1284	10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &
1285	- 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
1286	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
1287	) * adv_sca_5
1288
1289	ENDDO
1290
1291	ENDIF
1292
1293	flux_t(0) = 0.0_wp
1294	diss_t(0) = 0.0_wp
1295	flux_d = 0.0_wp
1296	diss_d = 0.0_wp
1297	!
1298	!-- Now compute the fluxes and tendency terms for the horizontal and
1299	!-- vertical parts up to the top of the highest topography.
1300	DO k = nzb+1, nzb_max
1301	!
1302	!-- Note: It is faster to conduct all multiplications explicitly, e.g.
1303	!-- * adv_sca_5 ... than to determine a factor and multiplicate the
1304	!-- flux at the end.
1305
1306	ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
1307	ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
1308	ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
1309
1310	u_comp = u(k,j,i+1) - u_gtrans
1311	flux_r(k) = u_comp * ( &
1312	( 37.0_wp * ibit2 * adv_sca_5 &
1313	+ 7.0_wp * ibit1 * adv_sca_3 &
1314	+ ibit0 * adv_sca_1 &
1315	) * &
1316	( sk(k,j,i+1) + sk(k,j,i) ) &
1317	- ( 8.0_wp * ibit2 * adv_sca_5 &
1318	+ ibit1 * adv_sca_3 &
1319	) * &
1320	( sk(k,j,i+2) + sk(k,j,i-1) ) &
1321	+ ( ibit2 * adv_sca_5 &
1322	) * &
1323	( sk(k,j,i+3) + sk(k,j,i-2) ) &
1324	)
1325
1326	diss_r(k) = -ABS( u_comp ) * ( &
1327	( 10.0_wp * ibit2 * adv_sca_5 &
1328	+ 3.0_wp * ibit1 * adv_sca_3 &
1329	+ ibit0 * adv_sca_1 &
1330	) * &
1331	( sk(k,j,i+1) - sk(k,j,i) ) &
1332	- ( 5.0_wp * ibit2 * adv_sca_5 &
1333	+ ibit1 * adv_sca_3 &
1334	) * &
1335	( sk(k,j,i+2) - sk(k,j,i-1) ) &
1336	+ ( ibit2 * adv_sca_5 &
1337	) * &
1338	( sk(k,j,i+3) - sk(k,j,i-2) ) &
1339	)
1340
1341	ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
1342	ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
1343	ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
1344
1345	v_comp = v(k,j+1,i) - v_gtrans
1346	flux_n(k) = v_comp * ( &
1347	( 37.0_wp * ibit5 * adv_sca_5 &
1348	+ 7.0_wp * ibit4 * adv_sca_3 &
1349	+ ibit3 * adv_sca_1 &
1350	) * &
1351	( sk(k,j+1,i) + sk(k,j,i) ) &
1352	- ( 8.0_wp * ibit5 * adv_sca_5 &
1353	+ ibit4 * adv_sca_3 &
1354	) * &
1355	( sk(k,j+2,i) + sk(k,j-1,i) ) &
1356	+ ( ibit5 * adv_sca_5 &
1357	) * &
1358	( sk(k,j+3,i) + sk(k,j-2,i) ) &
1359	)
1360
1361	diss_n(k) = -ABS( v_comp ) * ( &
1362	( 10.0_wp * ibit5 * adv_sca_5 &
1363	+ 3.0_wp * ibit4 * adv_sca_3 &
1364	+ ibit3 * adv_sca_1 &
1365	) * &
1366	( sk(k,j+1,i) - sk(k,j,i) ) &
1367	- ( 5.0_wp * ibit5 * adv_sca_5 &
1368	+ ibit4 * adv_sca_3 &
1369	) * &
1370	( sk(k,j+2,i) - sk(k,j-1,i) ) &
1371	+ ( ibit5 * adv_sca_5 &
1372	) * &
1373	( sk(k,j+3,i) - sk(k,j-2,i) ) &
1374	)
1375	!
1376	!-- k index has to be modified near bottom and top, else array
1377	!-- subscripts will be exceeded.
1378	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
1379	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
1380	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
1381
1382	k_ppp = k + 3 * ibit8
1383	k_pp = k + 2 * ( 1 - ibit6 )
1384	k_mm = k - 2 * ibit8
1385
1386
1387	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1388	( 37.0_wp * ibit8 * adv_sca_5 &
1389	+ 7.0_wp * ibit7 * adv_sca_3 &
1390	+ ibit6 * adv_sca_1 &
1391	) * &
1392	( sk(k+1,j,i) + sk(k,j,i) ) &
1393	- ( 8.0_wp * ibit8 * adv_sca_5 &
1394	+ ibit7 * adv_sca_3 &
1395	) * &
1396	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1397	+ ( ibit8 * adv_sca_5 &
1398	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1399	)
1400
1401	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1402	( 10.0_wp * ibit8 * adv_sca_5 &
1403	+ 3.0_wp * ibit7 * adv_sca_3 &
1404	+ ibit6 * adv_sca_1 &
1405	) * &
1406	( sk(k+1,j,i) - sk(k,j,i) ) &
1407	- ( 5.0_wp * ibit8 * adv_sca_5 &
1408	+ ibit7 * adv_sca_3 &
1409	) * &
1410	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1411	+ ( ibit8 * adv_sca_5 &
1412	) * &
1413	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1414	)
1415	!
1416	!-- Calculate the divergence of the velocity field. A respective
1417	!-- correction is needed to overcome numerical instabilities caused
1418	!-- by a not sufficient reduction of divergences near topography.
1419	div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &
1420	- u(k,j,i) * ( IBITS(advc_flags_1(k,j,i-1),0,1) &
1421	+ IBITS(advc_flags_1(k,j,i-1),1,1) &
1422	+ IBITS(advc_flags_1(k,j,i-1),2,1) &
1423	) &
1424	) * rho_air(k) * ddx &
1425	+ ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &
1426	- v(k,j,i) * ( IBITS(advc_flags_1(k,j-1,i),3,1) &
1427	+ IBITS(advc_flags_1(k,j-1,i),4,1) &
1428	+ IBITS(advc_flags_1(k,j-1,i),5,1) &
1429	) &
1430	) * rho_air(k) * ddy &
1431	+ ( w(k,j,i) * rho_air_zw(k) * &
1432	( ibit6 + ibit7 + ibit8 ) &
1433	- w(k-1,j,i) * rho_air_zw(k-1) * &
1434	( IBITS(advc_flags_1(k-1,j,i),6,1) &
1435	+ IBITS(advc_flags_1(k-1,j,i),7,1) &
1436	+ IBITS(advc_flags_1(k-1,j,i),8,1) &
1437	) &
1438	) * ddzw(k)
1439
1440
1441	tend(k,j,i) = tend(k,j,i) - ( &
1442	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1443	swap_diss_x_local(k,j,tn) ) * ddx &
1444	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1445	swap_diss_y_local(k,tn) ) * ddy &
1446	+ ( ( flux_t(k) + diss_t(k) ) - &
1447	( flux_d + diss_d ) &
1448	) * drho_air(k) * ddzw(k) &
1449	) + sk(k,j,i) * div
1450
1451	swap_flux_y_local(k,tn) = flux_n(k)
1452	swap_diss_y_local(k,tn) = diss_n(k)
1453	swap_flux_x_local(k,j,tn) = flux_r(k)
1454	swap_diss_x_local(k,j,tn) = diss_r(k)
1455	flux_d = flux_t(k)
1456	diss_d = diss_t(k)
1457
1458	ENDDO
1459	!
1460	!-- Now compute the fluxes and tendency terms for the horizontal and
1461	!-- vertical parts above the top of the highest topography. No degradation
1462	!-- for the horizontal parts, but for the vertical it is stell needed.
1463	DO k = nzb_max+1, nzt
1464
1465	u_comp = u(k,j,i+1) - u_gtrans
1466	flux_r(k) = u_comp * ( &
1467	37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) &
1468	- 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
1469	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
1470	diss_r(k) = -ABS( u_comp ) * ( &
1471	10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) &
1472	- 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
1473	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
1474
1475	v_comp = v(k,j+1,i) - v_gtrans
1476	flux_n(k) = v_comp * ( &
1477	37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) &
1478	- 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
1479	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
1480	diss_n(k) = -ABS( v_comp ) * ( &
1481	10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) &
1482	- 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
1483	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
1484	!
1485	!-- k index has to be modified near bottom and top, else array
1486	!-- subscripts will be exceeded.
1487	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
1488	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
1489	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
1490
1491	k_ppp = k + 3 * ibit8
1492	k_pp = k + 2 * ( 1 - ibit6 )
1493	k_mm = k - 2 * ibit8
1494
1495	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1496	( 37.0_wp * ibit8 * adv_sca_5 &
1497	+ 7.0_wp * ibit7 * adv_sca_3 &
1498	+ ibit6 * adv_sca_1 &
1499	) * &
1500	( sk(k+1,j,i) + sk(k,j,i) ) &
1501	- ( 8.0_wp * ibit8 * adv_sca_5 &
1502	+ ibit7 * adv_sca_3 &
1503	) * &
1504	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1505	+ ( ibit8 * adv_sca_5 &
1506	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1507	)
1508
1509	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1510	( 10.0_wp * ibit8 * adv_sca_5 &
1511	+ 3.0_wp * ibit7 * adv_sca_3 &
1512	+ ibit6 * adv_sca_1 &
1513	) * &
1514	( sk(k+1,j,i) - sk(k,j,i) ) &
1515	- ( 5.0_wp * ibit8 * adv_sca_5 &
1516	+ ibit7 * adv_sca_3 &
1517	) * &
1518	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1519	+ ( ibit8 * adv_sca_5 &
1520	) * &
1521	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1522	)
1523	!
1524	!-- Calculate the divergence of the velocity field. A respective
1525	!-- correction is needed to overcome numerical instabilities introduced
1526	!-- by a not sufficient reduction of divergences near topography.
1527	div = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx &
1528	+ ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy &
1529	+ ( w(k,j,i) * rho_air_zw(k) - &
1530	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
1531
1532	tend(k,j,i) = tend(k,j,i) - ( &
1533	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1534	swap_diss_x_local(k,j,tn) ) * ddx &
1535	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1536	swap_diss_y_local(k,tn) ) * ddy &
1537	+ ( ( flux_t(k) + diss_t(k) ) - &
1538	( flux_d + diss_d ) &
1539	) * drho_air(k) * ddzw(k) &
1540	) + sk(k,j,i) * div
1541
1542
1543	swap_flux_y_local(k,tn) = flux_n(k)
1544	swap_diss_y_local(k,tn) = diss_n(k)
1545	swap_flux_x_local(k,j,tn) = flux_r(k)
1546	swap_diss_x_local(k,j,tn) = diss_r(k)
1547	flux_d = flux_t(k)
1548	diss_d = diss_t(k)
1549
1550	ENDDO
1551
1552	!
1553	!-- Evaluation of statistics.
1554	SELECT CASE ( sk_char )
1555
1556	CASE ( 'pt' )
1557
1558	DO k = nzb, nzt
1559	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + &
1560	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1561	* ( w(k,j,i) - hom(k,1,3,0) ) &
1562	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1563	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1564	) * weight_substep(intermediate_timestep_count)
1565	ENDDO
1566
1567	CASE ( 'sa' )
1568
1569	DO k = nzb, nzt
1570	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + &
1571	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1572	* ( w(k,j,i) - hom(k,1,3,0) ) &
1573	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1574	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1575	) * weight_substep(intermediate_timestep_count)
1576	ENDDO
1577
1578	CASE ( 'q' )
1579
1580	DO k = nzb, nzt
1581	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + &
1582	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1583	* ( w(k,j,i) - hom(k,1,3,0) ) &
1584	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1585	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1586	) * weight_substep(intermediate_timestep_count)
1587	ENDDO
1588
1589	CASE ( 'qc' )
1590
1591	DO k = nzb, nzt
1592	sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) + &
1593	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1594	* ( w(k,j,i) - hom(k,1,3,0) ) &
1595	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1596	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1597	) * weight_substep(intermediate_timestep_count)
1598	ENDDO
1599
1600
1601	CASE ( 'qr' )
1602
1603	DO k = nzb, nzt
1604	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + &
1605	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1606	* ( w(k,j,i) - hom(k,1,3,0) ) &
1607	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1608	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1609	) * weight_substep(intermediate_timestep_count)
1610	ENDDO
1611
1612	CASE ( 'nc' )
1613
1614	DO k = nzb, nzt
1615	sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) + &
1616	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1617	* ( w(k,j,i) - hom(k,1,3,0) ) &
1618	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1619	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1620	) * weight_substep(intermediate_timestep_count)
1621	ENDDO
1622
1623	CASE ( 'nr' )
1624
1625	DO k = nzb, nzt
1626	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + &
1627	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1628	* ( w(k,j,i) - hom(k,1,3,0) ) &
1629	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1630	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1631	) * weight_substep(intermediate_timestep_count)
1632	ENDDO
1633
1634	CASE ( 's' )
1635
1636	DO k = nzb, nzt
1637	sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) + &
1638	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1639	* ( w(k,j,i) - hom(k,1,3,0) ) &
1640	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1641	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1642	) * weight_substep(intermediate_timestep_count)
1643	ENDDO
1644
1645	END SELECT
1646
1647	END SUBROUTINE advec_s_ws_ij
1648
1649
1650
1651
1652	!------------------------------------------------------------------------------!
1653	! Description:
1654	! ------------
1655	!> Advection of u-component - Call for grid point i,j
1656	!------------------------------------------------------------------------------!
1657	SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
1658
1659	USE arrays_3d, &
1660	ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,&
1661	drho_air, rho_air, rho_air_zw
1662
1663	USE constants, &
1664	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
1665
1666	USE control_parameters, &
1667	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
1668
1669	USE grid_variables, &
1670	ONLY: ddx, ddy
1671
1672	USE indices, &
1673	ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1
1674
1675	USE kinds
1676
1677	USE statistics, &
1678	ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
1679
1680	IMPLICIT NONE
1681
1682	INTEGER(iwp) :: i !<
1683	INTEGER(iwp) :: ibit9 !<
1684	INTEGER(iwp) :: ibit10 !<
1685	INTEGER(iwp) :: ibit11 !<
1686	INTEGER(iwp) :: ibit12 !<
1687	INTEGER(iwp) :: ibit13 !<
1688	INTEGER(iwp) :: ibit14 !<
1689	INTEGER(iwp) :: ibit15 !<
1690	INTEGER(iwp) :: ibit16 !<
1691	INTEGER(iwp) :: ibit17 !<
1692	INTEGER(iwp) :: i_omp !<
1693	INTEGER(iwp) :: j !<
1694	INTEGER(iwp) :: k !<
1695	INTEGER(iwp) :: k_mm !<
1696	INTEGER(iwp) :: k_pp !<
1697	INTEGER(iwp) :: k_ppp !<
1698	INTEGER(iwp) :: tn !<
1699
1700	REAL(wp) :: diss_d !<
1701	REAL(wp) :: div !<
1702	REAL(wp) :: flux_d !<
1703	REAL(wp) :: gu !<
1704	REAL(wp) :: gv !<
1705	REAL(wp) :: u_comp_l !<
1706	REAL(wp) :: v_comp !<
1707	REAL(wp) :: w_comp !<
1708
1709	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
1710	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
1711	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
1712	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
1713	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
1714	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
1715	REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !<
1716
1717	gu = 2.0_wp * u_gtrans
1718	gv = 2.0_wp * v_gtrans
1719	!
1720	!-- Compute southside fluxes for the respective boundary of PE
1721	IF ( j == nys ) THEN
1722
1723	DO k = nzb+1, nzb_max
1724
1725	ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
1726	ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
1727	ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
1728
1729	v_comp = v(k,j,i) + v(k,j,i-1) - gv
1730	flux_s_u(k,tn) = v_comp * ( &
1731	( 37.0_wp * ibit14 * adv_mom_5 &
1732	+ 7.0_wp * ibit13 * adv_mom_3 &
1733	+ ibit12 * adv_mom_1 &
1734	) * &
1735	( u(k,j,i) + u(k,j-1,i) ) &
1736	- ( 8.0_wp * ibit14 * adv_mom_5 &
1737	+ ibit13 * adv_mom_3 &
1738	) * &
1739	( u(k,j+1,i) + u(k,j-2,i) ) &
1740	+ ( ibit14 * adv_mom_5 &
1741	) * &
1742	( u(k,j+2,i) + u(k,j-3,i) ) &
1743	)
1744
1745	diss_s_u(k,tn) = - ABS ( v_comp ) * ( &
1746	( 10.0_wp * ibit14 * adv_mom_5 &
1747	+ 3.0_wp * ibit13 * adv_mom_3 &
1748	+ ibit12 * adv_mom_1 &
1749	) * &
1750	( u(k,j,i) - u(k,j-1,i) ) &
1751	- ( 5.0_wp * ibit14 * adv_mom_5 &
1752	+ ibit13 * adv_mom_3 &
1753	) * &
1754	( u(k,j+1,i) - u(k,j-2,i) ) &
1755	+ ( ibit14 * adv_mom_5 &
1756	) * &
1757	( u(k,j+2,i) - u(k,j-3,i) ) &
1758	)
1759
1760	ENDDO
1761
1762	DO k = nzb_max+1, nzt
1763
1764	v_comp = v(k,j,i) + v(k,j,i-1) - gv
1765	flux_s_u(k,tn) = v_comp * ( &
1766	37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) &
1767	- 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) &
1768	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
1769	diss_s_u(k,tn) = - ABS(v_comp) * ( &
1770	10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) &
1771	- 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) &
1772	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
1773
1774	ENDDO
1775
1776	ENDIF
1777	!
1778	!-- Compute leftside fluxes for the respective boundary of PE
1779	IF ( i == i_omp ) THEN
1780
1781	DO k = nzb+1, nzb_max
1782
1783	ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
1784	ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
1785	ibit9 = IBITS(advc_flags_1(k,j,i-1),9,1)
1786
1787	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
1788	flux_l_u(k,j,tn) = u_comp_l * ( &
1789	( 37.0_wp * ibit11 * adv_mom_5 &
1790	+ 7.0_wp * ibit10 * adv_mom_3 &
1791	+ ibit9 * adv_mom_1 &
1792	) * &
1793	( u(k,j,i) + u(k,j,i-1) ) &
1794	- ( 8.0_wp * ibit11 * adv_mom_5 &
1795	+ ibit10 * adv_mom_3 &
1796	) * &
1797	( u(k,j,i+1) + u(k,j,i-2) ) &
1798	+ ( ibit11 * adv_mom_5 &
1799	) * &
1800	( u(k,j,i+2) + u(k,j,i-3) ) &
1801	)
1802
1803	diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( &
1804	( 10.0_wp * ibit11 * adv_mom_5 &
1805	+ 3.0_wp * ibit10 * adv_mom_3 &
1806	+ ibit9 * adv_mom_1 &
1807	) * &
1808	( u(k,j,i) - u(k,j,i-1) ) &
1809	- ( 5.0_wp * ibit11 * adv_mom_5 &
1810	+ ibit10 * adv_mom_3 &
1811	) * &
1812	( u(k,j,i+1) - u(k,j,i-2) ) &
1813	+ ( ibit11 * adv_mom_5 &
1814	) * &
1815	( u(k,j,i+2) - u(k,j,i-3) ) &
1816	)
1817
1818	ENDDO
1819
1820	DO k = nzb_max+1, nzt
1821
1822	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
1823	flux_l_u(k,j,tn) = u_comp_l * ( &
1824	37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) &
1825	- 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) &
1826	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
1827	diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( &
1828	10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) &
1829	- 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) &
1830	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
1831
1832	ENDDO
1833
1834	ENDIF
1835
1836	flux_t(0) = 0.0_wp
1837	diss_t(0) = 0.0_wp
1838	flux_d = 0.0_wp
1839	diss_d = 0.0_wp
1840	!
1841	!-- Now compute the fluxes tendency terms for the horizontal and
1842	!-- vertical parts.
1843	DO k = nzb+1, nzb_max
1844
1845	ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
1846	ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
1847	ibit9 = IBITS(advc_flags_1(k,j,i),9,1)
1848
1849	u_comp(k) = u(k,j,i+1) + u(k,j,i)
1850	flux_r(k) = ( u_comp(k) - gu ) * ( &
1851	( 37.0_wp * ibit11 * adv_mom_5 &
1852	+ 7.0_wp * ibit10 * adv_mom_3 &
1853	+ ibit9 * adv_mom_1 &
1854	) * &
1855	( u(k,j,i+1) + u(k,j,i) ) &
1856	- ( 8.0_wp * ibit11 * adv_mom_5 &
1857	+ ibit10 * adv_mom_3 &
1858	) * &
1859	( u(k,j,i+2) + u(k,j,i-1) ) &
1860	+ ( ibit11 * adv_mom_5 &
1861	) * &
1862	( u(k,j,i+3) + u(k,j,i-2) ) &
1863	)
1864
1865	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
1866	( 10.0_wp * ibit11 * adv_mom_5 &
1867	+ 3.0_wp * ibit10 * adv_mom_3 &
1868	+ ibit9 * adv_mom_1 &
1869	) * &
1870	( u(k,j,i+1) - u(k,j,i) ) &
1871	- ( 5.0_wp * ibit11 * adv_mom_5 &
1872	+ ibit10 * adv_mom_3 &
1873	) * &
1874	( u(k,j,i+2) - u(k,j,i-1) ) &
1875	+ ( ibit11 * adv_mom_5 &
1876	) * &
1877	( u(k,j,i+3) - u(k,j,i-2) ) &
1878	)
1879
1880	ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
1881	ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
1882	ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
1883
1884	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
1885	flux_n(k) = v_comp * ( &
1886	( 37.0_wp * ibit14 * adv_mom_5 &
1887	+ 7.0_wp * ibit13 * adv_mom_3 &
1888	+ ibit12 * adv_mom_1 &
1889	) * &
1890	( u(k,j+1,i) + u(k,j,i) ) &
1891	- ( 8.0_wp * ibit14 * adv_mom_5 &
1892	+ ibit13 * adv_mom_3 &
1893	) * &
1894	( u(k,j+2,i) + u(k,j-1,i) ) &
1895	+ ( ibit14 * adv_mom_5 &
1896	) * &
1897	( u(k,j+3,i) + u(k,j-2,i) ) &
1898	)
1899
1900	diss_n(k) = - ABS ( v_comp ) * ( &
1901	( 10.0_wp * ibit14 * adv_mom_5 &
1902	+ 3.0_wp * ibit13 * adv_mom_3 &
1903	+ ibit12 * adv_mom_1 &
1904	) * &
1905	( u(k,j+1,i) - u(k,j,i) ) &
1906	- ( 5.0_wp * ibit14 * adv_mom_5 &
1907	+ ibit13 * adv_mom_3 &
1908	) * &
1909	( u(k,j+2,i) - u(k,j-1,i) ) &
1910	+ ( ibit14 * adv_mom_5 &
1911	) * &
1912	( u(k,j+3,i) - u(k,j-2,i) ) &
1913	)
1914	!
1915	!-- k index has to be modified near bottom and top, else array
1916	!-- subscripts will be exceeded.
1917	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
1918	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
1919	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
1920
1921	k_ppp = k + 3 * ibit17
1922	k_pp = k + 2 * ( 1 - ibit15 )
1923	k_mm = k - 2 * ibit17
1924
1925	w_comp = w(k,j,i) + w(k,j,i-1)
1926	flux_t(k) = w_comp * rho_air_zw(k) * ( &
1927	( 37.0_wp * ibit17 * adv_mom_5 &
1928	+ 7.0_wp * ibit16 * adv_mom_3 &
1929	+ ibit15 * adv_mom_1 &
1930	) * &
1931	( u(k+1,j,i) + u(k,j,i) ) &
1932	- ( 8.0_wp * ibit17 * adv_mom_5 &
1933	+ ibit16 * adv_mom_3 &
1934	) * &
1935	( u(k_pp,j,i) + u(k-1,j,i) ) &
1936	+ ( ibit17 * adv_mom_5 &
1937	) * &
1938	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
1939	)
1940
1941	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
1942	( 10.0_wp * ibit17 * adv_mom_5 &
1943	+ 3.0_wp * ibit16 * adv_mom_3 &
1944	+ ibit15 * adv_mom_1 &
1945	) * &
1946	( u(k+1,j,i) - u(k,j,i) ) &
1947	- ( 5.0_wp * ibit17 * adv_mom_5 &
1948	+ ibit16 * adv_mom_3 &
1949	) * &
1950	( u(k_pp,j,i) - u(k-1,j,i) ) &
1951	+ ( ibit17 * adv_mom_5 &
1952	) * &
1953	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
1954	)
1955	!
1956	!-- Calculate the divergence of the velocity field. A respective
1957	!-- correction is needed to overcome numerical instabilities introduced
1958	!-- by a not sufficient reduction of divergences near topography.
1959	div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) &
1960	- ( u(k,j,i) + u(k,j,i-1) ) &
1961	* ( IBITS(advc_flags_1(k,j,i-1),9,1) &
1962	+ IBITS(advc_flags_1(k,j,i-1),10,1) &
1963	+ IBITS(advc_flags_1(k,j,i-1),11,1) &
1964	) &
1965	) * rho_air(k) * ddx &
1966	+ ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) &
1967	- ( v(k,j,i) + v(k,j,i-1 ) ) &
1968	* ( IBITS(advc_flags_1(k,j-1,i),12,1) &
1969	+ IBITS(advc_flags_1(k,j-1,i),13,1) &
1970	+ IBITS(advc_flags_1(k,j-1,i),14,1) &
1971	) &
1972	) * rho_air(k) * ddy &
1973	+ ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) &
1974	- ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
1975	* ( IBITS(advc_flags_1(k-1,j,i),15,1) &
1976	+ IBITS(advc_flags_1(k-1,j,i),16,1) &
1977	+ IBITS(advc_flags_1(k-1,j,i),17,1) &
1978	) &
1979	) * ddzw(k) &
1980	) * 0.5_wp
1981
1982
1983	tend(k,j,i) = tend(k,j,i) - ( &
1984	( flux_r(k) + diss_r(k) &
1985	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
1986	+ ( flux_n(k) + diss_n(k) &
1987	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
1988	+ ( ( flux_t(k) + diss_t(k) ) &
1989	- ( flux_d + diss_d ) &
1990	) * drho_air(k) * ddzw(k) &
1991	) + div * u(k,j,i)
1992
1993	flux_l_u(k,j,tn) = flux_r(k)
1994	diss_l_u(k,j,tn) = diss_r(k)
1995	flux_s_u(k,tn) = flux_n(k)
1996	diss_s_u(k,tn) = diss_n(k)
1997	flux_d = flux_t(k)
1998	diss_d = diss_t(k)
1999	!
2000	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2001	!-- right evaluation when gallilei_trans = .T. .
2002	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2003	+ ( flux_r(k) &
2004	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2005	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2006	+ diss_r(k) &
2007	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2008	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2009	) * weight_substep(intermediate_timestep_count)
2010	!
2011	!-- Statistical Evaluation of w'u'.
2012	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2013	+ ( flux_t(k) &
2014	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2015	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2016	+ diss_t(k) &
2017	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2018	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2019	) * weight_substep(intermediate_timestep_count)
2020	ENDDO
2021
2022	DO k = nzb_max+1, nzt
2023
2024	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2025	flux_r(k) = ( u_comp(k) - gu ) * ( &
2026	37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &
2027	- 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &
2028	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2029	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2030	10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &
2031	- 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &
2032	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2033
2034	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2035	flux_n(k) = v_comp * ( &
2036	37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &
2037	- 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &
2038	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2039	diss_n(k) = - ABS( v_comp ) * ( &
2040	10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &
2041	- 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &
2042	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2043	!
2044	!-- k index has to be modified near bottom and top, else array
2045	!-- subscripts will be exceeded.
2046	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
2047	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
2048	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
2049
2050	k_ppp = k + 3 * ibit17
2051	k_pp = k + 2 * ( 1 - ibit15 )
2052	k_mm = k - 2 * ibit17
2053
2054	w_comp = w(k,j,i) + w(k,j,i-1)
2055	flux_t(k) = w_comp * rho_air_zw(k) * ( &
2056	( 37.0_wp * ibit17 * adv_mom_5 &
2057	+ 7.0_wp * ibit16 * adv_mom_3 &
2058	+ ibit15 * adv_mom_1 &
2059	) * &
2060	( u(k+1,j,i) + u(k,j,i) ) &
2061	- ( 8.0_wp * ibit17 * adv_mom_5 &
2062	+ ibit16 * adv_mom_3 &
2063	) * &
2064	( u(k_pp,j,i) + u(k-1,j,i) ) &
2065	+ ( ibit17 * adv_mom_5 &
2066	) * &
2067	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2068	)
2069
2070	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
2071	( 10.0_wp * ibit17 * adv_mom_5 &
2072	+ 3.0_wp * ibit16 * adv_mom_3 &
2073	+ ibit15 * adv_mom_1 &
2074	) * &
2075	( u(k+1,j,i) - u(k,j,i) ) &
2076	- ( 5.0_wp * ibit17 * adv_mom_5 &
2077	+ ibit16 * adv_mom_3 &
2078	) * &
2079	( u(k_pp,j,i) - u(k-1,j,i) ) &
2080	+ ( ibit17 * adv_mom_5 &
2081	) * &
2082	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2083	)
2084	!
2085	!-- Calculate the divergence of the velocity field. A respective
2086	!-- correction is needed to overcome numerical instabilities introduced
2087	!-- by a not sufficient reduction of divergences near topography.
2088	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
2089	* rho_air(k)&
2090	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
2091	* rho_air(k)&
2092	+ ( w_comp * rho_air_zw(k) - &
2093	( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
2094	) * ddzw(k) &
2095	) * 0.5_wp
2096
2097	tend(k,j,i) = tend(k,j,i) - ( &
2098	( flux_r(k) + diss_r(k) &
2099	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
2100	+ ( flux_n(k) + diss_n(k) &
2101	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
2102	+ ( ( flux_t(k) + diss_t(k) ) &
2103	- ( flux_d + diss_d ) &
2104	) * drho_air(k) * ddzw(k) &
2105	) + div * u(k,j,i)
2106
2107	flux_l_u(k,j,tn) = flux_r(k)
2108	diss_l_u(k,j,tn) = diss_r(k)
2109	flux_s_u(k,tn) = flux_n(k)
2110	diss_s_u(k,tn) = diss_n(k)
2111	flux_d = flux_t(k)
2112	diss_d = diss_t(k)
2113	!
2114	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2115	!-- right evaluation when gallilei_trans = .T. .
2116	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2117	+ ( flux_r(k) &
2118	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2119	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2120	+ diss_r(k) &
2121	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2122	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2123	) * weight_substep(intermediate_timestep_count)
2124	!
2125	!-- Statistical Evaluation of w'u'.
2126	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2127	+ ( flux_t(k) &
2128	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2129	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2130	+ diss_t(k) &
2131	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2132	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2133	) * weight_substep(intermediate_timestep_count)
2134	ENDDO
2135
2136	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
2137
2138
2139
2140	END SUBROUTINE advec_u_ws_ij
2141
2142
2143
2144	!-----------------------------------------------------------------------------!
2145	! Description:
2146	! ------------
2147	!> Advection of v-component - Call for grid point i,j
2148	!-----------------------------------------------------------------------------!
2149	SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
2150
2151	USE arrays_3d, &
2152	ONLY: ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w, &
2153	drho_air, rho_air, rho_air_zw
2154
2155	USE constants, &
2156	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
2157
2158	USE control_parameters, &
2159	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
2160
2161	USE grid_variables, &
2162	ONLY: ddx, ddy
2163
2164	USE indices, &
2165	ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_max, nzt, advc_flags_1
2166
2167	USE kinds
2168
2169	USE statistics, &
2170	ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep
2171
2172	IMPLICIT NONE
2173
2174	INTEGER(iwp) :: i !<
2175	INTEGER(iwp) :: ibit18 !<
2176	INTEGER(iwp) :: ibit19 !<
2177	INTEGER(iwp) :: ibit20 !<
2178	INTEGER(iwp) :: ibit21 !<
2179	INTEGER(iwp) :: ibit22 !<
2180	INTEGER(iwp) :: ibit23 !<
2181	INTEGER(iwp) :: ibit24 !<
2182	INTEGER(iwp) :: ibit25 !<
2183	INTEGER(iwp) :: ibit26 !<
2184	INTEGER(iwp) :: i_omp !<
2185	INTEGER(iwp) :: j !<
2186	INTEGER(iwp) :: k !<
2187	INTEGER(iwp) :: k_mm !<
2188	INTEGER(iwp) :: k_pp !<
2189	INTEGER(iwp) :: k_ppp !<
2190	INTEGER(iwp) :: tn !<
2191
2192	REAL(wp) :: diss_d !<
2193	REAL(wp) :: div !<
2194	REAL(wp) :: flux_d !<
2195	REAL(wp) :: gu !<
2196	REAL(wp) :: gv !<
2197	REAL(wp) :: u_comp !<
2198	REAL(wp) :: v_comp_l !<
2199	REAL(wp) :: w_comp !<
2200
2201	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
2202	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
2203	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
2204	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
2205	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
2206	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
2207	REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !<
2208
2209	gu = 2.0_wp * u_gtrans
2210	gv = 2.0_wp * v_gtrans
2211
2212	!
2213	!-- Compute leftside fluxes for the respective boundary.
2214	IF ( i == i_omp ) THEN
2215
2216	DO k = nzb+1, nzb_max
2217
2218	ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
2219	ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
2220	ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
2221
2222	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2223	flux_l_v(k,j,tn) = u_comp * ( &
2224	( 37.0_wp * ibit20 * adv_mom_5 &
2225	+ 7.0_wp * ibit19 * adv_mom_3 &
2226	+ ibit18 * adv_mom_1 &
2227	) * &
2228	( v(k,j,i) + v(k,j,i-1) ) &
2229	- ( 8.0_wp * ibit20 * adv_mom_5 &
2230	+ ibit19 * adv_mom_3 &
2231	) * &
2232	( v(k,j,i+1) + v(k,j,i-2) ) &
2233	+ ( ibit20 * adv_mom_5 &
2234	) * &
2235	( v(k,j,i+2) + v(k,j,i-3) ) &
2236	)
2237
2238	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
2239	( 10.0_wp * ibit20 * adv_mom_5 &
2240	+ 3.0_wp * ibit19 * adv_mom_3 &
2241	+ ibit18 * adv_mom_1 &
2242	) * &
2243	( v(k,j,i) - v(k,j,i-1) ) &
2244	- ( 5.0_wp * ibit20 * adv_mom_5 &
2245	+ ibit19 * adv_mom_3 &
2246	) * &
2247	( v(k,j,i+1) - v(k,j,i-2) ) &
2248	+ ( ibit20 * adv_mom_5 &
2249	) * &
2250	( v(k,j,i+2) - v(k,j,i-3) ) &
2251	)
2252
2253	ENDDO
2254
2255	DO k = nzb_max+1, nzt
2256
2257	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2258	flux_l_v(k,j,tn) = u_comp * ( &
2259	37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) &
2260	- 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) &
2261	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
2262	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
2263	10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) &
2264	- 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) &
2265	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
2266
2267	ENDDO
2268
2269	ENDIF
2270	!
2271	!-- Compute southside fluxes for the respective boundary.
2272	IF ( j == nysv ) THEN
2273
2274	DO k = nzb+1, nzb_max
2275
2276	ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
2277	ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
2278	ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
2279
2280	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2281	flux_s_v(k,tn) = v_comp_l * ( &
2282	( 37.0_wp * ibit23 * adv_mom_5 &
2283	+ 7.0_wp * ibit22 * adv_mom_3 &
2284	+ ibit21 * adv_mom_1 &
2285	) * &
2286	( v(k,j,i) + v(k,j-1,i) ) &
2287	- ( 8.0_wp * ibit23 * adv_mom_5 &
2288	+ ibit22 * adv_mom_3 &
2289	) * &
2290	( v(k,j+1,i) + v(k,j-2,i) ) &
2291	+ ( ibit23 * adv_mom_5 &
2292	) * &
2293	( v(k,j+2,i) + v(k,j-3,i) ) &
2294	)
2295
2296	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2297	( 10.0_wp * ibit23 * adv_mom_5 &
2298	+ 3.0_wp * ibit22 * adv_mom_3 &
2299	+ ibit21 * adv_mom_1 &
2300	) * &
2301	( v(k,j,i) - v(k,j-1,i) ) &
2302	- ( 5.0_wp * ibit23 * adv_mom_5 &
2303	+ ibit22 * adv_mom_3 &
2304	) * &
2305	( v(k,j+1,i) - v(k,j-2,i) ) &
2306	+ ( ibit23 * adv_mom_5 &
2307	) * &
2308	( v(k,j+2,i) - v(k,j-3,i) ) &
2309	)
2310
2311	ENDDO
2312
2313	DO k = nzb_max+1, nzt
2314
2315	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2316	flux_s_v(k,tn) = v_comp_l * ( &
2317	37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) &
2318	- 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) &
2319	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
2320	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2321	10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) &
2322	- 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) &
2323	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
2324
2325	ENDDO
2326
2327	ENDIF
2328
2329	flux_t(0) = 0.0_wp
2330	diss_t(0) = 0.0_wp
2331	flux_d = 0.0_wp
2332	diss_d = 0.0_wp
2333	!
2334	!-- Now compute the fluxes and tendency terms for the horizontal and
2335	!-- verical parts.
2336	DO k = nzb+1, nzb_max
2337
2338	ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
2339	ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
2340	ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
2341
2342	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2343	flux_r(k) = u_comp * ( &
2344	( 37.0_wp * ibit20 * adv_mom_5 &
2345	+ 7.0_wp * ibit19 * adv_mom_3 &
2346	+ ibit18 * adv_mom_1 &
2347	) * &
2348	( v(k,j,i+1) + v(k,j,i) ) &
2349	- ( 8.0_wp * ibit20 * adv_mom_5 &
2350	+ ibit19 * adv_mom_3 &
2351	) * &
2352	( v(k,j,i+2) + v(k,j,i-1) ) &
2353	+ ( ibit20 * adv_mom_5 &
2354	) * &
2355	( v(k,j,i+3) + v(k,j,i-2) ) &
2356	)
2357
2358	diss_r(k) = - ABS( u_comp ) * ( &
2359	( 10.0_wp * ibit20 * adv_mom_5 &
2360	+ 3.0_wp * ibit19 * adv_mom_3 &
2361	+ ibit18 * adv_mom_1 &
2362	) * &
2363	( v(k,j,i+1) - v(k,j,i) ) &
2364	- ( 5.0_wp * ibit20 * adv_mom_5 &
2365	+ ibit19 * adv_mom_3 &
2366	) * &
2367	( v(k,j,i+2) - v(k,j,i-1) ) &
2368	+ ( ibit20 * adv_mom_5 &
2369	) * &
2370	( v(k,j,i+3) - v(k,j,i-2) ) &
2371	)
2372
2373	ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
2374	ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
2375	ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
2376
2377
2378	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2379	flux_n(k) = ( v_comp(k) - gv ) * ( &
2380	( 37.0_wp * ibit23 * adv_mom_5 &
2381	+ 7.0_wp * ibit22 * adv_mom_3 &
2382	+ ibit21 * adv_mom_1 &
2383	) * &
2384	( v(k,j+1,i) + v(k,j,i) ) &
2385	- ( 8.0_wp * ibit23 * adv_mom_5 &
2386	+ ibit22 * adv_mom_3 &
2387	) * &
2388	( v(k,j+2,i) + v(k,j-1,i) ) &
2389	+ ( ibit23 * adv_mom_5 &
2390	) * &
2391	( v(k,j+3,i) + v(k,j-2,i) ) &
2392	)
2393
2394	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2395	( 10.0_wp * ibit23 * adv_mom_5 &
2396	+ 3.0_wp * ibit22 * adv_mom_3 &
2397	+ ibit21 * adv_mom_1 &
2398	) * &
2399	( v(k,j+1,i) - v(k,j,i) ) &
2400	- ( 5.0_wp * ibit23 * adv_mom_5 &
2401	+ ibit22 * adv_mom_3 &
2402	) * &
2403	( v(k,j+2,i) - v(k,j-1,i) ) &
2404	+ ( ibit23 * adv_mom_5 &
2405	) * &
2406	( v(k,j+3,i) - v(k,j-2,i) ) &
2407	)
2408	!
2409	!-- k index has to be modified near bottom and top, else array
2410	!-- subscripts will be exceeded.
2411	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
2412	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
2413	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
2414
2415	k_ppp = k + 3 * ibit26
2416	k_pp = k + 2 * ( 1 - ibit24 )
2417	k_mm = k - 2 * ibit26
2418
2419	w_comp = w(k,j-1,i) + w(k,j,i)
2420	flux_t(k) = w_comp * rho_air_zw(k) * ( &
2421	( 37.0_wp * ibit26 * adv_mom_5 &
2422	+ 7.0_wp * ibit25 * adv_mom_3 &
2423	+ ibit24 * adv_mom_1 &
2424	) * &
2425	( v(k+1,j,i) + v(k,j,i) ) &
2426	- ( 8.0_wp * ibit26 * adv_mom_5 &
2427	+ ibit25 * adv_mom_3 &
2428	) * &
2429	( v(k_pp,j,i) + v(k-1,j,i) ) &
2430	+ ( ibit26 * adv_mom_5 &
2431	) * &
2432	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2433	)
2434
2435	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
2436	( 10.0_wp * ibit26 * adv_mom_5 &
2437	+ 3.0_wp * ibit25 * adv_mom_3 &
2438	+ ibit24 * adv_mom_1 &
2439	) * &
2440	( v(k+1,j,i) - v(k,j,i) ) &
2441	- ( 5.0_wp * ibit26 * adv_mom_5 &
2442	+ ibit25 * adv_mom_3 &
2443	) * &
2444	( v(k_pp,j,i) - v(k-1,j,i) ) &
2445	+ ( ibit26 * adv_mom_5 &
2446	) * &
2447	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2448	)
2449	!
2450	!-- Calculate the divergence of the velocity field. A respective
2451	!-- correction is needed to overcome numerical instabilities introduced
2452	!-- by a not sufficient reduction of divergences near topography.
2453	div = ( ( ( u_comp + gu ) &
2454	* ( ibit18 + ibit19 + ibit20 ) &
2455	- ( u(k,j-1,i) + u(k,j,i) ) &
2456	* ( IBITS(advc_flags_1(k,j,i-1),18,1) &
2457	+ IBITS(advc_flags_1(k,j,i-1),19,1) &
2458	+ IBITS(advc_flags_1(k,j,i-1),20,1) &
2459	) &
2460	) * rho_air(k) * ddx &
2461	+ ( v_comp(k) &
2462	* ( ibit21 + ibit22 + ibit23 ) &
2463	- ( v(k,j,i) + v(k,j-1,i) ) &
2464	* ( IBITS(advc_flags_1(k,j-1,i),21,1) &
2465	+ IBITS(advc_flags_1(k,j-1,i),22,1) &
2466	+ IBITS(advc_flags_1(k,j-1,i),23,1) &
2467	) &
2468	) * rho_air(k) * ddy &
2469	+ ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 ) &
2470	- ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
2471	* ( IBITS(advc_flags_1(k-1,j,i),24,1) &
2472	+ IBITS(advc_flags_1(k-1,j,i),25,1) &
2473	+ IBITS(advc_flags_1(k-1,j,i),26,1) &
2474	) &
2475	) * ddzw(k) &
2476	) * 0.5_wp
2477
2478
2479	tend(k,j,i) = tend(k,j,i) - ( &
2480	( flux_r(k) + diss_r(k) &
2481	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2482	+ ( flux_n(k) + diss_n(k) &
2483	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2484	+ ( ( flux_t(k) + diss_t(k) ) &
2485	- ( flux_d + diss_d ) &
2486	) * drho_air(k) * ddzw(k) &
2487	) + v(k,j,i) * div
2488
2489	flux_l_v(k,j,tn) = flux_r(k)
2490	diss_l_v(k,j,tn) = diss_r(k)
2491	flux_s_v(k,tn) = flux_n(k)
2492	diss_s_v(k,tn) = diss_n(k)
2493	flux_d = flux_t(k)
2494	diss_d = diss_t(k)
2495
2496	!
2497	!-- Statistical Evaluation of v'v'. The factor has to be applied for
2498	!-- right evaluation when gallilei_trans = .T. .
2499	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
2500	+ ( flux_n(k) &
2501	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2502	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
2503	+ diss_n(k) &
2504	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2505	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
2506	) * weight_substep(intermediate_timestep_count)
2507	!
2508	!-- Statistical Evaluation of w'u'.
2509	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
2510	+ ( flux_t(k) &
2511	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2512	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2513	+ diss_t(k) &
2514	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2515	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2516	) * weight_substep(intermediate_timestep_count)
2517
2518	ENDDO
2519
2520	DO k = nzb_max+1, nzt
2521
2522	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2523	flux_r(k) = u_comp * ( &
2524	37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &
2525	- 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &
2526	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
2527
2528	diss_r(k) = - ABS( u_comp ) * ( &
2529	10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &
2530	- 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &
2531	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
2532
2533
2534	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2535	flux_n(k) = ( v_comp(k) - gv ) * ( &
2536	37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &
2537	- 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &
2538	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
2539
2540	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2541	10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &
2542	- 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &
2543	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
2544	!
2545	!-- k index has to be modified near bottom and top, else array
2546	!-- subscripts will be exceeded.
2547	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
2548	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
2549	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
2550
2551	k_ppp = k + 3 * ibit26
2552	k_pp = k + 2 * ( 1 - ibit24 )
2553	k_mm = k - 2 * ibit26
2554
2555	w_comp = w(k,j-1,i) + w(k,j,i)
2556	flux_t(k) = w_comp * rho_air_zw(k) * ( &
2557	( 37.0_wp * ibit26 * adv_mom_5 &
2558	+ 7.0_wp * ibit25 * adv_mom_3 &
2559	+ ibit24 * adv_mom_1 &
2560	) * &
2561	( v(k+1,j,i) + v(k,j,i) ) &
2562	- ( 8.0_wp * ibit26 * adv_mom_5 &
2563	+ ibit25 * adv_mom_3 &
2564	) * &
2565	( v(k_pp,j,i) + v(k-1,j,i) ) &
2566	+ ( ibit26 * adv_mom_5 &
2567	) * &
2568	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2569	)
2570
2571	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
2572	( 10.0_wp * ibit26 * adv_mom_5 &
2573	+ 3.0_wp * ibit25 * adv_mom_3 &
2574	+ ibit24 * adv_mom_1 &
2575	) * &
2576	( v(k+1,j,i) - v(k,j,i) ) &
2577	- ( 5.0_wp * ibit26 * adv_mom_5 &
2578	+ ibit25 * adv_mom_3 &
2579	) * &
2580	( v(k_pp,j,i) - v(k-1,j,i) ) &
2581	+ ( ibit26 * adv_mom_5 &
2582	) * &
2583	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2584	)
2585	!
2586	!-- Calculate the divergence of the velocity field. A respective
2587	!-- correction is needed to overcome numerical instabilities introduced
2588	!-- by a not sufficient reduction of divergences near topography.
2589	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
2590	* rho_air(k)&
2591	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
2592	* rho_air(k)&
2593	+ ( w_comp * rho_air_zw(k) - &
2594	( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
2595	) * ddzw(k) &
2596	) * 0.5_wp
2597
2598	tend(k,j,i) = tend(k,j,i) - ( &
2599	( flux_r(k) + diss_r(k) &
2600	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2601	+ ( flux_n(k) + diss_n(k) &
2602	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2603	+ ( ( flux_t(k) + diss_t(k) ) &
2604	- ( flux_d + diss_d ) &
2605	) * drho_air(k) * ddzw(k) &
2606	) + v(k,j,i) * div
2607
2608	flux_l_v(k,j,tn) = flux_r(k)
2609	diss_l_v(k,j,tn) = diss_r(k)
2610	flux_s_v(k,tn) = flux_n(k)
2611	diss_s_v(k,tn) = diss_n(k)
2612	flux_d = flux_t(k)
2613	diss_d = diss_t(k)
2614
2615	!
2616	!-- Statistical Evaluation of v'v'. The factor has to be applied for
2617	!-- right evaluation when gallilei_trans = .T. .
2618	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
2619	+ ( flux_n(k) &
2620	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2621	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
2622	+ diss_n(k) &
2623	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2624	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
2625	) * weight_substep(intermediate_timestep_count)
2626	!
2627	!-- Statistical Evaluation of w'u'.
2628	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
2629	+ ( flux_t(k) &
2630	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2631	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
2632	+ diss_t(k) &
2633	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
2634	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
2635	) * weight_substep(intermediate_timestep_count)
2636
2637	ENDDO
2638	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
2639
2640
2641	END SUBROUTINE advec_v_ws_ij
2642
2643
2644
2645	!------------------------------------------------------------------------------!
2646	! Description:
2647	! ------------
2648	!> Advection of w-component - Call for grid point i,j
2649	!------------------------------------------------------------------------------!
2650	SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
2651
2652	USE arrays_3d, &
2653	ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,&
2654	drho_air_zw, rho_air, rho_air_zw
2655
2656	USE constants, &
2657	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
2658
2659	USE control_parameters, &
2660	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
2661
2662	USE grid_variables, &
2663	ONLY: ddx, ddy
2664
2665	USE indices, &
2666	ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1, &
2667	advc_flags_2
2668
2669	USE kinds
2670
2671	USE statistics, &
2672	ONLY: hom, sums_ws2_ws_l, weight_substep
2673
2674	IMPLICIT NONE
2675
2676	INTEGER(iwp) :: i !<
2677	INTEGER(iwp) :: ibit27 !<
2678	INTEGER(iwp) :: ibit28 !<
2679	INTEGER(iwp) :: ibit29 !<
2680	INTEGER(iwp) :: ibit30 !<
2681	INTEGER(iwp) :: ibit31 !<
2682	INTEGER(iwp) :: ibit32 !<
2683	INTEGER(iwp) :: ibit33 !<
2684	INTEGER(iwp) :: ibit34 !<
2685	INTEGER(iwp) :: ibit35 !<
2686	INTEGER(iwp) :: i_omp !<
2687	INTEGER(iwp) :: j !<
2688	INTEGER(iwp) :: k !<
2689	INTEGER(iwp) :: k_mm !<
2690	INTEGER(iwp) :: k_pp !<
2691	INTEGER(iwp) :: k_ppp !<
2692	INTEGER(iwp) :: tn !<
2693
2694	REAL(wp) :: diss_d !<
2695	REAL(wp) :: div !<
2696	REAL(wp) :: flux_d !<
2697	REAL(wp) :: gu !<
2698	REAL(wp) :: gv !<
2699	REAL(wp) :: u_comp !<
2700	REAL(wp) :: v_comp !<
2701	REAL(wp) :: w_comp !<
2702
2703	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
2704	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
2705	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
2706	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
2707	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
2708	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
2709
2710	gu = 2.0_wp * u_gtrans
2711	gv = 2.0_wp * v_gtrans
2712
2713	!
2714	!-- Compute southside fluxes for the respective boundary.
2715	IF ( j == nys ) THEN
2716
2717	DO k = nzb+1, nzb_max
2718	ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
2719	ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
2720	ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
2721
2722	v_comp = v(k+1,j,i) + v(k,j,i) - gv
2723	flux_s_w(k,tn) = v_comp * ( &
2724	( 37.0_wp * ibit32 * adv_mom_5 &
2725	+ 7.0_wp * ibit31 * adv_mom_3 &
2726	+ ibit30 * adv_mom_1 &
2727	) * &
2728	( w(k,j,i) + w(k,j-1,i) ) &
2729	- ( 8.0_wp * ibit32 * adv_mom_5 &
2730	+ ibit31 * adv_mom_3 &
2731	) * &
2732	( w(k,j+1,i) + w(k,j-2,i) ) &
2733	+ ( ibit32 * adv_mom_5 &
2734	) * &
2735	( w(k,j+2,i) + w(k,j-3,i) ) &
2736	)
2737
2738	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
2739	( 10.0_wp * ibit32 * adv_mom_5 &
2740	+ 3.0_wp * ibit31 * adv_mom_3 &
2741	+ ibit30 * adv_mom_1 &
2742	) * &
2743	( w(k,j,i) - w(k,j-1,i) ) &
2744	- ( 5.0_wp * ibit32 * adv_mom_5 &
2745	+ ibit31 * adv_mom_3 &
2746	) * &
2747	( w(k,j+1,i) - w(k,j-2,i) ) &
2748	+ ( ibit32 * adv_mom_5 &
2749	) * &
2750	( w(k,j+2,i) - w(k,j-3,i) ) &
2751	)
2752
2753	ENDDO
2754
2755	DO k = nzb_max+1, nzt
2756
2757	v_comp = v(k+1,j,i) + v(k,j,i) - gv
2758	flux_s_w(k,tn) = v_comp * ( &
2759	37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) &
2760	- 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) &
2761	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
2762	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
2763	10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) &
2764	- 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) &
2765	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
2766
2767	ENDDO
2768
2769	ENDIF
2770	!
2771	!-- Compute leftside fluxes for the respective boundary.
2772	IF ( i == i_omp ) THEN
2773
2774	DO k = nzb+1, nzb_max
2775
2776	ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
2777	ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
2778	ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
2779
2780	u_comp = u(k+1,j,i) + u(k,j,i) - gu
2781	flux_l_w(k,j,tn) = u_comp * ( &
2782	( 37.0_wp * ibit29 * adv_mom_5 &
2783	+ 7.0_wp * ibit28 * adv_mom_3 &
2784	+ ibit27 * adv_mom_1 &
2785	) * &
2786	( w(k,j,i) + w(k,j,i-1) ) &
2787	- ( 8.0_wp * ibit29 * adv_mom_5 &
2788	+ ibit28 * adv_mom_3 &
2789	) * &
2790	( w(k,j,i+1) + w(k,j,i-2) ) &
2791	+ ( ibit29 * adv_mom_5 &
2792	) * &
2793	( w(k,j,i+2) + w(k,j,i-3) ) &
2794	)
2795
2796	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
2797	( 10.0_wp * ibit29 * adv_mom_5 &
2798	+ 3.0_wp * ibit28 * adv_mom_3 &
2799	+ ibit27 * adv_mom_1 &
2800	) * &
2801	( w(k,j,i) - w(k,j,i-1) ) &
2802	- ( 5.0_wp * ibit29 * adv_mom_5 &
2803	+ ibit28 * adv_mom_3 &
2804	) * &
2805	( w(k,j,i+1) - w(k,j,i-2) ) &
2806	+ ( ibit29 * adv_mom_5 &
2807	) * &
2808	( w(k,j,i+2) - w(k,j,i-3) ) &
2809	)
2810
2811	ENDDO
2812
2813	DO k = nzb_max+1, nzt
2814
2815	u_comp = u(k+1,j,i) + u(k,j,i) - gu
2816	flux_l_w(k,j,tn) = u_comp * ( &
2817	37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) &
2818	- 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) &
2819	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
2820	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
2821	10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) &
2822	- 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) &
2823	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
2824
2825	ENDDO
2826
2827	ENDIF
2828	!
2829	!-- The lower flux has to be calculated explicetely for the tendency at
2830	!-- the first w-level. For topography wall this is done implicitely by
2831	!-- advc_flags_1.
2832	k = nzb + 1
2833	w_comp = w(k,j,i) + w(k-1,j,i)
2834	flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
2835	diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
2836	flux_d = flux_t(0)
2837	diss_d = diss_t(0)
2838	!
2839	!-- Now compute the fluxes and tendency terms for the horizontal
2840	!-- and vertical parts.
2841	DO k = nzb+1, nzb_max
2842
2843	ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
2844	ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
2845	ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
2846
2847	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
2848	flux_r(k) = u_comp * ( &
2849	( 37.0_wp * ibit29 * adv_mom_5 &
2850	+ 7.0_wp * ibit28 * adv_mom_3 &
2851	+ ibit27 * adv_mom_1 &
2852	) * &
2853	( w(k,j,i+1) + w(k,j,i) ) &
2854	- ( 8.0_wp * ibit29 * adv_mom_5 &
2855	+ ibit28 * adv_mom_3 &
2856	) * &
2857	( w(k,j,i+2) + w(k,j,i-1) ) &
2858	+ ( ibit29 * adv_mom_5 &
2859	) * &
2860	( w(k,j,i+3) + w(k,j,i-2) ) &
2861	)
2862
2863	diss_r(k) = - ABS( u_comp ) * ( &
2864	( 10.0_wp * ibit29 * adv_mom_5 &
2865	+ 3.0_wp * ibit28 * adv_mom_3 &
2866	+ ibit27 * adv_mom_1 &
2867	) * &
2868	( w(k,j,i+1) - w(k,j,i) ) &
2869	- ( 5.0_wp * ibit29 * adv_mom_5 &
2870	+ ibit28 * adv_mom_3 &
2871	) * &
2872	( w(k,j,i+2) - w(k,j,i-1) ) &
2873	+ ( ibit29 * adv_mom_5 &
2874	) * &
2875	( w(k,j,i+3) - w(k,j,i-2) ) &
2876	)
2877
2878	ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
2879	ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
2880	ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
2881
2882	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
2883	flux_n(k) = v_comp * ( &
2884	( 37.0_wp * ibit32 * adv_mom_5 &
2885	+ 7.0_wp * ibit31 * adv_mom_3 &
2886	+ ibit30 * adv_mom_1 &
2887	) * &
2888	( w(k,j+1,i) + w(k,j,i) ) &
2889	- ( 8.0_wp * ibit32 * adv_mom_5 &
2890	+ ibit31 * adv_mom_3 &
2891	) * &
2892	( w(k,j+2,i) + w(k,j-1,i) ) &
2893	+ ( ibit32 * adv_mom_5 &
2894	) * &
2895	( w(k,j+3,i) + w(k,j-2,i) ) &
2896	)
2897
2898	diss_n(k) = - ABS( v_comp ) * ( &
2899	( 10.0_wp * ibit32 * adv_mom_5 &
2900	+ 3.0_wp * ibit31 * adv_mom_3 &
2901	+ ibit30 * adv_mom_1 &
2902	) * &
2903	( w(k,j+1,i) - w(k,j,i) ) &
2904	- ( 5.0_wp * ibit32 * adv_mom_5 &
2905	+ ibit31 * adv_mom_3 &
2906	) * &
2907	( w(k,j+2,i) - w(k,j-1,i) ) &
2908	+ ( ibit32 * adv_mom_5 &
2909	) * &
2910	( w(k,j+3,i) - w(k,j-2,i) ) &
2911	)
2912	!
2913	!-- k index has to be modified near bottom and top, else array
2914	!-- subscripts will be exceeded.
2915	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
2916	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
2917	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
2918
2919	k_ppp = k + 3 * ibit35
2920	k_pp = k + 2 * ( 1 - ibit33 )
2921	k_mm = k - 2 * ibit35
2922
2923	w_comp = w(k+1,j,i) + w(k,j,i)
2924	flux_t(k) = w_comp * rho_air(k+1) * ( &
2925	( 37.0_wp * ibit35 * adv_mom_5 &
2926	+ 7.0_wp * ibit34 * adv_mom_3 &
2927	+ ibit33 * adv_mom_1 &
2928	) * &
2929	( w(k+1,j,i) + w(k,j,i) ) &
2930	- ( 8.0_wp * ibit35 * adv_mom_5 &
2931	+ ibit34 * adv_mom_3 &
2932	) * &
2933	( w(k_pp,j,i) + w(k-1,j,i) ) &
2934	+ ( ibit35 * adv_mom_5 &
2935	) * &
2936	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
2937	)
2938
2939	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
2940	( 10.0_wp * ibit35 * adv_mom_5 &
2941	+ 3.0_wp * ibit34 * adv_mom_3 &
2942	+ ibit33 * adv_mom_1 &
2943	) * &
2944	( w(k+1,j,i) - w(k,j,i) ) &
2945	- ( 5.0_wp * ibit35 * adv_mom_5 &
2946	+ ibit34 * adv_mom_3 &
2947	) * &
2948	( w(k_pp,j,i) - w(k-1,j,i) ) &
2949	+ ( ibit35 * adv_mom_5 &
2950	) * &
2951	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
2952	)
2953
2954	!
2955	!-- Calculate the divergence of the velocity field. A respective
2956	!-- correction is needed to overcome numerical instabilities introduced
2957	!-- by a not sufficient reduction of divergences near topography.
2958	div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) &
2959	- ( u(k+1,j,i) + u(k,j,i) ) &
2960	* ( IBITS(advc_flags_1(k,j,i-1),27,1) &
2961	+ IBITS(advc_flags_1(k,j,i-1),28,1) &
2962	+ IBITS(advc_flags_1(k,j,i-1),29,1) &
2963	) &
2964	) * rho_air_zw(k) * ddx &
2965	+ ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) &
2966	- ( v(k+1,j,i) + v(k,j,i) ) &
2967	* ( IBITS(advc_flags_1(k,j-1,i),30,1) &
2968	+ IBITS(advc_flags_1(k,j-1,i),31,1) &
2969	+ IBITS(advc_flags_2(k,j-1,i),0,1) &
2970	) &
2971	) * rho_air_zw(k) * ddy &
2972	+ ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) &
2973	- ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
2974	* ( IBITS(advc_flags_2(k-1,j,i),1,1) &
2975	+ IBITS(advc_flags_2(k-1,j,i),2,1) &
2976	+ IBITS(advc_flags_2(k-1,j,i),3,1) &
2977	) &
2978	) * ddzu(k+1) &
2979	) * 0.5_wp
2980
2981
2982	tend(k,j,i) = tend(k,j,i) - ( &
2983	( flux_r(k) + diss_r(k) &
2984	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
2985	+ ( flux_n(k) + diss_n(k) &
2986	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
2987	+ ( ( flux_t(k) + diss_t(k) ) &
2988	- ( flux_d + diss_d ) &
2989	) * drho_air_zw(k) * ddzu(k+1) &
2990	) + div * w(k,j,i)
2991
2992	flux_l_w(k,j,tn) = flux_r(k)
2993	diss_l_w(k,j,tn) = diss_r(k)
2994	flux_s_w(k,tn) = flux_n(k)
2995	diss_s_w(k,tn) = diss_n(k)
2996	flux_d = flux_t(k)
2997	diss_d = diss_t(k)
2998	!
2999	!-- Statistical Evaluation of w'w'.
3000	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
3001	+ ( flux_t(k) &
3002	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3003	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
3004	+ diss_t(k) &
3005	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3006	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
3007	) * weight_substep(intermediate_timestep_count)
3008
3009	ENDDO
3010
3011	DO k = nzb_max+1, nzt
3012
3013	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
3014	flux_r(k) = u_comp * ( &
3015	37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) &
3016	- 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) &
3017	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
3018
3019	diss_r(k) = - ABS( u_comp ) * ( &
3020	10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) &
3021	- 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) &
3022	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
3023
3024	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
3025	flux_n(k) = v_comp * ( &
3026	37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) &
3027	- 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) &
3028	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
3029
3030	diss_n(k) = - ABS( v_comp ) * ( &
3031	10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) &
3032	- 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) &
3033	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
3034	!
3035	!-- k index has to be modified near bottom and top, else array
3036	!-- subscripts will be exceeded.
3037	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
3038	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
3039	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
3040
3041	k_ppp = k + 3 * ibit35
3042	k_pp = k + 2 * ( 1 - ibit33 )
3043	k_mm = k - 2 * ibit35
3044
3045	w_comp = w(k+1,j,i) + w(k,j,i)
3046	flux_t(k) = w_comp * rho_air(k+1) * ( &
3047	( 37.0_wp * ibit35 * adv_mom_5 &
3048	+ 7.0_wp * ibit34 * adv_mom_3 &
3049	+ ibit33 * adv_mom_1 &
3050	) * &
3051	( w(k+1,j,i) + w(k,j,i) ) &
3052	- ( 8.0_wp * ibit35 * adv_mom_5 &
3053	+ ibit34 * adv_mom_3 &
3054	) * &
3055	( w(k_pp,j,i) + w(k-1,j,i) ) &
3056	+ ( ibit35 * adv_mom_5 &
3057	) * &
3058	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
3059	)
3060
3061	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
3062	( 10.0_wp * ibit35 * adv_mom_5 &
3063	+ 3.0_wp * ibit34 * adv_mom_3 &
3064	+ ibit33 * adv_mom_1 &
3065	) * &
3066	( w(k+1,j,i) - w(k,j,i) ) &
3067	- ( 5.0_wp * ibit35 * adv_mom_5 &
3068	+ ibit34 * adv_mom_3 &
3069	) * &
3070	( w(k_pp,j,i) - w(k-1,j,i) ) &
3071	+ ( ibit35 * adv_mom_5 &
3072	) * &
3073	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
3074	)
3075	!
3076	!-- Calculate the divergence of the velocity field. A respective
3077	!-- correction is needed to overcome numerical instabilities introduced
3078	!-- by a not sufficient reduction of divergences near topography.
3079	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
3080	* rho_air_zw(k) &
3081	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
3082	* rho_air_zw(k) &
3083	+ ( w_comp * rho_air(k+1) - &
3084	( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
3085	) * ddzu(k+1) &
3086	) * 0.5_wp
3087
3088	tend(k,j,i) = tend(k,j,i) - ( &
3089	( flux_r(k) + diss_r(k) &
3090	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
3091	+ ( flux_n(k) + diss_n(k) &
3092	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
3093	+ ( ( flux_t(k) + diss_t(k) ) &
3094	- ( flux_d + diss_d ) &
3095	) * drho_air_zw(k) * ddzu(k+1) &
3096	) + div * w(k,j,i)
3097
3098	flux_l_w(k,j,tn) = flux_r(k)
3099	diss_l_w(k,j,tn) = diss_r(k)
3100	flux_s_w(k,tn) = flux_n(k)
3101	diss_s_w(k,tn) = diss_n(k)
3102	flux_d = flux_t(k)
3103	diss_d = diss_t(k)
3104	!
3105	!-- Statistical Evaluation of w'w'.
3106	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
3107	+ ( flux_t(k) &
3108	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3109	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
3110	+ diss_t(k) &
3111	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
3112	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
3113	) * weight_substep(intermediate_timestep_count)
3114
3115	ENDDO
3116
3117
3118	END SUBROUTINE advec_w_ws_ij
3119
3120
3121	!------------------------------------------------------------------------------!
3122	! Description:
3123	! ------------
3124	!> Scalar advection - Call for all grid points
3125	!------------------------------------------------------------------------------!
3126	SUBROUTINE advec_s_ws( sk, sk_char )
3127
3128	USE arrays_3d, &
3129	ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
3130
3131	USE constants, &
3132	ONLY: adv_sca_1, adv_sca_3, adv_sca_5
3133
3134	USE control_parameters, &
3135	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
3136
3137	USE grid_variables, &
3138	ONLY: ddx, ddy
3139
3140	USE indices, &
3141	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
3142	nzt, advc_flags_1
3143
3144	USE kinds
3145
3146	USE statistics, &
3147	ONLY: hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, &
3148	sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l, &
3149	sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep
3150
3151
3152
3153	IMPLICIT NONE
3154
3155	CHARACTER (LEN = *), INTENT(IN) :: sk_char !<
3156
3157	INTEGER(iwp) :: i !<
3158	INTEGER(iwp) :: ibit0 !<
3159	INTEGER(iwp) :: ibit1 !<
3160	INTEGER(iwp) :: ibit2 !<
3161	INTEGER(iwp) :: ibit3 !<
3162	INTEGER(iwp) :: ibit4 !<
3163	INTEGER(iwp) :: ibit5 !<
3164	INTEGER(iwp) :: ibit6 !<
3165	INTEGER(iwp) :: ibit7 !<
3166	INTEGER(iwp) :: ibit8 !<
3167	INTEGER(iwp) :: j !<
3168	INTEGER(iwp) :: k !<
3169	INTEGER(iwp) :: k_mm !<
3170	INTEGER(iwp) :: k_mmm !<
3171	INTEGER(iwp) :: k_pp !<
3172	INTEGER(iwp) :: k_ppp !<
3173	INTEGER(iwp) :: tn = 0 !<
3174
3175	#if defined( __nopointer )
3176	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !<
3177	#else
3178	REAL(wp), DIMENSION(:,:,:), POINTER :: sk !<
3179	#endif
3180
3181	REAL(wp) :: diss_d !<
3182	REAL(wp) :: div !<
3183	REAL(wp) :: flux_d !<
3184	REAL(wp) :: u_comp !<
3185	REAL(wp) :: v_comp !<
3186
3187	REAL(wp), DIMENSION(nzb:nzt) :: diss_n !<
3188	REAL(wp), DIMENSION(nzb:nzt) :: diss_r !<
3189	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
3190	REAL(wp), DIMENSION(nzb:nzt) :: flux_n !<
3191	REAL(wp), DIMENSION(nzb:nzt) :: flux_r !<
3192	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
3193
3194	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local !<
3195	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local !<
3196
3197	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local !<
3198	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local !<
3199
3200
3201	!
3202	!-- Compute the fluxes for the whole left boundary of the processor domain.
3203	i = nxl
3204	DO j = nys, nyn
3205
3206	DO k = nzb+1, nzb_max
3207
3208	ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
3209	ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
3210	ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
3211
3212	u_comp = u(k,j,i) - u_gtrans
3213	swap_flux_x_local(k,j) = u_comp * ( &
3214	( 37.0_wp * ibit2 * adv_sca_5 &
3215	+ 7.0_wp * ibit1 * adv_sca_3 &
3216	+ ibit0 * adv_sca_1 &
3217	) * &
3218	( sk(k,j,i) + sk(k,j,i-1) ) &
3219	- ( 8.0_wp * ibit2 * adv_sca_5 &
3220	+ ibit1 * adv_sca_3 &
3221	) * &
3222	( sk(k,j,i+1) + sk(k,j,i-2) ) &
3223	+ ( ibit2 * adv_sca_5 &
3224	) * &
3225	( sk(k,j,i+2) + sk(k,j,i-3) ) &
3226	)
3227
3228	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
3229	( 10.0_wp * ibit2 * adv_sca_5 &
3230	+ 3.0_wp * ibit1 * adv_sca_3 &
3231	+ ibit0 * adv_sca_1 &
3232	) * &
3233	( sk(k,j,i) - sk(k,j,i-1) ) &
3234	- ( 5.0_wp * ibit2 * adv_sca_5 &
3235	+ ibit1 * adv_sca_3 &
3236	) * &
3237	( sk(k,j,i+1) - sk(k,j,i-2) ) &
3238	+ ( ibit2 * adv_sca_5 &
3239	) * &
3240	( sk(k,j,i+2) - sk(k,j,i-3) ) &
3241	)
3242
3243	ENDDO
3244
3245	DO k = nzb_max+1, nzt
3246
3247	u_comp = u(k,j,i) - u_gtrans
3248	swap_flux_x_local(k,j) = u_comp * ( &
3249	37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &
3250	- 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
3251	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
3252	) * adv_sca_5
3253
3254	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
3255	10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &
3256	- 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
3257	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
3258	) * adv_sca_5
3259
3260	ENDDO
3261
3262	ENDDO
3263
3264	DO i = nxl, nxr
3265
3266	j = nys
3267	DO k = nzb+1, nzb_max
3268
3269	ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
3270	ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
3271	ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
3272
3273	v_comp = v(k,j,i) - v_gtrans
3274	swap_flux_y_local(k) = v_comp * ( &
3275	( 37.0_wp * ibit5 * adv_sca_5 &
3276	+ 7.0_wp * ibit4 * adv_sca_3 &
3277	+ ibit3 * adv_sca_1 &
3278	) * &
3279	( sk(k,j,i) + sk(k,j-1,i) ) &
3280	- ( 8.0_wp * ibit5 * adv_sca_5 &
3281	+ ibit4 * adv_sca_3 &
3282	) * &
3283	( sk(k,j+1,i) + sk(k,j-2,i) ) &
3284	+ ( ibit5 * adv_sca_5 &
3285	) * &
3286	( sk(k,j+2,i) + sk(k,j-3,i) ) &
3287	)
3288
3289	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
3290	( 10.0_wp * ibit5 * adv_sca_5 &
3291	+ 3.0_wp * ibit4 * adv_sca_3 &
3292	+ ibit3 * adv_sca_1 &
3293	) * &
3294	( sk(k,j,i) - sk(k,j-1,i) ) &
3295	- ( 5.0_wp * ibit5 * adv_sca_5 &
3296	+ ibit4 * adv_sca_3 &
3297	) * &
3298	( sk(k,j+1,i) - sk(k,j-2,i) ) &
3299	+ ( ibit5 * adv_sca_5 &
3300	) * &
3301	( sk(k,j+2,i) - sk(k,j-3,i) ) &
3302	)
3303
3304	ENDDO
3305	!
3306	!-- Above to the top of the highest topography. No degradation necessary.
3307	DO k = nzb_max+1, nzt
3308
3309	v_comp = v(k,j,i) - v_gtrans
3310	swap_flux_y_local(k) = v_comp * ( &
3311	37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &
3312	- 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
3313	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
3314	) * adv_sca_5
3315	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
3316	10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &
3317	- 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
3318	+ sk(k,j+2,i) - sk(k,j-3,i) &
3319	) * adv_sca_5
3320
3321	ENDDO
3322
3323	DO j = nys, nyn
3324
3325	flux_t(0) = 0.0_wp
3326	diss_t(0) = 0.0_wp
3327	flux_d = 0.0_wp
3328	diss_d = 0.0_wp
3329
3330	DO k = nzb+1, nzb_max
3331
3332	ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
3333	ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
3334	ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
3335
3336	u_comp = u(k,j,i+1) - u_gtrans
3337	flux_r(k) = u_comp * ( &
3338	( 37.0_wp * ibit2 * adv_sca_5 &
3339	+ 7.0_wp * ibit1 * adv_sca_3 &
3340	+ ibit0 * adv_sca_1 &
3341	) * &
3342	( sk(k,j,i+1) + sk(k,j,i) ) &
3343	- ( 8.0_wp * ibit2 * adv_sca_5 &
3344	+ ibit1 * adv_sca_3 &
3345	) * &
3346	( sk(k,j,i+2) + sk(k,j,i-1) ) &
3347	+ ( ibit2 * adv_sca_5 &
3348	) * &
3349	( sk(k,j,i+3) + sk(k,j,i-2) ) &
3350	)
3351
3352	diss_r(k) = -ABS( u_comp ) * ( &
3353	( 10.0_wp * ibit2 * adv_sca_5 &
3354	+ 3.0_wp * ibit1 * adv_sca_3 &
3355	+ ibit0 * adv_sca_1 &
3356	) * &
3357	( sk(k,j,i+1) - sk(k,j,i) ) &
3358	- ( 5.0_wp * ibit2 * adv_sca_5 &
3359	+ ibit1 * adv_sca_3 &
3360	) * &
3361	( sk(k,j,i+2) - sk(k,j,i-1) ) &
3362	+ ( ibit2 * adv_sca_5 &
3363	) * &
3364	( sk(k,j,i+3) - sk(k,j,i-2) ) &
3365	)
3366
3367	ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
3368	ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
3369	ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
3370
3371	v_comp = v(k,j+1,i) - v_gtrans
3372	flux_n(k) = v_comp * ( &
3373	( 37.0_wp * ibit5 * adv_sca_5 &
3374	+ 7.0_wp * ibit4 * adv_sca_3 &
3375	+ ibit3 * adv_sca_1 &
3376	) * &
3377	( sk(k,j+1,i) + sk(k,j,i) ) &
3378	- ( 8.0_wp * ibit5 * adv_sca_5 &
3379	+ ibit4 * adv_sca_3 &
3380	) * &
3381	( sk(k,j+2,i) + sk(k,j-1,i) ) &
3382	+ ( ibit5 * adv_sca_5 &
3383	) * &
3384	( sk(k,j+3,i) + sk(k,j-2,i) ) &
3385	)
3386
3387	diss_n(k) = -ABS( v_comp ) * ( &
3388	( 10.0_wp * ibit5 * adv_sca_5 &
3389	+ 3.0_wp * ibit4 * adv_sca_3 &
3390	+ ibit3 * adv_sca_1 &
3391	) * &
3392	( sk(k,j+1,i) - sk(k,j,i) ) &
3393	- ( 5.0_wp * ibit5 * adv_sca_5 &
3394	+ ibit4 * adv_sca_3 &
3395	) * &
3396	( sk(k,j+2,i) - sk(k,j-1,i) ) &
3397	+ ( ibit5 * adv_sca_5 &
3398	) * &
3399	( sk(k,j+3,i) - sk(k,j-2,i) ) &
3400	)
3401	!
3402	!-- k index has to be modified near bottom and top, else array
3403	!-- subscripts will be exceeded.
3404	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
3405	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
3406	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
3407
3408	k_ppp = k + 3 * ibit8
3409	k_pp = k + 2 * ( 1 - ibit6 )
3410	k_mm = k - 2 * ibit8
3411
3412
3413	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
3414	( 37.0_wp * ibit8 * adv_sca_5 &
3415	+ 7.0_wp * ibit7 * adv_sca_3 &
3416	+ ibit6 * adv_sca_1 &
3417	) * &
3418	( sk(k+1,j,i) + sk(k,j,i) ) &
3419	- ( 8.0_wp * ibit8 * adv_sca_5 &
3420	+ ibit7 * adv_sca_3 &
3421	) * &
3422	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
3423	+ ( ibit8 * adv_sca_5 &
3424	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
3425	)
3426
3427	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
3428	( 10.0_wp * ibit8 * adv_sca_5 &
3429	+ 3.0_wp * ibit7 * adv_sca_3 &
3430	+ ibit6 * adv_sca_1 &
3431	) * &
3432	( sk(k+1,j,i) - sk(k,j,i) ) &
3433	- ( 5.0_wp * ibit8 * adv_sca_5 &
3434	+ ibit7 * adv_sca_3 &
3435	) * &
3436	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
3437	+ ( ibit8 * adv_sca_5 &
3438	) * &
3439	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
3440	)
3441	!
3442	!-- Calculate the divergence of the velocity field. A respective
3443	!-- correction is needed to overcome numerical instabilities caused
3444	!-- by a not sufficient reduction of divergences near topography.
3445	div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &
3446	- u(k,j,i) * ( IBITS(advc_flags_1(k,j,i-1),0,1) &
3447	+ IBITS(advc_flags_1(k,j,i-1),1,1) &
3448	+ IBITS(advc_flags_1(k,j,i-1),2,1) &
3449	) &
3450	) * rho_air(k) * ddx &
3451	+ ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &
3452	- v(k,j,i) * ( IBITS(advc_flags_1(k,j-1,i),3,1) &
3453	+ IBITS(advc_flags_1(k,j-1,i),4,1) &
3454	+ IBITS(advc_flags_1(k,j-1,i),5,1) &
3455	) &
3456	) * rho_air(k) * ddy &
3457	+ ( w(k,j,i) * rho_air_zw(k) * &
3458	( ibit6 + ibit7 + ibit8 ) &
3459	- w(k-1,j,i) * rho_air_zw(k-1) * &
3460	( IBITS(advc_flags_1(k-1,j,i),6,1) &
3461	+ IBITS(advc_flags_1(k-1,j,i),7,1) &
3462	+ IBITS(advc_flags_1(k-1,j,i),8,1) &
3463	) &
3464	) * ddzw(k)
3465
3466
3467	tend(k,j,i) = tend(k,j,i) - ( &
3468	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
3469	swap_diss_x_local(k,j) ) * ddx &
3470	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
3471	swap_diss_y_local(k) ) * ddy &
3472	+ ( ( flux_t(k) + diss_t(k) ) - &
3473	( flux_d + diss_d ) &
3474	) * drho_air(k) * ddzw(k) &
3475	) + sk(k,j,i) * div
3476
3477	swap_flux_y_local(k) = flux_n(k)
3478	swap_diss_y_local(k) = diss_n(k)
3479	swap_flux_x_local(k,j) = flux_r(k)
3480	swap_diss_x_local(k,j) = diss_r(k)
3481	flux_d = flux_t(k)
3482	diss_d = diss_t(k)
3483
3484	ENDDO
3485
3486	DO k = nzb_max+1, nzt
3487
3488	u_comp = u(k,j,i+1) - u_gtrans
3489	flux_r(k) = u_comp * ( &
3490	37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) &
3491	- 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
3492	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
3493	diss_r(k) = -ABS( u_comp ) * ( &
3494	10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) &
3495	- 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
3496	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
3497
3498	v_comp = v(k,j+1,i) - v_gtrans
3499	flux_n(k) = v_comp * ( &
3500	37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) &
3501	- 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
3502	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
3503	diss_n(k) = -ABS( v_comp ) * ( &
3504	10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) &
3505	- 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
3506	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
3507	!
3508	!-- k index has to be modified near bottom and top, else array
3509	!-- subscripts will be exceeded.
3510	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
3511	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
3512	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
3513
3514	k_ppp = k + 3 * ibit8
3515	k_pp = k + 2 * ( 1 - ibit6 )
3516	k_mm = k - 2 * ibit8
3517
3518
3519	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
3520	( 37.0_wp * ibit8 * adv_sca_5 &
3521	+ 7.0_wp * ibit7 * adv_sca_3 &
3522	+ ibit6 * adv_sca_1 &
3523	) * &
3524	( sk(k+1,j,i) + sk(k,j,i) ) &
3525	- ( 8.0_wp * ibit8 * adv_sca_5 &
3526	+ ibit7 * adv_sca_3 &
3527	) * &
3528	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
3529	+ ( ibit8 * adv_sca_5 &
3530	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
3531	)
3532
3533	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
3534	( 10.0_wp * ibit8 * adv_sca_5 &
3535	+ 3.0_wp * ibit7 * adv_sca_3 &
3536	+ ibit6 * adv_sca_1 &
3537	) * &
3538	( sk(k+1,j,i) - sk(k,j,i) ) &
3539	- ( 5.0_wp * ibit8 * adv_sca_5 &
3540	+ ibit7 * adv_sca_3 &
3541	) * &
3542	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
3543	+ ( ibit8 * adv_sca_5 &
3544	) * &
3545	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
3546	)
3547	!
3548	!-- Calculate the divergence of the velocity field. A respective
3549	!-- correction is needed to overcome numerical instabilities introduced
3550	!-- by a not sufficient reduction of divergences near topography.
3551	div = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx &
3552	+ ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy &
3553	+ ( w(k,j,i) * rho_air_zw(k) - &
3554	w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
3555
3556	tend(k,j,i) = tend(k,j,i) - ( &
3557	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
3558	swap_diss_x_local(k,j) ) * ddx &
3559	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
3560	swap_diss_y_local(k) ) * ddy &
3561	+ ( ( flux_t(k) + diss_t(k) ) - &
3562	( flux_d + diss_d ) &
3563	) * drho_air(k) * ddzw(k) &
3564	) + sk(k,j,i) * div
3565
3566	swap_flux_y_local(k) = flux_n(k)
3567	swap_diss_y_local(k) = diss_n(k)
3568	swap_flux_x_local(k,j) = flux_r(k)
3569	swap_diss_x_local(k,j) = diss_r(k)
3570	flux_d = flux_t(k)
3571	diss_d = diss_t(k)
3572
3573	ENDDO
3574	!
3575	!-- Evaluation of statistics.
3576	SELECT CASE ( sk_char )
3577
3578	CASE ( 'pt' )
3579	DO k = nzb, nzt
3580	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
3581	+ ( flux_t(k) &
3582	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3583	* ( w(k,j,i) - hom(k,1,3,0) ) &
3584	+ diss_t(k) &
3585	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3586	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3587	) * weight_substep(intermediate_timestep_count)
3588	ENDDO
3589	CASE ( 'sa' )
3590	DO k = nzb, nzt
3591	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
3592	+ ( flux_t(k) &
3593	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3594	* ( w(k,j,i) - hom(k,1,3,0) ) &
3595	+ diss_t(k) &
3596	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3597	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3598	) * weight_substep(intermediate_timestep_count)
3599	ENDDO
3600	CASE ( 'q' )
3601	DO k = nzb, nzt
3602	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
3603	+ ( flux_t(k) &
3604	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3605	* ( w(k,j,i) - hom(k,1,3,0) ) &
3606	+ diss_t(k) &
3607	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3608	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3609	) * weight_substep(intermediate_timestep_count)
3610	ENDDO
3611	CASE ( 'qc' )
3612	DO k = nzb, nzt
3613	sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) &
3614	+ ( flux_t(k) &
3615	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3616	* ( w(k,j,i) - hom(k,1,3,0) ) &
3617	+ diss_t(k) &
3618	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3619	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3620	) * weight_substep(intermediate_timestep_count)
3621	ENDDO
3622	CASE ( 'qr' )
3623	DO k = nzb, nzt
3624	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) &
3625	+ ( flux_t(k) &
3626	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3627	* ( w(k,j,i) - hom(k,1,3,0) ) &
3628	+ diss_t(k) &
3629	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3630	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3631	) * weight_substep(intermediate_timestep_count)
3632	ENDDO
3633	CASE ( 'nc' )
3634	DO k = nzb, nzt
3635	sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) &
3636	+ ( flux_t(k) &
3637	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3638	* ( w(k,j,i) - hom(k,1,3,0) ) &
3639	+ diss_t(k) &
3640	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3641	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3642	) * weight_substep(intermediate_timestep_count)
3643	ENDDO
3644	CASE ( 'nr' )
3645	DO k = nzb, nzt
3646	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) &
3647	+ ( flux_t(k) &
3648	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3649	* ( w(k,j,i) - hom(k,1,3,0) ) &
3650	+ diss_t(k) &
3651	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3652	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3653	) * weight_substep(intermediate_timestep_count)
3654	ENDDO
3655	CASE ( 's' )
3656	DO k = nzb, nzt
3657	sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) &
3658	+ ( flux_t(k) &
3659	/ ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
3660	* ( w(k,j,i) - hom(k,1,3,0) ) &
3661	+ diss_t(k) &
3662	/ ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
3663	* ABS(w(k,j,i) - hom(k,1,3,0) ) &
3664	) * weight_substep(intermediate_timestep_count)
3665	ENDDO
3666
3667
3668	END SELECT
3669
3670	ENDDO
3671	ENDDO
3672
3673	END SUBROUTINE advec_s_ws
3674
3675
3676	!------------------------------------------------------------------------------!
3677	! Description:
3678	! ------------
3679	!> Advection of u - Call for all grid points
3680	!------------------------------------------------------------------------------!
3681	SUBROUTINE advec_u_ws
3682
3683	USE arrays_3d, &
3684	ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
3685
3686	USE constants, &
3687	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
3688
3689	USE control_parameters, &
3690	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
3691
3692	USE grid_variables, &
3693	ONLY: ddx, ddy
3694
3695	USE indices, &
3696	ONLY: nxl, nxlu, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1
3697
3698	USE kinds
3699
3700	USE statistics, &
3701	ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
3702
3703	IMPLICIT NONE
3704
3705	INTEGER(iwp) :: i !<
3706	INTEGER(iwp) :: ibit9 !<
3707	INTEGER(iwp) :: ibit10 !<
3708	INTEGER(iwp) :: ibit11 !<
3709	INTEGER(iwp) :: ibit12 !<
3710	INTEGER(iwp) :: ibit13 !<
3711	INTEGER(iwp) :: ibit14 !<
3712	INTEGER(iwp) :: ibit15 !<
3713	INTEGER(iwp) :: ibit16 !<
3714	INTEGER(iwp) :: ibit17 !<
3715	INTEGER(iwp) :: j !<
3716	INTEGER(iwp) :: k !<
3717	INTEGER(iwp) :: k_mm !<
3718	INTEGER(iwp) :: k_pp !<
3719	INTEGER(iwp) :: k_ppp !<
3720	INTEGER(iwp) :: tn = 0 !<
3721
3722	REAL(wp) :: diss_d !<
3723	REAL(wp) :: div !<
3724	REAL(wp) :: flux_d !<
3725	REAL(wp) :: gu !<
3726	REAL(wp) :: gv !<
3727	REAL(wp) :: v_comp !<
3728	REAL(wp) :: w_comp !<
3729
3730	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u !<
3731	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u !<
3732
3733	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u !<
3734	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u !<
3735
3736	REAL(wp), DIMENSION(nzb:nzt) :: diss_n !<
3737	REAL(wp), DIMENSION(nzb:nzt) :: diss_r !<
3738	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
3739	REAL(wp), DIMENSION(nzb:nzt) :: flux_n !<
3740	REAL(wp), DIMENSION(nzb:nzt) :: flux_r !<
3741	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
3742	REAL(wp), DIMENSION(nzb:nzt) :: u_comp !<
3743
3744	gu = 2.0_wp * u_gtrans
3745	gv = 2.0_wp * v_gtrans
3746
3747	!
3748	!-- Compute the fluxes for the whole left boundary of the processor domain.
3749	i = nxlu
3750	DO j = nys, nyn
3751	DO k = nzb+1, nzb_max
3752
3753	ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
3754	ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
3755	ibit9 = IBITS(advc_flags_1(k,j,i-1),9,1)
3756
3757	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
3758	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
3759	( 37.0_wp * ibit11 * adv_mom_5 &
3760	+ 7.0_wp * ibit10 * adv_mom_3 &
3761	+ ibit9 * adv_mom_1 &
3762	) * &
3763	( u(k,j,i) + u(k,j,i-1) ) &
3764	- ( 8.0_wp * ibit11 * adv_mom_5 &
3765	+ ibit10 * adv_mom_3 &
3766	) * &
3767	( u(k,j,i+1) + u(k,j,i-2) ) &
3768	+ ( ibit11 * adv_mom_5 &
3769	) * &
3770	( u(k,j,i+2) + u(k,j,i-3) ) &
3771	)
3772
3773	swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( &
3774	( 10.0_wp * ibit11 * adv_mom_5 &
3775	+ 3.0_wp * ibit10 * adv_mom_3 &
3776	+ ibit9 * adv_mom_1 &
3777	) * &
3778	( u(k,j,i) - u(k,j,i-1) ) &
3779	- ( 5.0_wp * ibit11 * adv_mom_5 &
3780	+ ibit10 * adv_mom_3 &
3781	) * &
3782	( u(k,j,i+1) - u(k,j,i-2) ) &
3783	+ ( ibit11 * adv_mom_5 &
3784	) * &
3785	( u(k,j,i+2) - u(k,j,i-3) ) &
3786	)
3787
3788	ENDDO
3789
3790	DO k = nzb_max+1, nzt
3791
3792	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
3793	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
3794	37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) &
3795	- 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) &
3796	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
3797	swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( &
3798	10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) &
3799	- 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) &
3800	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
3801
3802	ENDDO
3803	ENDDO
3804
3805	DO i = nxlu, nxr
3806	!
3807	!-- The following loop computes the fluxes for the south boundary points
3808	j = nys
3809	DO k = nzb+1, nzb_max
3810
3811	ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
3812	ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
3813	ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
3814
3815	v_comp = v(k,j,i) + v(k,j,i-1) - gv
3816	swap_flux_y_local_u(k) = v_comp * ( &
3817	( 37.0_wp * ibit14 * adv_mom_5 &
3818	+ 7.0_wp * ibit13 * adv_mom_3 &
3819	+ ibit12 * adv_mom_1 &
3820	) * &
3821	( u(k,j,i) + u(k,j-1,i) ) &
3822	- ( 8.0_wp * ibit14 * adv_mom_5 &
3823	+ ibit13 * adv_mom_3 &
3824	) * &
3825	( u(k,j+1,i) + u(k,j-2,i) ) &
3826	+ ( ibit14 * adv_mom_5 &
3827	) * &
3828	( u(k,j+2,i) + u(k,j-3,i) ) &
3829	)
3830
3831	swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( &
3832	( 10.0_wp * ibit14 * adv_mom_5 &
3833	+ 3.0_wp * ibit13 * adv_mom_3 &
3834	+ ibit12 * adv_mom_1 &
3835	) * &
3836	( u(k,j,i) - u(k,j-1,i) ) &
3837	- ( 5.0_wp * ibit14 * adv_mom_5 &
3838	+ ibit13 * adv_mom_3 &
3839	) * &
3840	( u(k,j+1,i) - u(k,j-2,i) ) &
3841	+ ( ibit14 * adv_mom_5 &
3842	) * &
3843	( u(k,j+2,i) - u(k,j-3,i) ) &
3844	)
3845
3846	ENDDO
3847
3848	DO k = nzb_max+1, nzt
3849
3850	v_comp = v(k,j,i) + v(k,j,i-1) - gv
3851	swap_flux_y_local_u(k) = v_comp * ( &
3852	37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) &
3853	- 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) &
3854	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
3855	swap_diss_y_local_u(k) = - ABS(v_comp) * ( &
3856	10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) &
3857	- 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) &
3858	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
3859
3860	ENDDO
3861	!
3862	!-- Computation of interior fluxes and tendency terms
3863	DO j = nys, nyn
3864
3865	flux_t(0) = 0.0_wp
3866	diss_t(0) = 0.0_wp
3867	flux_d = 0.0_wp
3868	diss_d = 0.0_wp
3869
3870	DO k = nzb+1, nzb_max
3871
3872	ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
3873	ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
3874	ibit9 = IBITS(advc_flags_1(k,j,i),9,1)
3875
3876	u_comp(k) = u(k,j,i+1) + u(k,j,i)
3877	flux_r(k) = ( u_comp(k) - gu ) * ( &
3878	( 37.0_wp * ibit11 * adv_mom_5 &
3879	+ 7.0_wp * ibit10 * adv_mom_3 &
3880	+ ibit9 * adv_mom_1 &
3881	) * &
3882	( u(k,j,i+1) + u(k,j,i) ) &
3883	- ( 8.0_wp * ibit11 * adv_mom_5 &
3884	+ ibit10 * adv_mom_3 &
3885	) * &
3886	( u(k,j,i+2) + u(k,j,i-1) ) &
3887	+ ( ibit11 * adv_mom_5 &
3888	) * &
3889	( u(k,j,i+3) + u(k,j,i-2) ) &
3890	)
3891
3892	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
3893	( 10.0_wp * ibit11 * adv_mom_5 &
3894	+ 3.0_wp * ibit10 * adv_mom_3 &
3895	+ ibit9 * adv_mom_1 &
3896	) * &
3897	( u(k,j,i+1) - u(k,j,i) ) &
3898	- ( 5.0_wp * ibit11 * adv_mom_5 &
3899	+ ibit10 * adv_mom_3 &
3900	) * &
3901	( u(k,j,i+2) - u(k,j,i-1) ) &
3902	+ ( ibit11 * adv_mom_5 &
3903	) * &
3904	( u(k,j,i+3) - u(k,j,i-2) ) &
3905	)
3906
3907	ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
3908	ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
3909	ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
3910
3911	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
3912	flux_n(k) = v_comp * ( &
3913	( 37.0_wp * ibit14 * adv_mom_5 &
3914	+ 7.0_wp * ibit13 * adv_mom_3 &
3915	+ ibit12 * adv_mom_1 &
3916	) * &
3917	( u(k,j+1,i) + u(k,j,i) ) &
3918	- ( 8.0_wp * ibit14 * adv_mom_5 &
3919	+ ibit13 * adv_mom_3 &
3920	) * &
3921	( u(k,j+2,i) + u(k,j-1,i) ) &
3922	+ ( ibit14 * adv_mom_5 &
3923	) * &
3924	( u(k,j+3,i) + u(k,j-2,i) ) &
3925	)
3926
3927	diss_n(k) = - ABS ( v_comp ) * ( &
3928	( 10.0_wp * ibit14 * adv_mom_5 &
3929	+ 3.0_wp * ibit13 * adv_mom_3 &
3930	+ ibit12 * adv_mom_1 &
3931	) * &
3932	( u(k,j+1,i) - u(k,j,i) ) &
3933	- ( 5.0_wp * ibit14 * adv_mom_5 &
3934	+ ibit13 * adv_mom_3 &
3935	) * &
3936	( u(k,j+2,i) - u(k,j-1,i) ) &
3937	+ ( ibit14 * adv_mom_5 &
3938	) * &
3939	( u(k,j+3,i) - u(k,j-2,i) ) &
3940	)
3941	!
3942	!-- k index has to be modified near bottom and top, else array
3943	!-- subscripts will be exceeded.
3944	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
3945	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
3946	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
3947
3948	k_ppp = k + 3 * ibit17
3949	k_pp = k + 2 * ( 1 - ibit15 )
3950	k_mm = k - 2 * ibit17
3951
3952	w_comp = w(k,j,i) + w(k,j,i-1)
3953	flux_t(k) = w_comp * rho_air_zw(k) * ( &
3954	( 37.0_wp * ibit17 * adv_mom_5 &
3955	+ 7.0_wp * ibit16 * adv_mom_3 &
3956	+ ibit15 * adv_mom_1 &
3957	) * &
3958	( u(k+1,j,i) + u(k,j,i) ) &
3959	- ( 8.0_wp * ibit17 * adv_mom_5 &
3960	+ ibit16 * adv_mom_3 &
3961	) * &
3962	( u(k_pp,j,i) + u(k-1,j,i) ) &
3963	+ ( ibit17 * adv_mom_5 &
3964	) * &
3965	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
3966	)
3967
3968	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
3969	( 10.0_wp * ibit17 * adv_mom_5 &
3970	+ 3.0_wp * ibit16 * adv_mom_3 &
3971	+ ibit15 * adv_mom_1 &
3972	) * &
3973	( u(k+1,j,i) - u(k,j,i) ) &
3974	- ( 5.0_wp * ibit17 * adv_mom_5 &
3975	+ ibit16 * adv_mom_3 &
3976	) * &
3977	( u(k_pp,j,i) - u(k-1,j,i) ) &
3978	+ ( ibit17 * adv_mom_5 &
3979	) * &
3980	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
3981	)
3982	!
3983	!-- Calculate the divergence of the velocity field. A respective
3984	!-- correction is needed to overcome numerical instabilities caused
3985	!-- by a not sufficient reduction of divergences near topography.
3986	div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) &
3987	- ( u(k,j,i) + u(k,j,i-1) ) &
3988	* ( IBITS(advc_flags_1(k,j,i-1),9,1) &
3989	+ IBITS(advc_flags_1(k,j,i-1),10,1) &
3990	+ IBITS(advc_flags_1(k,j,i-1),11,1) &
3991	) &
3992	) * rho_air(k) * ddx &
3993	+ ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) &
3994	- ( v(k,j,i) + v(k,j,i-1 ) ) &
3995	* ( IBITS(advc_flags_1(k,j-1,i),12,1) &
3996	+ IBITS(advc_flags_1(k,j-1,i),13,1) &
3997	+ IBITS(advc_flags_1(k,j-1,i),14,1) &
3998	) &
3999	) * rho_air(k) * ddy &
4000	+ ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) &
4001	- ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
4002	* ( IBITS(advc_flags_1(k-1,j,i),15,1) &
4003	+ IBITS(advc_flags_1(k-1,j,i),16,1) &
4004	+ IBITS(advc_flags_1(k-1,j,i),17,1) &
4005	) &
4006	) * ddzw(k) &
4007	) * 0.5_wp
4008
4009
4010
4011	tend(k,j,i) = tend(k,j,i) - ( &
4012	( flux_r(k) + diss_r(k) &
4013	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
4014	+ ( flux_n(k) + diss_n(k) &
4015	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
4016	+ ( ( flux_t(k) + diss_t(k) ) &
4017	- ( flux_d + diss_d ) &
4018	) * drho_air(k) * ddzw(k) &
4019	) + div * u(k,j,i)
4020
4021	swap_flux_x_local_u(k,j) = flux_r(k)
4022	swap_diss_x_local_u(k,j) = diss_r(k)
4023	swap_flux_y_local_u(k) = flux_n(k)
4024	swap_diss_y_local_u(k) = diss_n(k)
4025	flux_d = flux_t(k)
4026	diss_d = diss_t(k)
4027	!
4028	!-- Statistical Evaluation of u'u'. The factor has to be applied
4029	!-- for right evaluation when gallilei_trans = .T. .
4030	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
4031	+ ( flux_r(k) &
4032	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4033	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
4034	+ diss_r(k) &
4035	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4036	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
4037	) * weight_substep(intermediate_timestep_count)
4038	!
4039	!-- Statistical Evaluation of w'u'.
4040	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
4041	+ ( flux_t(k) &
4042	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4043	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4044	+ diss_t(k) &
4045	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4046	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4047	) * weight_substep(intermediate_timestep_count)
4048
4049	ENDDO
4050
4051	DO k = nzb_max+1, nzt
4052
4053	u_comp(k) = u(k,j,i+1) + u(k,j,i)
4054	flux_r(k) = ( u_comp(k) - gu ) * ( &
4055	37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &
4056	- 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &
4057	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
4058	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
4059	10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &
4060	- 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &
4061	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
4062
4063	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
4064	flux_n(k) = v_comp * ( &
4065	37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &
4066	- 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &
4067	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
4068	diss_n(k) = - ABS( v_comp ) * ( &
4069	10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &
4070	- 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &
4071	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
4072	!
4073	!-- k index has to be modified near bottom and top, else array
4074	!-- subscripts will be exceeded.
4075	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
4076	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
4077	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
4078
4079	k_ppp = k + 3 * ibit17
4080	k_pp = k + 2 * ( 1 - ibit15 )
4081	k_mm = k - 2 * ibit17
4082
4083	w_comp = w(k,j,i) + w(k,j,i-1)
4084	flux_t(k) = w_comp * rho_air_zw(k) * ( &
4085	( 37.0_wp * ibit17 * adv_mom_5 &
4086	+ 7.0_wp * ibit16 * adv_mom_3 &
4087	+ ibit15 * adv_mom_1 &
4088	) * &
4089	( u(k+1,j,i) + u(k,j,i) ) &
4090	- ( 8.0_wp * ibit17 * adv_mom_5 &
4091	+ ibit16 * adv_mom_3 &
4092	) * &
4093	( u(k_pp,j,i) + u(k-1,j,i) ) &
4094	+ ( ibit17 * adv_mom_5 &
4095	) * &
4096	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
4097	)
4098
4099	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
4100	( 10.0_wp * ibit17 * adv_mom_5 &
4101	+ 3.0_wp * ibit16 * adv_mom_3 &
4102	+ ibit15 * adv_mom_1 &
4103	) * &
4104	( u(k+1,j,i) - u(k,j,i) ) &
4105	- ( 5.0_wp * ibit17 * adv_mom_5 &
4106	+ ibit16 * adv_mom_3 &
4107	) * &
4108	( u(k_pp,j,i) - u(k-1,j,i) ) &
4109	+ ( ibit17 * adv_mom_5 &
4110	) * &
4111	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
4112	)
4113	!
4114	!-- Calculate the divergence of the velocity field. A respective
4115	!-- correction is needed to overcome numerical instabilities caused
4116	!-- by a not sufficient reduction of divergences near topography.
4117	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
4118	* rho_air(k)&
4119	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
4120	* rho_air(k)&
4121	+ ( w_comp * rho_air_zw(k) - &
4122	( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) &
4123	) * ddzw(k) &
4124	) * 0.5_wp
4125
4126	tend(k,j,i) = tend(k,j,i) - ( &
4127	( flux_r(k) + diss_r(k) &
4128	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
4129	+ ( flux_n(k) + diss_n(k) &
4130	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
4131	+ ( ( flux_t(k) + diss_t(k) ) &
4132	- ( flux_d + diss_d ) &
4133	) * drho_air(k) * ddzw(k) &
4134	) + div * u(k,j,i)
4135
4136	swap_flux_x_local_u(k,j) = flux_r(k)
4137	swap_diss_x_local_u(k,j) = diss_r(k)
4138	swap_flux_y_local_u(k) = flux_n(k)
4139	swap_diss_y_local_u(k) = diss_n(k)
4140	flux_d = flux_t(k)
4141	diss_d = diss_t(k)
4142	!
4143	!-- Statistical Evaluation of u'u'. The factor has to be applied
4144	!-- for right evaluation when gallilei_trans = .T. .
4145	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
4146	+ ( flux_r(k) &
4147	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4148	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
4149	+ diss_r(k) &
4150	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
4151	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
4152	) * weight_substep(intermediate_timestep_count)
4153	!
4154	!-- Statistical Evaluation of w'u'.
4155	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
4156	+ ( flux_t(k) &
4157	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4158	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4159	+ diss_t(k) &
4160	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4161	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4162	) * weight_substep(intermediate_timestep_count)
4163	ENDDO
4164	ENDDO
4165	ENDDO
4166	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
4167
4168
4169	END SUBROUTINE advec_u_ws
4170
4171
4172	!------------------------------------------------------------------------------!
4173	! Description:
4174	! ------------
4175	!> Advection of v - Call for all grid points
4176	!------------------------------------------------------------------------------!
4177	SUBROUTINE advec_v_ws
4178
4179	USE arrays_3d, &
4180	ONLY: ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
4181
4182	USE constants, &
4183	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
4184
4185	USE control_parameters, &
4186	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
4187
4188	USE grid_variables, &
4189	ONLY: ddx, ddy
4190
4191	USE indices, &
4192	ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_max, nzt, advc_flags_1
4193
4194	USE kinds
4195
4196	USE statistics, &
4197	ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep
4198
4199	IMPLICIT NONE
4200
4201
4202	INTEGER(iwp) :: i !<
4203	INTEGER(iwp) :: ibit18 !<
4204	INTEGER(iwp) :: ibit19 !<
4205	INTEGER(iwp) :: ibit20 !<
4206	INTEGER(iwp) :: ibit21 !<
4207	INTEGER(iwp) :: ibit22 !<
4208	INTEGER(iwp) :: ibit23 !<
4209	INTEGER(iwp) :: ibit24 !<
4210	INTEGER(iwp) :: ibit25 !<
4211	INTEGER(iwp) :: ibit26 !<
4212	INTEGER(iwp) :: j !<
4213	INTEGER(iwp) :: k !<
4214	INTEGER(iwp) :: k_mm !<
4215	INTEGER(iwp) :: k_pp !<
4216	INTEGER(iwp) :: k_ppp !<
4217	INTEGER(iwp) :: tn = 0 !<
4218
4219	REAL(wp) :: diss_d !<
4220	REAL(wp) :: div !<
4221	REAL(wp) :: flux_d !<
4222	REAL(wp) :: gu !<
4223	REAL(wp) :: gv !<
4224	REAL(wp) :: u_comp !<
4225	REAL(wp) :: w_comp !<
4226
4227	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v !<
4228	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v !<
4229
4230	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v !<
4231	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v !<
4232
4233	REAL(wp), DIMENSION(nzb:nzt) :: diss_n !<
4234	REAL(wp), DIMENSION(nzb:nzt) :: diss_r !<
4235	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
4236	REAL(wp), DIMENSION(nzb:nzt) :: flux_n !<
4237	REAL(wp), DIMENSION(nzb:nzt) :: flux_r !<
4238	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
4239	REAL(wp), DIMENSION(nzb:nzt) :: v_comp !<
4240
4241	gu = 2.0_wp * u_gtrans
4242	gv = 2.0_wp * v_gtrans
4243	!
4244	!-- First compute the whole left boundary of the processor domain
4245	i = nxl
4246	DO j = nysv, nyn
4247	DO k = nzb+1, nzb_max
4248
4249	ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
4250	ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
4251	ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
4252
4253	u_comp = u(k,j-1,i) + u(k,j,i) - gu
4254	swap_flux_x_local_v(k,j) = u_comp * ( &
4255	( 37.0_wp * ibit20 * adv_mom_5 &
4256	+ 7.0_wp * ibit19 * adv_mom_3 &
4257	+ ibit18 * adv_mom_1 &
4258	) * &
4259	( v(k,j,i) + v(k,j,i-1) ) &
4260	- ( 8.0_wp * ibit20 * adv_mom_5 &
4261	+ ibit19 * adv_mom_3 &
4262	) * &
4263	( v(k,j,i+1) + v(k,j,i-2) ) &
4264	+ ( ibit20 * adv_mom_5 &
4265	) * &
4266	( v(k,j,i+2) + v(k,j,i-3) ) &
4267	)
4268
4269	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
4270	( 10.0_wp * ibit20 * adv_mom_5 &
4271	+ 3.0_wp * ibit19 * adv_mom_3 &
4272	+ ibit18 * adv_mom_1 &
4273	) * &
4274	( v(k,j,i) - v(k,j,i-1) ) &
4275	- ( 5.0_wp * ibit20 * adv_mom_5 &
4276	+ ibit19 * adv_mom_3 &
4277	) * &
4278	( v(k,j,i+1) - v(k,j,i-2) ) &
4279	+ ( ibit20 * adv_mom_5 &
4280	) * &
4281	( v(k,j,i+2) - v(k,j,i-3) ) &
4282	)
4283
4284	ENDDO
4285
4286	DO k = nzb_max+1, nzt
4287
4288	u_comp = u(k,j-1,i) + u(k,j,i) - gu
4289	swap_flux_x_local_v(k,j) = u_comp * ( &
4290	37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) &
4291	- 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) &
4292	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
4293	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
4294	10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) &
4295	- 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) &
4296	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
4297
4298	ENDDO
4299
4300	ENDDO
4301
4302	DO i = nxl, nxr
4303
4304	j = nysv
4305	DO k = nzb+1, nzb_max
4306
4307	ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
4308	ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
4309	ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
4310
4311	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
4312	swap_flux_y_local_v(k) = v_comp(k) * ( &
4313	( 37.0_wp * ibit23 * adv_mom_5 &
4314	+ 7.0_wp * ibit22 * adv_mom_3 &
4315	+ ibit21 * adv_mom_1 &
4316	) * &
4317	( v(k,j,i) + v(k,j-1,i) ) &
4318	- ( 8.0_wp * ibit23 * adv_mom_5 &
4319	+ ibit22 * adv_mom_3 &
4320	) * &
4321	( v(k,j+1,i) + v(k,j-2,i) ) &
4322	+ ( ibit23 * adv_mom_5 &
4323	) * &
4324	( v(k,j+2,i) + v(k,j-3,i) ) &
4325	)
4326
4327	swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( &
4328	( 10.0_wp * ibit23 * adv_mom_5 &
4329	+ 3.0_wp * ibit22 * adv_mom_3 &
4330	+ ibit21 * adv_mom_1 &
4331	) * &
4332	( v(k,j,i) - v(k,j-1,i) ) &
4333	- ( 5.0_wp * ibit23 * adv_mom_5 &
4334	+ ibit22 * adv_mom_3 &
4335	) * &
4336	( v(k,j+1,i) - v(k,j-2,i) ) &
4337	+ ( ibit23 * adv_mom_5 &
4338	) * &
4339	( v(k,j+2,i) - v(k,j-3,i) ) &
4340	)
4341
4342	ENDDO
4343
4344	DO k = nzb_max+1, nzt
4345
4346	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
4347	swap_flux_y_local_v(k) = v_comp(k) * ( &
4348	37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) &
4349	- 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) &
4350	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
4351	swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( &
4352	10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) &
4353	- 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) &
4354	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
4355
4356	ENDDO
4357
4358	DO j = nysv, nyn
4359
4360	flux_t(0) = 0.0_wp
4361	diss_t(0) = 0.0_wp
4362	flux_d = 0.0_wp
4363	diss_d = 0.0_wp
4364
4365	DO k = nzb+1, nzb_max
4366
4367	ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
4368	ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
4369	ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
4370
4371	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
4372	flux_r(k) = u_comp * ( &
4373	( 37.0_wp * ibit20 * adv_mom_5 &
4374	+ 7.0_wp * ibit19 * adv_mom_3 &
4375	+ ibit18 * adv_mom_1 &
4376	) * &
4377	( v(k,j,i+1) + v(k,j,i) ) &
4378	- ( 8.0_wp * ibit20 * adv_mom_5 &
4379	+ ibit19 * adv_mom_3 &
4380	) * &
4381	( v(k,j,i+2) + v(k,j,i-1) ) &
4382	+ ( ibit20 * adv_mom_5 &
4383	) * &
4384	( v(k,j,i+3) + v(k,j,i-2) ) &
4385	)
4386
4387	diss_r(k) = - ABS( u_comp ) * ( &
4388	( 10.0_wp * ibit20 * adv_mom_5 &
4389	+ 3.0_wp * ibit19 * adv_mom_3 &
4390	+ ibit18 * adv_mom_1 &
4391	) * &
4392	( v(k,j,i+1) - v(k,j,i) ) &
4393	- ( 5.0_wp * ibit20 * adv_mom_5 &
4394	+ ibit19 * adv_mom_3 &
4395	) * &
4396	( v(k,j,i+2) - v(k,j,i-1) ) &
4397	+ ( ibit20 * adv_mom_5 &
4398	) * &
4399	( v(k,j,i+3) - v(k,j,i-2) ) &
4400	)
4401
4402	ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
4403	ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
4404	ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
4405
4406	v_comp(k) = v(k,j+1,i) + v(k,j,i)
4407	flux_n(k) = ( v_comp(k) - gv ) * ( &
4408	( 37.0_wp * ibit23 * adv_mom_5 &
4409	+ 7.0_wp * ibit22 * adv_mom_3 &
4410	+ ibit21 * adv_mom_1 &
4411	) * &
4412	( v(k,j+1,i) + v(k,j,i) ) &
4413	- ( 8.0_wp * ibit23 * adv_mom_5 &
4414	+ ibit22 * adv_mom_3 &
4415	) * &
4416	( v(k,j+2,i) + v(k,j-1,i) ) &
4417	+ ( ibit23 * adv_mom_5 &
4418	) * &
4419	( v(k,j+3,i) + v(k,j-2,i) ) &
4420	)
4421
4422	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
4423	( 10.0_wp * ibit23 * adv_mom_5 &
4424	+ 3.0_wp * ibit22 * adv_mom_3 &
4425	+ ibit21 * adv_mom_1 &
4426	) * &
4427	( v(k,j+1,i) - v(k,j,i) ) &
4428	- ( 5.0_wp * ibit23 * adv_mom_5 &
4429	+ ibit22 * adv_mom_3 &
4430	) * &
4431	( v(k,j+2,i) - v(k,j-1,i) ) &
4432	+ ( ibit23 * adv_mom_5 &
4433	) * &
4434	( v(k,j+3,i) - v(k,j-2,i) ) &
4435	)
4436	!
4437	!-- k index has to be modified near bottom and top, else array
4438	!-- subscripts will be exceeded.
4439	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
4440	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
4441	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
4442
4443	k_ppp = k + 3 * ibit26
4444	k_pp = k + 2 * ( 1 - ibit24 )
4445	k_mm = k - 2 * ibit26
4446
4447	w_comp = w(k,j-1,i) + w(k,j,i)
4448	flux_t(k) = w_comp * rho_air_zw(k) * ( &
4449	( 37.0_wp * ibit26 * adv_mom_5 &
4450	+ 7.0_wp * ibit25 * adv_mom_3 &
4451	+ ibit24 * adv_mom_1 &
4452	) * &
4453	( v(k+1,j,i) + v(k,j,i) ) &
4454	- ( 8.0_wp * ibit26 * adv_mom_5 &
4455	+ ibit25 * adv_mom_3 &
4456	) * &
4457	( v(k_pp,j,i) + v(k-1,j,i) ) &
4458	+ ( ibit26 * adv_mom_5 &
4459	) * &
4460	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
4461	)
4462
4463	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
4464	( 10.0_wp * ibit26 * adv_mom_5 &
4465	+ 3.0_wp * ibit25 * adv_mom_3 &
4466	+ ibit24 * adv_mom_1 &
4467	) * &
4468	( v(k+1,j,i) - v(k,j,i) ) &
4469	- ( 5.0_wp * ibit26 * adv_mom_5 &
4470	+ ibit25 * adv_mom_3 &
4471	) * &
4472	( v(k_pp,j,i) - v(k-1,j,i) ) &
4473	+ ( ibit26 * adv_mom_5 &
4474	) * &
4475	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
4476	)
4477	!
4478	!-- Calculate the divergence of the velocity field. A respective
4479	!-- correction is needed to overcome numerical instabilities caused
4480	!-- by a not sufficient reduction of divergences near topography.
4481	div = ( ( ( u_comp + gu ) &
4482	* ( ibit18 + ibit19 + ibit20 ) &
4483	- ( u(k,j-1,i) + u(k,j,i) ) &
4484	* ( IBITS(advc_flags_1(k,j,i-1),18,1) &
4485	+ IBITS(advc_flags_1(k,j,i-1),19,1) &
4486	+ IBITS(advc_flags_1(k,j,i-1),20,1) &
4487	) &
4488	) * rho_air(k) * ddx &
4489	+ ( v_comp(k) &
4490	* ( ibit21 + ibit22 + ibit23 ) &
4491	- ( v(k,j,i) + v(k,j-1,i) ) &
4492	* ( IBITS(advc_flags_1(k,j-1,i),21,1) &
4493	+ IBITS(advc_flags_1(k,j-1,i),22,1) &
4494	+ IBITS(advc_flags_1(k,j-1,i),23,1) &
4495	) &
4496	) * rho_air(k) * ddy &
4497	+ ( w_comp * rho_air_zw(k) &
4498	* ( ibit24 + ibit25 + ibit26 ) &
4499	- ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
4500	* ( IBITS(advc_flags_1(k-1,j,i),24,1) &
4501	+ IBITS(advc_flags_1(k-1,j,i),25,1) &
4502	+ IBITS(advc_flags_1(k-1,j,i),26,1) &
4503	) &
4504	) * ddzw(k) &
4505	) * 0.5_wp
4506
4507
4508	tend(k,j,i) = tend(k,j,i) - ( &
4509	( flux_r(k) + diss_r(k) &
4510	- swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) &
4511	) * ddx &
4512	+ ( flux_n(k) + diss_n(k) &
4513	- swap_flux_y_local_v(k) - swap_diss_y_local_v(k) &
4514	) * ddy &
4515	+ ( ( flux_t(k) + diss_t(k) ) &
4516	- ( flux_d + diss_d ) &
4517	) * drho_air(k) * ddzw(k) &
4518	) + v(k,j,i) * div
4519
4520	swap_flux_x_local_v(k,j) = flux_r(k)
4521	swap_diss_x_local_v(k,j) = diss_r(k)
4522	swap_flux_y_local_v(k) = flux_n(k)
4523	swap_diss_y_local_v(k) = diss_n(k)
4524	flux_d = flux_t(k)
4525	diss_d = diss_t(k)
4526
4527	!
4528	!-- Statistical Evaluation of v'v'. The factor has to be applied
4529	!-- for right evaluation when gallilei_trans = .T. .
4530	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
4531	+ ( flux_n(k) &
4532	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4533	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
4534	+ diss_n(k) &
4535	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4536	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
4537	) * weight_substep(intermediate_timestep_count)
4538	!
4539	!-- Statistical Evaluation of w'u'.
4540	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
4541	+ ( flux_t(k) &
4542	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4543	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4544	+ diss_t(k) &
4545	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4546	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4547	) * weight_substep(intermediate_timestep_count)
4548
4549	ENDDO
4550
4551	DO k = nzb_max+1, nzt
4552
4553	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
4554	flux_r(k) = u_comp * ( &
4555	37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &
4556	- 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &
4557	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
4558
4559	diss_r(k) = - ABS( u_comp ) * ( &
4560	10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &
4561	- 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &
4562	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
4563
4564
4565	v_comp(k) = v(k,j+1,i) + v(k,j,i)
4566	flux_n(k) = ( v_comp(k) - gv ) * ( &
4567	37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &
4568	- 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &
4569	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
4570
4571	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
4572	10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &
4573	- 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &
4574	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
4575	!
4576	!-- k index has to be modified near bottom and top, else array
4577	!-- subscripts will be exceeded.
4578	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
4579	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
4580	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
4581
4582	k_ppp = k + 3 * ibit26
4583	k_pp = k + 2 * ( 1 - ibit24 )
4584	k_mm = k - 2 * ibit26
4585
4586	w_comp = w(k,j-1,i) + w(k,j,i)
4587	flux_t(k) = w_comp * rho_air_zw(k) * ( &
4588	( 37.0_wp * ibit26 * adv_mom_5 &
4589	+ 7.0_wp * ibit25 * adv_mom_3 &
4590	+ ibit24 * adv_mom_1 &
4591	) * &
4592	( v(k+1,j,i) + v(k,j,i) ) &
4593	- ( 8.0_wp * ibit26 * adv_mom_5 &
4594	+ ibit25 * adv_mom_3 &
4595	) * &
4596	( v(k_pp,j,i) + v(k-1,j,i) ) &
4597	+ ( ibit26 * adv_mom_5 &
4598	) * &
4599	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
4600	)
4601
4602	diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( &
4603	( 10.0_wp * ibit26 * adv_mom_5 &
4604	+ 3.0_wp * ibit25 * adv_mom_3 &
4605	+ ibit24 * adv_mom_1 &
4606	) * &
4607	( v(k+1,j,i) - v(k,j,i) ) &
4608	- ( 5.0_wp * ibit26 * adv_mom_5 &
4609	+ ibit25 * adv_mom_3 &
4610	) * &
4611	( v(k_pp,j,i) - v(k-1,j,i) ) &
4612	+ ( ibit26 * adv_mom_5 &
4613	) * &
4614	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
4615	)
4616	!
4617	!-- Calculate the divergence of the velocity field. A respective
4618	!-- correction is needed to overcome numerical instabilities caused
4619	!-- by a not sufficient reduction of divergences near topography.
4620	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
4621	* rho_air(k)&
4622	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
4623	* rho_air(k)&
4624	+ ( w_comp * rho_air_zw(k) - &
4625	( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) &
4626	) * ddzw(k) &
4627	) * 0.5_wp
4628
4629	tend(k,j,i) = tend(k,j,i) - ( &
4630	( flux_r(k) + diss_r(k) &
4631	- swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) &
4632	) * ddx &
4633	+ ( flux_n(k) + diss_n(k) &
4634	- swap_flux_y_local_v(k) - swap_diss_y_local_v(k) &
4635	) * ddy &
4636	+ ( ( flux_t(k) + diss_t(k) ) &
4637	- ( flux_d + diss_d ) &
4638	) * drho_air(k) * ddzw(k) &
4639	) + v(k,j,i) * div
4640
4641	swap_flux_x_local_v(k,j) = flux_r(k)
4642	swap_diss_x_local_v(k,j) = diss_r(k)
4643	swap_flux_y_local_v(k) = flux_n(k)
4644	swap_diss_y_local_v(k) = diss_n(k)
4645	flux_d = flux_t(k)
4646	diss_d = diss_t(k)
4647
4648	!
4649	!-- Statistical Evaluation of v'v'. The factor has to be applied
4650	!-- for right evaluation when gallilei_trans = .T. .
4651	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
4652	+ ( flux_n(k) &
4653	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4654	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
4655	+ diss_n(k) &
4656	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
4657	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
4658	) * weight_substep(intermediate_timestep_count)
4659	!
4660	!-- Statistical Evaluation of w'u'.
4661	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
4662	+ ( flux_t(k) &
4663	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4664	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
4665	+ diss_t(k) &
4666	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
4667	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
4668	) * weight_substep(intermediate_timestep_count)
4669
4670	ENDDO
4671	ENDDO
4672	ENDDO
4673	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
4674
4675
4676	END SUBROUTINE advec_v_ws
4677
4678
4679	!------------------------------------------------------------------------------!
4680	! Description:
4681	! ------------
4682	!> Advection of w - Call for all grid points
4683	!------------------------------------------------------------------------------!
4684	SUBROUTINE advec_w_ws
4685
4686	USE arrays_3d, &
4687	ONLY: ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
4688
4689	USE constants, &
4690	ONLY: adv_mom_1, adv_mom_3, adv_mom_5
4691
4692	USE control_parameters, &
4693	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
4694
4695	USE grid_variables, &
4696	ONLY: ddx, ddy
4697
4698	USE indices, &
4699	ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, advc_flags_1, &
4700	advc_flags_2
4701
4702	USE kinds
4703
4704	USE statistics, &
4705	ONLY: hom, sums_ws2_ws_l, weight_substep
4706
4707	IMPLICIT NONE
4708
4709	INTEGER(iwp) :: i !<
4710	INTEGER(iwp) :: ibit27 !<
4711	INTEGER(iwp) :: ibit28 !<
4712	INTEGER(iwp) :: ibit29 !<
4713	INTEGER(iwp) :: ibit30 !<
4714	INTEGER(iwp) :: ibit31 !<
4715	INTEGER(iwp) :: ibit32 !<
4716	INTEGER(iwp) :: ibit33 !<
4717	INTEGER(iwp) :: ibit34 !<
4718	INTEGER(iwp) :: ibit35 !<
4719	INTEGER(iwp) :: j !<
4720	INTEGER(iwp) :: k !<
4721	INTEGER(iwp) :: k_mm !<
4722	INTEGER(iwp) :: k_pp !<
4723	INTEGER(iwp) :: k_ppp !<
4724	INTEGER(iwp) :: tn = 0 !<
4725
4726	REAL(wp) :: diss_d !<
4727	REAL(wp) :: div !<
4728	REAL(wp) :: flux_d !<
4729	REAL(wp) :: gu !<
4730	REAL(wp) :: gv !<
4731	REAL(wp) :: u_comp !<
4732	REAL(wp) :: v_comp !<
4733	REAL(wp) :: w_comp !<
4734
4735	REAL(wp), DIMENSION(nzb:nzt) :: diss_t !<
4736	REAL(wp), DIMENSION(nzb:nzt) :: flux_t !<
4737
4738	REAL(wp), DIMENSION(nzb+1:nzt) :: diss_n !<
4739	REAL(wp), DIMENSION(nzb+1:nzt) :: diss_r !<
4740	REAL(wp), DIMENSION(nzb+1:nzt) :: flux_n !<
4741	REAL(wp), DIMENSION(nzb+1:nzt) :: flux_r !<
4742	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_w !<
4743	REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_w !<
4744
4745	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_w !<
4746	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_w !<
4747
4748	gu = 2.0_wp * u_gtrans
4749	gv = 2.0_wp * v_gtrans
4750	!
4751	!-- compute the whole left boundary of the processor domain
4752	i = nxl
4753	DO j = nys, nyn
4754	DO k = nzb+1, nzb_max
4755
4756	ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
4757	ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
4758	ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
4759
4760	u_comp = u(k+1,j,i) + u(k,j,i) - gu
4761	swap_flux_x_local_w(k,j) = u_comp * ( &
4762	( 37.0_wp * ibit29 * adv_mom_5 &
4763	+ 7.0_wp * ibit28 * adv_mom_3 &
4764	+ ibit27 * adv_mom_1 &
4765	) * &
4766	( w(k,j,i) + w(k,j,i-1) ) &
4767	- ( 8.0_wp * ibit29 * adv_mom_5 &
4768	+ ibit28 * adv_mom_3 &
4769	) * &
4770	( w(k,j,i+1) + w(k,j,i-2) ) &
4771	+ ( ibit29 * adv_mom_5 &
4772	) * &
4773	( w(k,j,i+2) + w(k,j,i-3) ) &
4774	)
4775
4776	swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( &
4777	( 10.0_wp * ibit29 * adv_mom_5 &
4778	+ 3.0_wp * ibit28 * adv_mom_3 &
4779	+ ibit27 * adv_mom_1 &
4780	) * &
4781	( w(k,j,i) - w(k,j,i-1) ) &
4782	- ( 5.0_wp * ibit29 * adv_mom_5 &
4783	+ ibit28 * adv_mom_3 &
4784	) * &
4785	( w(k,j,i+1) - w(k,j,i-2) ) &
4786	+ ( ibit29 * adv_mom_5 &
4787	) * &
4788	( w(k,j,i+2) - w(k,j,i-3) ) &
4789	)
4790
4791	ENDDO
4792
4793	DO k = nzb_max+1, nzt
4794
4795	u_comp = u(k+1,j,i) + u(k,j,i) - gu
4796	swap_flux_x_local_w(k,j) = u_comp * ( &
4797	37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) &
4798	- 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) &
4799	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
4800	swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( &
4801	10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) &
4802	- 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) &
4803	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
4804
4805	ENDDO
4806
4807	ENDDO
4808
4809	DO i = nxl, nxr
4810
4811	j = nys
4812	DO k = nzb+1, nzb_max
4813
4814	ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
4815	ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
4816	ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
4817
4818	v_comp = v(k+1,j,i) + v(k,j,i) - gv
4819	swap_flux_y_local_w(k) = v_comp * ( &
4820	( 37.0_wp * ibit32 * adv_mom_5 &
4821	+ 7.0_wp * ibit31 * adv_mom_3 &
4822	+ ibit30 * adv_mom_1 &
4823	) * &
4824	( w(k,j,i) + w(k,j-1,i) ) &
4825	- ( 8.0_wp * ibit32 * adv_mom_5 &
4826	+ ibit31 * adv_mom_3 &
4827	) * &
4828	( w(k,j+1,i) + w(k,j-2,i) ) &
4829	+ ( ibit32 * adv_mom_5 &
4830	) * &
4831	( w(k,j+2,i) + w(k,j-3,i) ) &
4832	)
4833
4834	swap_diss_y_local_w(k) = - ABS( v_comp ) * ( &
4835	( 10.0_wp * ibit32 * adv_mom_5 &
4836	+ 3.0_wp * ibit31 * adv_mom_3 &
4837	+ ibit30 * adv_mom_1 &
4838	) * &
4839	( w(k,j,i) - w(k,j-1,i) ) &
4840	- ( 5.0_wp * ibit32 * adv_mom_5 &
4841	+ ibit31 * adv_mom_3 &
4842	) * &
4843	( w(k,j+1,i) - w(k,j-2,i) ) &
4844	+ ( ibit32 * adv_mom_5 &
4845	) * &
4846	( w(k,j+2,i) - w(k,j-3,i) ) &
4847	)
4848
4849	ENDDO
4850
4851	DO k = nzb_max+1, nzt
4852
4853	v_comp = v(k+1,j,i) + v(k,j,i) - gv
4854	swap_flux_y_local_w(k) = v_comp * ( &
4855	37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) &
4856	- 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) &
4857	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
4858	swap_diss_y_local_w(k) = - ABS( v_comp ) * ( &
4859	10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) &
4860	- 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) &
4861	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
4862
4863	ENDDO
4864
4865	DO j = nys, nyn
4866
4867	!
4868	!-- The lower flux has to be calculated explicetely for the tendency
4869	!-- at the first w-level. For topography wall this is done implicitely
4870	!-- by advc_flags_1.
4871	k = nzb + 1
4872	w_comp = w(k,j,i) + w(k-1,j,i)
4873	flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
4874	diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
4875	flux_d = flux_t(0)
4876	diss_d = diss_t(0)
4877
4878	DO k = nzb+1, nzb_max
4879
4880	ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
4881	ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
4882	ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
4883
4884	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
4885	flux_r(k) = u_comp * ( &
4886	( 37.0_wp * ibit29 * adv_mom_5 &
4887	+ 7.0_wp * ibit28 * adv_mom_3 &
4888	+ ibit27 * adv_mom_1 &
4889	) * &
4890	( w(k,j,i+1) + w(k,j,i) ) &
4891	- ( 8.0_wp * ibit29 * adv_mom_5 &
4892	+ ibit28 * adv_mom_3 &
4893	) * &
4894	( w(k,j,i+2) + w(k,j,i-1) ) &
4895	+ ( ibit29 * adv_mom_5 &
4896	) * &
4897	( w(k,j,i+3) + w(k,j,i-2) ) &
4898	)
4899
4900	diss_r(k) = - ABS( u_comp ) * ( &
4901	( 10.0_wp * ibit29 * adv_mom_5 &
4902	+ 3.0_wp * ibit28 * adv_mom_3 &
4903	+ ibit27 * adv_mom_1 &
4904	) * &
4905	( w(k,j,i+1) - w(k,j,i) ) &
4906	- ( 5.0_wp * ibit29 * adv_mom_5 &
4907	+ ibit28 * adv_mom_3 &
4908	) * &
4909	( w(k,j,i+2) - w(k,j,i-1) ) &
4910	+ ( ibit29 * adv_mom_5 &
4911	) * &
4912	( w(k,j,i+3) - w(k,j,i-2) ) &
4913	)
4914
4915	ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
4916	ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
4917	ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
4918
4919	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
4920	flux_n(k) = v_comp * ( &
4921	( 37.0_wp * ibit32 * adv_mom_5 &
4922	+ 7.0_wp * ibit31 * adv_mom_3 &
4923	+ ibit30 * adv_mom_1 &
4924	) * &
4925	( w(k,j+1,i) + w(k,j,i) ) &
4926	- ( 8.0_wp * ibit32 * adv_mom_5 &
4927	+ ibit31 * adv_mom_3 &
4928	) * &
4929	( w(k,j+2,i) + w(k,j-1,i) ) &
4930	+ ( ibit32 * adv_mom_5 &
4931	) * &
4932	( w(k,j+3,i) + w(k,j-2,i) ) &
4933	)
4934
4935	diss_n(k) = - ABS( v_comp ) * ( &
4936	( 10.0_wp * ibit32 * adv_mom_5 &
4937	+ 3.0_wp * ibit31 * adv_mom_3 &
4938	+ ibit30 * adv_mom_1 &
4939	) * &
4940	( w(k,j+1,i) - w(k,j,i) ) &
4941	- ( 5.0_wp * ibit32 * adv_mom_5 &
4942	+ ibit31 * adv_mom_3 &
4943	) * &
4944	( w(k,j+2,i) - w(k,j-1,i) ) &
4945	+ ( ibit32 * adv_mom_5 &
4946	) * &
4947	( w(k,j+3,i) - w(k,j-2,i) ) &
4948	)
4949	!
4950	!-- k index has to be modified near bottom and top, else array
4951	!-- subscripts will be exceeded.
4952	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
4953	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
4954	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
4955
4956	k_ppp = k + 3 * ibit35
4957	k_pp = k + 2 * ( 1 - ibit33 )
4958	k_mm = k - 2 * ibit35
4959
4960	w_comp = w(k+1,j,i) + w(k,j,i)
4961	flux_t(k) = w_comp * rho_air(k+1) * ( &
4962	( 37.0_wp * ibit35 * adv_mom_5 &
4963	+ 7.0_wp * ibit34 * adv_mom_3 &
4964	+ ibit33 * adv_mom_1 &
4965	) * &
4966	( w(k+1,j,i) + w(k,j,i) ) &
4967	- ( 8.0_wp * ibit35 * adv_mom_5 &
4968	+ ibit34 * adv_mom_3 &
4969	) * &
4970	( w(k_pp,j,i) + w(k-1,j,i) ) &
4971	+ ( ibit35 * adv_mom_5 &
4972	) * &
4973	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
4974	)
4975
4976	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
4977	( 10.0_wp * ibit35 * adv_mom_5 &
4978	+ 3.0_wp * ibit34 * adv_mom_3 &
4979	+ ibit33 * adv_mom_1 &
4980	) * &
4981	( w(k+1,j,i) - w(k,j,i) ) &
4982	- ( 5.0_wp * ibit35 * adv_mom_5 &
4983	+ ibit34 * adv_mom_3 &
4984	) * &
4985	( w(k_pp,j,i) - w(k-1,j,i) ) &
4986	+ ( ibit35 * adv_mom_5 &
4987	) * &
4988	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
4989	)
4990	!
4991	!-- Calculate the divergence of the velocity field. A respective
4992	!-- correction is needed to overcome numerical instabilities caused
4993	!-- by a not sufficient reduction of divergences near topography.
4994	div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) &
4995	- ( u(k+1,j,i) + u(k,j,i) ) &
4996	* ( IBITS(advc_flags_1(k,j,i-1),27,1) &
4997	+ IBITS(advc_flags_1(k,j,i-1),28,1) &
4998	+ IBITS(advc_flags_1(k,j,i-1),29,1) &
4999	) &
5000	) * rho_air_zw(k) * ddx &
5001	+ ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) &
5002	- ( v(k+1,j,i) + v(k,j,i) ) &
5003	* ( IBITS(advc_flags_1(k,j-1,i),30,1) &
5004	+ IBITS(advc_flags_1(k,j-1,i),31,1) &
5005	+ IBITS(advc_flags_2(k,j-1,i),0,1) &
5006	) &
5007	) * rho_air_zw(k) * ddy &
5008	+ ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) &
5009	- ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
5010	* ( IBITS(advc_flags_2(k-1,j,i),1,1) &
5011	+ IBITS(advc_flags_2(k-1,j,i),2,1) &
5012	+ IBITS(advc_flags_2(k-1,j,i),3,1) &
5013	) &
5014	) * ddzu(k+1) &
5015	) * 0.5_wp
5016
5017
5018
5019	tend(k,j,i) = tend(k,j,i) - ( &
5020	( flux_r(k) + diss_r(k) &
5021	- swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) &
5022	) * ddx &
5023	+ ( flux_n(k) + diss_n(k) &
5024	- swap_flux_y_local_w(k) - swap_diss_y_local_w(k) &
5025	) * ddy &
5026	+ ( ( flux_t(k) + diss_t(k) ) &
5027	- ( flux_d + diss_d ) &
5028	) * drho_air_zw(k) * ddzu(k+1) &
5029	) + div * w(k,j,i)
5030
5031	swap_flux_x_local_w(k,j) = flux_r(k)
5032	swap_diss_x_local_w(k,j) = diss_r(k)
5033	swap_flux_y_local_w(k) = flux_n(k)
5034	swap_diss_y_local_w(k) = diss_n(k)
5035	flux_d = flux_t(k)
5036	diss_d = diss_t(k)
5037
5038	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
5039	+ ( flux_t(k) &
5040	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5041	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
5042	+ diss_t(k) &
5043	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5044	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
5045	) * weight_substep(intermediate_timestep_count)
5046
5047	ENDDO
5048
5049	DO k = nzb_max+1, nzt
5050
5051	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
5052	flux_r(k) = u_comp * ( &
5053	37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) &
5054	- 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) &
5055	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
5056
5057	diss_r(k) = - ABS( u_comp ) * ( &
5058	10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) &
5059	- 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) &
5060	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
5061
5062	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
5063	flux_n(k) = v_comp * ( &
5064	37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) &
5065	- 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) &
5066	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
5067
5068	diss_n(k) = - ABS( v_comp ) * ( &
5069	10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) &
5070	- 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) &
5071	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
5072	!
5073	!-- k index has to be modified near bottom and top, else array
5074	!-- subscripts will be exceeded.
5075	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
5076	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
5077	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
5078
5079	k_ppp = k + 3 * ibit35
5080	k_pp = k + 2 * ( 1 - ibit33 )
5081	k_mm = k - 2 * ibit35
5082
5083	w_comp = w(k+1,j,i) + w(k,j,i)
5084	flux_t(k) = w_comp * rho_air(k+1) * ( &
5085	( 37.0_wp * ibit35 * adv_mom_5 &
5086	+ 7.0_wp * ibit34 * adv_mom_3 &
5087	+ ibit33 * adv_mom_1 &
5088	) * &
5089	( w(k+1,j,i) + w(k,j,i) ) &
5090	- ( 8.0_wp * ibit35 * adv_mom_5 &
5091	+ ibit34 * adv_mom_3 &
5092	) * &
5093	( w(k_pp,j,i) + w(k-1,j,i) ) &
5094	+ ( ibit35 * adv_mom_5 &
5095	) * &
5096	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
5097	)
5098
5099	diss_t(k) = - ABS( w_comp ) * rho_air(k+1) * ( &
5100	( 10.0_wp * ibit35 * adv_mom_5 &
5101	+ 3.0_wp * ibit34 * adv_mom_3 &
5102	+ ibit33 * adv_mom_1 &
5103	) * &
5104	( w(k+1,j,i) - w(k,j,i) ) &
5105	- ( 5.0_wp * ibit35 * adv_mom_5 &
5106	+ ibit34 * adv_mom_3 &
5107	) * &
5108	( w(k_pp,j,i) - w(k-1,j,i) ) &
5109	+ ( ibit35 * adv_mom_5 &
5110	) * &
5111	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
5112	)
5113	!
5114	!-- Calculate the divergence of the velocity field. A respective
5115	!-- correction is needed to overcome numerical instabilities caused
5116	!-- by a not sufficient reduction of divergences near topography.
5117	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
5118	* rho_air_zw(k) &
5119	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
5120	* rho_air_zw(k) &
5121	+ ( w_comp * rho_air(k+1) - &
5122	( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
5123	) * ddzu(k+1) &
5124	) * 0.5_wp
5125
5126	tend(k,j,i) = tend(k,j,i) - ( &
5127	( flux_r(k) + diss_r(k) &
5128	- swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) &
5129	) * ddx &
5130	+ ( flux_n(k) + diss_n(k) &
5131	- swap_flux_y_local_w(k) - swap_diss_y_local_w(k) &
5132	) * ddy &
5133	+ ( ( flux_t(k) + diss_t(k) ) &
5134	- ( flux_d + diss_d ) &
5135	) * drho_air_zw(k) * ddzu(k+1) &
5136	) + div * w(k,j,i)
5137
5138	swap_flux_x_local_w(k,j) = flux_r(k)
5139	swap_diss_x_local_w(k,j) = diss_r(k)
5140	swap_flux_y_local_w(k) = flux_n(k)
5141	swap_diss_y_local_w(k) = diss_n(k)
5142	flux_d = flux_t(k)
5143	diss_d = diss_t(k)
5144
5145	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
5146	+ ( flux_t(k) &
5147	* ( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5148	/ ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) &
5149	+ diss_t(k) &
5150	* ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) &
5151	/ ( ABS( w_comp ) + 1.0E-20_wp ) &
5152	) * weight_substep(intermediate_timestep_count)
5153
5154	ENDDO
5155	ENDDO
5156	ENDDO
5157
5158	END SUBROUTINE advec_w_ws
5159
5160	END MODULE advec_ws

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |