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

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

Moved ocean_mode specific code from advec_ws to ocean_mod + implemented ocean_actions

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