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

Last change on this file since 1557 was 1557, checked in by suehring, 7 years ago

Enable monotone advection for scalars in combination with fifth-order scheme using monotonic limiter.

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