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

Last change on this file since 4340 was 4340, checked in by Giersch, 4 years ago

Topography closed channel flow with symmetric boundaries implemented, ID tag in radiation module corrected

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