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

Last change on this file since 1374 was 1374, checked in by raasch, 7 years ago

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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