Home

Context Navigation

advec_ws.f90 @ 1567

Last change on this file since 1567 was 1567, checked in by suehring, 10 years ago
Bugfixes in monotonic limter.
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 342.4 KB

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