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

Last change on this file since 4502 was 4502, checked in by schwenkel, 15 months ago

Implementation of ice microphysics

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