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

Last change on this file since 4318 was 4318, checked in by Giersch, 17 months ago

Indirect indexing for calculating vertical fluxes close to boundaries is only used for loop indices where it is really necessary

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