Home

Context Navigation

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

Last change on this file since 1128 was 1128, checked in by raasch, 12 years ago
asynchronous transfer of ghost point data for acc-optimized version
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 271.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |