Home

Context Navigation

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

Last change on this file since 873 was 863, checked in by suehring, 13 years ago
last commit documented
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 191.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |