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

Last change on this file since 3870 was 3870, checked in by knoop, 2 years ago

Moving prognostic equations of bcm into bulk_cloud_model_mod

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