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

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

last commit documented

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