Home

Context Navigation

advec_ws.f90 @ 3298

Last change on this file since 3298 was 3298, checked in by kanani, 6 years ago
Merge chemistry branch at r3297 to trunk
Property svn:keywords set to `Id`
File size: 280.0 KB

Line
1	!> @file advec_ws.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: advec_ws.f90 3298 2018-10-02 12:21:11Z kanani $
27	! Add formatted note from ketelsen about statistics for chemistry
28	!
29	! 3294 2018-10-01 02:37:10Z raasch
30	! ocean renamed ocean_mode
31	!
32	! 3274 2018-09-24 15:42:55Z knoop
33	! Modularization of all bulk cloud physics code components
34	!
35	! 3241 2018-09-12 15:02:00Z raasch
36	! unused variables removed
37	!
38	! 3183 2018-07-27 14:25:55Z suehring
39	! Rename variables for boundary-condition flags and for offline nesting mode
40	!
41	! 3182 2018-07-27 13:36:03Z suehring
42	! Bugfix in calculation of left-sided fluxes for u-component in OpenMP case.
43	!
44	! 2731 2018-01-09 17:44:02Z suehring
45	! Enable loop vectorization by splitting the k-loop
46	!
47	! 2718 2018-01-02 08:49:38Z maronga
48	! Corrected "Former revisions" section
49	!
50	! 2696 2017-12-14 17:12:51Z kanani
51	! Change in file header (GPL part)
52	! Implement advection for TKE-dissipation in case of RANS-TKE-e closure (TG)
53	! Allocate advc_flags_1/2 within ws_init_flags instead of init_grid
54	! Change argument list for exchange_horiz_2d_int (MS)
55	!
56	! 2329 2017-08-03 14:24:56Z knoop
57	! Bugfix concerning density in divergence correction close to buildings
58	!
59	! 2292 2017-06-20 09:51:42Z schwenkel
60	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
61	! includes two more prognostic equations for cloud drop concentration (nc)
62	! and cloud water content (qc).
63	!
64	! 2233 2017-05-30 18:08:54Z suehring
65	!
66	! 2232 2017-05-30 17:47:52Z suehring
67	! Rename wall_flags_0 and wall_flags_00 into advc_flags_1 and advc_flags_2,
68	! respectively.
69	! Set advc_flags_1/2 on basis of wall_flags_0/00 instead of nzb_s/u/v/w_inner.
70	! Setting advc_flags_1/2 also for downward-facing walls
71	!
72	! 2200 2017-04-11 11:37:51Z suehring
73	! monotonic_adjustment removed
74	!
75	! 2118 2017-01-17 16:38:49Z raasch
76	! OpenACC version of subroutines removed
77	!
78	! 2037 2016-10-26 11:15:40Z knoop
79	! Anelastic approximation implemented
80	!
81	! 2000 2016-08-20 18:09:15Z knoop
82	! Forced header and separation lines into 80 columns
83	!
84	! 1996 2016-08-18 11:42:29Z suehring
85	! Bugfix concerning calculation of turbulent of turbulent fluxes
86	!
87	! 1960 2016-07-12 16:34:24Z suehring
88	! Separate humidity and passive scalar
89	!
90	! 1942 2016-06-14 12:18:18Z suehring
91	! Initialization of flags for ws-scheme moved from init_grid.
92	!
93	! 1873 2016-04-18 14:50:06Z maronga
94	! Module renamed (removed _mod)
95	!
96	!
97	! 1850 2016-04-08 13:29:27Z maronga
98	! Module renamed
99	!
100	!
101	! 1822 2016-04-07 07:49:42Z hoffmann
102	! icloud_scheme removed, microphysics_seifert added
103	!
104	! 1682 2015-10-07 23:56:08Z knoop
105	! Code annotations made doxygen readable
106	!
107	! 1630 2015-08-26 16:57:23Z suehring
108	!
109	!
110	! 1629 2015-08-26 16:56:11Z suehring
111	! Bugfix concerning wall_flags at left and south PE boundaries
112	!
113	! 1581 2015-04-10 13:45:59Z suehring
114	!
115	!
116	! 1580 2015-04-10 13:43:49Z suehring
117	! Bugfix: statistical evaluation of scalar fluxes in case of monotonic limiter
118	!
119	! 1567 2015-03-10 17:57:55Z suehring
120	! Bugfixes in monotonic limiter.
121	!
122	! 2015-03-09 13:10:37Z heinze
123	! Bugfix: REAL constants provided with KIND-attribute in call of
124	! intrinsic functions like MAX and MIN
125	!
126	! 1557 2015-03-05 16:43:04Z suehring
127	! Enable monotone advection for scalars using monotonic limiter
128	!
129	! 1374 2014-04-25 12:55:07Z raasch
130	! missing variables added to ONLY list
131	!
132	! 1361 2014-04-16 15:17:48Z hoffmann
133	! accelerator and vector version for qr and nr added
134	!
135	! 1353 2014-04-08 15:21:23Z heinze
136	! REAL constants provided with KIND-attribute,
137	! module kinds added
138	! some formatting adjustments
139	!
140	! 1322 2014-03-20 16:38:49Z raasch
141	! REAL constants defined as wp-kind
142	!
143	! 1320 2014-03-20 08:40:49Z raasch
144	! ONLY-attribute added to USE-statements,
145	! kind-parameters added to all INTEGER and REAL declaration statements,
146	! kinds are defined in new module kinds,
147	! old module precision_kind is removed,
148	! revision history before 2012 removed,
149	! comment fields (!:) to be used for variable explanations added to
150	! all variable declaration statements
151	!
152	! 1257 2013-11-08 15:18:40Z raasch
153	! accelerator loop directives removed
154	!
155	! 1221 2013-09-10 08:59:13Z raasch
156	! wall_flags_00 introduced, which holds bits 32-...
157	!
158	! 1128 2013-04-12 06:19:32Z raasch
159	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
160	! j_north
161	!
162	! 1115 2013-03-26 18:16:16Z hoffmann
163	! calculation of qr and nr is restricted to precipitation
164	!
165	! 1053 2012-11-13 17:11:03Z hoffmann
166	! necessary expansions according to the two new prognostic equations (nr, qr)
167	! of the two-moment cloud physics scheme:
168	! +flux_l_, flux_s_, diss_l_, diss_s_, sums_ws*s_ws_l
169	!
170	! 1036 2012-10-22 13:43:42Z raasch
171	! code put under GPL (PALM 3.9)
172	!
173	! 1027 2012-10-15 17:18:39Z suehring
174	! Bugfix in calculation indices k_mm, k_pp in accelerator version
175	!
176	! 1019 2012-09-28 06:46:45Z raasch
177	! small change in comment lines
178	!
179	! 1015 2012-09-27 09:23:24Z raasch
180	! accelerator versions (*_acc) added
181	!
182	! 1010 2012-09-20 07:59:54Z raasch
183	! cpp switch __nopointer added for pointer free version
184	!
185	! 888 2012-04-20 15:03:46Z suehring
186	! Number of IBITS() calls with identical arguments is reduced.
187	!
188	! 862 2012-03-26 14:21:38Z suehring
189	! ws-scheme also work with topography in combination with vector version.
190	! ws-scheme also work with outflow boundaries in combination with
191	! vector version.
192	! Degradation of the applied order of scheme is now steered by multiplying with
193	! Integer advc_flags_1. 2nd order scheme, WS3 and WS5 are calculated on each
194	! grid point and mulitplied with the appropriate flag.
195	! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
196	! term derived according to the 4th and 6th order terms is applied. It turns
197	! out that diss_2nd does not provide sufficient dissipation near walls.
198	! Therefore, the function diss_2nd is removed.
199	! Near walls a divergence correction is necessary to overcome numerical
200	! instabilities due to too less divergence reduction of the velocity field.
201	! boundary_flags and logicals steering the degradation are removed.
202	! Empty SUBROUTINE local_diss removed.
203	! Further formatting adjustments.
204	!
205	! 801 2012-01-10 17:30:36Z suehring
206	! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
207	! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
208	! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
209	! dimension.
210	!
211	! Initial revision
212	!
213	! 411 2009-12-11 12:31:43 Z suehring
214	!
215	! Description:
216	! ------------
217	!> Advection scheme for scalars and momentum using the flux formulation of
218	!> Wicker and Skamarock 5th order. Additionally the module contains of a
219	!> routine using for initialisation and steering of the statical evaluation.
220	!> The computation of turbulent fluxes takes place inside the advection
221	!> routines.
222	!> Near non-cyclic boundaries the order of the applied advection scheme is
223	!> degraded.
224	!> A divergence correction is applied. It is necessary for topography, since
225	!> the divergence is not sufficiently reduced, resulting in erroneous fluxes and
226	!> partly numerical instabilities.
227	!-----------------------------------------------------------------------------!
228	MODULE advec_ws
229
230	USE kinds
231
232	IMPLICIT NONE
233
234	REAL(wp) :: adv_mom_1 !< 1/4 - constant used in 5th-order advection scheme for momentum advection (1st-order part)
235	REAL(wp) :: adv_mom_3 !< 1/24 - constant used in 5th-order advection scheme for momentum advection (3rd-order part)
236	REAL(wp) :: adv_mom_5 !< 1/120 - constant used in 5th-order advection scheme for momentum advection (5th-order part)
237	REAL(wp) :: adv_sca_1 !< 1/2 - constant used in 5th-order advection scheme for scalar advection (1st-order part)
238	REAL(wp) :: adv_sca_3 !< 1/12 - constant used in 5th-order advection scheme for scalar advection (3rd-order part)
239	REAL(wp) :: adv_sca_5 !< 1/60 - constant used in 5th-order advection scheme for scalar advection (5th-order part)
240
241	PRIVATE
242	PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init, &
243	ws_init_flags, ws_statistics
244
245	INTERFACE ws_init
246	MODULE PROCEDURE ws_init
247	END INTERFACE ws_init
248
249	INTERFACE ws_init_flags
250	MODULE PROCEDURE ws_init_flags
251	END INTERFACE ws_init_flags
252
253	INTERFACE ws_statistics
254	MODULE PROCEDURE ws_statistics
255	END INTERFACE ws_statistics
256
257	INTERFACE advec_s_ws
258	MODULE PROCEDURE advec_s_ws
259	MODULE PROCEDURE advec_s_ws_ij
260	END INTERFACE advec_s_ws
261
262	INTERFACE advec_u_ws
263	MODULE PROCEDURE advec_u_ws
264	MODULE PROCEDURE advec_u_ws_ij
265	END INTERFACE advec_u_ws
266
267	INTERFACE advec_v_ws
268	MODULE PROCEDURE advec_v_ws
269	MODULE PROCEDURE advec_v_ws_ij
270	END INTERFACE advec_v_ws
271
272	INTERFACE advec_w_ws
273	MODULE PROCEDURE advec_w_ws
274	MODULE PROCEDURE advec_w_ws_ij
275	END INTERFACE advec_w_ws
276
277	CONTAINS
278
279
280	!------------------------------------------------------------------------------!
281	! Description:
282	! ------------
283	!> Initialization of WS-scheme
284	!------------------------------------------------------------------------------!
285	SUBROUTINE ws_init
286
287	USE arrays_3d, &
288	ONLY: diss_l_diss, diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, &
289	diss_l_q, diss_l_qc, diss_l_qr, diss_l_s, diss_l_sa, &
290	diss_l_u, diss_l_v, diss_l_w, flux_l_diss, flux_l_e, &
291	flux_l_nc, flux_l_nr, flux_l_pt, flux_l_q, flux_l_qc, &
292	flux_l_qr, flux_l_s, flux_l_sa, flux_l_u, flux_l_v, &
293	flux_l_w, diss_s_diss, diss_s_e, diss_s_nc, diss_s_nr, &
294	diss_s_pt, diss_s_q, diss_s_qc, diss_s_qr, diss_s_s, &
295	diss_s_sa, diss_s_u, diss_s_v, diss_s_w, flux_s_diss, &
296	flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, flux_s_q, &
297	flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, flux_s_u, &
298	flux_s_v, flux_s_w
299
300	USE control_parameters, &
301	ONLY: humidity, loop_optimization, passive_scalar, ocean_mode, &
302	rans_tke_e, ws_scheme_mom, ws_scheme_sca
303
304	USE indices, &
305	ONLY: nyn, nys, nzb, nzt
306
307	USE kinds
308
309	USE bulk_cloud_model_mod, &
310	ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
311
312	USE pegrid
313
314	USE statistics, &
315	ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
316	sums_wsnrs_ws_l,sums_wspts_ws_l, sums_wsqcs_ws_l, &
317	sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l, &
318	sums_wssas_ws_l, sums_wsss_ws_l, sums_wsus_ws_l, &
319	sums_wsvs_ws_l
320
321
322	!
323	!-- Set the appropriate factors for scalar and momentum advection.
324	adv_sca_5 = 1.0_wp / 60.0_wp
325	adv_sca_3 = 1.0_wp / 12.0_wp
326	adv_sca_1 = 1.0_wp / 2.0_wp
327	adv_mom_5 = 1.0_wp / 120.0_wp
328	adv_mom_3 = 1.0_wp / 24.0_wp
329	adv_mom_1 = 1.0_wp / 4.0_wp
330	!
331	!-- Arrays needed for statical evaluation of fluxes.
332	IF ( ws_scheme_mom ) THEN
333
334	ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), &
335	sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), &
336	sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
337	sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
338	sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
339
340	sums_wsus_ws_l = 0.0_wp
341	sums_wsvs_ws_l = 0.0_wp
342	sums_us2_ws_l = 0.0_wp
343	sums_vs2_ws_l = 0.0_wp
344	sums_ws2_ws_l = 0.0_wp
345
346	ENDIF
347
348	IF ( ws_scheme_sca ) THEN
349
350	ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
351	sums_wspts_ws_l = 0.0_wp
352
353	IF ( humidity ) THEN
354	ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
355	sums_wsqs_ws_l = 0.0_wp
356	ENDIF
357
358	IF ( passive_scalar ) THEN
359	ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
360	sums_wsss_ws_l = 0.0_wp
361	ENDIF
362
363	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
364	ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
365	ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
366	sums_wsqcs_ws_l = 0.0_wp
367	sums_wsncs_ws_l = 0.0_wp
368	ENDIF
369
370	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
371	ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
372	ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
373	sums_wsqrs_ws_l = 0.0_wp
374	sums_wsnrs_ws_l = 0.0_wp
375	ENDIF
376
377	IF ( ocean_mode ) THEN
378	ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
379	sums_wssas_ws_l = 0.0_wp
380	ENDIF
381
382	ENDIF
383
384	!
385	!-- Arrays needed for reasons of speed optimization for cache version.
386	!-- For the vector version the buffer arrays are not necessary,
387	!-- because the the fluxes can swapped directly inside the loops of the
388	!-- advection routines.
389	IF ( loop_optimization /= 'vector' ) THEN
390
391	IF ( ws_scheme_mom ) THEN
392
393	ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &
394	flux_s_v(nzb+1:nzt,0:threads_per_task-1), &
395	flux_s_w(nzb+1:nzt,0:threads_per_task-1), &
396	diss_s_u(nzb+1:nzt,0:threads_per_task-1), &
397	diss_s_v(nzb+1:nzt,0:threads_per_task-1), &
398	diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
399	ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
400	flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
401	flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
402	diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
403	diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
404	diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
405
406	ENDIF
407
408	IF ( ws_scheme_sca ) THEN
409
410	ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), &
411	flux_s_e(nzb+1:nzt,0:threads_per_task-1), &
412	diss_s_pt(nzb+1:nzt,0:threads_per_task-1), &
413	diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
414	ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
415	flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
416	diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
417	diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
418
419	IF ( rans_tke_e ) THEN
420	ALLOCATE( flux_s_diss(nzb+1:nzt,0:threads_per_task-1), &
421	diss_s_diss(nzb+1:nzt,0:threads_per_task-1) )
422	ALLOCATE( flux_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
423	diss_l_diss(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
424	ENDIF
425
426	IF ( humidity ) THEN
427	ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), &
428	diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
429	ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
430	diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
431	ENDIF
432
433	IF ( passive_scalar ) THEN
434	ALLOCATE( flux_s_s(nzb+1:nzt,0:threads_per_task-1), &
435	diss_s_s(nzb+1:nzt,0:threads_per_task-1) )
436	ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
437	diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
438	ENDIF
439
440	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
441	ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1), &
442	diss_s_qc(nzb+1:nzt,0:threads_per_task-1), &
443	flux_s_nc(nzb+1:nzt,0:threads_per_task-1), &
444	diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
445	ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
446	diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
447	flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
448	diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
449	ENDIF
450
451	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
452	ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), &
453	diss_s_qr(nzb+1:nzt,0:threads_per_task-1), &
454	flux_s_nr(nzb+1:nzt,0:threads_per_task-1), &
455	diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
456	ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
457	diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
458	flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
459	diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
460	ENDIF
461
462	IF ( ocean_mode ) THEN
463	ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), &
464	diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
465	ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
466	diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
467	ENDIF
468
469	ENDIF
470
471	ENDIF
472
473	END SUBROUTINE ws_init
474
475	!------------------------------------------------------------------------------!
476	! Description:
477	! ------------
478	!> Initialization of flags for WS-scheme used to degrade the order of the scheme
479	!> near walls.
480	!------------------------------------------------------------------------------!
481	SUBROUTINE ws_init_flags
482
483	USE control_parameters, &
484	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
485	bc_dirichlet_s, bc_radiation_l, bc_radiation_n, &
486	bc_radiation_r, bc_radiation_s, momentum_advec, scalar_advec
487
488	USE indices, &
489	ONLY: advc_flags_1, advc_flags_2, nbgp, nxl, nxlg, nxlu, nxr, nxrg,&
490	nyn, nyng, nys, nysg, nysv, nzb, nzt, wall_flags_0
491
492	USE kinds
493
494	IMPLICIT NONE
495
496	INTEGER(iwp) :: i !< index variable along x
497	INTEGER(iwp) :: j !< index variable along y
498	INTEGER(iwp) :: k !< index variable along z
499	INTEGER(iwp) :: k_mm !< dummy index along z
500	INTEGER(iwp) :: k_pp !< dummy index along z
501	INTEGER(iwp) :: k_ppp !< dummy index along z
502
503	LOGICAL :: flag_set !< steering variable for advection flags
504
505	ALLOCATE( advc_flags_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
506	advc_flags_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
507	advc_flags_1 = 0
508	advc_flags_2 = 0
509	IF ( scalar_advec == 'ws-scheme' ) THEN
510	!
511	!-- Set flags to steer the degradation of the advection scheme in advec_ws
512	!-- near topography, inflow- and outflow boundaries as well as bottom and
513	!-- top of model domain. advc_flags_1 remains zero for all non-prognostic
514	!-- grid points.
515	DO i = nxl, nxr
516	DO j = nys, nyn
517	DO k = nzb+1, nzt
518	!
519	!-- scalar - x-direction
520	!-- WS1 (0), WS3 (1), WS5 (2)
521	IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0) &
522	.OR. .NOT. BTEST(wall_flags_0(k,j,i+2),0) &
523	.OR. .NOT. BTEST(wall_flags_0(k,j,i-1),0) ) &
524	.OR. ( ( bc_dirichlet_l .OR. bc_radiation_l ) &
525	.AND. i == nxl ) &
526	.OR. ( ( bc_dirichlet_r .OR. bc_radiation_r ) &
527	.AND. i == nxr ) ) &
528	THEN
529	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 0 )
530	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0) &
531	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
532	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
533	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
534	) .OR. &
535	( .NOT. BTEST(wall_flags_0(k,j,i-2),0) &
536	.AND. BTEST(wall_flags_0(k,j,i+1),0) &
537	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
538	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
539	) &
540	.OR. &
541	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
542	.AND. i == nxr-1 ) .OR. &
543	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
544	.AND. i == nxlu ) ) & ! why not nxl+1
545	THEN
546	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 1 )
547	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0) &
548	.AND. BTEST(wall_flags_0(k,j,i+2),0) &
549	.AND. BTEST(wall_flags_0(k,j,i+3),0) &
550	.AND. BTEST(wall_flags_0(k,j,i-1),0) &
551	.AND. BTEST(wall_flags_0(k,j,i-2),0) ) &
552	THEN
553	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 2 )
554	ENDIF
555	!
556	!-- scalar - y-direction
557	!-- WS1 (3), WS3 (4), WS5 (5)
558	IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0) &
559	.OR. .NOT. BTEST(wall_flags_0(k,j+2,i),0) &
560	.OR. .NOT. BTEST(wall_flags_0(k,j-1,i),0)) &
561	.OR. ( ( bc_dirichlet_s .OR. bc_radiation_s ) &
562	.AND. j == nys ) &
563	.OR. ( ( bc_dirichlet_n .OR. bc_radiation_n ) &
564	.AND. j == nyn ) ) &
565	THEN
566	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 3 )
567	!
568	!-- WS3
569	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0) &
570	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
571	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
572	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
573	) .OR. &
574	( .NOT. BTEST(wall_flags_0(k,j-2,i),0) &
575	.AND. BTEST(wall_flags_0(k,j+1,i),0) &
576	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
577	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
578	) &
579	.OR. &
580	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
581	.AND. j == nysv ) .OR. & ! why not nys+1
582	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
583	.AND. j == nyn-1 ) ) &
584	THEN
585	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 4 )
586	!
587	!-- WS5
588	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0) &
589	.AND. BTEST(wall_flags_0(k,j+2,i),0) &
590	.AND. BTEST(wall_flags_0(k,j+3,i),0) &
591	.AND. BTEST(wall_flags_0(k,j-1,i),0) &
592	.AND. BTEST(wall_flags_0(k,j-2,i),0) ) &
593	THEN
594	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 5 )
595	ENDIF
596	!
597	!-- scalar - z-direction
598	!-- WS1 (6), WS3 (7), WS5 (8)
599	IF ( k == nzb+1 ) THEN
600	k_mm = nzb
601	ELSE
602	k_mm = k - 2
603	ENDIF
604	IF ( k > nzt-1 ) THEN
605	k_pp = nzt+1
606	ELSE
607	k_pp = k + 2
608	ENDIF
609	IF ( k > nzt-2 ) THEN
610	k_ppp = nzt+1
611	ELSE
612	k_ppp = k + 3
613	ENDIF
614
615	flag_set = .FALSE.
616	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0) .AND. &
617	BTEST(wall_flags_0(k,j,i),0) .OR. &
618	.NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND. &
619	BTEST(wall_flags_0(k,j,i),0) .OR. &
620	k == nzt ) &
621	THEN
622	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 6 )
623	flag_set = .TRUE.
624	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0) .OR. &
625	.NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND. &
626	BTEST(wall_flags_0(k-1,j,i),0) .AND. &
627	BTEST(wall_flags_0(k,j,i),0) .AND. &
628	BTEST(wall_flags_0(k+1,j,i),0) .AND. &
629	BTEST(wall_flags_0(k_pp,j,i),0) .OR. &
630	k == nzt - 1 ) &
631	THEN
632	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 7 )
633	flag_set = .TRUE.
634	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0) &
635	.AND. BTEST(wall_flags_0(k-1,j,i),0) &
636	.AND. BTEST(wall_flags_0(k,j,i),0) &
637	.AND. BTEST(wall_flags_0(k+1,j,i),0) &
638	.AND. BTEST(wall_flags_0(k_pp,j,i),0) &
639	.AND. BTEST(wall_flags_0(k_ppp,j,i),0) &
640	.AND. .NOT. flag_set ) &
641	THEN
642	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 8 )
643	ENDIF
644
645	ENDDO
646	ENDDO
647	ENDDO
648	ENDIF
649
650	IF ( momentum_advec == 'ws-scheme' ) THEN
651	!
652	!-- Set advc_flags_1 to steer the degradation of the advection scheme in advec_ws
653	!-- near topography, inflow- and outflow boundaries as well as bottom and
654	!-- top of model domain. advc_flags_1 remains zero for all non-prognostic
655	!-- grid points.
656	DO i = nxl, nxr
657	DO j = nys, nyn
658	DO k = nzb+1, nzt
659	!
660	!-- At first, set flags to WS1.
661	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
662	!-- in order to handle the left/south flux.
663	!-- near vertical walls.
664	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
665	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
666	!
667	!-- u component - x-direction
668	!-- WS1 (9), WS3 (10), WS5 (11)
669	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1) .OR. &
670	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
671	.AND. i <= nxlu ) .OR. &
672	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
673	.AND. i == nxr ) ) &
674	THEN
675	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 )
676	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1) .AND. &
677	BTEST(wall_flags_0(k,j,i+1),1) .OR. &
678	.NOT. BTEST(wall_flags_0(k,j,i-1),1) ) &
679	.OR. &
680	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
681	.AND. i == nxr-1 ) .OR. &
682	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
683	.AND. i == nxlu+1) ) &
684	THEN
685	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 10 )
686	!
687	!-- Clear flag for WS1
688	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
689	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1) .AND. &
690	BTEST(wall_flags_0(k,j,i+2),1) .AND. &
691	BTEST(wall_flags_0(k,j,i-1),1) ) &
692	THEN
693	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 11 )
694	!
695	!-- Clear flag for WS1
696	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 )
697	ENDIF
698	!
699	!-- u component - y-direction
700	!-- WS1 (12), WS3 (13), WS5 (14)
701	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1) .OR. &
702	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
703	.AND. j == nys ) .OR. &
704	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
705	.AND. j == nyn ) ) &
706	THEN
707	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 )
708	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1) .AND. &
709	BTEST(wall_flags_0(k,j+1,i),1) .OR. &
710	.NOT. BTEST(wall_flags_0(k,j-1,i),1) ) &
711	.OR. &
712	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
713	.AND. j == nysv ) .OR. &
714	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
715	.AND. j == nyn-1 ) ) &
716	THEN
717	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 13 )
718	!
719	!-- Clear flag for WS1
720	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
721	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1) .AND. &
722	BTEST(wall_flags_0(k,j+2,i),1) .AND. &
723	BTEST(wall_flags_0(k,j-1,i),1) ) &
724	THEN
725	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 14 )
726	!
727	!-- Clear flag for WS1
728	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 )
729	ENDIF
730	!
731	!-- u component - z-direction
732	!-- WS1 (15), WS3 (16), WS5 (17)
733	IF ( k == nzb+1 ) THEN
734	k_mm = nzb
735	ELSE
736	k_mm = k - 2
737	ENDIF
738	IF ( k > nzt-1 ) THEN
739	k_pp = nzt+1
740	ELSE
741	k_pp = k + 2
742	ENDIF
743	IF ( k > nzt-2 ) THEN
744	k_ppp = nzt+1
745	ELSE
746	k_ppp = k + 3
747	ENDIF
748
749	flag_set = .FALSE.
750	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1) .AND. &
751	BTEST(wall_flags_0(k,j,i),1) .OR. &
752	.NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
753	BTEST(wall_flags_0(k,j,i),1) .OR. &
754	k == nzt ) &
755	THEN
756	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 15 )
757	flag_set = .TRUE.
758	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1) .OR. &
759	.NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND. &
760	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
761	BTEST(wall_flags_0(k,j,i),1) .AND. &
762	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
763	BTEST(wall_flags_0(k_pp,j,i),1) .OR. &
764	k == nzt - 1 ) &
765	THEN
766	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 16 )
767	flag_set = .TRUE.
768	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1) .AND. &
769	BTEST(wall_flags_0(k-1,j,i),1) .AND. &
770	BTEST(wall_flags_0(k,j,i),1) .AND. &
771	BTEST(wall_flags_0(k+1,j,i),1) .AND. &
772	BTEST(wall_flags_0(k_pp,j,i),1) .AND. &
773	BTEST(wall_flags_0(k_ppp,j,i),1) .AND. &
774	.NOT. flag_set ) &
775	THEN
776	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 17 )
777	ENDIF
778
779	ENDDO
780	ENDDO
781	ENDDO
782
783	DO i = nxl, nxr
784	DO j = nys, nyn
785	DO k = nzb+1, nzt
786	!
787	!-- At first, set flags to WS1.
788	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
789	!-- in order to handle the left/south flux.
790	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
791	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
792	!
793	!-- v component - x-direction
794	!-- WS1 (18), WS3 (19), WS5 (20)
795	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2) .OR. &
796	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
797	.AND. i == nxl ) .OR. &
798	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
799	.AND. i == nxr ) ) &
800	THEN
801	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 )
802	!
803	!-- WS3
804	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2) .AND. &
805	BTEST(wall_flags_0(k,j,i+1),2) ) .OR. &
806	.NOT. BTEST(wall_flags_0(k,j,i-1),2) &
807	.OR. &
808	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
809	.AND. i == nxr-1 ) .OR. &
810	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
811	.AND. i == nxlu ) ) &
812	THEN
813	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 19 )
814	!
815	!-- Clear flag for WS1
816	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
817	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2) .AND. &
818	BTEST(wall_flags_0(k,j,i+2),2) .AND. &
819	BTEST(wall_flags_0(k,j,i-1),2) ) &
820	THEN
821	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 20 )
822	!
823	!-- Clear flag for WS1
824	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 )
825	ENDIF
826	!
827	!-- v component - y-direction
828	!-- WS1 (21), WS3 (22), WS5 (23)
829	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR. &
830	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
831	.AND. j <= nysv ) .OR. &
832	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
833	.AND. j == nyn ) ) &
834	THEN
835	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 )
836	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2) .AND. &
837	BTEST(wall_flags_0(k,j+1,i),2) .OR. &
838	.NOT. BTEST(wall_flags_0(k,j-1,i),2) ) &
839	.OR. &
840	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
841	.AND. j == nysv+1) .OR. &
842	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
843	.AND. j == nyn-1 ) ) &
844	THEN
845	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 22 )
846	!
847	!-- Clear flag for WS1
848	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
849	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2) .AND. &
850	BTEST(wall_flags_0(k,j+2,i),2) .AND. &
851	BTEST(wall_flags_0(k,j-1,i),2) ) &
852	THEN
853	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 23 )
854	!
855	!-- Clear flag for WS1
856	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 )
857	ENDIF
858	!
859	!-- v component - z-direction
860	!-- WS1 (24), WS3 (25), WS5 (26)
861	IF ( k == nzb+1 ) THEN
862	k_mm = nzb
863	ELSE
864	k_mm = k - 2
865	ENDIF
866	IF ( k > nzt-1 ) THEN
867	k_pp = nzt+1
868	ELSE
869	k_pp = k + 2
870	ENDIF
871	IF ( k > nzt-2 ) THEN
872	k_ppp = nzt+1
873	ELSE
874	k_ppp = k + 3
875	ENDIF
876
877	flag_set = .FALSE.
878	IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2) .AND. &
879	BTEST(wall_flags_0(k,j,i),2) .OR. &
880	.NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
881	BTEST(wall_flags_0(k,j,i),2) .OR. &
882	k == nzt ) &
883	THEN
884	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 24 )
885	flag_set = .TRUE.
886	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2) .OR. &
887	.NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND. &
888	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
889	BTEST(wall_flags_0(k,j,i),2) .AND. &
890	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
891	BTEST(wall_flags_0(k_pp,j,i),2) .OR. &
892	k == nzt - 1 ) &
893	THEN
894	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 25 )
895	flag_set = .TRUE.
896	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2) .AND. &
897	BTEST(wall_flags_0(k-1,j,i),2) .AND. &
898	BTEST(wall_flags_0(k,j,i),2) .AND. &
899	BTEST(wall_flags_0(k+1,j,i),2) .AND. &
900	BTEST(wall_flags_0(k_pp,j,i),2) .AND. &
901	BTEST(wall_flags_0(k_ppp,j,i),2) .AND. &
902	.NOT. flag_set ) &
903	THEN
904	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 26 )
905	ENDIF
906
907	ENDDO
908	ENDDO
909	ENDDO
910	DO i = nxl, nxr
911	DO j = nys, nyn
912	DO k = nzb+1, nzt
913	!
914	!-- At first, set flags to WS1.
915	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
916	!-- in order to handle the left/south flux.
917	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
918	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
919	!
920	!-- w component - x-direction
921	!-- WS1 (27), WS3 (28), WS5 (29)
922	IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR. &
923	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
924	.AND. i == nxl ) .OR. &
925	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
926	.AND. i == nxr ) ) &
927	THEN
928	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 )
929	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3) .AND. &
930	BTEST(wall_flags_0(k,j,i+1),3) .OR. &
931	.NOT. BTEST(wall_flags_0(k,j,i-1),3) ) &
932	.OR. &
933	( ( bc_dirichlet_r .OR. bc_radiation_r ) &
934	.AND. i == nxr-1 ) .OR. &
935	( ( bc_dirichlet_l .OR. bc_radiation_l ) &
936	.AND. i == nxlu ) ) &
937	THEN
938	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 28 )
939	!
940	!-- Clear flag for WS1
941	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
942	ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3) .AND. &
943	BTEST(wall_flags_0(k,j,i+2),3) .AND. &
944	BTEST(wall_flags_0(k,j,i-1),3) ) &
945	THEN
946	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i),29 )
947	!
948	!-- Clear flag for WS1
949	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 )
950	ENDIF
951	!
952	!-- w component - y-direction
953	!-- WS1 (30), WS3 (31), WS5 (32)
954	IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR. &
955	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
956	.AND. j == nys ) .OR. &
957	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
958	.AND. j == nyn ) ) &
959	THEN
960	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 )
961	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3) .AND. &
962	BTEST(wall_flags_0(k,j+1,i),3) .OR. &
963	.NOT. BTEST(wall_flags_0(k,j-1,i),3) ) &
964	.OR. &
965	( ( bc_dirichlet_s .OR. bc_radiation_s ) &
966	.AND. j == nysv ) .OR. &
967	( ( bc_dirichlet_n .OR. bc_radiation_n ) &
968	.AND. j == nyn-1 ) ) &
969	THEN
970	advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 31 )
971	!
972	!-- Clear flag for WS1
973	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
974	ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3) .AND. &
975	BTEST(wall_flags_0(k,j+2,i),3) .AND. &
976	BTEST(wall_flags_0(k,j-1,i),3) ) &
977	THEN
978	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 0 )
979	!
980	!-- Clear flag for WS1
981	advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 )
982	ENDIF
983	!
984	!-- w component - z-direction
985	!-- WS1 (33), WS3 (34), WS5 (35)
986	flag_set = .FALSE.
987	IF ( k == nzb+1 ) THEN
988	k_mm = nzb
989	ELSE
990	k_mm = k - 2
991	ENDIF
992	IF ( k > nzt-1 ) THEN
993	k_pp = nzt+1
994	ELSE
995	k_pp = k + 2
996	ENDIF
997	IF ( k > nzt-2 ) THEN
998	k_ppp = nzt+1
999	ELSE
1000	k_ppp = k + 3
1001	ENDIF
1002
1003	IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
1004	.NOT. BTEST(wall_flags_0(k,j,i),3) .AND. &
1005	BTEST(wall_flags_0(k+1,j,i),3) ) .OR. &
1006	( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. &
1007	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
1008	( .NOT. BTEST(wall_flags_0(k+1,j,i),3) .AND. &
1009	BTEST(wall_flags_0(k,j,i),3) ) .OR. &
1010	k == nzt ) &
1011	THEN
1012	!
1013	!-- Please note, at k == nzb_w_inner(j,i) a flag is explictely
1014	!-- set, although this is not a prognostic level. However,
1015	!-- contrary to the advection of u,v and s this is necessary
1016	!-- because flux_t(nzb_w_inner(j,i)) is used for the tendency
1017	!-- at k == nzb_w_inner(j,i)+1.
1018	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 1 )
1019	flag_set = .TRUE.
1020	ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3) .OR. &
1021	.NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND. &
1022	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
1023	BTEST(wall_flags_0(k,j,i),3) .AND. &
1024	BTEST(wall_flags_0(k+1,j,i),3) .OR. &
1025	k == nzt - 1 ) &
1026	THEN
1027	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 2 )
1028	flag_set = .TRUE.
1029	ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3) .AND. &
1030	BTEST(wall_flags_0(k-1,j,i),3) .AND. &
1031	BTEST(wall_flags_0(k,j,i),3) .AND. &
1032	BTEST(wall_flags_0(k+1,j,i),3) .AND. &
1033	BTEST(wall_flags_0(k_pp,j,i),3) .AND. &
1034	BTEST(wall_flags_0(k_ppp,j,i),3) .AND. &
1035	.NOT. flag_set ) &
1036	THEN
1037	advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 3 )
1038	ENDIF
1039
1040	ENDDO
1041	ENDDO
1042	ENDDO
1043
1044	ENDIF
1045
1046
1047	!
1048	!-- Exchange 3D integer wall_flags.
1049	IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' &
1050	) THEN
1051	!
1052	!-- Exchange ghost points for advection flags
1053	CALL exchange_horiz_int( advc_flags_1, nys, nyn, nxl, nxr, nzt, nbgp )
1054	CALL exchange_horiz_int( advc_flags_2, nys, nyn, nxl, nxr, nzt, nbgp )
1055	!
1056	!-- Set boundary flags at inflow and outflow boundary in case of
1057	!-- non-cyclic boundary conditions.
1058	IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN
1059	advc_flags_1(:,:,nxl-1) = advc_flags_1(:,:,nxl)
1060	advc_flags_2(:,:,nxl-1) = advc_flags_2(:,:,nxl)
1061	ENDIF
1062
1063	IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN
1064	advc_flags_1(:,:,nxr+1) = advc_flags_1(:,:,nxr)
1065	advc_flags_2(:,:,nxr+1) = advc_flags_2(:,:,nxr)
1066	ENDIF
1067
1068	IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN
1069	advc_flags_1(:,nyn+1,:) = advc_flags_1(:,nyn,:)
1070	advc_flags_2(:,nyn+1,:) = advc_flags_2(:,nyn,:)
1071	ENDIF
1072
1073	IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN
1074	advc_flags_1(:,nys-1,:) = advc_flags_1(:,nys,:)
1075	advc_flags_2(:,nys-1,:) = advc_flags_2(:,nys,:)
1076	ENDIF
1077
1078	ENDIF
1079
1080
1081	END SUBROUTINE ws_init_flags
1082
1083
1084	!------------------------------------------------------------------------------!
1085	! Description:
1086	! ------------
1087	!> Initialize variables used for storing statistic quantities (fluxes, variances)
1088	!------------------------------------------------------------------------------!
1089	SUBROUTINE ws_statistics
1090
1091	USE control_parameters, &
1092	ONLY: humidity, passive_scalar, ocean_mode, ws_scheme_mom, &
1093	ws_scheme_sca
1094
1095	USE kinds
1096
1097	USE bulk_cloud_model_mod, &
1098	ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
1099
1100	USE statistics, &
1101	ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
1102	sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqcs_ws_l, &
1103	sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l, &
1104	sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l
1105
1106
1107	IMPLICIT NONE
1108
1109	!
1110	!-- The arrays needed for statistical evaluation are set to to 0 at the
1111	!-- beginning of prognostic_equations.
1112	IF ( ws_scheme_mom ) THEN
1113	sums_wsus_ws_l = 0.0_wp
1114	sums_wsvs_ws_l = 0.0_wp
1115	sums_us2_ws_l = 0.0_wp
1116	sums_vs2_ws_l = 0.0_wp
1117	sums_ws2_ws_l = 0.0_wp
1118	ENDIF
1119
1120	IF ( ws_scheme_sca ) THEN
1121	sums_wspts_ws_l = 0.0_wp
1122	IF ( humidity ) sums_wsqs_ws_l = 0.0_wp
1123	IF ( passive_scalar ) sums_wsss_ws_l = 0.0_wp
1124	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
1125	sums_wsqcs_ws_l = 0.0_wp
1126	sums_wsncs_ws_l = 0.0_wp
1127	ENDIF
1128	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
1129	sums_wsqrs_ws_l = 0.0_wp
1130	sums_wsnrs_ws_l = 0.0_wp
1131	ENDIF
1132	IF ( ocean_mode ) sums_wssas_ws_l = 0.0_wp
1133
1134	ENDIF
1135
1136	END SUBROUTINE ws_statistics
1137
1138
1139	!------------------------------------------------------------------------------!
1140	! Description:
1141	! ------------
1142	!> Scalar advection - Call for grid point i,j
1143	!------------------------------------------------------------------------------!
1144	SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local, &
1145	swap_diss_y_local, swap_flux_x_local, &
1146	swap_diss_x_local, i_omp, tn )
1147
1148	USE arrays_3d, &
1149	ONLY: ddzw, drho_air, tend, u, v, w, rho_air_zw
1150
1151	USE control_parameters, &
1152	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
1153
1154	USE grid_variables, &
1155	ONLY: ddx, ddy
1156
1157	USE indices, &
1158	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
1159	nzt, advc_flags_1
1160
1161	USE kinds
1162
1163	USE pegrid
1164
1165	USE statistics, &
1166	ONLY: hom, sums_wsncs_ws_l, sums_wsnrs_ws_l, sums_wspts_ws_l, &
1167	sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l, &
1168	sums_wssas_ws_l, sums_wsss_ws_l, weight_substep
1169
1170
1171	IMPLICIT NONE
1172
1173	CHARACTER (LEN = *), INTENT(IN) :: sk_char !<
1174
1175	INTEGER(iwp) :: i !<
1176	INTEGER(iwp) :: ibit0 !<
1177	INTEGER(iwp) :: ibit1 !<
1178	INTEGER(iwp) :: ibit2 !<
1179	INTEGER(iwp) :: ibit3 !<
1180	INTEGER(iwp) :: ibit4 !<
1181	INTEGER(iwp) :: ibit5 !<
1182	INTEGER(iwp) :: ibit6 !<
1183	INTEGER(iwp) :: ibit7 !<
1184	INTEGER(iwp) :: ibit8 !<
1185	INTEGER(iwp) :: i_omp !<
1186	INTEGER(iwp) :: j !<
1187	INTEGER(iwp) :: k !<
1188	INTEGER(iwp) :: k_mm !<
1189	INTEGER(iwp) :: k_pp !<
1190	INTEGER(iwp) :: k_ppp !<
1191	INTEGER(iwp) :: tn !<
1192
1193	REAL(wp) :: diss_d !<
1194	REAL(wp) :: div !<
1195	REAL(wp) :: flux_d !<
1196	REAL(wp) :: u_comp !<
1197	REAL(wp) :: v_comp !<
1198
1199	#if defined( __nopointer )
1200	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !<
1201	#else
1202	REAL(wp), DIMENSION(:,:,:), POINTER :: sk !<
1203	#endif
1204	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
1205	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
1206	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
1207	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
1208	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
1209	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
1210
1211	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !<
1212	REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !<
1213
1214	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !<
1215	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !<
1216
1217
1218	!
1219	!-- Compute southside fluxes of the respective PE bounds.
1220	IF ( j == nys ) THEN
1221	!
1222	!-- Up to the top of the highest topography.
1223	DO k = nzb+1, nzb_max
1224
1225	ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)
1226	ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)
1227	ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
1228
1229	v_comp = v(k,j,i) - v_gtrans
1230	swap_flux_y_local(k,tn) = v_comp * ( &
1231	( 37.0_wp * ibit5 * adv_sca_5 &
1232	+ 7.0_wp * ibit4 * adv_sca_3 &
1233	+ ibit3 * adv_sca_1 &
1234	) * &
1235	( sk(k,j,i) + sk(k,j-1,i) ) &
1236	- ( 8.0_wp * ibit5 * adv_sca_5 &
1237	+ ibit4 * adv_sca_3 &
1238	) * &
1239	( sk(k,j+1,i) + sk(k,j-2,i) ) &
1240	+ ( ibit5 * adv_sca_5 &
1241	) * &
1242	( sk(k,j+2,i) + sk(k,j-3,i) ) &
1243	)
1244
1245	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1246	( 10.0_wp * ibit5 * adv_sca_5 &
1247	+ 3.0_wp * ibit4 * adv_sca_3 &
1248	+ ibit3 * adv_sca_1 &
1249	) * &
1250	( sk(k,j,i) - sk(k,j-1,i) ) &
1251	- ( 5.0_wp * ibit5 * adv_sca_5 &
1252	+ ibit4 * adv_sca_3 &
1253	) * &
1254	( sk(k,j+1,i) - sk(k,j-2,i) ) &
1255	+ ( ibit5 * adv_sca_5 &
1256	) * &
1257	( sk(k,j+2,i) - sk(k,j-3,i) ) &
1258	)
1259
1260	ENDDO
1261	!
1262	!-- Above to the top of the highest topography. No degradation necessary.
1263	DO k = nzb_max+1, nzt
1264
1265	v_comp = v(k,j,i) - v_gtrans
1266	swap_flux_y_local(k,tn) = v_comp * ( &
1267	37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &
1268	- 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
1269	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
1270	) * adv_sca_5
1271	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
1272	10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &
1273	- 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
1274	+ sk(k,j+2,i) - sk(k,j-3,i) &
1275	) * adv_sca_5
1276
1277	ENDDO
1278
1279	ENDIF
1280	!
1281	!-- Compute leftside fluxes of the respective PE bounds.
1282	IF ( i == i_omp ) THEN
1283
1284	DO k = nzb+1, nzb_max
1285
1286	ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)
1287	ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)
1288	ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
1289
1290	u_comp = u(k,j,i) - u_gtrans
1291	swap_flux_x_local(k,j,tn) = u_comp * ( &
1292	( 37.0_wp * ibit2 * adv_sca_5 &
1293	+ 7.0_wp * ibit1 * adv_sca_3 &
1294	+ ibit0 * adv_sca_1 &
1295	) * &
1296	( sk(k,j,i) + sk(k,j,i-1) ) &
1297	- ( 8.0_wp * ibit2 * adv_sca_5 &
1298	+ ibit1 * adv_sca_3 &
1299	) * &
1300	( sk(k,j,i+1) + sk(k,j,i-2) ) &
1301	+ ( ibit2 * adv_sca_5 &
1302	) * &
1303	( sk(k,j,i+2) + sk(k,j,i-3) ) &
1304	)
1305
1306	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1307	( 10.0_wp * ibit2 * adv_sca_5 &
1308	+ 3.0_wp * ibit1 * adv_sca_3 &
1309	+ ibit0 * adv_sca_1 &
1310	) * &
1311	( sk(k,j,i) - sk(k,j,i-1) ) &
1312	- ( 5.0_wp * ibit2 * adv_sca_5 &
1313	+ ibit1 * adv_sca_3 &
1314	) * &
1315	( sk(k,j,i+1) - sk(k,j,i-2) ) &
1316	+ ( ibit2 * adv_sca_5 &
1317	) * &
1318	( sk(k,j,i+2) - sk(k,j,i-3) ) &
1319	)
1320
1321	ENDDO
1322
1323	DO k = nzb_max+1, nzt
1324
1325	u_comp = u(k,j,i) - u_gtrans
1326	swap_flux_x_local(k,j,tn) = u_comp * ( &
1327	37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &
1328	- 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
1329	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
1330	) * adv_sca_5
1331
1332	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
1333	10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &
1334	- 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
1335	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
1336	) * adv_sca_5
1337
1338	ENDDO
1339
1340	ENDIF
1341
1342	flux_t(0) = 0.0_wp
1343	diss_t(0) = 0.0_wp
1344	!
1345	!-- Now compute the fluxes and tendency terms for the horizontal and
1346	!-- vertical parts up to the top of the highest topography.
1347	DO k = nzb+1, nzb_max
1348	!
1349	!-- Note: It is faster to conduct all multiplications explicitly, e.g.
1350	!-- * adv_sca_5 ... than to determine a factor and multiplicate the
1351	!-- flux at the end.
1352
1353	ibit2 = IBITS(advc_flags_1(k,j,i),2,1)
1354	ibit1 = IBITS(advc_flags_1(k,j,i),1,1)
1355	ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
1356
1357	u_comp = u(k,j,i+1) - u_gtrans
1358	flux_r(k) = u_comp * ( &
1359	( 37.0_wp * ibit2 * adv_sca_5 &
1360	+ 7.0_wp * ibit1 * adv_sca_3 &
1361	+ ibit0 * adv_sca_1 &
1362	) * &
1363	( sk(k,j,i+1) + sk(k,j,i) ) &
1364	- ( 8.0_wp * ibit2 * adv_sca_5 &
1365	+ ibit1 * adv_sca_3 &
1366	) * &
1367	( sk(k,j,i+2) + sk(k,j,i-1) ) &
1368	+ ( ibit2 * adv_sca_5 &
1369	) * &
1370	( sk(k,j,i+3) + sk(k,j,i-2) ) &
1371	)
1372
1373	diss_r(k) = -ABS( u_comp ) * ( &
1374	( 10.0_wp * ibit2 * adv_sca_5 &
1375	+ 3.0_wp * ibit1 * adv_sca_3 &
1376	+ ibit0 * adv_sca_1 &
1377	) * &
1378	( sk(k,j,i+1) - sk(k,j,i) ) &
1379	- ( 5.0_wp * ibit2 * adv_sca_5 &
1380	+ ibit1 * adv_sca_3 &
1381	) * &
1382	( sk(k,j,i+2) - sk(k,j,i-1) ) &
1383	+ ( ibit2 * adv_sca_5 &
1384	) * &
1385	( sk(k,j,i+3) - sk(k,j,i-2) ) &
1386	)
1387
1388	ibit5 = IBITS(advc_flags_1(k,j,i),5,1)
1389	ibit4 = IBITS(advc_flags_1(k,j,i),4,1)
1390	ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
1391
1392	v_comp = v(k,j+1,i) - v_gtrans
1393	flux_n(k) = v_comp * ( &
1394	( 37.0_wp * ibit5 * adv_sca_5 &
1395	+ 7.0_wp * ibit4 * adv_sca_3 &
1396	+ ibit3 * adv_sca_1 &
1397	) * &
1398	( sk(k,j+1,i) + sk(k,j,i) ) &
1399	- ( 8.0_wp * ibit5 * adv_sca_5 &
1400	+ ibit4 * adv_sca_3 &
1401	) * &
1402	( sk(k,j+2,i) + sk(k,j-1,i) ) &
1403	+ ( ibit5 * adv_sca_5 &
1404	) * &
1405	( sk(k,j+3,i) + sk(k,j-2,i) ) &
1406	)
1407
1408	diss_n(k) = -ABS( v_comp ) * ( &
1409	( 10.0_wp * ibit5 * adv_sca_5 &
1410	+ 3.0_wp * ibit4 * adv_sca_3 &
1411	+ ibit3 * adv_sca_1 &
1412	) * &
1413	( sk(k,j+1,i) - sk(k,j,i) ) &
1414	- ( 5.0_wp * ibit5 * adv_sca_5 &
1415	+ ibit4 * adv_sca_3 &
1416	) * &
1417	( sk(k,j+2,i) - sk(k,j-1,i) ) &
1418	+ ( ibit5 * adv_sca_5 &
1419	) * &
1420	( sk(k,j+3,i) - sk(k,j-2,i) ) &
1421	)
1422	!
1423	!-- k index has to be modified near bottom and top, else array
1424	!-- subscripts will be exceeded.
1425	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
1426	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
1427	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
1428
1429	k_ppp = k + 3 * ibit8
1430	k_pp = k + 2 * ( 1 - ibit6 )
1431	k_mm = k - 2 * ibit8
1432
1433
1434	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1435	( 37.0_wp * ibit8 * adv_sca_5 &
1436	+ 7.0_wp * ibit7 * adv_sca_3 &
1437	+ ibit6 * adv_sca_1 &
1438	) * &
1439	( sk(k+1,j,i) + sk(k,j,i) ) &
1440	- ( 8.0_wp * ibit8 * adv_sca_5 &
1441	+ ibit7 * adv_sca_3 &
1442	) * &
1443	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1444	+ ( ibit8 * adv_sca_5 &
1445	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1446	)
1447
1448	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1449	( 10.0_wp * ibit8 * adv_sca_5 &
1450	+ 3.0_wp * ibit7 * adv_sca_3 &
1451	+ ibit6 * adv_sca_1 &
1452	) * &
1453	( sk(k+1,j,i) - sk(k,j,i) ) &
1454	- ( 5.0_wp * ibit8 * adv_sca_5 &
1455	+ ibit7 * adv_sca_3 &
1456	) * &
1457	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1458	+ ( ibit8 * adv_sca_5 &
1459	) * &
1460	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1461	)
1462	ENDDO
1463
1464	DO k = nzb+1, nzb_max
1465
1466	flux_d = flux_t(k-1)
1467	diss_d = diss_t(k-1)
1468	!
1469	!-- Calculate the divergence of the velocity field. A respective
1470	!-- correction is needed to overcome numerical instabilities caused
1471	!-- by a not sufficient reduction of divergences near topography.
1472	div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &
1473	- u(k,j,i) * ( IBITS(advc_flags_1(k,j,i-1),0,1) &
1474	+ IBITS(advc_flags_1(k,j,i-1),1,1) &
1475	+ IBITS(advc_flags_1(k,j,i-1),2,1) &
1476	) &
1477	) * ddx &
1478	+ ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &
1479	- v(k,j,i) * ( IBITS(advc_flags_1(k,j-1,i),3,1) &
1480	+ IBITS(advc_flags_1(k,j-1,i),4,1) &
1481	+ IBITS(advc_flags_1(k,j-1,i),5,1) &
1482	) &
1483	) * ddy &
1484	+ ( w(k,j,i) * rho_air_zw(k) * &
1485	( ibit6 + ibit7 + ibit8 ) &
1486	- w(k-1,j,i) * rho_air_zw(k-1) * &
1487	( IBITS(advc_flags_1(k-1,j,i),6,1) &
1488	+ IBITS(advc_flags_1(k-1,j,i),7,1) &
1489	+ IBITS(advc_flags_1(k-1,j,i),8,1) &
1490	) &
1491	) * drho_air(k) * ddzw(k)
1492
1493
1494	tend(k,j,i) = tend(k,j,i) - ( &
1495	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1496	swap_diss_x_local(k,j,tn) ) * ddx &
1497	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1498	swap_diss_y_local(k,tn) ) * ddy &
1499	+ ( ( flux_t(k) + diss_t(k) ) - &
1500	( flux_d + diss_d ) &
1501	) * drho_air(k) * ddzw(k) &
1502	) + sk(k,j,i) * div
1503
1504	swap_flux_y_local(k,tn) = flux_n(k)
1505	swap_diss_y_local(k,tn) = diss_n(k)
1506	swap_flux_x_local(k,j,tn) = flux_r(k)
1507	swap_diss_x_local(k,j,tn) = diss_r(k)
1508
1509	ENDDO
1510	!
1511	!-- Now compute the fluxes and tendency terms for the horizontal and
1512	!-- vertical parts above the top of the highest topography. No degradation
1513	!-- for the horizontal parts, but for the vertical it is stell needed.
1514	DO k = nzb_max+1, nzt
1515
1516	u_comp = u(k,j,i+1) - u_gtrans
1517	flux_r(k) = u_comp * ( &
1518	37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) &
1519	- 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
1520	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
1521	diss_r(k) = -ABS( u_comp ) * ( &
1522	10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) &
1523	- 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
1524	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
1525
1526	v_comp = v(k,j+1,i) - v_gtrans
1527	flux_n(k) = v_comp * ( &
1528	37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) &
1529	- 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
1530	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
1531	diss_n(k) = -ABS( v_comp ) * ( &
1532	10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) &
1533	- 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
1534	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
1535	!
1536	!-- k index has to be modified near bottom and top, else array
1537	!-- subscripts will be exceeded.
1538	ibit8 = IBITS(advc_flags_1(k,j,i),8,1)
1539	ibit7 = IBITS(advc_flags_1(k,j,i),7,1)
1540	ibit6 = IBITS(advc_flags_1(k,j,i),6,1)
1541
1542	k_ppp = k + 3 * ibit8
1543	k_pp = k + 2 * ( 1 - ibit6 )
1544	k_mm = k - 2 * ibit8
1545
1546	flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &
1547	( 37.0_wp * ibit8 * adv_sca_5 &
1548	+ 7.0_wp * ibit7 * adv_sca_3 &
1549	+ ibit6 * adv_sca_1 &
1550	) * &
1551	( sk(k+1,j,i) + sk(k,j,i) ) &
1552	- ( 8.0_wp * ibit8 * adv_sca_5 &
1553	+ ibit7 * adv_sca_3 &
1554	) * &
1555	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
1556	+ ( ibit8 * adv_sca_5 &
1557	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
1558	)
1559
1560	diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &
1561	( 10.0_wp * ibit8 * adv_sca_5 &
1562	+ 3.0_wp * ibit7 * adv_sca_3 &
1563	+ ibit6 * adv_sca_1 &
1564	) * &
1565	( sk(k+1,j,i) - sk(k,j,i) ) &
1566	- ( 5.0_wp * ibit8 * adv_sca_5 &
1567	+ ibit7 * adv_sca_3 &
1568	) * &
1569	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
1570	+ ( ibit8 * adv_sca_5 &
1571	) * &
1572	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
1573	)
1574
1575	ENDDO
1576
1577	DO k = nzb_max+1, nzt
1578
1579	flux_d = flux_t(k-1)
1580	diss_d = diss_t(k-1)
1581
1582	!
1583	!-- Calculate the divergence of the velocity field. A respective
1584	!-- correction is needed to overcome numerical instabilities introduced
1585	!-- by a not sufficient reduction of divergences near topography.
1586	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
1587	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
1588	+ ( w(k,j,i) * rho_air_zw(k) - &
1589	w(k-1,j,i) * rho_air_zw(k-1) &
1590	) * drho_air(k) * ddzw(k)
1591
1592	tend(k,j,i) = tend(k,j,i) - ( &
1593	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
1594	swap_diss_x_local(k,j,tn) ) * ddx &
1595	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
1596	swap_diss_y_local(k,tn) ) * ddy &
1597	+ ( ( flux_t(k) + diss_t(k) ) - &
1598	( flux_d + diss_d ) &
1599	) * drho_air(k) * ddzw(k) &
1600	) + sk(k,j,i) * div
1601
1602
1603	swap_flux_y_local(k,tn) = flux_n(k)
1604	swap_diss_y_local(k,tn) = diss_n(k)
1605	swap_flux_x_local(k,j,tn) = flux_r(k)
1606	swap_diss_x_local(k,j,tn) = diss_r(k)
1607
1608	ENDDO
1609
1610	!
1611	!-- Evaluation of statistics.
1612	SELECT CASE ( sk_char )
1613
1614	CASE ( 'pt' )
1615
1616	DO k = nzb, nzt
1617	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + &
1618	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1619	* ( w(k,j,i) - hom(k,1,3,0) ) &
1620	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1621	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1622	) * weight_substep(intermediate_timestep_count)
1623	ENDDO
1624
1625	CASE ( 'sa' )
1626
1627	DO k = nzb, nzt
1628	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + &
1629	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1630	* ( w(k,j,i) - hom(k,1,3,0) ) &
1631	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1632	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1633	) * weight_substep(intermediate_timestep_count)
1634	ENDDO
1635
1636	CASE ( 'q' )
1637
1638	DO k = nzb, nzt
1639	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + &
1640	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1641	* ( w(k,j,i) - hom(k,1,3,0) ) &
1642	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1643	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1644	) * weight_substep(intermediate_timestep_count)
1645	ENDDO
1646
1647	CASE ( 'qc' )
1648
1649	DO k = nzb, nzt
1650	sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) + &
1651	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1652	* ( w(k,j,i) - hom(k,1,3,0) ) &
1653	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1654	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1655	) * weight_substep(intermediate_timestep_count)
1656	ENDDO
1657
1658
1659	CASE ( 'qr' )
1660
1661	DO k = nzb, nzt
1662	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + &
1663	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1664	* ( w(k,j,i) - hom(k,1,3,0) ) &
1665	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1666	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1667	) * weight_substep(intermediate_timestep_count)
1668	ENDDO
1669
1670	CASE ( 'nc' )
1671
1672	DO k = nzb, nzt
1673	sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) + &
1674	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1675	* ( w(k,j,i) - hom(k,1,3,0) ) &
1676	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1677	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1678	) * weight_substep(intermediate_timestep_count)
1679	ENDDO
1680
1681	CASE ( 'nr' )
1682
1683	DO k = nzb, nzt
1684	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + &
1685	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1686	* ( w(k,j,i) - hom(k,1,3,0) ) &
1687	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1688	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1689	) * weight_substep(intermediate_timestep_count)
1690	ENDDO
1691
1692	CASE ( 's' )
1693
1694	DO k = nzb, nzt
1695	sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) + &
1696	( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) &
1697	* ( w(k,j,i) - hom(k,1,3,0) ) &
1698	+ diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) &
1699	* ABS( w(k,j,i) - hom(k,1,3,0) ) &
1700	) * weight_substep(intermediate_timestep_count)
1701	ENDDO
1702
1703	! CASE ( 'kc' )
1704	!kk Has to be implemented for kpp chemistry
1705
1706
1707	END SELECT
1708
1709	END SUBROUTINE advec_s_ws_ij
1710
1711
1712
1713
1714	!------------------------------------------------------------------------------!
1715	! Description:
1716	! ------------
1717	!> Advection of u-component - Call for grid point i,j
1718	!------------------------------------------------------------------------------!
1719	SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
1720
1721	USE arrays_3d, &
1722	ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,&
1723	drho_air, rho_air_zw
1724
1725	USE control_parameters, &
1726	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
1727
1728	USE grid_variables, &
1729	ONLY: ddx, ddy
1730
1731	USE indices, &
1732	ONLY: nxlu, nys, nzb, nzb_max, nzt, advc_flags_1
1733
1734	USE kinds
1735
1736	USE statistics, &
1737	ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep
1738
1739	IMPLICIT NONE
1740
1741	INTEGER(iwp) :: i !<
1742	INTEGER(iwp) :: ibit9 !<
1743	INTEGER(iwp) :: ibit10 !<
1744	INTEGER(iwp) :: ibit11 !<
1745	INTEGER(iwp) :: ibit12 !<
1746	INTEGER(iwp) :: ibit13 !<
1747	INTEGER(iwp) :: ibit14 !<
1748	INTEGER(iwp) :: ibit15 !<
1749	INTEGER(iwp) :: ibit16 !<
1750	INTEGER(iwp) :: ibit17 !<
1751	INTEGER(iwp) :: i_omp !<
1752	INTEGER(iwp) :: j !<
1753	INTEGER(iwp) :: k !<
1754	INTEGER(iwp) :: k_mm !<
1755	INTEGER(iwp) :: k_pp !<
1756	INTEGER(iwp) :: k_ppp !<
1757	INTEGER(iwp) :: tn !<
1758
1759	REAL(wp) :: diss_d !<
1760	REAL(wp) :: div !<
1761	REAL(wp) :: flux_d !<
1762	REAL(wp) :: gu !<
1763	REAL(wp) :: gv !<
1764	REAL(wp) :: u_comp_l !<
1765
1766	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
1767	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
1768	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
1769	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
1770	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
1771	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
1772	REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !<
1773	REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !<
1774	REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !<
1775
1776	gu = 2.0_wp * u_gtrans
1777	gv = 2.0_wp * v_gtrans
1778	!
1779	!-- Compute southside fluxes for the respective boundary of PE
1780	IF ( j == nys ) THEN
1781
1782	DO k = nzb+1, nzb_max
1783
1784	ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)
1785	ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)
1786	ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)
1787
1788	v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv
1789	flux_s_u(k,tn) = v_comp(k) * ( &
1790	( 37.0_wp * ibit14 * adv_mom_5 &
1791	+ 7.0_wp * ibit13 * adv_mom_3 &
1792	+ ibit12 * adv_mom_1 &
1793	) * &
1794	( u(k,j,i) + u(k,j-1,i) ) &
1795	- ( 8.0_wp * ibit14 * adv_mom_5 &
1796	+ ibit13 * adv_mom_3 &
1797	) * &
1798	( u(k,j+1,i) + u(k,j-2,i) ) &
1799	+ ( ibit14 * adv_mom_5 &
1800	) * &
1801	( u(k,j+2,i) + u(k,j-3,i) ) &
1802	)
1803
1804	diss_s_u(k,tn) = - ABS ( v_comp(k) ) * ( &
1805	( 10.0_wp * ibit14 * adv_mom_5 &
1806	+ 3.0_wp * ibit13 * adv_mom_3 &
1807	+ ibit12 * adv_mom_1 &
1808	) * &
1809	( u(k,j,i) - u(k,j-1,i) ) &
1810	- ( 5.0_wp * ibit14 * adv_mom_5 &
1811	+ ibit13 * adv_mom_3 &
1812	) * &
1813	( u(k,j+1,i) - u(k,j-2,i) ) &
1814	+ ( ibit14 * adv_mom_5 &
1815	) * &
1816	( u(k,j+2,i) - u(k,j-3,i) ) &
1817	)
1818
1819	ENDDO
1820
1821	DO k = nzb_max+1, nzt
1822
1823	v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv
1824	flux_s_u(k,tn) = v_comp(k) * ( &
1825	37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) &
1826	- 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) &
1827	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
1828	diss_s_u(k,tn) = - ABS(v_comp(k)) * ( &
1829	10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) &
1830	- 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) &
1831	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
1832
1833	ENDDO
1834
1835	ENDIF
1836	!
1837	!-- Compute leftside fluxes for the respective boundary of PE
1838	IF ( i == i_omp .OR. i == nxlu ) THEN
1839
1840	DO k = nzb+1, nzb_max
1841
1842	ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)
1843	ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)
1844	ibit9 = IBITS(advc_flags_1(k,j,i-1),9,1)
1845
1846	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
1847	flux_l_u(k,j,tn) = u_comp_l * ( &
1848	( 37.0_wp * ibit11 * adv_mom_5 &
1849	+ 7.0_wp * ibit10 * adv_mom_3 &
1850	+ ibit9 * adv_mom_1 &
1851	) * &
1852	( u(k,j,i) + u(k,j,i-1) ) &
1853	- ( 8.0_wp * ibit11 * adv_mom_5 &
1854	+ ibit10 * adv_mom_3 &
1855	) * &
1856	( u(k,j,i+1) + u(k,j,i-2) ) &
1857	+ ( ibit11 * adv_mom_5 &
1858	) * &
1859	( u(k,j,i+2) + u(k,j,i-3) ) &
1860	)
1861
1862	diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( &
1863	( 10.0_wp * ibit11 * adv_mom_5 &
1864	+ 3.0_wp * ibit10 * adv_mom_3 &
1865	+ ibit9 * adv_mom_1 &
1866	) * &
1867	( u(k,j,i) - u(k,j,i-1) ) &
1868	- ( 5.0_wp * ibit11 * adv_mom_5 &
1869	+ ibit10 * adv_mom_3 &
1870	) * &
1871	( u(k,j,i+1) - u(k,j,i-2) ) &
1872	+ ( ibit11 * adv_mom_5 &
1873	) * &
1874	( u(k,j,i+2) - u(k,j,i-3) ) &
1875	)
1876
1877	ENDDO
1878
1879	DO k = nzb_max+1, nzt
1880
1881	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
1882	flux_l_u(k,j,tn) = u_comp_l * ( &
1883	37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) &
1884	- 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) &
1885	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
1886	diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( &
1887	10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) &
1888	- 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) &
1889	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
1890
1891	ENDDO
1892
1893	ENDIF
1894
1895	flux_t(0) = 0.0_wp
1896	diss_t(0) = 0.0_wp
1897	w_comp(0) = 0.0_wp
1898	!
1899	!-- Now compute the fluxes tendency terms for the horizontal and
1900	!-- vertical parts.
1901	DO k = nzb+1, nzb_max
1902
1903	ibit11 = IBITS(advc_flags_1(k,j,i),11,1)
1904	ibit10 = IBITS(advc_flags_1(k,j,i),10,1)
1905	ibit9 = IBITS(advc_flags_1(k,j,i),9,1)
1906
1907	u_comp(k) = u(k,j,i+1) + u(k,j,i)
1908	flux_r(k) = ( u_comp(k) - gu ) * ( &
1909	( 37.0_wp * ibit11 * adv_mom_5 &
1910	+ 7.0_wp * ibit10 * adv_mom_3 &
1911	+ ibit9 * adv_mom_1 &
1912	) * &
1913	( u(k,j,i+1) + u(k,j,i) ) &
1914	- ( 8.0_wp * ibit11 * adv_mom_5 &
1915	+ ibit10 * adv_mom_3 &
1916	) * &
1917	( u(k,j,i+2) + u(k,j,i-1) ) &
1918	+ ( ibit11 * adv_mom_5 &
1919	) * &
1920	( u(k,j,i+3) + u(k,j,i-2) ) &
1921	)
1922
1923	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
1924	( 10.0_wp * ibit11 * adv_mom_5 &
1925	+ 3.0_wp * ibit10 * adv_mom_3 &
1926	+ ibit9 * adv_mom_1 &
1927	) * &
1928	( u(k,j,i+1) - u(k,j,i) ) &
1929	- ( 5.0_wp * ibit11 * adv_mom_5 &
1930	+ ibit10 * adv_mom_3 &
1931	) * &
1932	( u(k,j,i+2) - u(k,j,i-1) ) &
1933	+ ( ibit11 * adv_mom_5 &
1934	) * &
1935	( u(k,j,i+3) - u(k,j,i-2) ) &
1936	)
1937
1938	ibit14 = IBITS(advc_flags_1(k,j,i),14,1)
1939	ibit13 = IBITS(advc_flags_1(k,j,i),13,1)
1940	ibit12 = IBITS(advc_flags_1(k,j,i),12,1)
1941
1942	v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
1943	flux_n(k) = v_comp(k) * ( &
1944	( 37.0_wp * ibit14 * adv_mom_5 &
1945	+ 7.0_wp * ibit13 * adv_mom_3 &
1946	+ ibit12 * adv_mom_1 &
1947	) * &
1948	( u(k,j+1,i) + u(k,j,i) ) &
1949	- ( 8.0_wp * ibit14 * adv_mom_5 &
1950	+ ibit13 * adv_mom_3 &
1951	) * &
1952	( u(k,j+2,i) + u(k,j-1,i) ) &
1953	+ ( ibit14 * adv_mom_5 &
1954	) * &
1955	( u(k,j+3,i) + u(k,j-2,i) ) &
1956	)
1957
1958	diss_n(k) = - ABS ( v_comp(k) ) * ( &
1959	( 10.0_wp * ibit14 * adv_mom_5 &
1960	+ 3.0_wp * ibit13 * adv_mom_3 &
1961	+ ibit12 * adv_mom_1 &
1962	) * &
1963	( u(k,j+1,i) - u(k,j,i) ) &
1964	- ( 5.0_wp * ibit14 * adv_mom_5 &
1965	+ ibit13 * adv_mom_3 &
1966	) * &
1967	( u(k,j+2,i) - u(k,j-1,i) ) &
1968	+ ( ibit14 * adv_mom_5 &
1969	) * &
1970	( u(k,j+3,i) - u(k,j-2,i) ) &
1971	)
1972	!
1973	!-- k index has to be modified near bottom and top, else array
1974	!-- subscripts will be exceeded.
1975	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
1976	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
1977	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
1978
1979	k_ppp = k + 3 * ibit17
1980	k_pp = k + 2 * ( 1 - ibit15 )
1981	k_mm = k - 2 * ibit17
1982
1983	w_comp(k) = w(k,j,i) + w(k,j,i-1)
1984	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
1985	( 37.0_wp * ibit17 * adv_mom_5 &
1986	+ 7.0_wp * ibit16 * adv_mom_3 &
1987	+ ibit15 * adv_mom_1 &
1988	) * &
1989	( u(k+1,j,i) + u(k,j,i) ) &
1990	- ( 8.0_wp * ibit17 * adv_mom_5 &
1991	+ ibit16 * adv_mom_3 &
1992	) * &
1993	( u(k_pp,j,i) + u(k-1,j,i) ) &
1994	+ ( ibit17 * adv_mom_5 &
1995	) * &
1996	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
1997	)
1998
1999	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2000	( 10.0_wp * ibit17 * adv_mom_5 &
2001	+ 3.0_wp * ibit16 * adv_mom_3 &
2002	+ ibit15 * adv_mom_1 &
2003	) * &
2004	( u(k+1,j,i) - u(k,j,i) ) &
2005	- ( 5.0_wp * ibit17 * adv_mom_5 &
2006	+ ibit16 * adv_mom_3 &
2007	) * &
2008	( u(k_pp,j,i) - u(k-1,j,i) ) &
2009	+ ( ibit17 * adv_mom_5 &
2010	) * &
2011	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2012	)
2013	ENDDO
2014
2015	DO k = nzb+1, nzb_max
2016
2017	flux_d = flux_t(k-1)
2018	diss_d = diss_t(k-1)
2019	!
2020	!-- Calculate the divergence of the velocity field. A respective
2021	!-- correction is needed to overcome numerical instabilities introduced
2022	!-- by a not sufficient reduction of divergences near topography.
2023	div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) &
2024	- ( u(k,j,i) + u(k,j,i-1) ) &
2025	* ( IBITS(advc_flags_1(k,j,i-1),9,1) &
2026	+ IBITS(advc_flags_1(k,j,i-1),10,1) &
2027	+ IBITS(advc_flags_1(k,j,i-1),11,1) &
2028	) &
2029	) * ddx &
2030	+ ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 ) &
2031	- ( v(k,j,i) + v(k,j,i-1 ) ) &
2032	* ( IBITS(advc_flags_1(k,j-1,i),12,1) &
2033	+ IBITS(advc_flags_1(k,j-1,i),13,1) &
2034	+ IBITS(advc_flags_1(k,j-1,i),14,1) &
2035	) &
2036	) * ddy &
2037	+ ( w_comp(k) * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )&
2038	- w_comp(k-1) * rho_air_zw(k-1) &
2039	* ( IBITS(advc_flags_1(k-1,j,i),15,1) &
2040	+ IBITS(advc_flags_1(k-1,j,i),16,1) &
2041	+ IBITS(advc_flags_1(k-1,j,i),17,1) &
2042	) &
2043	) * drho_air(k) * ddzw(k) &
2044	) * 0.5_wp
2045
2046
2047	tend(k,j,i) = tend(k,j,i) - ( &
2048	( flux_r(k) + diss_r(k) &
2049	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
2050	+ ( flux_n(k) + diss_n(k) &
2051	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
2052	+ ( ( flux_t(k) + diss_t(k) ) &
2053	- ( flux_d + diss_d ) &
2054	) * drho_air(k) * ddzw(k) &
2055	) + div * u(k,j,i)
2056
2057	flux_l_u(k,j,tn) = flux_r(k)
2058	diss_l_u(k,j,tn) = diss_r(k)
2059	flux_s_u(k,tn) = flux_n(k)
2060	diss_s_u(k,tn) = diss_n(k)
2061	!
2062	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2063	!-- right evaluation when gallilei_trans = .T. .
2064	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2065	+ ( flux_r(k) &
2066	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2067	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2068	+ diss_r(k) &
2069	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2070	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2071	) * weight_substep(intermediate_timestep_count)
2072	!
2073	!-- Statistical Evaluation of w'u'.
2074	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2075	+ ( flux_t(k) &
2076	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2077	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
2078	+ diss_t(k) &
2079	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2080	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
2081	) * weight_substep(intermediate_timestep_count)
2082	ENDDO
2083
2084	DO k = nzb_max+1, nzt
2085
2086	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2087	flux_r(k) = ( u_comp(k) - gu ) * ( &
2088	37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &
2089	- 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &
2090	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2091	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2092	10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &
2093	- 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &
2094	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2095
2096	v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv
2097	flux_n(k) = v_comp(k) * ( &
2098	37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &
2099	- 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &
2100	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2101	diss_n(k) = - ABS( v_comp(k) ) * ( &
2102	10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &
2103	- 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &
2104	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2105	!
2106	!-- k index has to be modified near bottom and top, else array
2107	!-- subscripts will be exceeded.
2108	ibit17 = IBITS(advc_flags_1(k,j,i),17,1)
2109	ibit16 = IBITS(advc_flags_1(k,j,i),16,1)
2110	ibit15 = IBITS(advc_flags_1(k,j,i),15,1)
2111
2112	k_ppp = k + 3 * ibit17
2113	k_pp = k + 2 * ( 1 - ibit15 )
2114	k_mm = k - 2 * ibit17
2115
2116	w_comp(k) = w(k,j,i) + w(k,j,i-1)
2117	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2118	( 37.0_wp * ibit17 * adv_mom_5 &
2119	+ 7.0_wp * ibit16 * adv_mom_3 &
2120	+ ibit15 * adv_mom_1 &
2121	) * &
2122	( u(k+1,j,i) + u(k,j,i) ) &
2123	- ( 8.0_wp * ibit17 * adv_mom_5 &
2124	+ ibit16 * adv_mom_3 &
2125	) * &
2126	( u(k_pp,j,i) + u(k-1,j,i) ) &
2127	+ ( ibit17 * adv_mom_5 &
2128	) * &
2129	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2130	)
2131
2132	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2133	( 10.0_wp * ibit17 * adv_mom_5 &
2134	+ 3.0_wp * ibit16 * adv_mom_3 &
2135	+ ibit15 * adv_mom_1 &
2136	) * &
2137	( u(k+1,j,i) - u(k,j,i) ) &
2138	- ( 5.0_wp * ibit17 * adv_mom_5 &
2139	+ ibit16 * adv_mom_3 &
2140	) * &
2141	( u(k_pp,j,i) - u(k-1,j,i) ) &
2142	+ ( ibit17 * adv_mom_5 &
2143	) * &
2144	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2145	)
2146
2147	ENDDO
2148
2149	DO k = nzb_max+1, nzt
2150
2151	flux_d = flux_t(k-1)
2152	diss_d = diss_t(k-1)
2153	!
2154	!-- Calculate the divergence of the velocity field. A respective
2155	!-- correction is needed to overcome numerical instabilities introduced
2156	!-- by a not sufficient reduction of divergences near topography.
2157	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
2158	+ ( v_comp(k) + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
2159	+ ( w_comp(k) * rho_air_zw(k) - &
2160	w_comp(k-1) * rho_air_zw(k-1) &
2161	) * drho_air(k) * ddzw(k) &
2162	) * 0.5_wp
2163
2164	tend(k,j,i) = tend(k,j,i) - ( &
2165	( flux_r(k) + diss_r(k) &
2166	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
2167	+ ( flux_n(k) + diss_n(k) &
2168	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
2169	+ ( ( flux_t(k) + diss_t(k) ) &
2170	- ( flux_d + diss_d ) &
2171	) * drho_air(k) * ddzw(k) &
2172	) + div * u(k,j,i)
2173
2174	flux_l_u(k,j,tn) = flux_r(k)
2175	diss_l_u(k,j,tn) = diss_r(k)
2176	flux_s_u(k,tn) = flux_n(k)
2177	diss_s_u(k,tn) = diss_n(k)
2178	!
2179	!-- Statistical Evaluation of u'u'. The factor has to be applied for
2180	!-- right evaluation when gallilei_trans = .T. .
2181	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2182	+ ( flux_r(k) &
2183	* ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2184	/ ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &
2185	+ diss_r(k) &
2186	* ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &
2187	/ ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &
2188	) * weight_substep(intermediate_timestep_count)
2189	!
2190	!-- Statistical Evaluation of w'u'.
2191	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2192	+ ( flux_t(k) &
2193	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2194	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
2195	+ diss_t(k) &
2196	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2197	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
2198	) * weight_substep(intermediate_timestep_count)
2199	ENDDO
2200
2201	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
2202
2203
2204
2205	END SUBROUTINE advec_u_ws_ij
2206
2207
2208
2209	!-----------------------------------------------------------------------------!
2210	! Description:
2211	! ------------
2212	!> Advection of v-component - Call for grid point i,j
2213	!-----------------------------------------------------------------------------!
2214	SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
2215
2216	USE arrays_3d, &
2217	ONLY: ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w, &
2218	drho_air, rho_air_zw
2219
2220	USE control_parameters, &
2221	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
2222
2223	USE grid_variables, &
2224	ONLY: ddx, ddy
2225
2226	USE indices, &
2227	ONLY: nysv, nzb, nzb_max, nzt, advc_flags_1
2228
2229	USE kinds
2230
2231	USE statistics, &
2232	ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep
2233
2234	IMPLICIT NONE
2235
2236	INTEGER(iwp) :: i !<
2237	INTEGER(iwp) :: ibit18 !<
2238	INTEGER(iwp) :: ibit19 !<
2239	INTEGER(iwp) :: ibit20 !<
2240	INTEGER(iwp) :: ibit21 !<
2241	INTEGER(iwp) :: ibit22 !<
2242	INTEGER(iwp) :: ibit23 !<
2243	INTEGER(iwp) :: ibit24 !<
2244	INTEGER(iwp) :: ibit25 !<
2245	INTEGER(iwp) :: ibit26 !<
2246	INTEGER(iwp) :: i_omp !<
2247	INTEGER(iwp) :: j !<
2248	INTEGER(iwp) :: k !<
2249	INTEGER(iwp) :: k_mm !<
2250	INTEGER(iwp) :: k_pp !<
2251	INTEGER(iwp) :: k_ppp !<
2252	INTEGER(iwp) :: tn !<
2253
2254	REAL(wp) :: diss_d !<
2255	REAL(wp) :: div !<
2256	REAL(wp) :: flux_d !<
2257	REAL(wp) :: gu !<
2258	REAL(wp) :: gv !<
2259	REAL(wp) :: v_comp_l !<
2260
2261	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
2262	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
2263	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
2264	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
2265	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
2266	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
2267	REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !<
2268	REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !<
2269	REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !<
2270
2271	gu = 2.0_wp * u_gtrans
2272	gv = 2.0_wp * v_gtrans
2273
2274	!
2275	!-- Compute leftside fluxes for the respective boundary.
2276	IF ( i == i_omp ) THEN
2277
2278	DO k = nzb+1, nzb_max
2279
2280	ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)
2281	ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)
2282	ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)
2283
2284	u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu
2285	flux_l_v(k,j,tn) = u_comp(k) * ( &
2286	( 37.0_wp * ibit20 * adv_mom_5 &
2287	+ 7.0_wp * ibit19 * adv_mom_3 &
2288	+ ibit18 * adv_mom_1 &
2289	) * &
2290	( v(k,j,i) + v(k,j,i-1) ) &
2291	- ( 8.0_wp * ibit20 * adv_mom_5 &
2292	+ ibit19 * adv_mom_3 &
2293	) * &
2294	( v(k,j,i+1) + v(k,j,i-2) ) &
2295	+ ( ibit20 * adv_mom_5 &
2296	) * &
2297	( v(k,j,i+2) + v(k,j,i-3) ) &
2298	)
2299
2300	diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( &
2301	( 10.0_wp * ibit20 * adv_mom_5 &
2302	+ 3.0_wp * ibit19 * adv_mom_3 &
2303	+ ibit18 * adv_mom_1 &
2304	) * &
2305	( v(k,j,i) - v(k,j,i-1) ) &
2306	- ( 5.0_wp * ibit20 * adv_mom_5 &
2307	+ ibit19 * adv_mom_3 &
2308	) * &
2309	( v(k,j,i+1) - v(k,j,i-2) ) &
2310	+ ( ibit20 * adv_mom_5 &
2311	) * &
2312	( v(k,j,i+2) - v(k,j,i-3) ) &
2313	)
2314
2315	ENDDO
2316
2317	DO k = nzb_max+1, nzt
2318
2319	u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu
2320	flux_l_v(k,j,tn) = u_comp(k) * ( &
2321	37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) &
2322	- 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) &
2323	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
2324	diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( &
2325	10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) &
2326	- 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) &
2327	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
2328
2329	ENDDO
2330
2331	ENDIF
2332	!
2333	!-- Compute southside fluxes for the respective boundary.
2334	IF ( j == nysv ) THEN
2335
2336	DO k = nzb+1, nzb_max
2337
2338	ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)
2339	ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)
2340	ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)
2341
2342	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2343	flux_s_v(k,tn) = v_comp_l * ( &
2344	( 37.0_wp * ibit23 * adv_mom_5 &
2345	+ 7.0_wp * ibit22 * adv_mom_3 &
2346	+ ibit21 * adv_mom_1 &
2347	) * &
2348	( v(k,j,i) + v(k,j-1,i) ) &
2349	- ( 8.0_wp * ibit23 * adv_mom_5 &
2350	+ ibit22 * adv_mom_3 &
2351	) * &
2352	( v(k,j+1,i) + v(k,j-2,i) ) &
2353	+ ( ibit23 * adv_mom_5 &
2354	) * &
2355	( v(k,j+2,i) + v(k,j-3,i) ) &
2356	)
2357
2358	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2359	( 10.0_wp * ibit23 * adv_mom_5 &
2360	+ 3.0_wp * ibit22 * adv_mom_3 &
2361	+ ibit21 * adv_mom_1 &
2362	) * &
2363	( v(k,j,i) - v(k,j-1,i) ) &
2364	- ( 5.0_wp * ibit23 * adv_mom_5 &
2365	+ ibit22 * adv_mom_3 &
2366	) * &
2367	( v(k,j+1,i) - v(k,j-2,i) ) &
2368	+ ( ibit23 * adv_mom_5 &
2369	) * &
2370	( v(k,j+2,i) - v(k,j-3,i) ) &
2371	)
2372
2373	ENDDO
2374
2375	DO k = nzb_max+1, nzt
2376
2377	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
2378	flux_s_v(k,tn) = v_comp_l * ( &
2379	37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) &
2380	- 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) &
2381	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
2382	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
2383	10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) &
2384	- 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) &
2385	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
2386
2387	ENDDO
2388
2389	ENDIF
2390
2391	flux_t(0) = 0.0_wp
2392	diss_t(0) = 0.0_wp
2393	w_comp(0) = 0.0_wp
2394	!
2395	!-- Now compute the fluxes and tendency terms for the horizontal and
2396	!-- verical parts.
2397	DO k = nzb+1, nzb_max
2398
2399	ibit20 = IBITS(advc_flags_1(k,j,i),20,1)
2400	ibit19 = IBITS(advc_flags_1(k,j,i),19,1)
2401	ibit18 = IBITS(advc_flags_1(k,j,i),18,1)
2402
2403	u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
2404	flux_r(k) = u_comp(k) * ( &
2405	( 37.0_wp * ibit20 * adv_mom_5 &
2406	+ 7.0_wp * ibit19 * adv_mom_3 &
2407	+ ibit18 * adv_mom_1 &
2408	) * &
2409	( v(k,j,i+1) + v(k,j,i) ) &
2410	- ( 8.0_wp * ibit20 * adv_mom_5 &
2411	+ ibit19 * adv_mom_3 &
2412	) * &
2413	( v(k,j,i+2) + v(k,j,i-1) ) &
2414	+ ( ibit20 * adv_mom_5 &
2415	) * &
2416	( v(k,j,i+3) + v(k,j,i-2) ) &
2417	)
2418
2419	diss_r(k) = - ABS( u_comp(k) ) * ( &
2420	( 10.0_wp * ibit20 * adv_mom_5 &
2421	+ 3.0_wp * ibit19 * adv_mom_3 &
2422	+ ibit18 * adv_mom_1 &
2423	) * &
2424	( v(k,j,i+1) - v(k,j,i) ) &
2425	- ( 5.0_wp * ibit20 * adv_mom_5 &
2426	+ ibit19 * adv_mom_3 &
2427	) * &
2428	( v(k,j,i+2) - v(k,j,i-1) ) &
2429	+ ( ibit20 * adv_mom_5 &
2430	) * &
2431	( v(k,j,i+3) - v(k,j,i-2) ) &
2432	)
2433
2434	ibit23 = IBITS(advc_flags_1(k,j,i),23,1)
2435	ibit22 = IBITS(advc_flags_1(k,j,i),22,1)
2436	ibit21 = IBITS(advc_flags_1(k,j,i),21,1)
2437
2438
2439	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2440	flux_n(k) = ( v_comp(k) - gv ) * ( &
2441	( 37.0_wp * ibit23 * adv_mom_5 &
2442	+ 7.0_wp * ibit22 * adv_mom_3 &
2443	+ ibit21 * adv_mom_1 &
2444	) * &
2445	( v(k,j+1,i) + v(k,j,i) ) &
2446	- ( 8.0_wp * ibit23 * adv_mom_5 &
2447	+ ibit22 * adv_mom_3 &
2448	) * &
2449	( v(k,j+2,i) + v(k,j-1,i) ) &
2450	+ ( ibit23 * adv_mom_5 &
2451	) * &
2452	( v(k,j+3,i) + v(k,j-2,i) ) &
2453	)
2454
2455	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2456	( 10.0_wp * ibit23 * adv_mom_5 &
2457	+ 3.0_wp * ibit22 * adv_mom_3 &
2458	+ ibit21 * adv_mom_1 &
2459	) * &
2460	( v(k,j+1,i) - v(k,j,i) ) &
2461	- ( 5.0_wp * ibit23 * adv_mom_5 &
2462	+ ibit22 * adv_mom_3 &
2463	) * &
2464	( v(k,j+2,i) - v(k,j-1,i) ) &
2465	+ ( ibit23 * adv_mom_5 &
2466	) * &
2467	( v(k,j+3,i) - v(k,j-2,i) ) &
2468	)
2469	!
2470	!-- k index has to be modified near bottom and top, else array
2471	!-- subscripts will be exceeded.
2472	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
2473	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
2474	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
2475
2476	k_ppp = k + 3 * ibit26
2477	k_pp = k + 2 * ( 1 - ibit24 )
2478	k_mm = k - 2 * ibit26
2479
2480	w_comp(k) = w(k,j-1,i) + w(k,j,i)
2481	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2482	( 37.0_wp * ibit26 * adv_mom_5 &
2483	+ 7.0_wp * ibit25 * adv_mom_3 &
2484	+ ibit24 * adv_mom_1 &
2485	) * &
2486	( v(k+1,j,i) + v(k,j,i) ) &
2487	- ( 8.0_wp * ibit26 * adv_mom_5 &
2488	+ ibit25 * adv_mom_3 &
2489	) * &
2490	( v(k_pp,j,i) + v(k-1,j,i) ) &
2491	+ ( ibit26 * adv_mom_5 &
2492	) * &
2493	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2494	)
2495
2496	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2497	( 10.0_wp * ibit26 * adv_mom_5 &
2498	+ 3.0_wp * ibit25 * adv_mom_3 &
2499	+ ibit24 * adv_mom_1 &
2500	) * &
2501	( v(k+1,j,i) - v(k,j,i) ) &
2502	- ( 5.0_wp * ibit26 * adv_mom_5 &
2503	+ ibit25 * adv_mom_3 &
2504	) * &
2505	( v(k_pp,j,i) - v(k-1,j,i) ) &
2506	+ ( ibit26 * adv_mom_5 &
2507	) * &
2508	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2509	)
2510	ENDDO
2511
2512	DO k = nzb+1, nzb_max
2513
2514	flux_d = flux_t(k-1)
2515	diss_d = diss_t(k-1)
2516	!
2517	!-- Calculate the divergence of the velocity field. A respective
2518	!-- correction is needed to overcome numerical instabilities introduced
2519	!-- by a not sufficient reduction of divergences near topography.
2520	div = ( ( ( u_comp(k) + gu ) &
2521	* ( ibit18 + ibit19 + ibit20 ) &
2522	- ( u(k,j-1,i) + u(k,j,i) ) &
2523	* ( IBITS(advc_flags_1(k,j,i-1),18,1) &
2524	+ IBITS(advc_flags_1(k,j,i-1),19,1) &
2525	+ IBITS(advc_flags_1(k,j,i-1),20,1) &
2526	) &
2527	) * ddx &
2528	+ ( v_comp(k) &
2529	* ( ibit21 + ibit22 + ibit23 ) &
2530	- ( v(k,j,i) + v(k,j-1,i) ) &
2531	* ( IBITS(advc_flags_1(k,j-1,i),21,1) &
2532	+ IBITS(advc_flags_1(k,j-1,i),22,1) &
2533	+ IBITS(advc_flags_1(k,j-1,i),23,1) &
2534	) &
2535	) * ddy &
2536	+ ( w_comp(k) * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )&
2537	- w_comp(k-1) * rho_air_zw(k-1) &
2538	* ( IBITS(advc_flags_1(k-1,j,i),24,1) &
2539	+ IBITS(advc_flags_1(k-1,j,i),25,1) &
2540	+ IBITS(advc_flags_1(k-1,j,i),26,1) &
2541	) &
2542	) * drho_air(k) * ddzw(k) &
2543	) * 0.5_wp
2544
2545
2546	tend(k,j,i) = tend(k,j,i) - ( &
2547	( flux_r(k) + diss_r(k) &
2548	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2549	+ ( flux_n(k) + diss_n(k) &
2550	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2551	+ ( ( flux_t(k) + diss_t(k) ) &
2552	- ( flux_d + diss_d ) &
2553	) * drho_air(k) * ddzw(k) &
2554	) + v(k,j,i) * div
2555
2556	flux_l_v(k,j,tn) = flux_r(k)
2557	diss_l_v(k,j,tn) = diss_r(k)
2558	flux_s_v(k,tn) = flux_n(k)
2559	diss_s_v(k,tn) = diss_n(k)
2560
2561	!
2562	!-- Statistical Evaluation of v'v'. The factor has to be applied for
2563	!-- right evaluation when gallilei_trans = .T. .
2564	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
2565	+ ( flux_n(k) &
2566	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2567	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
2568	+ diss_n(k) &
2569	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2570	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
2571	) * weight_substep(intermediate_timestep_count)
2572	!
2573	!-- Statistical Evaluation of w'u'.
2574	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
2575	+ ( flux_t(k) &
2576	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2577	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
2578	+ diss_t(k) &
2579	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2580	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
2581	) * weight_substep(intermediate_timestep_count)
2582
2583	ENDDO
2584
2585	DO k = nzb_max+1, nzt
2586
2587	u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu
2588	flux_r(k) = u_comp(k) * ( &
2589	37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &
2590	- 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &
2591	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
2592
2593	diss_r(k) = - ABS( u_comp(k) ) * ( &
2594	10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &
2595	- 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &
2596	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
2597
2598
2599	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2600	flux_n(k) = ( v_comp(k) - gv ) * ( &
2601	37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &
2602	- 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &
2603	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
2604
2605	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2606	10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &
2607	- 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &
2608	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
2609	!
2610	!-- k index has to be modified near bottom and top, else array
2611	!-- subscripts will be exceeded.
2612	ibit26 = IBITS(advc_flags_1(k,j,i),26,1)
2613	ibit25 = IBITS(advc_flags_1(k,j,i),25,1)
2614	ibit24 = IBITS(advc_flags_1(k,j,i),24,1)
2615
2616	k_ppp = k + 3 * ibit26
2617	k_pp = k + 2 * ( 1 - ibit24 )
2618	k_mm = k - 2 * ibit26
2619
2620	w_comp(k) = w(k,j-1,i) + w(k,j,i)
2621	flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &
2622	( 37.0_wp * ibit26 * adv_mom_5 &
2623	+ 7.0_wp * ibit25 * adv_mom_3 &
2624	+ ibit24 * adv_mom_1 &
2625	) * &
2626	( v(k+1,j,i) + v(k,j,i) ) &
2627	- ( 8.0_wp * ibit26 * adv_mom_5 &
2628	+ ibit25 * adv_mom_3 &
2629	) * &
2630	( v(k_pp,j,i) + v(k-1,j,i) ) &
2631	+ ( ibit26 * adv_mom_5 &
2632	) * &
2633	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
2634	)
2635
2636	diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &
2637	( 10.0_wp * ibit26 * adv_mom_5 &
2638	+ 3.0_wp * ibit25 * adv_mom_3 &
2639	+ ibit24 * adv_mom_1 &
2640	) * &
2641	( v(k+1,j,i) - v(k,j,i) ) &
2642	- ( 5.0_wp * ibit26 * adv_mom_5 &
2643	+ ibit25 * adv_mom_3 &
2644	) * &
2645	( v(k_pp,j,i) - v(k-1,j,i) ) &
2646	+ ( ibit26 * adv_mom_5 &
2647	) * &
2648	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
2649	)
2650	ENDDO
2651
2652	DO k = nzb_max+1, nzt
2653
2654	flux_d = flux_t(k-1)
2655	diss_d = diss_t(k-1)
2656	!
2657	!-- Calculate the divergence of the velocity field. A respective
2658	!-- correction is needed to overcome numerical instabilities introduced
2659	!-- by a not sufficient reduction of divergences near topography.
2660	div = ( ( u_comp(k) + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
2661	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
2662	+ ( w_comp(k) * rho_air_zw(k) - &
2663	w_comp(k-1) * rho_air_zw(k-1) &
2664	) * drho_air(k) * ddzw(k) &
2665	) * 0.5_wp
2666
2667	tend(k,j,i) = tend(k,j,i) - ( &
2668	( flux_r(k) + diss_r(k) &
2669	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
2670	+ ( flux_n(k) + diss_n(k) &
2671	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
2672	+ ( ( flux_t(k) + diss_t(k) ) &
2673	- ( flux_d + diss_d ) &
2674	) * drho_air(k) * ddzw(k) &
2675	) + v(k,j,i) * div
2676
2677	flux_l_v(k,j,tn) = flux_r(k)
2678	diss_l_v(k,j,tn) = diss_r(k)
2679	flux_s_v(k,tn) = flux_n(k)
2680	diss_s_v(k,tn) = diss_n(k)
2681	!
2682	!-- Statistical Evaluation of v'v'. The factor has to be applied for
2683	!-- right evaluation when gallilei_trans = .T. .
2684	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
2685	+ ( flux_n(k) &
2686	* ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2687	/ ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &
2688	+ diss_n(k) &
2689	* ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &
2690	/ ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &
2691	) * weight_substep(intermediate_timestep_count)
2692	!
2693	!-- Statistical Evaluation of w'u'.
2694	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
2695	+ ( flux_t(k) &
2696	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2697	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
2698	+ diss_t(k) &
2699	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
2700	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
2701	) * weight_substep(intermediate_timestep_count)
2702
2703	ENDDO
2704	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
2705
2706
2707	END SUBROUTINE advec_v_ws_ij
2708
2709
2710
2711	!------------------------------------------------------------------------------!
2712	! Description:
2713	! ------------
2714	!> Advection of w-component - Call for grid point i,j
2715	!------------------------------------------------------------------------------!
2716	SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
2717
2718	USE arrays_3d, &
2719	ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,&
2720	drho_air_zw, rho_air
2721
2722	USE control_parameters, &
2723	ONLY: intermediate_timestep_count, u_gtrans, v_gtrans
2724
2725	USE grid_variables, &
2726	ONLY: ddx, ddy
2727
2728	USE indices, &
2729	ONLY: nys, nzb, nzb_max, nzt, advc_flags_1, advc_flags_2
2730
2731	USE kinds
2732
2733	USE statistics, &
2734	ONLY: hom, sums_ws2_ws_l, weight_substep
2735
2736	IMPLICIT NONE
2737
2738	INTEGER(iwp) :: i !<
2739	INTEGER(iwp) :: ibit27 !<
2740	INTEGER(iwp) :: ibit28 !<
2741	INTEGER(iwp) :: ibit29 !<
2742	INTEGER(iwp) :: ibit30 !<
2743	INTEGER(iwp) :: ibit31 !<
2744	INTEGER(iwp) :: ibit32 !<
2745	INTEGER(iwp) :: ibit33 !<
2746	INTEGER(iwp) :: ibit34 !<
2747	INTEGER(iwp) :: ibit35 !<
2748	INTEGER(iwp) :: i_omp !<
2749	INTEGER(iwp) :: j !<
2750	INTEGER(iwp) :: k !<
2751	INTEGER(iwp) :: k_mm !<
2752	INTEGER(iwp) :: k_pp !<
2753	INTEGER(iwp) :: k_ppp !<
2754	INTEGER(iwp) :: tn !<
2755
2756	REAL(wp) :: diss_d !<
2757	REAL(wp) :: div !<
2758	REAL(wp) :: flux_d !<
2759	REAL(wp) :: gu !<
2760	REAL(wp) :: gv !<
2761
2762	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !<
2763	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !<
2764	REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !<
2765	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !<
2766	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !<
2767	REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !<
2768	REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !<
2769	REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !<
2770	REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !<
2771
2772	gu = 2.0_wp * u_gtrans
2773	gv = 2.0_wp * v_gtrans
2774
2775	!
2776	!-- Compute southside fluxes for the respective boundary.
2777	IF ( j == nys ) THEN
2778
2779	DO k = nzb+1, nzb_max
2780	ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)
2781	ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)
2782	ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)
2783
2784	v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv
2785	flux_s_w(k,tn) = v_comp(k) * ( &
2786	( 37.0_wp * ibit32 * adv_mom_5 &
2787	+ 7.0_wp * ibit31 * adv_mom_3 &
2788	+ ibit30 * adv_mom_1 &
2789	) * &
2790	( w(k,j,i) + w(k,j-1,i) ) &
2791	- ( 8.0_wp * ibit32 * adv_mom_5 &
2792	+ ibit31 * adv_mom_3 &
2793	) * &
2794	( w(k,j+1,i) + w(k,j-2,i) ) &
2795	+ ( ibit32 * adv_mom_5 &
2796	) * &
2797	( w(k,j+2,i) + w(k,j-3,i) ) &
2798	)
2799
2800	diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( &
2801	( 10.0_wp * ibit32 * adv_mom_5 &
2802	+ 3.0_wp * ibit31 * adv_mom_3 &
2803	+ ibit30 * adv_mom_1 &
2804	) * &
2805	( w(k,j,i) - w(k,j-1,i) ) &
2806	- ( 5.0_wp * ibit32 * adv_mom_5 &
2807	+ ibit31 * adv_mom_3 &
2808	) * &
2809	( w(k,j+1,i) - w(k,j-2,i) ) &
2810	+ ( ibit32 * adv_mom_5 &
2811	) * &
2812	( w(k,j+2,i) - w(k,j-3,i) ) &
2813	)
2814
2815	ENDDO
2816
2817	DO k = nzb_max+1, nzt
2818
2819	v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv
2820	flux_s_w(k,tn) = v_comp(k) * ( &
2821	37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) &
2822	- 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) &
2823	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
2824	diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( &
2825	10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) &
2826	- 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) &
2827	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
2828
2829	ENDDO
2830
2831	ENDIF
2832	!
2833	!-- Compute leftside fluxes for the respective boundary.
2834	IF ( i == i_omp ) THEN
2835
2836	DO k = nzb+1, nzb_max
2837
2838	ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)
2839	ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)
2840	ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)
2841
2842	u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu
2843	flux_l_w(k,j,tn) = u_comp(k) * ( &
2844	( 37.0_wp * ibit29 * adv_mom_5 &
2845	+ 7.0_wp * ibit28 * adv_mom_3 &
2846	+ ibit27 * adv_mom_1 &
2847	) * &
2848	( w(k,j,i) + w(k,j,i-1) ) &
2849	- ( 8.0_wp * ibit29 * adv_mom_5 &
2850	+ ibit28 * adv_mom_3 &
2851	) * &
2852	( w(k,j,i+1) + w(k,j,i-2) ) &
2853	+ ( ibit29 * adv_mom_5 &
2854	) * &
2855	( w(k,j,i+2) + w(k,j,i-3) ) &
2856	)
2857
2858	diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( &
2859	( 10.0_wp * ibit29 * adv_mom_5 &
2860	+ 3.0_wp * ibit28 * adv_mom_3 &
2861	+ ibit27 * adv_mom_1 &
2862	) * &
2863	( w(k,j,i) - w(k,j,i-1) ) &
2864	- ( 5.0_wp * ibit29 * adv_mom_5 &
2865	+ ibit28 * adv_mom_3 &
2866	) * &
2867	( w(k,j,i+1) - w(k,j,i-2) ) &
2868	+ ( ibit29 * adv_mom_5 &
2869	) * &
2870	( w(k,j,i+2) - w(k,j,i-3) ) &
2871	)
2872
2873	ENDDO
2874
2875	DO k = nzb_max+1, nzt
2876
2877	u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu
2878	flux_l_w(k,j,tn) = u_comp(k) * ( &
2879	37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) &
2880	- 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) &
2881	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
2882	diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( &
2883	10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) &
2884	- 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) &
2885	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
2886
2887	ENDDO
2888
2889	ENDIF
2890	!
2891	!-- The lower flux has to be calculated explicetely for the tendency at
2892	!-- the first w-level. For topography wall this is done implicitely by
2893	!-- advc_flags_1.
2894	k = nzb + 1
2895	w_comp(k) = w(k,j,i) + w(k-1,j,i)
2896	flux_t(0) = w_comp(k) * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
2897	diss_t(0) = -ABS(w_comp(k)) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
2898	!
2899	!-- Now compute the fluxes and tendency terms for the horizontal
2900	!-- and vertical parts.
2901	DO k = nzb+1, nzb_max
2902
2903	ibit29 = IBITS(advc_flags_1(k,j,i),29,1)
2904	ibit28 = IBITS(advc_flags_1(k,j,i),28,1)
2905	ibit27 = IBITS(advc_flags_1(k,j,i),27,1)
2906
2907	u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
2908	flux_r(k) = u_comp(k) * ( &
2909	( 37.0_wp * ibit29 * adv_mom_5 &
2910	+ 7.0_wp * ibit28 * adv_mom_3 &
2911	+ ibit27 * adv_mom_1 &
2912	) * &
2913	( w(k,j,i+1) + w(k,j,i) ) &
2914	- ( 8.0_wp * ibit29 * adv_mom_5 &
2915	+ ibit28 * adv_mom_3 &
2916	) * &
2917	( w(k,j,i+2) + w(k,j,i-1) ) &
2918	+ ( ibit29 * adv_mom_5 &
2919	) * &
2920	( w(k,j,i+3) + w(k,j,i-2) ) &
2921	)
2922
2923	diss_r(k) = - ABS( u_comp(k) ) * ( &
2924	( 10.0_wp * ibit29 * adv_mom_5 &
2925	+ 3.0_wp * ibit28 * adv_mom_3 &
2926	+ ibit27 * adv_mom_1 &
2927	) * &
2928	( w(k,j,i+1) - w(k,j,i) ) &
2929	- ( 5.0_wp * ibit29 * adv_mom_5 &
2930	+ ibit28 * adv_mom_3 &
2931	) * &
2932	( w(k,j,i+2) - w(k,j,i-1) ) &
2933	+ ( ibit29 * adv_mom_5 &
2934	) * &
2935	( w(k,j,i+3) - w(k,j,i-2) ) &
2936	)
2937
2938	ibit32 = IBITS(advc_flags_2(k,j,i),0,1)
2939	ibit31 = IBITS(advc_flags_1(k,j,i),31,1)
2940	ibit30 = IBITS(advc_flags_1(k,j,i),30,1)
2941
2942	v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
2943	flux_n(k) = v_comp(k) * ( &
2944	( 37.0_wp * ibit32 * adv_mom_5 &
2945	+ 7.0_wp * ibit31 * adv_mom_3 &
2946	+ ibit30 * adv_mom_1 &
2947	) * &
2948	( w(k,j+1,i) + w(k,j,i) ) &
2949	- ( 8.0_wp * ibit32 * adv_mom_5 &
2950	+ ibit31 * adv_mom_3 &
2951	) * &
2952	( w(k,j+2,i) + w(k,j-1,i) ) &
2953	+ ( ibit32 * adv_mom_5 &
2954	) * &
2955	( w(k,j+3,i) + w(k,j-2,i) ) &
2956	)
2957
2958	diss_n(k) = - ABS( v_comp(k) ) * ( &
2959	( 10.0_wp * ibit32 * adv_mom_5 &
2960	+ 3.0_wp * ibit31 * adv_mom_3 &
2961	+ ibit30 * adv_mom_1 &
2962	) * &
2963	( w(k,j+1,i) - w(k,j,i) ) &
2964	- ( 5.0_wp * ibit32 * adv_mom_5 &
2965	+ ibit31 * adv_mom_3 &
2966	) * &
2967	( w(k,j+2,i) - w(k,j-1,i) ) &
2968	+ ( ibit32 * adv_mom_5 &
2969	) * &
2970	( w(k,j+3,i) - w(k,j-2,i) ) &
2971	)
2972	!
2973	!-- k index has to be modified near bottom and top, else array
2974	!-- subscripts will be exceeded.
2975	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
2976	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
2977	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
2978
2979	k_ppp = k + 3 * ibit35
2980	k_pp = k + 2 * ( 1 - ibit33 )
2981	k_mm = k - 2 * ibit35
2982
2983	w_comp(k) = w(k+1,j,i) + w(k,j,i)
2984	flux_t(k) = w_comp(k) * rho_air(k+1) * ( &
2985	( 37.0_wp * ibit35 * adv_mom_5 &
2986	+ 7.0_wp * ibit34 * adv_mom_3 &
2987	+ ibit33 * adv_mom_1 &
2988	) * &
2989	( w(k+1,j,i) + w(k,j,i) ) &
2990	- ( 8.0_wp * ibit35 * adv_mom_5 &
2991	+ ibit34 * adv_mom_3 &
2992	) * &
2993	( w(k_pp,j,i) + w(k-1,j,i) ) &
2994	+ ( ibit35 * adv_mom_5 &
2995	) * &
2996	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
2997	)
2998
2999	diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( &
3000	( 10.0_wp * ibit35 * adv_mom_5 &
3001	+ 3.0_wp * ibit34 * adv_mom_3 &
3002	+ ibit33 * adv_mom_1 &
3003	) * &
3004	( w(k+1,j,i) - w(k,j,i) ) &
3005	- ( 5.0_wp * ibit35 * adv_mom_5 &
3006	+ ibit34 * adv_mom_3 &
3007	) * &
3008	( w(k_pp,j,i) - w(k-1,j,i) ) &
3009	+ ( ibit35 * adv_mom_5 &
3010	) * &
3011	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
3012	)
3013	ENDDO
3014
3015	DO k = nzb+1, nzb_max
3016
3017	flux_d = flux_t(k-1)
3018	diss_d = diss_t(k-1)
3019	!
3020	!-- Calculate the divergence of the velocity field. A respective
3021	!-- correction is needed to overcome numerical instabilities introduced
3022	!-- by a not sufficient reduction of divergences near topography.
3023	div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 ) &
3024	- ( u(k+1,j,i) + u(k,j,i) ) &
3025	* ( IBITS(advc_flags_1(k,j,i-1),27,1) &
3026	+ IBITS(advc_flags_1(k,j,i-1),28,1) &
3027	+ IBITS(advc_flags_1(k,j,i-1),29,1) &
3028	) &
3029	) * ddx &
3030	+ ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 ) &
3031	- ( v(k+1,j,i) + v(k,j,i) ) &
3032	* ( IBITS(advc_flags_1(k,j-1,i),30,1) &
3033	+ IBITS(advc_flags_1(k,j-1,i),31,1) &
3034	+ IBITS(advc_flags_2(k,j-1,i),0,1) &
3035	) &
3036	) * ddy &
3037	+ ( w_comp(k) * rho_air(k+1) &
3038	* ( ibit33 + ibit34 + ibit35 ) &
3039	- ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &
3040	* ( IBITS(advc_flags_2(k-1,j,i),1,1) &
3041	+ IBITS(advc_flags_2(k-1,j,i),2,1) &
3042	+ IBITS(advc_flags_2(k-1,j,i),3,1) &
3043	) &
3044	) * drho_air_zw(k) * ddzu(k+1) &
3045	) * 0.5_wp
3046
3047
3048	tend(k,j,i) = tend(k,j,i) - ( &
3049	( flux_r(k) + diss_r(k) &
3050	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
3051	+ ( flux_n(k) + diss_n(k) &
3052	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
3053	+ ( ( flux_t(k) + diss_t(k) ) &
3054	- ( flux_d + diss_d ) &
3055	) * drho_air_zw(k) * ddzu(k+1) &
3056	) + div * w(k,j,i)
3057
3058	flux_l_w(k,j,tn) = flux_r(k)
3059	diss_l_w(k,j,tn) = diss_r(k)
3060	flux_s_w(k,tn) = flux_n(k)
3061	diss_s_w(k,tn) = diss_n(k)
3062	!
3063	!-- Statistical Evaluation of w'w'.
3064	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
3065	+ ( flux_t(k) &
3066	* ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
3067	/ ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &
3068	+ diss_t(k) &
3069	* ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &
3070	/ ( ABS( w_comp(k) ) + 1.0E-20_wp ) &
3071	) * weight_substep(intermediate_timestep_count)
3072
3073	ENDDO
3074
3075	DO k = nzb_max+1, nzt
3076
3077	u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
3078	flux_r(k) = u_comp(k) * ( &
3079	37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) &
3080	- 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) &
3081	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
3082
3083	diss_r(k) = - ABS( u_comp(k) ) * ( &
3084	10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) &
3085	- 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) &
3086	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
3087
3088	v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv
3089	flux_n(k) = v_comp(k) * ( &
3090	37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) &
3091	- 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) &
3092	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
3093
3094	diss_n(k) = - ABS( v_comp(k) ) * ( &
3095	10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) &
3096	- 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) &
3097	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
3098	!
3099	!-- k index has to be modified near bottom and top, else array
3100	!-- subscripts will be exceeded.
3101	ibit35 = IBITS(advc_flags_2(k,j,i),3,1)
3102	ibit34 = IBITS(advc_flags_2(k,j,i),2,1)
3103	ibit33 = IBITS(advc_flags_2(k,j,i),1,1)
3104
3105	k_ppp = k + 3 * ibit35
3106	k_pp = k + 2 * ( 1 - ibit33 )
3107	k_mm = k - 2 * ibit35
3108
3109	w_comp(k) = w(k+1,j,i) + w(k,j,i)
3110	flux_t(k) = w_comp(k) * rho_air(k+1) * ( &
3111	( 37.0_wp * ibit35 * adv_mom_5 &
3112	+ 7.0_wp * ibit34 *