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

Last change on this file since 4509 was 4509, checked in by raasch, 14 months ago

files re-formatted to follow the PALM coding standard

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