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

Last change on this file since 4324 was 4324, checked in by Giersch, 16 months ago

Vertical fluxes of w are now set to zero at nzt and nzt+1, setting of advection flags for fluxes in z-direction revised, comments extended

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