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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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