Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |