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

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

last commit documented

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