Home

Context Navigation

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

Last change on this file since 1391 was 1375, checked in by raasch, 11 years ago
last commit documented
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 296.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |