source: palm/trunk/SOURCE/advec_ws.f90

Last change on this file was 4828, checked in by Giersch, 7 months ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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