Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |