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

Last change on this file since 1823 was 1823, checked in by hoffmann, 5 years ago

last commit documented

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