Home

Context Navigation

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

Last change on this file since 1027 was 1027, checked in by suehring, 13 years ago
Bugfix in calculation indices k_mm, k_pp in accelerator version
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 268.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |