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

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

Including last commit, salsa dependency for advec_ws removed

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