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

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