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

Last change on this file since 3661 was 3661, checked in by suehring, 3 years ago

Split loops in advec_ws in order to reduce bit queries; Introduce new parameter to better control order degradation of advection scheme at non-cyclic boundaries; Remove setting of Neumann conditions for horizontal velocity variances; Minor bugfix in divergence correction in advection scheme (only has implications at downward-facing wall surfaces)

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