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

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

optimization of advec_ws

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