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

Last change on this file since 1877 was 1877, checked in by raasch, 5 years ago

file added without execute permits

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