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

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

last commit documented

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