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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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