Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |