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

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

Implementation of a monotonic flux limiter for vertical advection term in Wicker-Skamarock scheme. The flux limiter is currently only applied for passive scalars (passive scalar, chemical species, aerosols) within the region up to the highest topography, in order to avoid the built-up of large concentrations within poorly resolved cavities in urban environments. To enable the limiter monotonic_limiter_z = .T. must be set. Note, the limiter is currently only implemented for the cache-optimized version of advec_ws. Further changes in offline nesting: Set boundary condition for w at nzt+1 at all lateral boundaries (even though these won't enter the numerical solution), in order to avoid high vertical velocities in the run-control file which might built-up due to the mass-conservation; bugfix in offline nesting for chemical species

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