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

Last change on this file since 4178 was 4110, checked in by suehring, 2 years ago

last changes documented

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