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

Last change on this file since 3049 was 3022, checked in by suehring, 3 years ago

Revise recent bugfix in nested runs at left and south boundary; bugfix in advection of u in case of OpenMP parallelization; bugfix in plant transpiration

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