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

Last change on this file since 4317 was 4317, checked in by Giersch, 17 months ago

Comments revised/added, formatting improved, fluxes for u,v, and scalars are explicitly set to zero at nzt+1, fluxes of w-component are now calculated only until nzt-1 (Prognostic equation for w-velocity component ends at nzt-1)

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