Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |