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

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

last commit documented

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