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

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

last commit documented

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