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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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