Home

Context Navigation

advec_ws.f90 @ 1115

Last change on this file since 1115 was 1115, checked in by hoffmann, 12 years ago
optimization of two-moments cloud physics
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 271.7 KB

Line
1	MODULE advec_ws
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	! calculation of qr and nr is restricted to precipitation
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: advec_ws.f90 1115 2013-03-26 18:16:16Z hoffmann $
27	!
28	! 1053 2012-11-13 17:11:03Z hoffmann
29	! necessary expansions according to the two new prognostic equations (nr, qr)
30	! of the two-moment cloud physics scheme:
31	! +flux_l_, flux_s_, diss_l_, diss_s_, sums_ws*s_ws_l
32	!
33	! 1036 2012-10-22 13:43:42Z raasch
34	! code put under GPL (PALM 3.9)
35	!
36	! 1027 2012-10-15 17:18:39Z suehring
37	! Bugfix in calculation indices k_mm, k_pp in accelerator version
38	!
39	! 1019 2012-09-28 06:46:45Z raasch
40	! small change in comment lines
41	!
42	! 1015 2012-09-27 09:23:24Z raasch
43	! accelerator versions (*_acc) added
44	!
45	! 1010 2012-09-20 07:59:54Z raasch
46	! cpp switch __nopointer added for pointer free version
47	!
48	! 888 2012-04-20 15:03:46Z suehring
49	! Number of IBITS() calls with identical arguments is reduced.
50	!
51	! 862 2012-03-26 14:21:38Z suehring
52	! ws-scheme also work with topography in combination with vector version.
53	! ws-scheme also work with outflow boundaries in combination with
54	! vector version.
55	! Degradation of the applied order of scheme is now steered by multiplying with
56	! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each
57	! grid point and mulitplied with the appropriate flag.
58	! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
59	! term derived according to the 4th and 6th order terms is applied. It turns
60	! out that diss_2nd does not provide sufficient dissipation near walls.
61	! Therefore, the function diss_2nd is removed.
62	! Near walls a divergence correction is necessary to overcome numerical
63	! instabilities due to too less divergence reduction of the velocity field.
64	! boundary_flags and logicals steering the degradation are removed.
65	! Empty SUBROUTINE local_diss removed.
66	! Further formatting adjustments.
67	!
68	! 801 2012-01-10 17:30:36Z suehring
69	! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
70	! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
71	! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
72	! dimension.
73	!
74	! 743 2011-08-18 16:10:16Z suehring
75	! Evaluation of turbulent fluxes with WS-scheme only for the whole model
76	! domain. Therefor dimension of arrays needed for statistical evaluation
77	! decreased.
78	!
79	! 736 2011-08-17 14:13:26Z suehring
80	! Bugfix concerning OpenMP parallelization. i_omp introduced, because first
81	! index where fluxes on left side have to be calculated explicitly is
82	! different on each thread. Furthermore the swapping of fluxes is now
83	! thread-safe by adding an additional dimension.
84	!
85	! 713 2011-03-30 14:21:21Z suehring
86	! File reformatted.
87	! Bugfix in vertical advection of w concerning the optimized version for
88	! vector architecture.
89	! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as
90	! broken numbers.
91	!
92	! 709 2011-03-30 09:31:40Z raasch
93	! formatting adjustments
94	!
95	! 705 2011-03-25 11:21:43 Z suehring $
96	! Bugfix in declaration of logicals concerning outflow boundaries.
97	!
98	! 411 2009-12-11 12:31:43 Z suehring
99	! Allocation of weight_substep moved to init_3d_model.
100	! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters.
101	! Setting bc for the horizontal velocity variances added (moved from
102	! flow_statistics).
103	!
104	! Initial revision
105	!
106	! 411 2009-12-11 12:31:43 Z suehring
107	!
108	! Description:
109	! ------------
110	! Advection scheme for scalars and momentum using the flux formulation of
111	! Wicker and Skamarock 5th order. Additionally the module contains of a
112	! routine using for initialisation and steering of the statical evaluation.
113	! The computation of turbulent fluxes takes place inside the advection
114	! routines.
115	! Near non-cyclic boundaries the order of the applied advection scheme is
116	! degraded.
117	! A divergence correction is applied. It is necessary for topography, since
118	! the divergence is not sufficiently reduced, resulting in erroneous fluxes and
119	! partly numerical instabilities.
120	!-----------------------------------------------------------------------------!
121
122	PRIVATE
123	PUBLIC advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, &
124	advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc, &
125	ws_init, ws_statistics
126
127	INTERFACE ws_init
128	MODULE PROCEDURE ws_init
129	END INTERFACE ws_init
130
131	INTERFACE ws_statistics
132	MODULE PROCEDURE ws_statistics
133	END INTERFACE ws_statistics
134
135	INTERFACE advec_s_ws
136	MODULE PROCEDURE advec_s_ws
137	MODULE PROCEDURE advec_s_ws_ij
138	END INTERFACE advec_s_ws
139
140	INTERFACE advec_u_ws
141	MODULE PROCEDURE advec_u_ws
142	MODULE PROCEDURE advec_u_ws_ij
143	END INTERFACE advec_u_ws
144
145	INTERFACE advec_u_ws_acc
146	MODULE PROCEDURE advec_u_ws_acc
147	END INTERFACE advec_u_ws_acc
148
149	INTERFACE advec_v_ws
150	MODULE PROCEDURE advec_v_ws
151	MODULE PROCEDURE advec_v_ws_ij
152	END INTERFACE advec_v_ws
153
154	INTERFACE advec_v_ws_acc
155	MODULE PROCEDURE advec_v_ws_acc
156	END INTERFACE advec_v_ws_acc
157
158	INTERFACE advec_w_ws
159	MODULE PROCEDURE advec_w_ws
160	MODULE PROCEDURE advec_w_ws_ij
161	END INTERFACE advec_w_ws
162
163	INTERFACE advec_w_ws_acc
164	MODULE PROCEDURE advec_w_ws_acc
165	END INTERFACE advec_w_ws_acc
166
167	CONTAINS
168
169
170	!------------------------------------------------------------------------------!
171	! Initialization of WS-scheme
172	!------------------------------------------------------------------------------!
173	SUBROUTINE ws_init
174
175	USE arrays_3d
176	USE constants
177	USE control_parameters
178	USE indices
179	USE pegrid
180	USE statistics
181
182	!
183	!-- Set the appropriate factors for scalar and momentum advection.
184	adv_sca_5 = 1./60.
185	adv_sca_3 = 1./12.
186	adv_sca_1 = 1./2.
187	adv_mom_5 = 1./120.
188	adv_mom_3 = 1./24.
189	adv_mom_1 = 1./4.
190	!
191	!-- Arrays needed for statical evaluation of fluxes.
192	IF ( ws_scheme_mom ) THEN
193
194	ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), &
195	sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), &
196	sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
197	sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
198	sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
199
200	sums_wsus_ws_l = 0.0
201	sums_wsvs_ws_l = 0.0
202	sums_us2_ws_l = 0.0
203	sums_vs2_ws_l = 0.0
204	sums_ws2_ws_l = 0.0
205
206	ENDIF
207
208	IF ( ws_scheme_sca ) THEN
209
210	ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
211	sums_wspts_ws_l = 0.0
212
213	IF ( humidity .OR. passive_scalar ) THEN
214	ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
215	sums_wsqs_ws_l = 0.0
216	ENDIF
217
218	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
219	precipitation ) THEN
220	ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
221	ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
222	sums_wsqrs_ws_l = 0.0
223	sums_wsnrs_ws_l = 0.0
224	ENDIF
225
226	IF ( ocean ) THEN
227	ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
228	sums_wssas_ws_l = 0.0
229	ENDIF
230
231	ENDIF
232
233	!
234	!-- Arrays needed for reasons of speed optimization for cache version.
235	!-- For the vector version the buffer arrays are not necessary,
236	!-- because the the fluxes can swapped directly inside the loops of the
237	!-- advection routines.
238	IF ( loop_optimization /= 'vector' ) THEN
239
240	IF ( ws_scheme_mom ) THEN
241
242	ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &
243	flux_s_v(nzb+1:nzt,0:threads_per_task-1), &
244	flux_s_w(nzb+1:nzt,0:threads_per_task-1), &
245	diss_s_u(nzb+1:nzt,0:threads_per_task-1), &
246	diss_s_v(nzb+1:nzt,0:threads_per_task-1), &
247	diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
248	ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
249	flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
250	flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
251	diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
252	diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
253	diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
254
255	ENDIF
256
257	IF ( ws_scheme_sca ) THEN
258
259	ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), &
260	flux_s_e(nzb+1:nzt,0:threads_per_task-1), &
261	diss_s_pt(nzb+1:nzt,0:threads_per_task-1), &
262	diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
263	ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
264	flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
265	diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
266	diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
267
268	IF ( humidity .OR. passive_scalar ) THEN
269	ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), &
270	diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
271	ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
272	diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
273	ENDIF
274
275	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
276	precipitation ) THEN
277	ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), &
278	diss_s_qr(nzb+1:nzt,0:threads_per_task-1), &
279	flux_s_nr(nzb+1:nzt,0:threads_per_task-1), &
280	diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
281	ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
282	diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
283	flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
284	diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
285	ENDIF
286
287	IF ( ocean ) THEN
288	ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), &
289	diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
290	ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
291	diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
292	ENDIF
293
294	ENDIF
295
296	ENDIF
297
298	END SUBROUTINE ws_init
299
300
301	!------------------------------------------------------------------------------!
302	! Initialize variables used for storing statistic quantities (fluxes, variances)
303	!------------------------------------------------------------------------------!
304	SUBROUTINE ws_statistics
305
306	USE control_parameters
307	USE statistics
308
309	IMPLICIT NONE
310
311	!
312	!-- The arrays needed for statistical evaluation are set to to 0 at the
313	!-- beginning of prognostic_equations.
314	IF ( ws_scheme_mom ) THEN
315	sums_wsus_ws_l = 0.0
316	sums_wsvs_ws_l = 0.0
317	sums_us2_ws_l = 0.0
318	sums_vs2_ws_l = 0.0
319	sums_ws2_ws_l = 0.0
320	ENDIF
321
322	IF ( ws_scheme_sca ) THEN
323	sums_wspts_ws_l = 0.0
324	IF ( humidity .OR. passive_scalar ) sums_wsqs_ws_l = 0.0
325	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
326	precipitation ) THEN
327	sums_wsqrs_ws_l = 0.0
328	sums_wsnrs_ws_l = 0.0
329	ENDIF
330	IF ( ocean ) sums_wssas_ws_l = 0.0
331
332	ENDIF
333
334	END SUBROUTINE ws_statistics
335
336
337	!------------------------------------------------------------------------------!
338	! Scalar advection - Call for grid point i,j
339	!------------------------------------------------------------------------------!
340	SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local, &
341	swap_diss_y_local, swap_flux_x_local, &
342	swap_diss_x_local, i_omp, tn )
343
344	USE arrays_3d
345	USE constants
346	USE control_parameters
347	USE grid_variables
348	USE indices
349	USE pegrid
350	USE statistics
351
352	IMPLICIT NONE
353
354	INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, &
355	ibit7, ibit8, i_omp, j, k, k_mm, k_pp, k_ppp, tn
356	REAL :: diss_d, div, flux_d, u_comp, v_comp
357	#if defined( __nopointer )
358	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk
359	#else
360	REAL, DIMENSION(:,:,:), POINTER :: sk
361	#endif
362	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, &
363	flux_r, flux_t
364	REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local, &
365	swap_flux_y_local
366	REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: &
367	swap_diss_x_local, &
368	swap_flux_x_local
369	CHARACTER (LEN = *), INTENT(IN) :: sk_char
370
371	!
372	!-- Compute southside fluxes of the respective PE bounds.
373	IF ( j == nys ) THEN
374	!
375	!-- Up to the top of the highest topography.
376	DO k = nzb+1, nzb_max
377
378	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
379	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
380	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
381
382	v_comp = v(k,j,i) - v_gtrans
383	swap_flux_y_local(k,tn) = v_comp * ( &
384	( 37.0 * ibit5 * adv_sca_5 &
385	+ 7.0 * ibit4 * adv_sca_3 &
386	+ ibit3 * adv_sca_1 &
387	) * &
388	( sk(k,j,i) + sk(k,j-1,i) ) &
389	- ( 8.0 * ibit5 * adv_sca_5 &
390	+ ibit4 * adv_sca_3 &
391	) * &
392	( sk(k,j+1,i) + sk(k,j-2,i) ) &
393	+ ( ibit5 * adv_sca_5 &
394	) * &
395	( sk(k,j+2,i) + sk(k,j-3,i) ) &
396	)
397
398	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
399	( 10.0 * ibit5 * adv_sca_5 &
400	+ 3.0 * ibit4 * adv_sca_3 &
401	+ ibit3 * adv_sca_1 &
402	) * &
403	( sk(k,j,i) - sk(k,j-1,i) ) &
404	- ( 5.0 * ibit5 * adv_sca_5 &
405	+ ibit4 * adv_sca_3 &
406	) * &
407	( sk(k,j+1,i) - sk(k,j-2,i) ) &
408	+ ( ibit5 * adv_sca_5 &
409	) * &
410	( sk(k,j+2,i) - sk(k,j-3,i) ) &
411	)
412
413	ENDDO
414	!
415	!-- Above to the top of the highest topography. No degradation necessary.
416	DO k = nzb_max+1, nzt
417
418	v_comp = v(k,j,i) - v_gtrans
419	swap_flux_y_local(k,tn) = v_comp * ( &
420	37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) &
421	- 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
422	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
423	) * adv_sca_5
424	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
425	10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) &
426	- 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
427	+ sk(k,j+2,i) - sk(k,j-3,i) &
428	) * adv_sca_5
429
430	ENDDO
431
432	ENDIF
433	!
434	!-- Compute leftside fluxes of the respective PE bounds.
435	IF ( i == i_omp ) THEN
436
437	DO k = nzb+1, nzb_max
438
439	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
440	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
441	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
442
443	u_comp = u(k,j,i) - u_gtrans
444	swap_flux_x_local(k,j,tn) = u_comp * ( &
445	( 37.0 * ibit2 * adv_sca_5 &
446	+ 7.0 * ibit1 * adv_sca_3 &
447	+ ibit0 * adv_sca_1 &
448	) * &
449	( sk(k,j,i) + sk(k,j,i-1) ) &
450	- ( 8.0 * ibit2 * adv_sca_5 &
451	+ ibit1 * adv_sca_3 &
452	) * &
453	( sk(k,j,i+1) + sk(k,j,i-2) ) &
454	+ ( ibit2 * adv_sca_5 &
455	) * &
456	( sk(k,j,i+2) + sk(k,j,i-3) ) &
457	)
458
459	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
460	( 10.0 * ibit2 * adv_sca_5 &
461	+ 3.0 * ibit1 * adv_sca_3 &
462	+ ibit0 * adv_sca_1 &
463	) * &
464	( sk(k,j,i) - sk(k,j,i-1) ) &
465	- ( 5.0 * ibit2 * adv_sca_5 &
466	+ ibit1 * adv_sca_3 &
467	) * &
468	( sk(k,j,i+1) - sk(k,j,i-2) ) &
469	+ ( ibit2 * adv_sca_5 &
470	) * &
471	( sk(k,j,i+2) - sk(k,j,i-3) ) &
472	)
473
474	ENDDO
475
476	DO k = nzb_max+1, nzt
477
478	u_comp = u(k,j,i) - u_gtrans
479	swap_flux_x_local(k,j,tn) = u_comp * ( &
480	37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) &
481	- 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
482	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
483	) * adv_sca_5
484
485	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
486	10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) &
487	- 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
488	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
489	) * adv_sca_5
490
491	ENDDO
492
493	ENDIF
494
495	flux_t(0) = 0.0
496	diss_t(0) = 0.0
497	flux_d = 0.0
498	diss_d = 0.0
499	!
500	!-- Now compute the fluxes and tendency terms for the horizontal and
501	!-- vertical parts up to the top of the highest topography.
502	DO k = nzb+1, nzb_max
503	!
504	!-- Note: It is faster to conduct all multiplications explicitly, e.g.
505	!-- * adv_sca_5 ... than to determine a factor and multiplicate the
506	!-- flux at the end.
507
508	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
509	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
510	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
511
512	u_comp = u(k,j,i+1) - u_gtrans
513	flux_r(k) = u_comp * ( &
514	( 37.0 * ibit2 * adv_sca_5 &
515	+ 7.0 * ibit1 * adv_sca_3 &
516	+ ibit0 * adv_sca_1 &
517	) * &
518	( sk(k,j,i+1) + sk(k,j,i) ) &
519	- ( 8.0 * ibit2 * adv_sca_5 &
520	+ ibit1 * adv_sca_3 &
521	) * &
522	( sk(k,j,i+2) + sk(k,j,i-1) ) &
523	+ ( ibit2 * adv_sca_5 &
524	) * &
525	( sk(k,j,i+3) + sk(k,j,i-2) ) &
526	)
527
528	diss_r(k) = -ABS( u_comp ) * ( &
529	( 10.0 * ibit2 * adv_sca_5 &
530	+ 3.0 * ibit1 * adv_sca_3 &
531	+ ibit0 * adv_sca_1 &
532	) * &
533	( sk(k,j,i+1) - sk(k,j,i) ) &
534	- ( 5.0 * ibit2 * adv_sca_5 &
535	+ ibit1 * adv_sca_3 &
536	) * &
537	( sk(k,j,i+2) - sk(k,j,i-1) ) &
538	+ ( ibit2 * adv_sca_5 &
539	) * &
540	( sk(k,j,i+3) - sk(k,j,i-2) ) &
541	)
542
543	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
544	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
545	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
546
547	v_comp = v(k,j+1,i) - v_gtrans
548	flux_n(k) = v_comp * ( &
549	( 37.0 * ibit5 * adv_sca_5 &
550	+ 7.0 * ibit4 * adv_sca_3 &
551	+ ibit3 * adv_sca_1 &
552	) * &
553	( sk(k,j+1,i) + sk(k,j,i) ) &
554	- ( 8.0 * ibit5 * adv_sca_5 &
555	+ ibit4 * adv_sca_3 &
556	) * &
557	( sk(k,j+2,i) + sk(k,j-1,i) ) &
558	+ ( ibit5 * adv_sca_5 &
559	) * &
560	( sk(k,j+3,i) + sk(k,j-2,i) ) &
561	)
562
563	diss_n(k) = -ABS( v_comp ) * ( &
564	( 10.0 * ibit5 * adv_sca_5 &
565	+ 3.0 * ibit4 * adv_sca_3 &
566	+ ibit3 * adv_sca_1 &
567	) * &
568	( sk(k,j+1,i) - sk(k,j,i) ) &
569	- ( 5.0 * ibit5 * adv_sca_5 &
570	+ ibit4 * adv_sca_3 &
571	) * &
572	( sk(k,j+2,i) - sk(k,j-1,i) ) &
573	+ ( ibit5 * adv_sca_5 &
574	) * &
575	( sk(k,j+3,i) - sk(k,j-2,i) ) &
576	)
577	!
578	!-- k index has to be modified near bottom and top, else array
579	!-- subscripts will be exceeded.
580	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
581	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
582	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
583
584	k_ppp = k + 3 * ibit8
585	k_pp = k + 2 * ( 1 - ibit6 )
586	k_mm = k - 2 * ibit8
587
588
589	flux_t(k) = w(k,j,i) * ( &
590	( 37.0 * ibit8 * adv_sca_5 &
591	+ 7.0 * ibit7 * adv_sca_3 &
592	+ ibit6 * adv_sca_1 &
593	) * &
594	( sk(k+1,j,i) + sk(k,j,i) ) &
595	- ( 8.0 * ibit8 * adv_sca_5 &
596	+ ibit7 * adv_sca_3 &
597	) * &
598	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
599	+ ( ibit8 * adv_sca_5 &
600	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
601	)
602
603	diss_t(k) = -ABS( w(k,j,i) ) * ( &
604	( 10.0 * ibit8 * adv_sca_5 &
605	+ 3.0 * ibit7 * adv_sca_3 &
606	+ ibit6 * adv_sca_1 &
607	) * &
608	( sk(k+1,j,i) - sk(k,j,i) ) &
609	- ( 5.0 * ibit8 * adv_sca_5 &
610	+ ibit7 * adv_sca_3 &
611	) * &
612	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
613	+ ( ibit8 * adv_sca_5 &
614	) * &
615	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
616	)
617	!
618	!-- Calculate the divergence of the velocity field. A respective
619	!-- correction is needed to overcome numerical instabilities caused
620	!-- by a not sufficient reduction of divergences near topography.
621	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
622	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
623	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
624
625	tend(k,j,i) = tend(k,j,i) - ( &
626	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
627	swap_diss_x_local(k,j,tn) ) * ddx &
628	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
629	swap_diss_y_local(k,tn) ) * ddy &
630	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
631	) * ddzw(k) &
632	) + sk(k,j,i) * div
633
634	swap_flux_y_local(k,tn) = flux_n(k)
635	swap_diss_y_local(k,tn) = diss_n(k)
636	swap_flux_x_local(k,j,tn) = flux_r(k)
637	swap_diss_x_local(k,j,tn) = diss_r(k)
638	flux_d = flux_t(k)
639	diss_d = diss_t(k)
640
641	ENDDO
642	!
643	!-- Now compute the fluxes and tendency terms for the horizontal and
644	!-- vertical parts above the top of the highest topography. No degradation
645	!-- for the horizontal parts, but for the vertical it is stell needed.
646	DO k = nzb_max+1, nzt
647
648	u_comp = u(k,j,i+1) - u_gtrans
649	flux_r(k) = u_comp * ( &
650	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
651	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
652	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
653	diss_r(k) = -ABS( u_comp ) * ( &
654	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
655	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
656	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
657
658	v_comp = v(k,j+1,i) - v_gtrans
659	flux_n(k) = v_comp * ( &
660	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
661	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
662	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
663	diss_n(k) = -ABS( v_comp ) * ( &
664	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
665	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
666	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
667	!
668	!-- k index has to be modified near bottom and top, else array
669	!-- subscripts will be exceeded.
670	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
671	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
672	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
673
674	k_ppp = k + 3 * ibit8
675	k_pp = k + 2 * ( 1 - ibit6 )
676	k_mm = k - 2 * ibit8
677
678
679	flux_t(k) = w(k,j,i) * ( &
680	( 37.0 * ibit8 * adv_sca_5 &
681	+ 7.0 * ibit7 * adv_sca_3 &
682	+ ibit6 * adv_sca_1 &
683	) * &
684	( sk(k+1,j,i) + sk(k,j,i) ) &
685	- ( 8.0 * ibit8 * adv_sca_5 &
686	+ ibit7 * adv_sca_3 &
687	) * &
688	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
689	+ ( ibit8 * adv_sca_5 &
690	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
691	)
692
693	diss_t(k) = -ABS( w(k,j,i) ) * ( &
694	( 10.0 * ibit8 * adv_sca_5 &
695	+ 3.0 * ibit7 * adv_sca_3 &
696	+ ibit6 * adv_sca_1 &
697	) * &
698	( sk(k+1,j,i) - sk(k,j,i) ) &
699	- ( 5.0 * ibit8 * adv_sca_5 &
700	+ ibit7 * adv_sca_3 &
701	) * &
702	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
703	+ ( ibit8 * adv_sca_5 &
704	) * &
705	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
706	)
707	!
708	!-- Calculate the divergence of the velocity field. A respective
709	!-- correction is needed to overcome numerical instabilities introduced
710	!-- by a not sufficient reduction of divergences near topography.
711	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
712	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
713	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
714
715	tend(k,j,i) = tend(k,j,i) - ( &
716	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
717	swap_diss_x_local(k,j,tn) ) * ddx &
718	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
719	swap_diss_y_local(k,tn) ) * ddy &
720	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
721	) * ddzw(k) &
722	) + sk(k,j,i) * div
723
724	swap_flux_y_local(k,tn) = flux_n(k)
725	swap_diss_y_local(k,tn) = diss_n(k)
726	swap_flux_x_local(k,j,tn) = flux_r(k)
727	swap_diss_x_local(k,j,tn) = diss_r(k)
728	flux_d = flux_t(k)
729	diss_d = diss_t(k)
730
731	ENDDO
732
733	!
734	!-- Evaluation of statistics
735	SELECT CASE ( sk_char )
736
737	CASE ( 'pt' )
738
739	DO k = nzb, nzt
740	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + &
741	( flux_t(k) + diss_t(k) ) &
742	* weight_substep(intermediate_timestep_count)
743	ENDDO
744
745	CASE ( 'sa' )
746
747	DO k = nzb, nzt
748	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + &
749	( flux_t(k) + diss_t(k) ) &
750	* weight_substep(intermediate_timestep_count)
751	ENDDO
752
753	CASE ( 'q' )
754
755	DO k = nzb, nzt
756	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + &
757	( flux_t(k) + diss_t(k) ) &
758	* weight_substep(intermediate_timestep_count)
759	ENDDO
760
761	CASE ( 'qr' )
762
763	DO k = nzb, nzt
764	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + &
765	( flux_t(k) + diss_t(k) ) &
766	* weight_substep(intermediate_timestep_count)
767	ENDDO
768
769	CASE ( 'nr' )
770
771	DO k = nzb, nzt
772	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + &
773	( flux_t(k) + diss_t(k) ) &
774	* weight_substep(intermediate_timestep_count)
775	ENDDO
776
777	END SELECT
778
779	END SUBROUTINE advec_s_ws_ij
780
781
782
783
784	!------------------------------------------------------------------------------!
785	! Advection of u-component - Call for grid point i,j
786	!------------------------------------------------------------------------------!
787	SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
788
789	USE arrays_3d
790	USE constants
791	USE control_parameters
792	USE grid_variables
793	USE indices
794	USE statistics
795
796	IMPLICIT NONE
797
798	INTEGER :: i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15, &
799	ibit16, ibit17, i_omp, j, k, k_mm, k_pp, k_ppp, tn
800	REAL :: diss_d, div, flux_d, gu, gv, u_comp_l, v_comp, w_comp
801	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
802	flux_t, u_comp
803
804	gu = 2.0 * u_gtrans
805	gv = 2.0 * v_gtrans
806	!
807	!-- Compute southside fluxes for the respective boundary of PE
808	IF ( j == nys ) THEN
809
810	DO k = nzb+1, nzb_max
811
812	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
813	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
814	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
815
816	v_comp = v(k,j,i) + v(k,j,i-1) - gv
817	flux_s_u(k,tn) = v_comp * ( &
818	( 37.0 * ibit14 * adv_mom_5 &
819	+ 7.0 * ibit13 * adv_mom_3 &
820	+ ibit12 * adv_mom_1 &
821	) * &
822	( u(k,j,i) + u(k,j-1,i) ) &
823	- ( 8.0 * ibit14 * adv_mom_5 &
824	+ ibit13 * adv_mom_3 &
825	) * &
826	( u(k,j+1,i) + u(k,j-2,i) ) &
827	+ ( ibit14 * adv_mom_5 &
828	) * &
829	( u(k,j+2,i) + u(k,j-3,i) ) &
830	)
831
832	diss_s_u(k,tn) = - ABS ( v_comp ) * ( &
833	( 10.0 * ibit14 * adv_mom_5 &
834	+ 3.0 * ibit13 * adv_mom_3 &
835	+ ibit12 * adv_mom_1 &
836	) * &
837	( u(k,j,i) - u(k,j-1,i) ) &
838	- ( 5.0 * ibit14 * adv_mom_5 &
839	+ ibit13 * adv_mom_3 &
840	) * &
841	( u(k,j+1,i) - u(k,j-2,i) ) &
842	+ ( ibit14 * adv_mom_5 &
843	) * &
844	( u(k,j+2,i) - u(k,j-3,i) ) &
845	)
846
847	ENDDO
848
849	DO k = nzb_max+1, nzt
850
851	v_comp = v(k,j,i) + v(k,j,i-1) - gv
852	flux_s_u(k,tn) = v_comp * ( &
853	37.0 * ( u(k,j,i) + u(k,j-1,i) ) &
854	- 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) &
855	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
856	diss_s_u(k,tn) = - ABS(v_comp) * ( &
857	10.0 * ( u(k,j,i) - u(k,j-1,i) ) &
858	- 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) &
859	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
860
861	ENDDO
862
863	ENDIF
864	!
865	!-- Compute leftside fluxes for the respective boundary of PE
866	IF ( i == i_omp ) THEN
867
868	DO k = nzb+1, nzb_max
869
870	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
871	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
872	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
873
874	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
875	flux_l_u(k,j,tn) = u_comp_l * ( &
876	( 37.0 * ibit11 * adv_mom_5 &
877	+ 7.0 * ibit10 * adv_mom_3 &
878	+ ibit9 * adv_mom_1 &
879	) * &
880	( u(k,j,i) + u(k,j,i-1) ) &
881	- ( 8.0 * ibit11 * adv_mom_5 &
882	+ ibit10 * adv_mom_3 &
883	) * &
884	( u(k,j,i+1) + u(k,j,i-2) ) &
885	+ ( ibit11 * adv_mom_5 &
886	) * &
887	( u(k,j,i+2) + u(k,j,i-3) ) &
888	)
889
890	diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( &
891	( 10.0 * ibit11 * adv_mom_5 &
892	+ 3.0 * ibit10 * adv_mom_3 &
893	+ ibit9 * adv_mom_1 &
894	) * &
895	( u(k,j,i) - u(k,j,i-1) ) &
896	- ( 5.0 * ibit11 * adv_mom_5 &
897	+ ibit10 * adv_mom_3 &
898	) * &
899	( u(k,j,i+1) - u(k,j,i-2) ) &
900	+ ( ibit11 * adv_mom_5 &
901	) * &
902	( u(k,j,i+2) - u(k,j,i-3) ) &
903	)
904
905	ENDDO
906
907	DO k = nzb_max+1, nzt
908
909	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
910	flux_l_u(k,j,tn) = u_comp_l * ( &
911	37.0 * ( u(k,j,i) + u(k,j,i-1) ) &
912	- 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) &
913	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
914	diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( &
915	10.0 * ( u(k,j,i) - u(k,j,i-1) ) &
916	- 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) &
917	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
918
919	ENDDO
920
921	ENDIF
922
923	flux_t(0) = 0.0
924	diss_t(0) = 0.0
925	flux_d = 0.0
926	diss_d = 0.0
927	!
928	!-- Now compute the fluxes tendency terms for the horizontal and
929	!-- vertical parts.
930	DO k = nzb+1, nzb_max
931
932	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
933	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
934	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
935
936	u_comp(k) = u(k,j,i+1) + u(k,j,i)
937	flux_r(k) = ( u_comp(k) - gu ) * ( &
938	( 37.0 * ibit11 * adv_mom_5 &
939	+ 7.0 * ibit10 * adv_mom_3 &
940	+ ibit9 * adv_mom_1 &
941	) * &
942	( u(k,j,i+1) + u(k,j,i) ) &
943	- ( 8.0 * ibit11 * adv_mom_5 &
944	+ ibit10 * adv_mom_3 &
945	) * &
946	( u(k,j,i+2) + u(k,j,i-1) ) &
947	+ ( ibit11 * adv_mom_5 &
948	) * &
949	( u(k,j,i+3) + u(k,j,i-2) ) &
950	)
951
952	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
953	( 10.0 * ibit11 * adv_mom_5 &
954	+ 3.0 * ibit10 * adv_mom_3 &
955	+ ibit9 * adv_mom_1 &
956	) * &
957	( u(k,j,i+1) - u(k,j,i) ) &
958	- ( 5.0 * ibit11 * adv_mom_5 &
959	+ ibit10 * adv_mom_3 &
960	) * &
961	( u(k,j,i+2) - u(k,j,i-1) ) &
962	+ ( ibit11 * adv_mom_5 &
963	) * &
964	( u(k,j,i+3) - u(k,j,i-2) ) &
965	)
966
967	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
968	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
969	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
970
971	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
972	flux_n(k) = v_comp * ( &
973	( 37.0 * ibit14 * adv_mom_5 &
974	+ 7.0 * ibit13 * adv_mom_3 &
975	+ ibit12 * adv_mom_1 &
976	) * &
977	( u(k,j+1,i) + u(k,j,i) ) &
978	- ( 8.0 * ibit14 * adv_mom_5 &
979	+ ibit13 * adv_mom_3 &
980	) * &
981	( u(k,j+2,i) + u(k,j-1,i) ) &
982	+ ( ibit14 * adv_mom_5 &
983	) * &
984	( u(k,j+3,i) + u(k,j-2,i) ) &
985	)
986
987	diss_n(k) = - ABS ( v_comp ) * ( &
988	( 10.0 * ibit14 * adv_mom_5 &
989	+ 3.0 * ibit13 * adv_mom_3 &
990	+ ibit12 * adv_mom_1 &
991	) * &
992	( u(k,j+1,i) - u(k,j,i) ) &
993	- ( 5.0 * ibit14 * adv_mom_5 &
994	+ ibit13 * adv_mom_3 &
995	) * &
996	( u(k,j+2,i) - u(k,j-1,i) ) &
997	+ ( ibit14 * adv_mom_5 &
998	) * &
999	( u(k,j+3,i) - u(k,j-2,i) ) &
1000	)
1001	!
1002	!-- k index has to be modified near bottom and top, else array
1003	!-- subscripts will be exceeded.
1004	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
1005	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
1006	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
1007
1008	k_ppp = k + 3 * ibit17
1009	k_pp = k + 2 * ( 1 - ibit15 )
1010	k_mm = k - 2 * ibit17
1011
1012	w_comp = w(k,j,i) + w(k,j,i-1)
1013	flux_t(k) = w_comp * ( &
1014	( 37.0 * ibit17 * adv_mom_5 &
1015	+ 7.0 * ibit16 * adv_mom_3 &
1016	+ ibit15 * adv_mom_1 &
1017	) * &
1018	( u(k+1,j,i) + u(k,j,i) ) &
1019	- ( 8.0 * ibit17 * adv_mom_5 &
1020	+ ibit16 * adv_mom_3 &
1021	) * &
1022	( u(k_pp,j,i) + u(k-1,j,i) ) &
1023	+ ( ibit17 * adv_mom_5 &
1024	) * &
1025	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
1026	)
1027
1028	diss_t(k) = - ABS( w_comp ) * ( &
1029	( 10.0 * ibit17 * adv_mom_5 &
1030	+ 3.0 * ibit16 * adv_mom_3 &
1031	+ ibit15 * adv_mom_1 &
1032	) * &
1033	( u(k+1,j,i) - u(k,j,i) ) &
1034	- ( 5.0 * ibit17 * adv_mom_5 &
1035	+ ibit16 * adv_mom_3 &
1036	) * &
1037	( u(k_pp,j,i) - u(k-1,j,i) ) &
1038	+ ( ibit17 * adv_mom_5 &
1039	) * &
1040	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
1041	)
1042	!
1043	!-- Calculate the divergence of the velocity field. A respective
1044	!-- correction is needed to overcome numerical instabilities introduced
1045	!-- by a not sufficient reduction of divergences near topography.
1046	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
1047	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
1048	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) &
1049	) * 0.5
1050
1051	tend(k,j,i) = tend(k,j,i) - ( &
1052	( flux_r(k) + diss_r(k) &
1053	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
1054	+ ( flux_n(k) + diss_n(k) &
1055	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
1056	+ ( flux_t(k) + diss_t(k) &
1057	- flux_d - diss_d ) * ddzw(k) &
1058	) + div * u(k,j,i)
1059
1060	flux_l_u(k,j,tn) = flux_r(k)
1061	diss_l_u(k,j,tn) = diss_r(k)
1062	flux_s_u(k,tn) = flux_n(k)
1063	diss_s_u(k,tn) = diss_n(k)
1064	flux_d = flux_t(k)
1065	diss_d = diss_t(k)
1066	!
1067	!-- Statistical Evaluation of u'u'. The factor has to be applied for
1068	!-- right evaluation when gallilei_trans = .T. .
1069	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
1070	+ ( flux_r(k) * &
1071	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
1072	/ ( u_comp(k) - gu + 1.0E-20 ) &
1073	+ diss_r(k) * &
1074	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
1075	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
1076	* weight_substep(intermediate_timestep_count)
1077	!
1078	!-- Statistical Evaluation of w'u'.
1079	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
1080	+ ( flux_t(k) + diss_t(k) ) &
1081	* weight_substep(intermediate_timestep_count)
1082	ENDDO
1083
1084	DO k = nzb_max+1, nzt
1085
1086	u_comp(k) = u(k,j,i+1) + u(k,j,i)
1087	flux_r(k) = ( u_comp(k) - gu ) * ( &
1088	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
1089	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
1090	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
1091	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
1092	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
1093	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
1094	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
1095
1096	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
1097	flux_n(k) = v_comp * ( &
1098	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
1099	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
1100	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
1101	diss_n(k) = - ABS( v_comp ) * ( &
1102	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
1103	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
1104	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
1105	!
1106	!-- k index has to be modified near bottom and top, else array
1107	!-- subscripts will be exceeded.
1108	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
1109	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
1110	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
1111
1112	k_ppp = k + 3 * ibit17
1113	k_pp = k + 2 * ( 1 - ibit15 )
1114	k_mm = k - 2 * ibit17
1115
1116	w_comp = w(k,j,i) + w(k,j,i-1)
1117	flux_t(k) = w_comp * ( &
1118	( 37.0 * ibit17 * adv_mom_5 &
1119	+ 7.0 * ibit16 * adv_mom_3 &
1120	+ ibit15 * adv_mom_1 &
1121	) * &
1122	( u(k+1,j,i) + u(k,j,i) ) &
1123	- ( 8.0 * ibit17 * adv_mom_5 &
1124	+ ibit16 * adv_mom_3 &
1125	) * &
1126	( u(k_pp,j,i) + u(k-1,j,i) ) &
1127	+ ( ibit17 * adv_mom_5 &
1128	) * &
1129	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
1130	)
1131
1132	diss_t(k) = - ABS( w_comp ) * ( &
1133	( 10.0 * ibit17 * adv_mom_5 &
1134	+ 3.0 * ibit16 * adv_mom_3 &
1135	+ ibit15 * adv_mom_1 &
1136	) * &
1137	( u(k+1,j,i) - u(k,j,i) ) &
1138	- ( 5.0 * ibit17 * adv_mom_5 &
1139	+ ibit16 * adv_mom_3 &
1140	) * &
1141	( u(k_pp,j,i) - u(k-1,j,i) ) &
1142	+ ( ibit17 * adv_mom_5 &
1143	) * &
1144	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
1145	)
1146	!
1147	!-- Calculate the divergence of the velocity field. A respective
1148	!-- correction is needed to overcome numerical instabilities introduced
1149	!-- by a not sufficient reduction of divergences near topography.
1150	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
1151	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
1152	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) &
1153	) * 0.5
1154
1155	tend(k,j,i) = tend(k,j,i) - ( &
1156	( flux_r(k) + diss_r(k) &
1157	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
1158	+ ( flux_n(k) + diss_n(k) &
1159	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
1160	+ ( flux_t(k) + diss_t(k) &
1161	- flux_d - diss_d ) * ddzw(k) &
1162	) + div * u(k,j,i)
1163
1164	flux_l_u(k,j,tn) = flux_r(k)
1165	diss_l_u(k,j,tn) = diss_r(k)
1166	flux_s_u(k,tn) = flux_n(k)
1167	diss_s_u(k,tn) = diss_n(k)
1168	flux_d = flux_t(k)
1169	diss_d = diss_t(k)
1170	!
1171	!-- Statistical Evaluation of u'u'. The factor has to be applied for
1172	!-- right evaluation when gallilei_trans = .T. .
1173	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
1174	+ ( flux_r(k) * &
1175	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
1176	/ ( u_comp(k) - gu + 1.0E-20 ) &
1177	+ diss_r(k) * &
1178	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
1179	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
1180	* weight_substep(intermediate_timestep_count)
1181	!
1182	!-- Statistical Evaluation of w'u'.
1183	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
1184	+ ( flux_t(k) + diss_t(k) ) &
1185	* weight_substep(intermediate_timestep_count)
1186	ENDDO
1187
1188	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
1189
1190
1191
1192	END SUBROUTINE advec_u_ws_ij
1193
1194
1195
1196	!-----------------------------------------------------------------------------!
1197	! Advection of v-component - Call for grid point i,j
1198	!-----------------------------------------------------------------------------!
1199	SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
1200
1201	USE arrays_3d
1202	USE constants
1203	USE control_parameters
1204	USE grid_variables
1205	USE indices
1206	USE statistics
1207
1208	IMPLICIT NONE
1209
1210	INTEGER :: i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, &
1211	ibit25, ibit26, i_omp, j, k, k_mm, k_pp, k_ppp, tn
1212	REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp_l, w_comp
1213	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
1214	flux_t, v_comp
1215
1216	gu = 2.0 * u_gtrans
1217	gv = 2.0 * v_gtrans
1218
1219	!
1220	!-- Compute leftside fluxes for the respective boundary.
1221	IF ( i == i_omp ) THEN
1222
1223	DO k = nzb+1, nzb_max
1224
1225	ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
1226	ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
1227	ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
1228
1229	u_comp = u(k,j-1,i) + u(k,j,i) - gu
1230	flux_l_v(k,j,tn) = u_comp * ( &
1231	( 37.0 * ibit20 * adv_mom_5 &
1232	+ 7.0 * ibit19 * adv_mom_3 &
1233	+ ibit18 * adv_mom_1 &
1234	) * &
1235	( v(k,j,i) + v(k,j,i-1) ) &
1236	- ( 8.0 * ibit20 * adv_mom_5 &
1237	+ ibit19 * adv_mom_3 &
1238	) * &
1239	( v(k,j,i+1) + v(k,j,i-2) ) &
1240	+ ( ibit20 * adv_mom_5 &
1241	) * &
1242	( v(k,j,i+2) + v(k,j,i-3) ) &
1243	)
1244
1245	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
1246	( 10.0 * ibit20 * adv_mom_5 &
1247	+ 3.0 * ibit19 * adv_mom_3 &
1248	+ ibit18 * adv_mom_1 &
1249	) * &
1250	( v(k,j,i) - v(k,j,i-1) ) &
1251	- ( 5.0 * ibit20 * adv_mom_5 &
1252	+ ibit19 * adv_mom_3 &
1253	) * &
1254	( v(k,j,i+1) - v(k,j,i-2) ) &
1255	+ ( ibit20 * adv_mom_5 &
1256	) * &
1257	( v(k,j,i+2) - v(k,j,i-3) ) &
1258	)
1259
1260	ENDDO
1261
1262	DO k = nzb_max+1, nzt
1263
1264	u_comp = u(k,j-1,i) + u(k,j,i) - gu
1265	flux_l_v(k,j,tn) = u_comp * ( &
1266	37.0 * ( v(k,j,i) + v(k,j,i-1) ) &
1267	- 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) &
1268	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
1269	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
1270	10.0 * ( v(k,j,i) - v(k,j,i-1) ) &
1271	- 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) &
1272	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
1273
1274	ENDDO
1275
1276	ENDIF
1277	!
1278	!-- Compute southside fluxes for the respective boundary.
1279	IF ( j == nysv ) THEN
1280
1281	DO k = nzb+1, nzb_max
1282
1283	ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
1284	ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
1285	ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
1286
1287	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
1288	flux_s_v(k,tn) = v_comp_l * ( &
1289	( 37.0 * ibit23 * adv_mom_5 &
1290	+ 7.0 * ibit22 * adv_mom_3 &
1291	+ ibit21 * adv_mom_1 &
1292	) * &
1293	( v(k,j,i) + v(k,j-1,i) ) &
1294	- ( 8.0 * ibit23 * adv_mom_5 &
1295	+ ibit22 * adv_mom_3 &
1296	) * &
1297	( v(k,j+1,i) + v(k,j-2,i) ) &
1298	+ ( ibit23 * adv_mom_5 &
1299	) * &
1300	( v(k,j+2,i) + v(k,j-3,i) ) &
1301	)
1302
1303	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
1304	( 10.0 * ibit23 * adv_mom_5 &
1305	+ 3.0 * ibit22 * adv_mom_3 &
1306	+ ibit21 * adv_mom_1 &
1307	) * &
1308	( v(k,j,i) - v(k,j-1,i) ) &
1309	- ( 5.0 * ibit23 * adv_mom_5 &
1310	+ ibit22 * adv_mom_3 &
1311	) * &
1312	( v(k,j+1,i) - v(k,j-2,i) ) &
1313	+ ( ibit23 * adv_mom_5 &
1314	) * &
1315	( v(k,j+2,i) - v(k,j-3,i) ) &
1316	)
1317
1318	ENDDO
1319
1320	DO k = nzb_max+1, nzt
1321
1322	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
1323	flux_s_v(k,tn) = v_comp_l * ( &
1324	37.0 * ( v(k,j,i) + v(k,j-1,i) ) &
1325	- 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) &
1326	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
1327	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
1328	10.0 * ( v(k,j,i) - v(k,j-1,i) ) &
1329	- 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) &
1330	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
1331
1332	ENDDO
1333
1334	ENDIF
1335
1336	flux_t(0) = 0.0
1337	diss_t(0) = 0.0
1338	flux_d = 0.0
1339	diss_d = 0.0
1340	!
1341	!-- Now compute the fluxes and tendency terms for the horizontal and
1342	!-- verical parts.
1343	DO k = nzb+1, nzb_max
1344
1345	ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
1346	ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
1347	ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
1348
1349	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1350	flux_r(k) = u_comp * ( &
1351	( 37.0 * ibit20 * adv_mom_5 &
1352	+ 7.0 * ibit19 * adv_mom_3 &
1353	+ ibit18 * adv_mom_1 &
1354	) * &
1355	( v(k,j,i+1) + v(k,j,i) ) &
1356	- ( 8.0 * ibit20 * adv_mom_5 &
1357	+ ibit19 * adv_mom_3 &
1358	) * &
1359	( v(k,j,i+2) + v(k,j,i-1) ) &
1360	+ ( ibit20 * adv_mom_5 &
1361	) * &
1362	( v(k,j,i+3) + v(k,j,i-2) ) &
1363	)
1364
1365	diss_r(k) = - ABS( u_comp ) * ( &
1366	( 10.0 * ibit20 * adv_mom_5 &
1367	+ 3.0 * ibit19 * adv_mom_3 &
1368	+ ibit18 * adv_mom_1 &
1369	) * &
1370	( v(k,j,i+1) - v(k,j,i) ) &
1371	- ( 5.0 * ibit20 * adv_mom_5 &
1372	+ ibit19 * adv_mom_3 &
1373	) * &
1374	( v(k,j,i+2) - v(k,j,i-1) ) &
1375	+ ( ibit20 * adv_mom_5 &
1376	) * &
1377	( v(k,j,i+3) - v(k,j,i-2) ) &
1378	)
1379
1380	ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
1381	ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
1382	ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
1383
1384
1385	v_comp(k) = v(k,j+1,i) + v(k,j,i)
1386	flux_n(k) = ( v_comp(k) - gv ) * ( &
1387	( 37.0 * ibit23 * adv_mom_5 &
1388	+ 7.0 * ibit22 * adv_mom_3 &
1389	+ ibit21 * adv_mom_1 &
1390	) * &
1391	( v(k,j+1,i) + v(k,j,i) ) &
1392	- ( 8.0 * ibit23 * adv_mom_5 &
1393	+ ibit22 * adv_mom_3 &
1394	) * &
1395	( v(k,j+2,i) + v(k,j-1,i) ) &
1396	+ ( ibit23 * adv_mom_5 &
1397	) * &
1398	( v(k,j+3,i) + v(k,j-2,i) ) &
1399	)
1400
1401	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
1402	( 10.0 * ibit23 * adv_mom_5 &
1403	+ 3.0 * ibit22 * adv_mom_3 &
1404	+ ibit21 * adv_mom_1 &
1405	) * &
1406	( v(k,j+1,i) - v(k,j,i) ) &
1407	- ( 5.0 * ibit23 * adv_mom_5 &
1408	+ ibit22 * adv_mom_3 &
1409	) * &
1410	( v(k,j+2,i) - v(k,j-1,i) ) &
1411	+ ( ibit23 * adv_mom_5 &
1412	) * &
1413	( v(k,j+3,i) - v(k,j-2,i) ) &
1414	)
1415	!
1416	!-- k index has to be modified near bottom and top, else array
1417	!-- subscripts will be exceeded.
1418	ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
1419	ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
1420	ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
1421
1422	k_ppp = k + 3 * ibit26
1423	k_pp = k + 2 * ( 1 - ibit24 )
1424	k_mm = k - 2 * ibit26
1425
1426	w_comp = w(k,j-1,i) + w(k,j,i)
1427	flux_t(k) = w_comp * ( &
1428	( 37.0 * ibit26 * adv_mom_5 &
1429	+ 7.0 * ibit25 * adv_mom_3 &
1430	+ ibit24 * adv_mom_1 &
1431	) * &
1432	( v(k+1,j,i) + v(k,j,i) ) &
1433	- ( 8.0 * ibit26 * adv_mom_5 &
1434	+ ibit25 * adv_mom_3 &
1435	) * &
1436	( v(k_pp,j,i) + v(k-1,j,i) ) &
1437	+ ( ibit26 * adv_mom_5 &
1438	) * &
1439	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
1440	)
1441
1442	diss_t(k) = - ABS( w_comp ) * ( &
1443	( 10.0 * ibit26 * adv_mom_5 &
1444	+ 3.0 * ibit25 * adv_mom_3 &
1445	+ ibit24 * adv_mom_1 &
1446	) * &
1447	( v(k+1,j,i) - v(k,j,i) ) &
1448	- ( 5.0 * ibit26 * adv_mom_5 &
1449	+ ibit25 * adv_mom_3 &
1450	) * &
1451	( v(k_pp,j,i) - v(k-1,j,i) ) &
1452	+ ( ibit26 * adv_mom_5 &
1453	) * &
1454	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
1455	)
1456	!
1457	!-- Calculate the divergence of the velocity field. A respective
1458	!-- correction is needed to overcome numerical instabilities introduced
1459	!-- by a not sufficient reduction of divergences near topography.
1460	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
1461	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
1462	+ ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) &
1463	) * 0.5
1464
1465	tend(k,j,i) = tend(k,j,i) - ( &
1466	( flux_r(k) + diss_r(k) &
1467	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
1468	+ ( flux_n(k) + diss_n(k) &
1469	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
1470	+ ( flux_t(k) + diss_t(k) &
1471	- flux_d - diss_d ) * ddzw(k) &
1472	) + v(k,j,i) * div
1473
1474	flux_l_v(k,j,tn) = flux_r(k)
1475	diss_l_v(k,j,tn) = diss_r(k)
1476	flux_s_v(k,tn) = flux_n(k)
1477	diss_s_v(k,tn) = diss_n(k)
1478	flux_d = flux_t(k)
1479	diss_d = diss_t(k)
1480
1481	!
1482	!-- Statistical Evaluation of v'v'. The factor has to be applied for
1483	!-- right evaluation when gallilei_trans = .T. .
1484	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
1485	+ ( flux_n(k) &
1486	* ( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
1487	/ ( v_comp(k) - gv + 1.0E-20 ) &
1488	+ diss_n(k) &
1489	* ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
1490	/ ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) &
1491	* weight_substep(intermediate_timestep_count)
1492	!
1493	!-- Statistical Evaluation of w'v'.
1494	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
1495	+ ( flux_t(k) + diss_t(k) ) &
1496	* weight_substep(intermediate_timestep_count)
1497
1498	ENDDO
1499
1500	DO k = nzb_max+1, nzt
1501
1502	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1503	flux_r(k) = u_comp * ( &
1504	37.0 * ( v(k,j,i+1) + v(k,j,i) ) &
1505	- 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) &
1506	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1507
1508	diss_r(k) = - ABS( u_comp ) * ( &
1509	10.0 * ( v(k,j,i+1) - v(k,j,i) ) &
1510	- 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) &
1511	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1512
1513
1514	v_comp(k) = v(k,j+1,i) + v(k,j,i)
1515	flux_n(k) = ( v_comp(k) - gv ) * ( &
1516	37.0 * ( v(k,j+1,i) + v(k,j,i) ) &
1517	- 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) &
1518	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1519
1520	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
1521	10.0 * ( v(k,j+1,i) - v(k,j,i) ) &
1522	- 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) &
1523	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
1524	!
1525	!-- k index has to be modified near bottom and top, else array
1526	!-- subscripts will be exceeded.
1527	ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
1528	ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
1529	ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
1530
1531	k_ppp = k + 3 * ibit26
1532	k_pp = k + 2 * ( 1 - ibit24 )
1533	k_mm = k - 2 * ibit26
1534
1535	w_comp = w(k,j-1,i) + w(k,j,i)
1536	flux_t(k) = w_comp * ( &
1537	( 37.0 * ibit26 * adv_mom_5 &
1538	+ 7.0 * ibit25 * adv_mom_3 &
1539	+ ibit24 * adv_mom_1 &
1540	) * &
1541	( v(k+1,j,i) + v(k,j,i) ) &
1542	- ( 8.0 * ibit26 * adv_mom_5 &
1543	+ ibit25 * adv_mom_3 &
1544	) * &
1545	( v(k_pp,j,i) + v(k-1,j,i) ) &
1546	+ ( ibit26 * adv_mom_5 &
1547	) * &
1548	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
1549	)
1550
1551	diss_t(k) = - ABS( w_comp ) * ( &
1552	( 10.0 * ibit26 * adv_mom_5 &
1553	+ 3.0 * ibit25 * adv_mom_3 &
1554	+ ibit24 * adv_mom_1 &
1555	) * &
1556	( v(k+1,j,i) - v(k,j,i) ) &
1557	- ( 5.0 * ibit26 * adv_mom_5 &
1558	+ ibit25 * adv_mom_3 &
1559	) * &
1560	( v(k_pp,j,i) - v(k-1,j,i) ) &
1561	+ ( ibit26 * adv_mom_5 &
1562	) * &
1563	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
1564	)
1565	!
1566	!-- Calculate the divergence of the velocity field. A respective
1567	!-- correction is needed to overcome numerical instabilities introduced
1568	!-- by a not sufficient reduction of divergences near topography.
1569	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
1570	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
1571	+ ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) &
1572	) * 0.5
1573
1574	tend(k,j,i) = tend(k,j,i) - ( &
1575	( flux_r(k) + diss_r(k) &
1576	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
1577	+ ( flux_n(k) + diss_n(k) &
1578	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
1579	+ ( flux_t(k) + diss_t(k) &
1580	- flux_d - diss_d ) * ddzw(k) &
1581	) + v(k,j,i) * div
1582
1583	flux_l_v(k,j,tn) = flux_r(k)
1584	diss_l_v(k,j,tn) = diss_r(k)
1585	flux_s_v(k,tn) = flux_n(k)
1586	diss_s_v(k,tn) = diss_n(k)
1587	flux_d = flux_t(k)
1588	diss_d = diss_t(k)
1589
1590	!
1591	!-- Statistical Evaluation of v'v'. The factor has to be applied for
1592	!-- right evaluation when gallilei_trans = .T. .
1593	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
1594	+ ( flux_n(k) &
1595	* ( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
1596	/ ( v_comp(k) - gv + 1.0E-20 ) &
1597	+ diss_n(k) &
1598	* ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
1599	/ ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) &
1600	* weight_substep(intermediate_timestep_count)
1601	!
1602	!-- Statistical Evaluation of w'v'.
1603	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
1604	+ ( flux_t(k) + diss_t(k) ) &
1605	* weight_substep(intermediate_timestep_count)
1606
1607	ENDDO
1608	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
1609
1610
1611	END SUBROUTINE advec_v_ws_ij
1612
1613
1614
1615	!------------------------------------------------------------------------------!
1616	! Advection of w-component - Call for grid point i,j
1617	!------------------------------------------------------------------------------!
1618	SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
1619
1620	USE arrays_3d
1621	USE constants
1622	USE control_parameters
1623	USE grid_variables
1624	USE indices
1625	USE statistics
1626
1627	IMPLICIT NONE
1628
1629	INTEGER :: i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, &
1630	ibit34, ibit35, i_omp, j, k, k_mm, k_pp, k_ppp, tn
1631	REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp
1632	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
1633	flux_t
1634
1635	gu = 2.0 * u_gtrans
1636	gv = 2.0 * v_gtrans
1637
1638	!
1639	!-- Compute southside fluxes for the respective boundary.
1640	IF ( j == nys ) THEN
1641
1642	DO k = nzb+1, nzb_max
1643	ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
1644	ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
1645	ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
1646
1647	v_comp = v(k+1,j,i) + v(k,j,i) - gv
1648	flux_s_w(k,tn) = v_comp * ( &
1649	( 37.0 * ibit32 * adv_mom_5 &
1650	+ 7.0 * ibit31 * adv_mom_3 &
1651	+ ibit30 * adv_mom_1 &
1652	) * &
1653	( w(k,j,i) + w(k,j-1,i) ) &
1654	- ( 8.0 * ibit32 * adv_mom_5 &
1655	+ ibit31 * adv_mom_3 &
1656	) * &
1657	( w(k,j+1,i) + w(k,j-2,i) ) &
1658	+ ( ibit32 * adv_mom_5 &
1659	) * &
1660	( w(k,j+2,i) + w(k,j-3,i) ) &
1661	)
1662
1663	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
1664	( 10.0 * ibit32 * adv_mom_5 &
1665	+ 3.0 * ibit31 * adv_mom_3 &
1666	+ ibit30 * adv_mom_1 &
1667	) * &
1668	( w(k,j,i) - w(k,j-1,i) ) &
1669	- ( 5.0 * ibit32 * adv_mom_5 &
1670	+ ibit31 * adv_mom_3 &
1671	) * &
1672	( w(k,j+1,i) - w(k,j-2,i) ) &
1673	+ ( ibit32 * adv_mom_5 &
1674	) * &
1675	( w(k,j+2,i) - w(k,j-3,i) ) &
1676	)
1677
1678	ENDDO
1679
1680	DO k = nzb_max+1, nzt
1681
1682	v_comp = v(k+1,j,i) + v(k,j,i) - gv
1683	flux_s_w(k,tn) = v_comp * ( &
1684	37.0 * ( w(k,j,i) + w(k,j-1,i) ) &
1685	- 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) &
1686	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
1687	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
1688	10.0 * ( w(k,j,i) - w(k,j-1,i) ) &
1689	- 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) &
1690	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
1691
1692	ENDDO
1693
1694	ENDIF
1695	!
1696	!-- Compute leftside fluxes for the respective boundary.
1697	IF ( i == i_omp ) THEN
1698
1699	DO k = nzb+1, nzb_max
1700
1701	ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
1702	ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
1703	ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
1704
1705	u_comp = u(k+1,j,i) + u(k,j,i) - gu
1706	flux_l_w(k,j,tn) = u_comp * ( &
1707	( 37.0 * ibit29 * adv_mom_5 &
1708	+ 7.0 * ibit28 * adv_mom_3 &
1709	+ ibit27 * adv_mom_1 &
1710	) * &
1711	( w(k,j,i) + w(k,j,i-1) ) &
1712	- ( 8.0 * ibit29 * adv_mom_5 &
1713	+ ibit28 * adv_mom_3 &
1714	) * &
1715	( w(k,j,i+1) + w(k,j,i-2) ) &
1716	+ ( ibit29 * adv_mom_5 &
1717	) * &
1718	( w(k,j,i+2) + w(k,j,i-3) ) &
1719	)
1720
1721	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
1722	( 10.0 * ibit29 * adv_mom_5 &
1723	+ 3.0 * ibit28 * adv_mom_3 &
1724	+ ibit27 * adv_mom_1 &
1725	) * &
1726	( w(k,j,i) - w(k,j,i-1) ) &
1727	- ( 5.0 * ibit29 * adv_mom_5 &
1728	+ ibit28 * adv_mom_3 &
1729	) * &
1730	( w(k,j,i+1) - w(k,j,i-2) ) &
1731	+ ( ibit29 * adv_mom_5 &
1732	) * &
1733	( w(k,j,i+2) - w(k,j,i-3) ) &
1734	)
1735
1736	ENDDO
1737
1738	DO k = nzb_max+1, nzt
1739
1740	u_comp = u(k+1,j,i) + u(k,j,i) - gu
1741	flux_l_w(k,j,tn) = u_comp * ( &
1742	37.0 * ( w(k,j,i) + w(k,j,i-1) ) &
1743	- 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) &
1744	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
1745	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
1746	10.0 * ( w(k,j,i) - w(k,j,i-1) ) &
1747	- 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) &
1748	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
1749
1750	ENDDO
1751
1752	ENDIF
1753	!
1754	!-- The lower flux has to be calculated explicetely for the tendency at
1755	!-- the first w-level. For topography wall this is done implicitely by
1756	!-- wall_flags_0.
1757	k = nzb + 1
1758	w_comp = w(k,j,i) + w(k-1,j,i)
1759	flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
1760	diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
1761	flux_d = flux_t(0)
1762	diss_d = diss_t(0)
1763	!
1764	!-- Now compute the fluxes and tendency terms for the horizontal
1765	!-- and vertical parts.
1766	DO k = nzb+1, nzb_max
1767
1768	ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
1769	ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
1770	ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
1771
1772	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1773	flux_r(k) = u_comp * ( &
1774	( 37.0 * ibit29 * adv_mom_5 &
1775	+ 7.0 * ibit28 * adv_mom_3 &
1776	+ ibit27 * adv_mom_1 &
1777	) * &
1778	( w(k,j,i+1) + w(k,j,i) ) &
1779	- ( 8.0 * ibit29 * adv_mom_5 &
1780	+ ibit28 * adv_mom_3 &
1781	) * &
1782	( w(k,j,i+2) + w(k,j,i-1) ) &
1783	+ ( ibit29 * adv_mom_5 &
1784	) * &
1785	( w(k,j,i+3) + w(k,j,i-2) ) &
1786	)
1787
1788	diss_r(k) = - ABS( u_comp ) * ( &
1789	( 10.0 * ibit29 * adv_mom_5 &
1790	+ 3.0 * ibit28 * adv_mom_3 &
1791	+ ibit27 * adv_mom_1 &
1792	) * &
1793	( w(k,j,i+1) - w(k,j,i) ) &
1794	- ( 5.0 * ibit29 * adv_mom_5 &
1795	+ ibit28 * adv_mom_3 &
1796	) * &
1797	( w(k,j,i+2) - w(k,j,i-1) ) &
1798	+ ( ibit29 * adv_mom_5 &
1799	) * &
1800	( w(k,j,i+3) - w(k,j,i-2) ) &
1801	)
1802
1803	ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
1804	ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
1805	ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
1806
1807	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1808	flux_n(k) = v_comp * ( &
1809	( 37.0 * ibit32 * adv_mom_5 &
1810	+ 7.0 * ibit31 * adv_mom_3 &
1811	+ ibit30 * adv_mom_1 &
1812	) * &
1813	( w(k,j+1,i) + w(k,j,i) ) &
1814	- ( 8.0 * ibit32 * adv_mom_5 &
1815	+ ibit31 * adv_mom_3 &
1816	) * &
1817	( w(k,j+2,i) + w(k,j-1,i) ) &
1818	+ ( ibit32 * adv_mom_5 &
1819	) * &
1820	( w(k,j+3,i) + w(k,j-2,i) ) &
1821	)
1822
1823	diss_n(k) = - ABS( v_comp ) * ( &
1824	( 10.0 * ibit32 * adv_mom_5 &
1825	+ 3.0 * ibit31 * adv_mom_3 &
1826	+ ibit30 * adv_mom_1 &
1827	) * &
1828	( w(k,j+1,i) - w(k,j,i) ) &
1829	- ( 5.0 * ibit32 * adv_mom_5 &
1830	+ ibit31 * adv_mom_3 &
1831	) * &
1832	( w(k,j+2,i) - w(k,j-1,i) ) &
1833	+ ( ibit32 * adv_mom_5 &
1834	) * &
1835	( w(k,j+3,i) - w(k,j-2,i) ) &
1836	)
1837	!
1838	!-- k index has to be modified near bottom and top, else array
1839	!-- subscripts will be exceeded.
1840	ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
1841	ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
1842	ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
1843
1844	k_ppp = k + 3 * ibit35
1845	k_pp = k + 2 * ( 1 - ibit33 )
1846	k_mm = k - 2 * ibit35
1847
1848	w_comp = w(k+1,j,i) + w(k,j,i)
1849	flux_t(k) = w_comp * ( &
1850	( 37.0 * ibit35 * adv_mom_5 &
1851	+ 7.0 * ibit34 * adv_mom_3 &
1852	+ ibit33 * adv_mom_1 &
1853	) * &
1854	( w(k+1,j,i) + w(k,j,i) ) &
1855	- ( 8.0 * ibit35 * adv_mom_5 &
1856	+ ibit34 * adv_mom_3 &
1857	) * &
1858	( w(k_pp,j,i) + w(k-1,j,i) ) &
1859	+ ( ibit35 * adv_mom_5 &
1860	) * &
1861	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
1862	)
1863
1864	diss_t(k) = - ABS( w_comp ) * ( &
1865	( 10.0 * ibit35 * adv_mom_5 &
1866	+ 3.0 * ibit34 * adv_mom_3 &
1867	+ ibit33 * adv_mom_1 &
1868	) * &
1869	( w(k+1,j,i) - w(k,j,i) ) &
1870	- ( 5.0 * ibit35 * adv_mom_5 &
1871	+ ibit34 * adv_mom_3 &
1872	) * &
1873	( w(k_pp,j,i) - w(k-1,j,i) ) &
1874	+ ( ibit35 * adv_mom_5 &
1875	) * &
1876	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
1877	)
1878
1879	!
1880	!-- Calculate the divergence of the velocity field. A respective
1881	!-- correction is needed to overcome numerical instabilities introduced
1882	!-- by a not sufficient reduction of divergences near topography.
1883	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
1884	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
1885	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) &
1886	) * 0.5
1887
1888	tend(k,j,i) = tend(k,j,i) - ( &
1889	( flux_r(k) + diss_r(k) &
1890	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
1891	+ ( flux_n(k) + diss_n(k) &
1892	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
1893	+ ( flux_t(k) + diss_t(k) &
1894	- flux_d - diss_d ) * ddzu(k+1) &
1895	) + div * w(k,j,i)
1896
1897	flux_l_w(k,j,tn) = flux_r(k)
1898	diss_l_w(k,j,tn) = diss_r(k)
1899	flux_s_w(k,tn) = flux_n(k)
1900	diss_s_w(k,tn) = diss_n(k)
1901	flux_d = flux_t(k)
1902	diss_d = diss_t(k)
1903	!
1904	!-- Statistical Evaluation of w'w'.
1905	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
1906	+ ( flux_t(k) + diss_t(k) ) &
1907	* weight_substep(intermediate_timestep_count)
1908
1909	ENDDO
1910
1911	DO k = nzb_max+1, nzt
1912
1913	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1914	flux_r(k) = u_comp * ( &
1915	37.0 * ( w(k,j,i+1) + w(k,j,i) ) &
1916	- 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) &
1917	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1918
1919	diss_r(k) = - ABS( u_comp ) * ( &
1920	10.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1921	- 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) &
1922	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1923
1924	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1925	flux_n(k) = v_comp * ( &
1926	37.0 * ( w(k,j+1,i) + w(k,j,i) ) &
1927	- 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) &
1928	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1929
1930	diss_n(k) = - ABS( v_comp ) * ( &
1931	10.0 * ( w(k,j+1,i) - w(k,j,i) ) &
1932	- 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) &
1933	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1934	!
1935	!-- k index has to be modified near bottom and top, else array
1936	!-- subscripts will be exceeded.
1937	ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
1938	ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
1939	ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
1940
1941	k_ppp = k + 3 * ibit35
1942	k_pp = k + 2 * ( 1 - ibit33 )
1943	k_mm = k - 2 * ibit35
1944
1945	w_comp = w(k+1,j,i) + w(k,j,i)
1946	flux_t(k) = w_comp * ( &
1947	( 37.0 * ibit35 * adv_mom_5 &
1948	+ 7.0 * ibit34 * adv_mom_3 &
1949	+ ibit33 * adv_mom_1 &
1950	) * &
1951	( w(k+1,j,i) + w(k,j,i) ) &
1952	- ( 8.0 * ibit35 * adv_mom_5 &
1953	+ ibit34 * adv_mom_3 &
1954	) * &
1955	( w(k_pp,j,i) + w(k-1,j,i) ) &
1956	+ ( ibit35 * adv_mom_5 &
1957	) * &
1958	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
1959	)
1960
1961	diss_t(k) = - ABS( w_comp ) * ( &
1962	( 10.0 * ibit35 * adv_mom_5 &
1963	+ 3.0 * ibit34 * adv_mom_3 &
1964	+ ibit33 * adv_mom_1 &
1965	) * &
1966	( w(k+1,j,i) - w(k,j,i) ) &
1967	- ( 5.0 * ibit35 * adv_mom_5 &
1968	+ ibit34 * adv_mom_3 &
1969	) * &
1970	( w(k_pp,j,i) - w(k-1,j,i) ) &
1971	+ ( ibit35 * adv_mom_5 &
1972	) * &
1973	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
1974	)
1975	!
1976	!-- Calculate the divergence of the velocity field. A respective
1977	!-- correction is needed to overcome numerical instabilities introduced
1978	!-- by a not sufficient reduction of divergences near topography.
1979	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
1980	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
1981	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) &
1982	) * 0.5
1983
1984	tend(k,j,i) = tend(k,j,i) - ( &
1985	( flux_r(k) + diss_r(k) &
1986	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
1987	+ ( flux_n(k) + diss_n(k) &
1988	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
1989	+ ( flux_t(k) + diss_t(k) &
1990	- flux_d - diss_d ) * ddzu(k+1) &
1991	) + div * w(k,j,i)
1992
1993	flux_l_w(k,j,tn) = flux_r(k)
1994	diss_l_w(k,j,tn) = diss_r(k)
1995	flux_s_w(k,tn) = flux_n(k)
1996	diss_s_w(k,tn) = diss_n(k)
1997	flux_d = flux_t(k)
1998	diss_d = diss_t(k)
1999	!
2000	!-- Statistical Evaluation of w'w'.
2001	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
2002	+ ( flux_t(k) + diss_t(k) ) &
2003	* weight_substep(intermediate_timestep_count)
2004
2005	ENDDO
2006
2007
2008	END SUBROUTINE advec_w_ws_ij
2009
2010
2011	!------------------------------------------------------------------------------!
2012	! Scalar advection - Call for all grid points
2013	!------------------------------------------------------------------------------!
2014	SUBROUTINE advec_s_ws( sk, sk_char )
2015
2016	USE arrays_3d
2017	USE constants
2018	USE control_parameters
2019	USE grid_variables
2020	USE indices
2021	USE statistics
2022
2023	IMPLICIT NONE
2024
2025	INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, &
2026	ibit7, ibit8, j, k, k_mm, k_pp, k_ppp, tn = 0
2027	#if defined( __nopointer )
2028	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk
2029	#else
2030	REAL, DIMENSION(:,:,:), POINTER :: sk
2031	#endif
2032	REAL :: diss_d, div, flux_d, u_comp, v_comp
2033	REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
2034	flux_t
2035	REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local, swap_flux_y_local
2036	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local, &
2037	swap_flux_x_local
2038	CHARACTER (LEN = *), INTENT(IN) :: sk_char
2039
2040	!
2041	!-- Compute the fluxes for the whole left boundary of the processor domain.
2042	i = nxl
2043	DO j = nys, nyn
2044
2045	DO k = nzb+1, nzb_max
2046
2047	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
2048	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
2049	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
2050
2051	u_comp = u(k,j,i) - u_gtrans
2052	swap_flux_x_local(k,j) = u_comp * ( &
2053	( 37.0 * ibit2 * adv_sca_5 &
2054	+ 7.0 * ibit1 * adv_sca_3 &
2055	+ ibit0 * adv_sca_1 &
2056	) * &
2057	( sk(k,j,i) + sk(k,j,i-1) ) &
2058	- ( 8.0 * ibit2 * adv_sca_5 &
2059	+ ibit1 * adv_sca_3 &
2060	) * &
2061	( sk(k,j,i+1) + sk(k,j,i-2) ) &
2062	+ ( ibit2 * adv_sca_5 &
2063	) * &
2064	( sk(k,j,i+2) + sk(k,j,i-3) ) &
2065	)
2066
2067	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
2068	( 10.0 * ibit2 * adv_sca_5 &
2069	+ 3.0 * ibit1 * adv_sca_3 &
2070	+ ibit0 * adv_sca_1 &
2071	) * &
2072	( sk(k,j,i) - sk(k,j,i-1) ) &
2073	- ( 5.0 * ibit2 * adv_sca_5 &
2074	+ ibit1 * adv_sca_3 &
2075	) * &
2076	( sk(k,j,i+1) - sk(k,j,i-2) ) &
2077	+ ( ibit2 * adv_sca_5 &
2078	) * &
2079	( sk(k,j,i+2) - sk(k,j,i-3) ) &
2080	)
2081
2082	ENDDO
2083
2084	DO k = nzb_max+1, nzt
2085
2086	u_comp = u(k,j,i) - u_gtrans
2087	swap_flux_x_local(k,j) = u_comp * ( &
2088	37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) &
2089	- 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
2090	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
2091	) * adv_sca_5
2092
2093	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
2094	10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) &
2095	- 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
2096	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
2097	) * adv_sca_5
2098
2099	ENDDO
2100
2101	ENDDO
2102
2103	DO i = nxl, nxr
2104
2105	j = nys
2106	DO k = nzb+1, nzb_max
2107
2108	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
2109	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
2110	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
2111
2112	v_comp = v(k,j,i) - v_gtrans
2113	swap_flux_y_local(k) = v_comp * ( &
2114	( 37.0 * ibit5 * adv_sca_5 &
2115	+ 7.0 * ibit4 * adv_sca_3 &
2116	+ ibit3 * adv_sca_1 &
2117	) * &
2118	( sk(k,j,i) + sk(k,j-1,i) ) &
2119	- ( 8.0 * ibit5 * adv_sca_5 &
2120	+ ibit4 * adv_sca_3 &
2121	) * &
2122	( sk(k,j+1,i) + sk(k,j-2,i) ) &
2123	+ ( ibit5 * adv_sca_5 &
2124	) * &
2125	( sk(k,j+2,i) + sk(k,j-3,i) ) &
2126	)
2127
2128	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
2129	( 10.0 * ibit5 * adv_sca_5 &
2130	+ 3.0 * ibit4 * adv_sca_3 &
2131	+ ibit3 * adv_sca_1 &
2132	) * &
2133	( sk(k,j,i) - sk(k,j-1,i) ) &
2134	- ( 5.0 * ibit5 * adv_sca_5 &
2135	+ ibit4 * adv_sca_3 &
2136	) * &
2137	( sk(k,j+1,i) - sk(k,j-2,i) ) &
2138	+ ( ibit5 * adv_sca_5 &
2139	) * &
2140	( sk(k,j+2,i) - sk(k,j-3,i) ) &
2141	)
2142
2143	ENDDO
2144	!
2145	!-- Above to the top of the highest topography. No degradation necessary.
2146	DO k = nzb_max+1, nzt
2147
2148	v_comp = v(k,j,i) - v_gtrans
2149	swap_flux_y_local(k) = v_comp * ( &
2150	37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) &
2151	- 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
2152	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
2153	) * adv_sca_5
2154	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
2155	10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) &
2156	- 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
2157	+ sk(k,j+2,i) - sk(k,j-3,i) &
2158	) * adv_sca_5
2159
2160	ENDDO
2161
2162	DO j = nys, nyn
2163
2164	flux_t(0) = 0.0
2165	diss_t(0) = 0.0
2166	flux_d = 0.0
2167	diss_d = 0.0
2168
2169	DO k = nzb+1, nzb_max
2170
2171	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
2172	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
2173	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
2174
2175	u_comp = u(k,j,i+1) - u_gtrans
2176	flux_r(k) = u_comp * ( &
2177	( 37.0 * ibit2 * adv_sca_5 &
2178	+ 7.0 * ibit1 * adv_sca_3 &
2179	+ ibit0 * adv_sca_1 &
2180	) * &
2181	( sk(k,j,i+1) + sk(k,j,i) ) &
2182	- ( 8.0 * ibit2 * adv_sca_5 &
2183	+ ibit1 * adv_sca_3 &
2184	) * &
2185	( sk(k,j,i+2) + sk(k,j,i-1) ) &
2186	+ ( ibit2 * adv_sca_5 &
2187	) * &
2188	( sk(k,j,i+3) + sk(k,j,i-2) ) &
2189	)
2190
2191	diss_r(k) = -ABS( u_comp ) * ( &
2192	( 10.0 * ibit2 * adv_sca_5 &
2193	+ 3.0 * ibit1 * adv_sca_3 &
2194	+ ibit0 * adv_sca_1 &
2195	) * &
2196	( sk(k,j,i+1) - sk(k,j,i) ) &
2197	- ( 5.0 * ibit2 * adv_sca_5 &
2198	+ ibit1 * adv_sca_3 &
2199	) * &
2200	( sk(k,j,i+2) - sk(k,j,i-1) ) &
2201	+ ( ibit2 * adv_sca_5 &
2202	) * &
2203	( sk(k,j,i+3) - sk(k,j,i-2) ) &
2204	)
2205
2206	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
2207	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
2208	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
2209
2210	v_comp = v(k,j+1,i) - v_gtrans
2211	flux_n(k) = v_comp * ( &
2212	( 37.0 * ibit5 * adv_sca_5 &
2213	+ 7.0 * ibit4 * adv_sca_3 &
2214	+ ibit3 * adv_sca_1 &
2215	) * &
2216	( sk(k,j+1,i) + sk(k,j,i) ) &
2217	- ( 8.0 * ibit5 * adv_sca_5 &
2218	+ ibit4 * adv_sca_3 &
2219	) * &
2220	( sk(k,j+2,i) + sk(k,j-1,i) ) &
2221	+ ( ibit5 * adv_sca_5 &
2222	) * &
2223	( sk(k,j+3,i) + sk(k,j-2,i) ) &
2224	)
2225
2226	diss_n(k) = -ABS( v_comp ) * ( &
2227	( 10.0 * ibit5 * adv_sca_5 &
2228	+ 3.0 * ibit4 * adv_sca_3 &
2229	+ ibit3 * adv_sca_1 &
2230	) * &
2231	( sk(k,j+1,i) - sk(k,j,i) ) &
2232	- ( 5.0 * ibit5 * adv_sca_5 &
2233	+ ibit4 * adv_sca_3 &
2234	) * &
2235	( sk(k,j+2,i) - sk(k,j-1,i) ) &
2236	+ ( ibit5 * adv_sca_5 &
2237	) * &
2238	( sk(k,j+3,i) - sk(k,j-2,i) ) &
2239	)
2240	!
2241	!-- k index has to be modified near bottom and top, else array
2242	!-- subscripts will be exceeded.
2243	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
2244	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
2245	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
2246
2247	k_ppp = k + 3 * ibit8
2248	k_pp = k + 2 * ( 1 - ibit6 )
2249	k_mm = k - 2 * ibit8
2250
2251
2252	flux_t(k) = w(k,j,i) * ( &
2253	( 37.0 * ibit8 * adv_sca_5 &
2254	+ 7.0 * ibit7 * adv_sca_3 &
2255	+ ibit6 * adv_sca_1 &
2256	) * &
2257	( sk(k+1,j,i) + sk(k,j,i) ) &
2258	- ( 8.0 * ibit8 * adv_sca_5 &
2259	+ ibit7 * adv_sca_3 &
2260	) * &
2261	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
2262	+ ( ibit8 * adv_sca_5 &
2263	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
2264	)
2265
2266	diss_t(k) = -ABS( w(k,j,i) ) * ( &
2267	( 10.0 * ibit8 * adv_sca_5 &
2268	+ 3.0 * ibit7 * adv_sca_3 &
2269	+ ibit6 * adv_sca_1 &
2270	) * &
2271	( sk(k+1,j,i) - sk(k,j,i) ) &
2272	- ( 5.0 * ibit8 * adv_sca_5 &
2273	+ ibit7 * adv_sca_3 &
2274	) * &
2275	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
2276	+ ( ibit8 * adv_sca_5 &
2277	) * &
2278	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
2279	)
2280	!
2281	!-- Calculate the divergence of the velocity field. A respective
2282	!-- correction is needed to overcome numerical instabilities caused
2283	!-- by a not sufficient reduction of divergences near topography.
2284	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
2285	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
2286	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
2287
2288	tend(k,j,i) = tend(k,j,i) - ( &
2289	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
2290	swap_diss_x_local(k,j) ) * ddx &
2291	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
2292	swap_diss_y_local(k) ) * ddy &
2293	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
2294	) * ddzw(k) &
2295	) + sk(k,j,i) * div
2296
2297	swap_flux_y_local(k) = flux_n(k)
2298	swap_diss_y_local(k) = diss_n(k)
2299	swap_flux_x_local(k,j) = flux_r(k)
2300	swap_diss_x_local(k,j) = diss_r(k)
2301	flux_d = flux_t(k)
2302	diss_d = diss_t(k)
2303
2304	ENDDO
2305
2306	DO k = nzb_max+1, nzt
2307
2308	u_comp = u(k,j,i+1) - u_gtrans
2309	flux_r(k) = u_comp * ( &
2310	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
2311	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
2312	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
2313	diss_r(k) = -ABS( u_comp ) * ( &
2314	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
2315	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
2316	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
2317
2318	v_comp = v(k,j+1,i) - v_gtrans
2319	flux_n(k) = v_comp * ( &
2320	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
2321	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
2322	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
2323	diss_n(k) = -ABS( v_comp ) * ( &
2324	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
2325	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
2326	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
2327	!
2328	!-- k index has to be modified near bottom and top, else array
2329	!-- subscripts will be exceeded.
2330	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
2331	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
2332	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
2333
2334	k_ppp = k + 3 * ibit8
2335	k_pp = k + 2 * ( 1 - ibit6 )
2336	k_mm = k - 2 * ibit8
2337
2338
2339	flux_t(k) = w(k,j,i) * ( &
2340	( 37.0 * ibit8 * adv_sca_5 &
2341	+ 7.0 * ibit7 * adv_sca_3 &
2342	+ ibit6 * adv_sca_1 &
2343	) * &
2344	( sk(k+1,j,i) + sk(k,j,i) ) &
2345	- ( 8.0 * ibit8 * adv_sca_5 &
2346	+ ibit7 * adv_sca_3 &
2347	) * &
2348	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
2349	+ ( ibit8 * adv_sca_5 &
2350	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
2351	)
2352
2353	diss_t(k) = -ABS( w(k,j,i) ) * ( &
2354	( 10.0 * ibit8 * adv_sca_5 &
2355	+ 3.0 * ibit7 * adv_sca_3 &
2356	+ ibit6 * adv_sca_1 &
2357	) * &
2358	( sk(k+1,j,i) - sk(k,j,i) ) &
2359	- ( 5.0 * ibit8 * adv_sca_5 &
2360	+ ibit7 * adv_sca_3 &
2361	) * &
2362	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
2363	+ ( ibit8 * adv_sca_5 &
2364	) * &
2365	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
2366	)
2367	!
2368	!-- Calculate the divergence of the velocity field. A respective
2369	!-- correction is needed to overcome numerical instabilities introduced
2370	!-- by a not sufficient reduction of divergences near topography.
2371	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
2372	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
2373	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
2374
2375	tend(k,j,i) = tend(k,j,i) - ( &
2376	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
2377	swap_diss_x_local(k,j) ) * ddx &
2378	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
2379	swap_diss_y_local(k) ) * ddy &
2380	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
2381	) * ddzw(k) &
2382	) + sk(k,j,i) * div
2383
2384	swap_flux_y_local(k) = flux_n(k)
2385	swap_diss_y_local(k) = diss_n(k)
2386	swap_flux_x_local(k,j) = flux_r(k)
2387	swap_diss_x_local(k,j) = diss_r(k)
2388	flux_d = flux_t(k)
2389	diss_d = diss_t(k)
2390
2391	ENDDO
2392	!
2393	!-- evaluation of statistics
2394	SELECT CASE ( sk_char )
2395
2396	CASE ( 'pt' )
2397	DO k = nzb, nzt
2398	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
2399	+ ( flux_t(k) + diss_t(k) ) &
2400	* weight_substep(intermediate_timestep_count)
2401	ENDDO
2402	CASE ( 'sa' )
2403	DO k = nzb, nzt
2404	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
2405	+ ( flux_t(k) + diss_t(k) ) &
2406	* weight_substep(intermediate_timestep_count)
2407	ENDDO
2408	CASE ( 'q' )
2409	DO k = nzb, nzt
2410	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
2411	+ ( flux_t(k) + diss_t(k) ) &
2412	* weight_substep(intermediate_timestep_count)
2413	ENDDO
2414
2415	END SELECT
2416
2417	ENDDO
2418	ENDDO
2419
2420	END SUBROUTINE advec_s_ws
2421
2422
2423	!------------------------------------------------------------------------------!
2424	! Scalar advection - Call for all grid points - accelerator version
2425	!------------------------------------------------------------------------------!
2426	SUBROUTINE advec_s_ws_acc ( sk, sk_char )
2427
2428	USE arrays_3d
2429	USE constants
2430	USE control_parameters
2431	USE grid_variables
2432	USE indices
2433	USE statistics
2434
2435	IMPLICIT NONE
2436
2437	CHARACTER (LEN = *), INTENT(IN) :: sk_char
2438
2439	INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, &
2440	ibit7, ibit8, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
2441
2442	REAL :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, flux_d, &
2443	flux_l, flux_n, flux_r, flux_s, flux_t, u_comp, v_comp
2444
2445	REAL, INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk
2446
2447	!
2448	!-- Computation of fluxes and tendency terms
2449	!$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 )
2450	!$acc loop
2451	DO i = nxl, nxr
2452	DO j = nys, nyn
2453	!$acc loop vector( 32 )
2454	DO k = nzb+1, nzt
2455
2456	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
2457	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
2458	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
2459
2460	u_comp = u(k,j,i) - u_gtrans
2461	flux_l = u_comp * ( &
2462	( 37.0 * ibit2 * adv_sca_5 &
2463	+ 7.0 * ibit1 * adv_sca_3 &
2464	+ ibit0 * adv_sca_1 &
2465	) * &
2466	( sk(k,j,i) + sk(k,j,i-1) ) &
2467	- ( 8.0 * ibit2 * adv_sca_5 &
2468	+ ibit1 * adv_sca_3 &
2469	) * &
2470	( sk(k,j,i+1) + sk(k,j,i-2) ) &
2471	+ ( ibit2 * adv_sca_5 &
2472	) * &
2473	( sk(k,j,i+2) + sk(k,j,i-3) ) &
2474	)
2475
2476	diss_l = -ABS( u_comp ) * ( &
2477	( 10.0 * ibit2 * adv_sca_5 &
2478	+ 3.0 * ibit1 * adv_sca_3 &
2479	+ ibit0 * adv_sca_1 &
2480	) * &
2481	( sk(k,j,i) - sk(k,j,i-1) ) &
2482	- ( 5.0 * ibit2 * adv_sca_5 &
2483	+ ibit1 * adv_sca_3 &
2484	) * &
2485	( sk(k,j,i+1) - sk(k,j,i-2) ) &
2486	+ ( ibit2 * adv_sca_5 &
2487	) * &
2488	( sk(k,j,i+2) - sk(k,j,i-3) ) &
2489	)
2490
2491	u_comp = u(k,j,i+1) - u_gtrans
2492	flux_r = u_comp * ( &
2493	( 37.0 * ibit2 * adv_sca_5 &
2494	+ 7.0 * ibit1 * adv_sca_3 &
2495	+ ibit0 * adv_sca_1 &
2496	) * &
2497	( sk(k,j,i+1) + sk(k,j,i) ) &
2498	- ( 8.0 * ibit2 * adv_sca_5 &
2499	+ ibit1 * adv_sca_3 &
2500	) * &
2501	( sk(k,j,i+2) + sk(k,j,i-1) ) &
2502	+ ( ibit2 * adv_sca_5 &
2503	) * &
2504	( sk(k,j,i+3) + sk(k,j,i-2) ) &
2505	)
2506
2507	diss_r = -ABS( u_comp ) * ( &
2508	( 10.0 * ibit2 * adv_sca_5 &
2509	+ 3.0 * ibit1 * adv_sca_3 &
2510	+ ibit0 * adv_sca_1 &
2511	) * &
2512	( sk(k,j,i+1) - sk(k,j,i) ) &
2513	- ( 5.0 * ibit2 * adv_sca_5 &
2514	+ ibit1 * adv_sca_3 &
2515	) * &
2516	( sk(k,j,i+2) - sk(k,j,i-1) ) &
2517	+ ( ibit2 * adv_sca_5 &
2518	) * &
2519	( sk(k,j,i+3) - sk(k,j,i-2) ) &
2520	)
2521
2522	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
2523	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
2524	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
2525
2526	v_comp = v(k,j,i) - v_gtrans
2527	flux_s = v_comp * ( &
2528	( 37.0 * ibit5 * adv_sca_5 &
2529	+ 7.0 * ibit4 * adv_sca_3 &
2530	+ ibit3 * adv_sca_1 &
2531	) * &
2532	( sk(k,j,i) + sk(k,j-1,i) ) &
2533	- ( 8.0 * ibit5 * adv_sca_5 &
2534	+ ibit4 * adv_sca_3 &
2535	) * &
2536	( sk(k,j+1,i) + sk(k,j-2,i) ) &
2537	+ ( ibit5 * adv_sca_5 &
2538	) * &
2539	( sk(k,j+2,i) + sk(k,j-3,i) ) &
2540	)
2541
2542	diss_s = -ABS( v_comp ) * ( &
2543	( 10.0 * ibit5 * adv_sca_5 &
2544	+ 3.0 * ibit4 * adv_sca_3 &
2545	+ ibit3 * adv_sca_1 &
2546	) * &
2547	( sk(k,j,i) - sk(k,j-1,i) ) &
2548	- ( 5.0 * ibit5 * adv_sca_5 &
2549	+ ibit4 * adv_sca_3 &
2550	) * &
2551	( sk(k,j+1,i) - sk(k,j-2,i) ) &
2552	+ ( ibit5 * adv_sca_5 &
2553	) * &
2554	( sk(k,j+2,i) - sk(k,j-3,i) ) &
2555	)
2556
2557
2558	v_comp = v(k,j+1,i) - v_gtrans
2559	flux_n = v_comp * ( &
2560	( 37.0 * ibit5 * adv_sca_5 &
2561	+ 7.0 * ibit4 * adv_sca_3 &
2562	+ ibit3 * adv_sca_1 &
2563	) * &
2564	( sk(k,j+1,i) + sk(k,j,i) ) &
2565	- ( 8.0 * ibit5 * adv_sca_5 &
2566	+ ibit4 * adv_sca_3 &
2567	) * &
2568	( sk(k,j+2,i) + sk(k,j-1,i) ) &
2569	+ ( ibit5 * adv_sca_5 &
2570	) * &
2571	( sk(k,j+3,i) + sk(k,j-2,i) ) &
2572	)
2573
2574	diss_n = -ABS( v_comp ) * ( &
2575	( 10.0 * ibit5 * adv_sca_5 &
2576	+ 3.0 * ibit4 * adv_sca_3 &
2577	+ ibit3 * adv_sca_1 &
2578	) * &
2579	( sk(k,j+1,i) - sk(k,j,i) ) &
2580	- ( 5.0 * ibit5 * adv_sca_5 &
2581	+ ibit4 * adv_sca_3 &
2582	) * &
2583	( sk(k,j+2,i) - sk(k,j-1,i) ) &
2584	+ ( ibit5 * adv_sca_5 &
2585	) * &
2586	( sk(k,j+3,i) - sk(k,j-2,i) ) &
2587	)
2588
2589	!
2590	!-- indizes k_m, k_mm, ... should be known at these point
2591	ibit8 = IBITS(wall_flags_0(k-1,j,i),8,1)
2592	ibit7 = IBITS(wall_flags_0(k-1,j,i),7,1)
2593	ibit6 = IBITS(wall_flags_0(k-1,j,i),6,1)
2594
2595	k_pp = k + 2 * ibit8
2596	k_mm = k - 2 * ( ibit7 + ibit8 )
2597	k_mmm = k - 3 * ibit8
2598
2599	flux_d = w(k-1,j,i) * ( &
2600	( 37.0 * ibit8 * adv_sca_5 &
2601	+ 7.0 * ibit7 * adv_sca_3 &
2602	+ ibit6 * adv_sca_1 &
2603	) * &
2604	( sk(k,j,i) + sk(k-1,j,i) ) &
2605	- ( 8.0 * ibit8 * adv_sca_5 &
2606	+ ibit7 * adv_sca_3 &
2607	) * &
2608	( sk(k+1,j,i) + sk(k_mm,j,i) ) &
2609	+ ( ibit8 * adv_sca_5 &
2610	) * ( sk(k_pp,j,i)+ sk(k_mmm,j,i) ) &
2611	)
2612
2613	diss_d = -ABS( w(k-1,j,i) ) * ( &
2614	( 10.0 * ibit8 * adv_sca_5 &
2615	+ 3.0 * ibit7 * adv_sca_3 &
2616	+ ibit6 * adv_sca_1 &
2617	) * &
2618	( sk(k,j,i) - sk(k-1,j,i) ) &
2619	- ( 5.0 * ibit8 * adv_sca_5 &
2620	+ ibit7 * adv_sca_3 &
2621	) * &
2622	( sk(k+1,j,i) - sk(k_mm,j,i) ) &
2623	+ ( ibit8 * adv_sca_5 &
2624	) * &
2625	( sk(k_pp,j,i) - sk(k_mmm,j,i) ) &
2626	)
2627
2628	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
2629	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
2630	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
2631
2632	k_ppp = k + 3 * ibit8
2633	k_pp = k + 2 * ( 1 - ibit6 )
2634	k_mm = k - 2 * ibit8
2635
2636	flux_t = w(k,j,i) * ( &
2637	( 37.0 * ibit8 * adv_sca_5 &
2638	+ 7.0 * ibit7 * adv_sca_3 &
2639	+ ibit6 * adv_sca_1 &
2640	) * &
2641	( sk(k+1,j,i) + sk(k,j,i) ) &
2642	- ( 8.0 * ibit8 * adv_sca_5 &
2643	+ ibit7 * adv_sca_3 &
2644	) * &
2645	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
2646	+ ( ibit8 * adv_sca_5 &
2647	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
2648	)
2649
2650	diss_t = -ABS( w(k,j,i) ) * ( &
2651	( 10.0 * ibit8 * adv_sca_5 &
2652	+ 3.0 * ibit7 * adv_sca_3 &
2653	+ ibit6 * adv_sca_1 &
2654	) * &
2655	( sk(k+1,j,i) - sk(k,j,i) ) &
2656	- ( 5.0 * ibit8 * adv_sca_5 &
2657	+ ibit7 * adv_sca_3 &
2658	) * &
2659	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
2660	+ ( ibit8 * adv_sca_5 &
2661	) * &
2662	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
2663	)
2664	!
2665	!-- Calculate the divergence of the velocity field. A respective
2666	!-- correction is needed to overcome numerical instabilities caused
2667	!-- by a not sufficient reduction of divergences near topography.
2668	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
2669	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
2670	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
2671
2672	tend(k,j,i) = - ( &
2673	( flux_r + diss_r - flux_l - diss_l ) * ddx &
2674	+ ( flux_n + diss_n - flux_s - diss_s ) * ddy &
2675	+ ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)&
2676	) + div * sk(k,j,i)
2677
2678	!++
2679	!-- Evaluation of statistics
2680	! SELECT CASE ( sk_char )
2681	!
2682	! CASE ( 'pt' )
2683	! sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
2684	! + ( flux_t + diss_t ) &
2685	! * weight_substep(intermediate_timestep_count)
2686	! CASE ( 'sa' )
2687	! sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
2688	! + ( flux_t + diss_t ) &
2689	! * weight_substep(intermediate_timestep_count)
2690	! CASE ( 'q' )
2691	! sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
2692	! + ( flux_t + diss_t ) &
2693	! * weight_substep(intermediate_timestep_count)
2694	!
2695	! END SELECT
2696
2697	ENDDO
2698	ENDDO
2699	ENDDO
2700	!$acc end kernels
2701
2702	END SUBROUTINE advec_s_ws_acc
2703
2704
2705	!------------------------------------------------------------------------------!
2706	! Advection of u - Call for all grid points
2707	!------------------------------------------------------------------------------!
2708	SUBROUTINE advec_u_ws
2709
2710	USE arrays_3d
2711	USE constants
2712	USE control_parameters
2713	USE grid_variables
2714	USE indices
2715	USE statistics
2716
2717	IMPLICIT NONE
2718
2719	INTEGER :: i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15, &
2720	ibit16, ibit17, j, k, k_mm, k_pp, k_ppp, tn = 0
2721	REAL :: diss_d, div, flux_d, gu, gv, v_comp, w_comp
2722	REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u, swap_flux_y_local_u
2723	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u, &
2724	swap_flux_x_local_u
2725	REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
2726	flux_t, u_comp
2727
2728	gu = 2.0 * u_gtrans
2729	gv = 2.0 * v_gtrans
2730
2731	!
2732	!-- Compute the fluxes for the whole left boundary of the processor domain.
2733	i = nxlu
2734	DO j = nys, nyn
2735	DO k = nzb+1, nzb_max
2736
2737	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
2738	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
2739	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
2740
2741	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
2742	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
2743	( 37.0 * ibit11 * adv_mom_5 &
2744	+ 7.0 * ibit10 * adv_mom_3 &
2745	+ ibit9 * adv_mom_1 &
2746	) * &
2747	( u(k,j,i) + u(k,j,i-1) ) &
2748	- ( 8.0 * ibit11 * adv_mom_5 &
2749	+ ibit10 * adv_mom_3 &
2750	) * &
2751	( u(k,j,i+1) + u(k,j,i-2) ) &
2752	+ ( ibit11 * adv_mom_5 &
2753	) * &
2754	( u(k,j,i+2) + u(k,j,i-3) ) &
2755	)
2756
2757	swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( &
2758	( 10.0 * ibit11 * adv_mom_5 &
2759	+ 3.0 * ibit10 * adv_mom_3 &
2760	+ ibit9 * adv_mom_1 &
2761	) * &
2762	( u(k,j,i) - u(k,j,i-1) ) &
2763	- ( 5.0 * ibit11 * adv_mom_5 &
2764	+ ibit10 * adv_mom_3 &
2765	) * &
2766	( u(k,j,i+1) - u(k,j,i-2) ) &
2767	+ ( ibit11 * adv_mom_5 &
2768	) * &
2769	( u(k,j,i+2) - u(k,j,i-3) ) &
2770	)
2771
2772	ENDDO
2773
2774	DO k = nzb_max+1, nzt
2775
2776	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
2777	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
2778	37.0 * ( u(k,j,i) + u(k,j,i-1) ) &
2779	- 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) &
2780	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
2781	swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( &
2782	10.0 * ( u(k,j,i) - u(k,j,i-1) ) &
2783	- 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) &
2784	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
2785
2786	ENDDO
2787	ENDDO
2788
2789	DO i = nxlu, nxr
2790	!
2791	!-- The following loop computes the fluxes for the south boundary points
2792	j = nys
2793	DO k = nzb+1, nzb_max
2794
2795	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
2796	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
2797	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
2798
2799	v_comp = v(k,j,i) + v(k,j,i-1) - gv
2800	swap_flux_y_local_u(k) = v_comp * ( &
2801	( 37.0 * ibit14 * adv_mom_5 &
2802	+ 7.0 * ibit13 * adv_mom_3 &
2803	+ ibit12 * adv_mom_1 &
2804	) * &
2805	( u(k,j,i) + u(k,j-1,i) ) &
2806	- ( 8.0 * ibit14 * adv_mom_5 &
2807	+ ibit13 * adv_mom_3 &
2808	) * &
2809	( u(k,j+1,i) + u(k,j-2,i) ) &
2810	+ ( ibit14 * adv_mom_5 &
2811	) * &
2812	( u(k,j+2,i) + u(k,j-3,i) ) &
2813	)
2814
2815	swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( &
2816	( 10.0 * ibit14 * adv_mom_5 &
2817	+ 3.0 * ibit13 * adv_mom_3 &
2818	+ ibit12 * adv_mom_1 &
2819	) * &
2820	( u(k,j,i) - u(k,j-1,i) ) &
2821	- ( 5.0 * ibit14 * adv_mom_5 &
2822	+ ibit13 * adv_mom_3 &
2823	) * &
2824	( u(k,j+1,i) - u(k,j-2,i) ) &
2825	+ ( ibit14 * adv_mom_5 &
2826	) * &
2827	( u(k,j+2,i) - u(k,j-3,i) ) &
2828	)
2829
2830	ENDDO
2831
2832	DO k = nzb_max+1, nzt
2833
2834	v_comp = v(k,j,i) + v(k,j,i-1) - gv
2835	swap_flux_y_local_u(k) = v_comp * ( &
2836	37.0 * ( u(k,j,i) + u(k,j-1,i) ) &
2837	- 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) &
2838	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
2839	swap_diss_y_local_u(k) = - ABS(v_comp) * ( &
2840	10.0 * ( u(k,j,i) - u(k,j-1,i) ) &
2841	- 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) &
2842	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
2843
2844	ENDDO
2845	!
2846	!-- Computation of interior fluxes and tendency terms
2847	DO j = nys, nyn
2848
2849	flux_t(0) = 0.0
2850	diss_t(0) = 0.0
2851	flux_d = 0.0
2852	diss_d = 0.0
2853
2854	DO k = nzb+1, nzb_max
2855
2856	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
2857	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
2858	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
2859
2860	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2861	flux_r(k) = ( u_comp(k) - gu ) * ( &
2862	( 37.0 * ibit11 * adv_mom_5 &
2863	+ 7.0 * ibit10 * adv_mom_3 &
2864	+ ibit9 * adv_mom_1 &
2865	) * &
2866	( u(k,j,i+1) + u(k,j,i) ) &
2867	- ( 8.0 * ibit11 * adv_mom_5 &
2868	+ ibit10 * adv_mom_3 &
2869	) * &
2870	( u(k,j,i+2) + u(k,j,i-1) ) &
2871	+ ( ibit11 * adv_mom_5 &
2872	) * &
2873	( u(k,j,i+3) + u(k,j,i-2) ) &
2874	)
2875
2876	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2877	( 10.0 * ibit11 * adv_mom_5 &
2878	+ 3.0 * ibit10 * adv_mom_3 &
2879	+ ibit9 * adv_mom_1 &
2880	) * &
2881	( u(k,j,i+1) - u(k,j,i) ) &
2882	- ( 5.0 * ibit11 * adv_mom_5 &
2883	+ ibit10 * adv_mom_3 &
2884	) * &
2885	( u(k,j,i+2) - u(k,j,i-1) ) &
2886	+ ( ibit11 * adv_mom_5 &
2887	) * &
2888	( u(k,j,i+3) - u(k,j,i-2) ) &
2889	)
2890
2891	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
2892	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
2893	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
2894
2895	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2896	flux_n(k) = v_comp * ( &
2897	( 37.0 * ibit14 * adv_mom_5 &
2898	+ 7.0 * ibit13 * adv_mom_3 &
2899	+ ibit12 * adv_mom_1 &
2900	) * &
2901	( u(k,j+1,i) + u(k,j,i) ) &
2902	- ( 8.0 * ibit14 * adv_mom_5 &
2903	+ ibit13 * adv_mom_3 &
2904	) * &
2905	( u(k,j+2,i) + u(k,j-1,i) ) &
2906	+ ( ibit14 * adv_mom_5 &
2907	) * &
2908	( u(k,j+3,i) + u(k,j-2,i) ) &
2909	)
2910
2911	diss_n(k) = - ABS ( v_comp ) * ( &
2912	( 10.0 * ibit14 * adv_mom_5 &
2913	+ 3.0 * ibit13 * adv_mom_3 &
2914	+ ibit12 * adv_mom_1 &
2915	) * &
2916	( u(k,j+1,i) - u(k,j,i) ) &
2917	- ( 5.0 * ibit14 * adv_mom_5 &
2918	+ ibit13 * adv_mom_3 &
2919	) * &
2920	( u(k,j+2,i) - u(k,j-1,i) ) &
2921	+ ( ibit14 * adv_mom_5 &
2922	) * &
2923	( u(k,j+3,i) - u(k,j-2,i) ) &
2924	)
2925	!
2926	!-- k index has to be modified near bottom and top, else array
2927	!-- subscripts will be exceeded.
2928	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
2929	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
2930	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
2931
2932	k_ppp = k + 3 * ibit17
2933	k_pp = k + 2 * ( 1 - ibit15 )
2934	k_mm = k - 2 * ibit17
2935
2936	w_comp = w(k,j,i) + w(k,j,i-1)
2937	flux_t(k) = w_comp * ( &
2938	( 37.0 * ibit17 * adv_mom_5 &
2939	+ 7.0 * ibit16 * adv_mom_3 &
2940	+ ibit15 * adv_mom_1 &
2941	) * &
2942	( u(k+1,j,i) + u(k,j,i) ) &
2943	- ( 8.0 * ibit17 * adv_mom_5 &
2944	+ ibit16 * adv_mom_3 &
2945	) * &
2946	( u(k_pp,j,i) + u(k-1,j,i) ) &
2947	+ ( ibit17 * adv_mom_5 &
2948	) * &
2949	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2950	)
2951
2952	diss_t(k) = - ABS( w_comp ) * ( &
2953	( 10.0 * ibit17 * adv_mom_5 &
2954	+ 3.0 * ibit16 * adv_mom_3 &
2955	+ ibit15 * adv_mom_1 &
2956	) * &
2957	( u(k+1,j,i) - u(k,j,i) ) &
2958	- ( 5.0 * ibit17 * adv_mom_5 &
2959	+ ibit16 * adv_mom_3 &
2960	) * &
2961	( u(k_pp,j,i) - u(k-1,j,i) ) &
2962	+ ( ibit17 * adv_mom_5 &
2963	) * &
2964	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2965	)
2966	!
2967	!-- Calculate the divergence of the velocity field. A respective
2968	!-- correction is needed to overcome numerical instabilities caused
2969	!-- by a not sufficient reduction of divergences near topography.
2970	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
2971	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
2972	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) &
2973	* ddzw(k) &
2974	) * 0.5
2975
2976	tend(k,j,i) = tend(k,j,i) - ( &
2977	( flux_r(k) + diss_r(k) &
2978	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
2979	+ ( flux_n(k) + diss_n(k) &
2980	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
2981	+ ( flux_t(k) + diss_t(k) &
2982	- flux_d - diss_d &
2983	) * ddzw(k) &
2984	) + div * u(k,j,i)
2985
2986	swap_flux_x_local_u(k,j) = flux_r(k)
2987	swap_diss_x_local_u(k,j) = diss_r(k)
2988	swap_flux_y_local_u(k) = flux_n(k)
2989	swap_diss_y_local_u(k) = diss_n(k)
2990	flux_d = flux_t(k)
2991	diss_d = diss_t(k)
2992	!
2993	!-- Statistical Evaluation of u'u'. The factor has to be applied
2994	!-- for right evaluation when gallilei_trans = .T. .
2995	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2996	+ ( flux_r(k) * &
2997	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
2998	/ ( u_comp(k) - gu + 1.0E-20 ) &
2999	+ diss_r(k) * &
3000	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
3001	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
3002	* weight_substep(intermediate_timestep_count)
3003	!
3004	!-- Statistical Evaluation of w'u'.
3005	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
3006	+ ( flux_t(k) + diss_t(k) ) &
3007	* weight_substep(intermediate_timestep_count)
3008	ENDDO
3009
3010	DO k = nzb_max+1, nzt
3011
3012	u_comp(k) = u(k,j,i+1) + u(k,j,i)
3013	flux_r(k) = ( u_comp(k) - gu ) * ( &
3014	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
3015	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
3016	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
3017	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
3018	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
3019	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
3020	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
3021
3022	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
3023	flux_n(k) = v_comp * ( &
3024	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
3025	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
3026	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
3027	diss_n(k) = - ABS( v_comp ) * ( &
3028	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
3029	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
3030	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
3031	!
3032	!-- k index has to be modified near bottom and top, else array
3033	!-- subscripts will be exceeded.
3034	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
3035	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
3036	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
3037
3038	k_ppp = k + 3 * ibit17
3039	k_pp = k + 2 * ( 1 - ibit15 )
3040	k_mm = k - 2 * ibit17
3041
3042	w_comp = w(k,j,i) + w(k,j,i-1)
3043	flux_t(k) = w_comp * ( &
3044	( 37.0 * ibit17 * adv_mom_5 &
3045	+ 7.0 * ibit16 * adv_mom_3 &
3046	+ ibit15 * adv_mom_1 &
3047	) * &
3048	( u(k+1,j,i) + u(k,j,i) ) &
3049	- ( 8.0 * ibit17 * adv_mom_5 &
3050	+ ibit16 * adv_mom_3 &
3051	) * &
3052	( u(k_pp,j,i) + u(k-1,j,i) ) &
3053	+ ( ibit17 * adv_mom_5 &
3054	) * &
3055	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
3056	)
3057
3058	diss_t(k) = - ABS( w_comp ) * ( &
3059	( 10.0 * ibit17 * adv_mom_5 &
3060	+ 3.0 * ibit16 * adv_mom_3 &
3061	+ ibit15 * adv_mom_1 &
3062	) * &
3063	( u(k+1,j,i) - u(k,j,i) ) &
3064	- ( 5.0 * ibit17 * adv_mom_5 &
3065	+ ibit16 * adv_mom_3 &
3066	) * &
3067	( u(k_pp,j,i) - u(k-1,j,i) ) &
3068	+ ( ibit17 * adv_mom_5 &
3069	) * &
3070	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
3071	)
3072	!
3073	!-- Calculate the divergence of