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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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