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

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

Setting of advection flags for vertical fluxes of w revised, Bugfix: air density for vertical flux calculation of w at k=1 is considered now

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