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

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

Adjustments according new topography and surface-modelling concept implemented

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