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

Last change on this file since 1323 was 1323, checked in by raasch, 11 years ago

last commit documented

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