Home

Context Navigation

advec_ws.f90 @ 861

Last change on this file since 861 was 861, checked in by suehring, 12 years ago
WS5 is available in combination with topography. Version number changed from 3.8 to 3.8a. Modification in init_pt_anomaly.
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 190.6 KB

Line
1	MODULE advec_ws
2
3	!-----------------------------------------------------------------------------!
4	! Current revisions:
5	! ------------------
6	! ws-scheme also work with topography in combination with vector version.
7	! ws-scheme also work with outflow boundaries in combination with
8	! vector version.
9
10	! Degradation of the applied order of scheme is now steered by multiplying with
11	! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each
12	! grid point and mulitplied with the appropriate flag.
13	! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
14	! term derived according to the 4th and 6th order terms is applied. It turns
15	! out that diss_2nd does not provide sufficient dissipation near walls.
16	! Therefore, the function diss_2nd is removed.
17
18	! Near walls a divergence correction is necessary to overcome numerical
19	! instabilities due to too less divergence reduction of the velocity field.
20	!
21	! boundary_flags and logicals steering the degradation are removed.
22	! Empty SUBROUTINE local_diss removed.
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: advec_ws.f90 861 2012-03-26 14:18:34Z suehring $
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, tn, k_ppp, k_pp, k_mm
277	REAL :: flux_d, diss_d, u_comp, v_comp, div
278	REAL, DIMENSION(:,:,:), POINTER :: sk
279	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, &
280	flux_n, diss_n
281	REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local, &
282	swap_diss_y_local
283	REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: &
284	swap_flux_x_local, &
285	swap_diss_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, tn, k_ppp, k_pp, k_mm
674	REAL :: gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp, div
675	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, &
676	flux_n, diss_n, 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, tn, k_ppp, k_pp, k_mm
1061	REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp, div
1062	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_n, &
1063	diss_n, flux_r, diss_r, 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, tn, k_ppp, k_pp, k_mm
1454	REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp, div
1455	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n, &
1456	diss_n
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	flux_t(0) = 0.0
1571	diss_t(0) = 0.0
1572	flux_d = 0.0
1573	diss_d = 0.0
1574	!
1575	!-- Now compute the fluxes and tendency terms for the horizontal
1576	!-- and vertical parts.
1577	DO k = nzb+1, nzb_max
1578
1579	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1580	flux_r(k) = u_comp * ( &
1581	( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
1582	+ 7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
1583	+ IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 &
1584	) * &
1585	( w(k,j,i+1) + w(k,j,i) ) &
1586	- ( 8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
1587	+ IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
1588	) * &
1589	( w(k,j,i+2) + w(k,j,i-1) ) &
1590	+ ( IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
1591	) * &
1592	( w(k,j,i+3) + w(k,j,i-2) ) &
1593	)
1594
1595	diss_r(k) = - ABS( u_comp ) * ( &
1596	( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
1597	+ 3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
1598	+ IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 &
1599	) * &
1600	( w(k,j,i+1) - w(k,j,i) ) &
1601	- ( 5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
1602	+ IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
1603	) * &
1604	( w(k,j,i+2) - w(k,j,i-1) ) &
1605	+ ( IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
1606	) * &
1607	( w(k,j,i+3) - w(k,j,i-2) ) &
1608	)
1609
1610	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1611	flux_n(k) = v_comp * ( &
1612	( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
1613	+ 7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
1614	+ IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 &
1615	) * &
1616	( w(k,j+1,i) + w(k,j,i) ) &
1617	- ( 8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
1618	+ IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
1619	) * &
1620	( w(k,j+2,i) + w(k,j-1,i) ) &
1621	+ ( IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
1622	) * &
1623	( w(k,j+3,i) + w(k,j-2,i) ) &
1624	)
1625
1626	diss_n(k) = - ABS( v_comp ) * ( &
1627	( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
1628	+ 3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
1629	+ IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 &
1630	) * &
1631	( w(k,j+1,i) - w(k,j,i) ) &
1632	- ( 5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
1633	+ IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
1634	) * &
1635	( w(k,j+2,i) - w(k,j-1,i) ) &
1636	+ ( IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
1637	) * &
1638	( w(k,j+3,i) - w(k,j-2,i) ) &
1639	)
1640	!
1641	!-- k index has to be modified near bottom and top, else array
1642	!-- subscripts will be exceeded.
1643	k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1)
1644	k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1) )
1645	k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),35,1)
1646
1647	w_comp = w(k+1,j,i) + w(k,j,i)
1648	flux_t(k) = w_comp * ( &
1649	( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1650	+ 7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1651	+ IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 &
1652	) * &
1653	( w(k+1,j,i) + w(k,j,i) ) &
1654	- ( 8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1655	+ IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1656	) * &
1657	( w(k_pp,j,i) + w(k-1,j,i) ) &
1658	+ ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1659	) * &
1660	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
1661	)
1662
1663	diss_t(k) = - ABS( w_comp ) * ( &
1664	( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1665	+ 3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1666	+ IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 &
1667	) * &
1668	( w(k+1,j,i) - w(k,j,i) ) &
1669	- ( 5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1670	+ IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1671	) * &
1672	( w(k_pp,j,i) - w(k-1,j,i) ) &
1673	+ ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1674	) * &
1675	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
1676	)
1677	!
1678	!-- Calculate the divergence of the velocity field. A respective
1679	!-- correction is needed to overcome numerical instabilities introduced
1680	!-- by a not sufficient reduction of divergences near topography.
1681	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
1682	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
1683	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) &
1684	) * 0.5
1685
1686	tend(k,j,i) = tend(k,j,i) - ( &
1687	( flux_r(k) + diss_r(k) &
1688	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
1689	+ ( flux_n(k) + diss_n(k) &
1690	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
1691	+ ( flux_t(k) + diss_t(k) &
1692	- flux_d - diss_d ) * ddzu(k+1) &
1693	) + div * w(k,j,i)
1694
1695	flux_l_w(k,j,tn) = flux_r(k)
1696	diss_l_w(k,j,tn) = diss_r(k)
1697	flux_s_w(k,tn) = flux_n(k)
1698	diss_s_w(k,tn) = diss_n(k)
1699	flux_d = flux_t(k)
1700	diss_d = diss_t(k)
1701	!
1702	!-- Statistical Evaluation of w'w'.
1703	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
1704	+ ( flux_t(k) + diss_t(k) ) &
1705	* weight_substep(intermediate_timestep_count)
1706
1707	ENDDO
1708
1709	DO k = nzb_max+1, nzt
1710
1711	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1712	flux_r(k) = u_comp * ( &
1713	37.0 * ( w(k,j,i+1) + w(k,j,i) ) &
1714	- 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) &
1715	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1716
1717	diss_r(k) = - ABS( u_comp ) * ( &
1718	10.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1719	- 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) &
1720	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1721
1722	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1723	flux_n(k) = v_comp * ( &
1724	37.0 * ( w(k,j+1,i) + w(k,j,i) ) &
1725	- 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) &
1726	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1727
1728	diss_n(k) = - ABS( v_comp ) * ( &
1729	10.0 * ( w(k,j+1,i) - w(k,j,i) ) &
1730	- 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) &
1731	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1732	!
1733	!-- k index has to be modified near bottom and top, else array
1734	!-- subscripts will be exceeded.
1735	k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1)
1736	k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1) )
1737	k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),35,1)
1738
1739	w_comp = w(k+1,j,i) + w(k,j,i)
1740	flux_t(k) = w_comp * ( &
1741	( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1742	+ 7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1743	+ IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 &
1744	) * &
1745	( w(k+1,j,i) + w(k,j,i) ) &
1746	- ( 8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1747	+ IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1748	) * &
1749	( w(k_pp,j,i) + w(k-1,j,i) ) &
1750	+ ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1751	) * &
1752	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
1753	)
1754
1755	diss_t(k) = - ABS( w_comp ) * ( &
1756	( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1757	+ 3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1758	+ IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1 &
1759	) * &
1760	( w(k+1,j,i) - w(k,j,i) ) &
1761	- ( 5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1762	+ IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3 &
1763	) * &
1764	( w(k_pp,j,i) - w(k-1,j,i) ) &
1765	+ ( IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5 &
1766	) * &
1767	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
1768	)
1769	!
1770	!-- Calculate the divergence of the velocity field. A respective
1771	!-- correction is needed to overcome numerical instabilities introduced
1772	!-- by a not sufficient reduction of divergences near topography.
1773	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
1774	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
1775	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) &
1776	) * 0.5
1777
1778	tend(k,j,i) = tend(k,j,i) - ( &
1779	( flux_r(k) + diss_r(k) &
1780	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
1781	+ ( flux_n(k) + diss_n(k) &
1782	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
1783	+ ( flux_t(k) + diss_t(k) &
1784	- flux_d - diss_d ) * ddzu(k+1) &
1785	) + div * w(k,j,i)
1786
1787	flux_l_w(k,j,tn) = flux_r(k)
1788	diss_l_w(k,j,tn) = diss_r(k)
1789	flux_s_w(k,tn) = flux_n(k)
1790	diss_s_w(k,tn) = diss_n(k)
1791	flux_d = flux_t(k)
1792	diss_d = diss_t(k)
1793	!
1794	!-- Statistical Evaluation of w'w'.
1795	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
1796	+ ( flux_t(k) + diss_t(k) ) &
1797	* weight_substep(intermediate_timestep_count)
1798
1799	ENDDO
1800
1801
1802	END SUBROUTINE advec_w_ws_ij
1803
1804
1805	!------------------------------------------------------------------------------!
1806	! Scalar advection - Call for all grid points
1807	!------------------------------------------------------------------------------!
1808	SUBROUTINE advec_s_ws( sk, sk_char )
1809
1810	USE arrays_3d
1811	USE constants
1812	USE control_parameters
1813	USE grid_variables
1814	USE indices
1815	USE statistics
1816
1817	IMPLICIT NONE
1818
1819	INTEGER :: i, j, k, tn = 0, k_ppp, k_pp, k_mm
1820	REAL, DIMENSION(:,:,:), POINTER :: sk
1821	REAL :: flux_d, diss_d, u_comp, v_comp, div
1822	REAL, DIMENSION(nzb:nzt) :: flux_r, diss_r, flux_n, diss_n, flux_t, &
1823	diss_t
1824	REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, swap_diss_y_local
1825	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local, &
1826	swap_diss_x_local
1827	CHARACTER (LEN = *), INTENT(IN) :: sk_char
1828
1829	!
1830	!-- Compute the fluxes for the whole left boundary of the processor domain.
1831	i = nxl
1832	DO j = nys, nyn
1833
1834	DO k = nzb+1, nzb_max
1835
1836	u_comp = u(k,j,i) - u_gtrans
1837	swap_flux_x_local(k,j) = u_comp * ( &
1838	( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1839	+ 7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1840	+ IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 &
1841	) * &
1842	( sk(k,j,i) + sk(k,j,i-1) ) &
1843	- ( 8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1844	+ IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1845	) * &
1846	( sk(k,j,i+1) + sk(k,j,i-2) ) &
1847	+ ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1848	) * ( sk(k,j,i+2) + sk(k,j,i-3) ) &
1849	)
1850
1851	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
1852	( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1853	+ 3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1854	+ IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 &
1855	) * &
1856	( sk(k,j,i) - sk(k,j,i-1) ) &
1857	- ( 5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1858	+ IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1859	) * &
1860	( sk(k,j,i+1) - sk(k,j,i-2) ) &
1861	+ ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1862	) * &
1863	( sk(k,j,i+2) - sk(k,j,i-3) ) &
1864	)
1865
1866	ENDDO
1867
1868	DO k = nzb_max+1, nzt
1869
1870	u_comp = u(k,j,i) - u_gtrans
1871	swap_flux_x_local(k,j) = u_comp * ( &
1872	37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) &
1873	- 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
1874	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
1875	) * adv_sca_5
1876
1877	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
1878	10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) &
1879	- 5.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	ENDDO
1884
1885	ENDDO
1886
1887	DO i = nxl, nxr
1888
1889	j = nys
1890	DO k = nzb+1, nzb_max
1891
1892	v_comp = v(k,j,i) - v_gtrans
1893	swap_flux_y_local(k) = v_comp * ( &
1894	( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1895	+ 7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
1896	+ IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 &
1897	) * &
1898	( sk(k,j,i) + sk(k,j-1,i) ) &
1899	- ( 8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1900	+ IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
1901	) * &
1902	( sk(k,j+1,i) + sk(k,j-2,i) ) &
1903	+ ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1904	) * &
1905	( sk(k,j+2,i) + sk(k,j-3,i) ) &
1906	)
1907
1908	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
1909	( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1910	+ 3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
1911	+ IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 &
1912	) * &
1913	( sk(k,j,i) - sk(k,j-1,i) ) &
1914	- ( 5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1915	+ IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
1916	) * &
1917	( sk(k,j+1,i) - sk(k,j-2,i) ) &
1918	+ ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1919	) * &
1920	( sk(k,j+2,i) - sk(k,j-3,i) ) &
1921	)
1922
1923	ENDDO
1924	!
1925	!-- Above to the top of the highest topography. No degradation necessary.
1926	DO k = nzb_max+1, nzt
1927
1928	v_comp = v(k,j,i) - v_gtrans
1929	swap_flux_y_local(k) = v_comp * ( &
1930	37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) &
1931	- 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
1932	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
1933	) * adv_sca_5
1934	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
1935	10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) &
1936	- 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
1937	+ sk(k,j+2,i) - sk(k,j-3,i) &
1938	) * adv_sca_5
1939
1940	ENDDO
1941
1942	DO j = nys, nyn
1943
1944	flux_t(0) = 0.0
1945	diss_t(0) = 0.0
1946	flux_d = 0.0
1947	diss_d = 0.0
1948
1949	DO k = nzb+1, nzb_max
1950
1951	u_comp = u(k,j,i+1) - u_gtrans
1952	flux_r(k) = u_comp * ( &
1953	( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1954	+ 7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1955	+ IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 &
1956	) * &
1957	( sk(k,j,i+1) + sk(k,j,i) ) &
1958	- ( 8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1959	+ IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1960	) * &
1961	( sk(k,j,i+2) + sk(k,j,i-1) ) &
1962	+ ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1963	) * ( sk(k,j,i+3) + sk(k,j,i-2) ) &
1964	)
1965
1966	diss_r(k) = -ABS( u_comp ) * ( &
1967	( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1968	+ 3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1969	+ IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1 &
1970	) * &
1971	( sk(k,j,i+1) - sk(k,j,i) ) &
1972	- ( 5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1973	+ IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3 &
1974	) * &
1975	( sk(k,j,i+2) - sk(k,j,i-1) ) &
1976	+ ( IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5 &
1977	) * &
1978	( sk(k,j,i+3) - sk(k,j,i-2) ) &
1979	)
1980
1981	v_comp = v(k,j+1,i) - v_gtrans
1982	flux_n(k) = v_comp * ( &
1983	( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1984	+ 7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
1985	+ IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 &
1986	) * &
1987	( sk(k,j+1,i) + sk(k,j,i) ) &
1988	- ( 8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1989	+ IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
1990	) * &
1991	( sk(k,j+2,i) + sk(k,j-1,i) ) &
1992	+ ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1993	) * &
1994	( sk(k,j+3,i) + sk(k,j-2,i) ) &
1995	)
1996
1997	diss_n(k) = -ABS( v_comp ) * ( &
1998	( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
1999	+ 3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
2000	+ IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1 &
2001	) * &
2002	( sk(k,j+1,i) - sk(k,j,i) ) &
2003	- ( 5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
2004	+ IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3 &
2005	) * &
2006	( sk(k,j+2,i) - sk(k,j-1,i) ) &
2007	+ ( IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5 &
2008	) * &
2009	( sk(k,j+3,i) - sk(k,j-2,i) ) &
2010	)
2011	!
2012	!-- k index has to be modified near bottom and top, else array
2013	!-- subscripts will be exceeded.
2014	k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
2015	k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1) )
2016	k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
2017
2018
2019	flux_t(k) = w(k,j,i) * ( &
2020	( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2021	+ 7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2022	+ IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 &
2023	) * &
2024	( sk(k+1,j,i) + sk(k,j,i) ) &
2025	- ( 8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2026	+ IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2027	) * &
2028	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
2029	+ ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2030	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
2031	)
2032
2033	diss_t(k) = -ABS( w(k,j,i) ) * ( &
2034	( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2035	+ 3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2036	+ IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 &
2037	) * &
2038	( sk(k+1,j,i) - sk(k,j,i) ) &
2039	- ( 5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2040	+ IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2041	) * &
2042	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
2043	+ ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2044	) * &
2045	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
2046	)
2047	!
2048	!-- Calculate the divergence of the velocity field. A respective
2049	!-- correction is needed to overcome numerical instabilities caused
2050	!-- by a not sufficient reduction of divergences near topography.
2051	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
2052	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
2053	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
2054
2055	tend(k,j,i) = tend(k,j,i) - ( &
2056	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
2057	swap_diss_x_local(k,j) ) * ddx &
2058	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
2059	swap_diss_y_local(k) ) * ddy &
2060	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
2061	) * ddzw(k) &
2062	) + sk(k,j,i) * div
2063
2064	swap_flux_y_local(k) = flux_n(k)
2065	swap_diss_y_local(k) = diss_n(k)
2066	swap_flux_x_local(k,j) = flux_r(k)
2067	swap_diss_x_local(k,j) = diss_r(k)
2068	flux_d = flux_t(k)
2069	diss_d = diss_t(k)
2070
2071	ENDDO
2072
2073	DO k = nzb_max+1, nzt
2074
2075	u_comp = u(k,j,i+1) - u_gtrans
2076	flux_r(k) = u_comp * ( &
2077	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
2078	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
2079	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
2080	diss_r(k) = -ABS( u_comp ) * ( &
2081	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
2082	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
2083	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
2084
2085	v_comp = v(k,j+1,i) - v_gtrans
2086	flux_n(k) = v_comp * ( &
2087	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
2088	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
2089	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
2090	diss_n(k) = -ABS( v_comp ) * ( &
2091	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
2092	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
2093	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
2094	!
2095	!-- k index has to be modified near bottom and top, else array
2096	!-- subscripts will be exceeded.
2097	k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
2098	k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1) )
2099	k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
2100
2101	flux_t(k) = w(k,j,i) * ( &
2102	( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2103	+ 7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2104	+ IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 &
2105	) * &
2106	( sk(k+1,j,i) + sk(k,j,i) ) &
2107	- ( 8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2108	+ IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2109	) * &
2110	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
2111	+ ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2112	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
2113	)
2114
2115	diss_t(k) = -ABS( w(k,j,i) ) * ( &
2116	( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2117	+ 3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2118	+ IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1 &
2119	) * &
2120	( sk(k+1,j,i) - sk(k,j,i) ) &
2121	- ( 5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2122	+ IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3 &
2123	) * &
2124	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
2125	+ ( IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5 &
2126	) * &
2127	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
2128	)
2129	!
2130	!-- Calculate the divergence of the velocity field. A respective
2131	!-- correction is needed to overcome numerical instabilities introduced
2132	!-- by a not sufficient reduction of divergences near topography.
2133	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
2134	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
2135	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
2136
2137	tend(k,j,i) = tend(k,j,i) - ( &
2138	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
2139	swap_diss_x_local(k,j) ) * ddx &
2140	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
2141	swap_diss_y_local(k) ) * ddy &
2142	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
2143	) * ddzw(k) &
2144	) + sk(k,j,i) * div
2145
2146	swap_flux_y_local(k) = flux_n(k)
2147	swap_diss_y_local(k) = diss_n(k)
2148	swap_flux_x_local(k,j) = flux_r(k)
2149	swap_diss_x_local(k,j) = diss_r(k)
2150	flux_d = flux_t(k)
2151	diss_d = diss_t(k)
2152
2153	ENDDO
2154	!
2155	!-- evaluation of statistics
2156	SELECT CASE ( sk_char )
2157
2158	CASE ( 'pt' )
2159	DO k = nzb, nzt
2160	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
2161	+ ( flux_t(k) + diss_t(k) ) &
2162	* weight_substep(intermediate_timestep_count)
2163	ENDDO
2164	CASE ( 'sa' )
2165	DO k = nzb, nzt
2166	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
2167	+ ( flux_t(k) + diss_t(k) ) &
2168	* weight_substep(intermediate_timestep_count)
2169	ENDDO
2170	CASE ( 'q' )
2171	DO k = nzb, nzt
2172	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
2173	+ ( flux_t(k) + diss_t(k) ) &
2174	* weight_substep(intermediate_timestep_count)
2175	ENDDO
2176
2177	END SELECT
2178
2179	ENDDO
2180	ENDDO
2181
2182	END SUBROUTINE advec_s_ws
2183
2184
2185	!------------------------------------------------------------------------------!
2186	! Advection of u - Call for all grid points
2187	!------------------------------------------------------------------------------!
2188	SUBROUTINE advec_u_ws
2189
2190	USE arrays_3d
2191	USE constants
2192	USE control_parameters
2193	USE grid_variables
2194	USE indices
2195	USE statistics
2196
2197	IMPLICIT NONE
2198
2199	INTEGER :: i, j, k, tn = 0, k_ppp, k_pp, k_mm
2200	REAL :: gu, gv, flux_d, diss_d, v_comp, w_comp, div
2201	REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u, swap_diss_y_local_u
2202	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u, &
2203	swap_diss_x_local_u
2204	REAL, DIMENSION(nzb:nzt) :: flux_t, diss_t, flux_r, diss_r, flux_n, &
2205	diss_n, u_comp
2206
2207	gu = 2.0 * u_gtrans
2208	gv = 2.0 * v_gtrans
2209
2210	!
2211	!-- Compute the fluxes for the whole left boundary of the processor domain.
2212	i = nxlu
2213	DO j = nys, nyn
2214	DO k = nzb+1, nzb_max
2215
2216	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
2217	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
2218	( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2219	+ 7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2220	+ IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 &
2221	) * &
2222	( u(k,j,i) + u(k,j,i-1) ) &
2223	- ( 8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2224	+ IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2225	) * &
2226	( u(k,j,i+1) + u(k,j,i-2) ) &
2227	+ ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2228	) * &
2229	( u(k,j,i+2) + u(k,j,i-3) ) &
2230	)
2231
2232	swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( &
2233	( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2234	+ 3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2235	+ IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 &
2236	) * &
2237	( u(k,j,i) - u(k,j,i-1) ) &
2238	- ( 5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2239	+ IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2240	) * &
2241	( u(k,j,i+1) - u(k,j,i-2) ) &
2242	+ ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2243	) * &
2244	( u(k,j,i+2) - u(k,j,i-3) ) &
2245	)
2246
2247	ENDDO
2248
2249	DO k = nzb_max+1, nzt
2250
2251	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
2252	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
2253	37.0 * ( u(k,j,i) + u(k,j,i-1) ) &
2254	- 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) &
2255	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
2256	swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( &
2257	10.0 * ( u(k,j,i) - u(k,j,i-1) ) &
2258	- 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) &
2259	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
2260
2261	ENDDO
2262	ENDDO
2263
2264	DO i = nxlu, nxr
2265	!
2266	!-- The following loop computes the fluxes for the south boundary points
2267	j = nys
2268	DO k = nzb+1, nzb_max
2269
2270	v_comp = v(k,j,i) + v(k,j,i-1) - gv
2271	swap_flux_y_local_u(k) = v_comp * ( &
2272	( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2273	+ 7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2274	+ IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
2275	) * &
2276	( u(k,j,i) + u(k,j-1,i) ) &
2277	- ( 8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2278	+ IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2279	) * &
2280	( u(k,j+1,i) + u(k,j-2,i) ) &
2281	+ ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2282	) * &
2283	( u(k,j+2,i) + u(k,j-3,i) ) &
2284	)
2285
2286	swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( &
2287	( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2288	+ 3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2289	+ IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
2290	) * &
2291	( u(k,j,i) - u(k,j-1,i) ) &
2292	- ( 5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2293	+ IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2294	) * &
2295	( u(k,j+1,i) - u(k,j-2,i) ) &
2296	+ ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2297	) * &
2298	( u(k,j+2,i) - u(k,j-3,i) ) &
2299	)
2300
2301	ENDDO
2302
2303	DO k = nzb_max+1, nzt
2304
2305	v_comp = v(k,j,i) + v(k,j,i-1) - gv
2306	swap_flux_y_local_u(k) = v_comp * ( &
2307	37.0 * ( u(k,j,i) + u(k,j-1,i) ) &
2308	- 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) &
2309	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
2310	swap_diss_y_local_u(k) = - ABS(v_comp) * ( &
2311	10.0 * ( u(k,j,i) - u(k,j-1,i) ) &
2312	- 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) &
2313	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
2314
2315	ENDDO
2316	!
2317	!-- Computation of interior fluxes and tendency terms
2318	DO j = nys, nyn
2319
2320	flux_t(0) = 0.0
2321	diss_t(0) = 0.0
2322	flux_d = 0.0
2323	diss_d = 0.0
2324
2325	DO k = nzb+1, nzb_max
2326
2327	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2328	flux_r(k) = ( u_comp(k) - gu ) * ( &
2329	( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2330	+ 7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2331	+ IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 &
2332	) * &
2333	( u(k,j,i+1) + u(k,j,i) ) &
2334	- ( 8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2335	+ IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2336	) * &
2337	( u(k,j,i+2) + u(k,j,i-1) ) &
2338	+ ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2339	) * &
2340	( u(k,j,i+3) + u(k,j,i-2) ) &
2341	)
2342
2343	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2344	( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2345	+ 3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2346	+ IBITS(wall_flags_0(k,j,i),9,1) * adv_mom_1 &
2347	) * &
2348	( u(k,j,i+1) - u(k,j,i) ) &
2349	- ( 5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2350	+ IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2351	) * &
2352	( u(k,j,i+2) - u(k,j,i-1) ) &
2353	+ ( IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2354	) * &
2355	( u(k,j,i+3) - u(k,j,i-2) ) &
2356	)
2357
2358	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2359	flux_n(k) = v_comp * ( &
2360	( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2361	+ 7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2362	+ IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
2363	) * &
2364	( u(k,j+1,i) + u(k,j,i) ) &
2365	- ( 8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2366	+ IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2367	) * &
2368	( u(k,j+2,i) + u(k,j-1,i) ) &
2369	+ ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2370	) * &
2371	( u(k,j+3,i) + u(k,j-2,i) ) &
2372	)
2373
2374	diss_n(k) = - ABS ( v_comp ) * ( &
2375	( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2376	+ 3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2377	+ IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
2378	) * &
2379	( u(k,j+1,i) - u(k,j,i) ) &
2380	- ( 5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2381	+ IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2382	) * &
2383	( u(k,j+2,i) - u(k,j-1,i) ) &
2384	+ ( IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2385	) * &
2386	( u(k,j+3,i) - u(k,j-2,i) ) &
2387	)
2388	!
2389	!-- k index has to be modified near bottom and top, else array
2390	!-- subscripts will be exceeded.
2391	k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
2392	k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1) )
2393	k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
2394
2395	w_comp = w(k,j,i) + w(k,j,i-1)
2396	flux_t(k) = w_comp * ( &
2397	( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2398	+ 7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2399	+ IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 &
2400	) * &
2401	( u(k+1,j,i) + u(k,j,i) ) &
2402	- ( 8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2403	+ IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2404	) * &
2405	( u(k_pp,j,i) + u(k-1,j,i) ) &
2406	+ ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2407	) * &
2408	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2409	)
2410
2411	diss_t(k) = - ABS( w_comp ) * ( &
2412	( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2413	+ 3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2414	+ IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 &
2415	) * &
2416	( u(k+1,j,i) - u(k,j,i) ) &
2417	- ( 5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2418	+ IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2419	) * &
2420	( u(k_pp,j,i) - u(k-1,j,i) ) &
2421	+ ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2422	) * &
2423	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2424	)
2425	!
2426	!-- Calculate the divergence of the velocity field. A respective
2427	!-- correction is needed to overcome numerical instabilities caused
2428	!-- by a not sufficient reduction of divergences near topography.
2429	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
2430	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
2431	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) &
2432	* ddzw(k) &
2433	) * 0.5
2434
2435	tend(k,j,i) = tend(k,j,i) - ( &
2436	( flux_r(k) + diss_r(k) &
2437	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
2438	+ ( flux_n(k) + diss_n(k) &
2439	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
2440	+ ( flux_t(k) + diss_t(k) &
2441	- flux_d - diss_d &
2442	) * ddzw(k) &
2443	) + div * u(k,j,i)
2444
2445	swap_flux_x_local_u(k,j) = flux_r(k)
2446	swap_diss_x_local_u(k,j) = diss_r(k)
2447	swap_flux_y_local_u(k) = flux_n(k)
2448	swap_diss_y_local_u(k) = diss_n(k)
2449	flux_d = flux_t(k)
2450	diss_d = diss_t(k)
2451	!
2452	!-- Statistical Evaluation of u'u'. The factor has to be applied
2453	!-- for right evaluation when gallilei_trans = .T. .
2454	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2455	+ ( flux_r(k) * &
2456	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
2457	/ ( u_comp(k) - gu + 1.0E-20 ) &
2458	+ diss_r(k) * &
2459	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
2460	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
2461	* weight_substep(intermediate_timestep_count)
2462	!
2463	!-- Statistical Evaluation of w'u'.
2464	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2465	+ ( flux_t(k) + diss_t(k) ) &
2466	* weight_substep(intermediate_timestep_count)
2467	ENDDO
2468
2469	DO k = nzb_max+1, nzt
2470
2471	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2472	flux_r(k) = ( u_comp(k) - gu ) * ( &
2473	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
2474	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
2475	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2476	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2477	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
2478	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
2479	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2480
2481	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2482	flux_n(k) = v_comp * ( &
2483	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
2484	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
2485	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2486	diss_n(k) = - ABS( v_comp ) * ( &
2487	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
2488	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
2489	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2490	!
2491	!-- k index has to be modified near bottom and top, else array
2492	!-- subscripts will be exceeded.
2493	k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
2494	k_pp = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1) )
2495	k_mm = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
2496
2497	w_comp = w(k,j,i) + w(k,j,i-1)
2498	flux_t(k) = w_comp * ( &
2499	( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2500	+ 7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2501	+ IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 &
2502	) * &
2503	( u(k+1,j,i) + u(k,j,i) ) &
2504	- ( 8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2505	+ IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2506	) * &
2507	( u(k_pp,j,i) + u(k-1,j,i) ) &
2508	+ ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2509	) * &
2510	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
2511	)
2512
2513	diss_t(k) = - ABS( w_comp ) * ( &
2514	( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2515	+ 3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2516	+ IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1 &
2517	) * &
2518	( u(k+1,j,i) - u(k,j,i) ) &
2519	- ( 5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2520	+ IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3 &
2521	) * &
2522	( u(k_pp,j,i) - u(k-1,j,i) ) &
2523	+ ( IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5 &
2524	) * &
2525	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
2526	)
2527	!
2528	!-- Calculate the divergence of the velocity field. A respective
2529	!-- correction is needed to overcome numerical instabilities caused
2530	!-- by a not sufficient reduction of divergences near topography.
2531	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
2532	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
2533	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) &
2534	* ddzw(k) &
2535	) * 0.5
2536
2537	tend(k,j,i) = tend(k,j,i) - ( &
2538	( flux_r(k) + diss_r(k) &
2539	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
2540	+ ( flux_n(k) + diss_n(k) &
2541	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
2542	+ ( flux_t(k) + diss_t(k) &
2543	- flux_d - diss_d &
2544	) * ddzw(k) &
2545	) + div * u(k,j,i)
2546
2547	swap_flux_x_local_u(k,j) = flux_r(k)
2548	swap_diss_x_local_u(k,j) = diss_r(k)
2549	swap_flux_y_local_u(k) = flux_n(k)
2550	swap_diss_y_local_u(k) = diss_n(k)
2551	flux_d = flux_t(k)
2552	diss_d = diss_t(k)
2553	!
2554	!-- Statistical Evaluation of u'u'. The factor has to be applied
2555	!-- for right evaluation when gallilei_trans = .T. .
2556	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2557	+ ( flux_r(k) * &
2558	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
2559	/ ( u_comp(k) - gu + 1.0E-20 ) &
2560	+ diss_r(k) * &
2561	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
2562	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
2563	* weight_substep(intermediate_timestep_count)
2564	!
2565	!-- Statistical Evaluation of w'u'.
2566	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2567	+ ( flux_t(k) + diss_t(k) ) &
2568	* weight_substep(intermediate_timestep_count)
2569	ENDDO
2570	ENDDO
2571	ENDDO
2572	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
2573
2574
2575	END SUBROUTINE advec_u_ws
2576
2577
2578	!------------------------------------------------------------------------------!
2579	! Advection of v - Call for all grid points
2580	!------------------------------------------------------------------------------!
2581	SUBROUTINE advec_v_ws
2582
2583	USE arrays_3d
2584	USE constants
2585	USE control_parameters
2586	USE grid_variables
2587	USE indices
2588	USE statistics
2589
2590	IMPLICIT NONE
2591
2592
2593	INTEGER :: i, j, k, tn = 0, k_ppp, k_pp, k_mm
2594	REAL :: gu, gv, flux_d, diss_d, u_comp, w_comp, div
2595	REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v, swap_diss_y_local_v
2596	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v, &
2597	swap_diss_x_local_v
2598	REAL, DIMENSION(nzb:nzt) :: flux_t, diss_t, flux_n, diss_n, flux_r, &
2599	diss_r, v_comp
2600
2601	gu = 2.0 * u_gtrans
2602	gv = 2.0 * v_gtrans
2603	!
2604	!-- First compute the whole left boundary of the processor domain
2605	i = nxl
2606	DO j = nysv, nyn
2607	DO k = nzb+1, nzb_max
2608
2609	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2610	swap_flux_x_local_v(k,j) = u_comp * ( &
2611	( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2612	+ 7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2613	+ IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 &
2614	) * &
2615	( v(k,j,i) + v(k,j,i-1) ) &
2616	- ( 8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2617	+ IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2618	) * &
2619	( v(k,j,i+1) + v(k,j,i-2) ) &
2620	+ ( IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2621	) * &
2622	( v(k,j,i+2) + v(k,j,i-3) ) &
2623	)
2624
2625	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
2626	( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2627	+ 3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2628	+ IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 &
2629	) * &
2630	( v(k,j,i) - v(k,j,i-1) ) &
2631	- ( 5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2632	+ IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2633	) * &
2634	( v(k,j,i+1) - v(k,j,i-2) ) &
2635	+ ( IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2636	) * &
2637	( v(k,j,i+2) - v(k,j,i-3) ) &
2638	)
2639	</