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

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

last commit documented

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