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

Last change on this file since 3864 was 3864, checked in by monakurppa, 2 years ago

major changes in salsa: data input, format and performance

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