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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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