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

Last change on this file since 2119 was 2119, checked in by raasch, 5 years ago

last commit documented

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