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

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

Bugfixes in radiation_model; monotonic limiter in 5th order advection scheme removed

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