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

Last change on this file since 1822 was 1822, checked in by hoffmann, 5 years ago

changes in LPM and bulk cloud microphysics

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