Home

Context Navigation

advec_ws.f90 @ 1362

Last change on this file since 1362 was 1362, checked in by hoffmann, 10 years ago
last commit documented
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 296.8 KB

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