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

Last change on this file since 4329 was 4329, checked in by motisi, 16 months ago

Renamed wall_flags_0 to wall_flags_static_0

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