Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |