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

Last change on this file since 3294 was 3294, checked in by raasch, 3 years ago

modularization of the ocean code

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