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

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

Bugfixes in monotonic limter.

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