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

Last change on this file since 3589 was 3589, checked in by suehring, 3 years ago

Remove erroneous UTF encoding; last commit documented

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