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

Last change on this file since 1873 was 1873, checked in by maronga, 5 years ago

revised renaming of modules

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