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

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