Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |