source: palm/trunk/SOURCE/advec_ws_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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