Home

Context Navigation

advec_ws.f90 @ 1361

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