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

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

last commit documented

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