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

Last change on this file since 2329 was 2329, checked in by knoop, 4 years ago

Bugfix for topography usage with anelastic approximation and boussinesq approximation with air density != 1

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