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

Last change on this file since 1353 was 1353, checked in by heinze, 11 years ago

REAL constants provided with KIND-attribute

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