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

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

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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