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

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

Bugfix in degrading the order of the advection scheme at non-cyclic lateral domain boundaries

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