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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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