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

Last change on this file since 4330 was 4330, checked in by knoop, 16 months ago

Bugix: removed syntax error introduced by last commit

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