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

Last change on this file since 3871 was 3871, checked in by knoop, 2 years ago

Moving prognostic equations of bcm into bulk_cloud_model_mod

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