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

Last change on this file since 2292 was 2292, checked in by schwenkel, 4 years ago

implementation of new bulk microphysics scheme

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