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

Last change on this file since 3547 was 3547, checked in by suehring, 4 years ago

variable description added some routines

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