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

Last change on this file since 4313 was 4204, checked in by knoop, 22 months ago

Bugfix in advec_s_ws: Changed sk_num initialization default to avoid implicit SAVE-attribut.

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