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

Last change on this file since 1362 was 1362, checked in by hoffmann, 7 years ago

last commit documented

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