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

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

Merge branch salsa with trunk

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