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

Last change on this file since 4581 was 4581, checked in by suehring, 12 months ago

mesoscale nesting: omit explicit pressure forcing via geostrophic wind components; chemistry: enable profile output of vertical fluxes; urban-surface: bugfix in initialization in case of cyclic_fill

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