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

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

Changes documented

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