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

Last change on this file since 3634 was 3634, checked in by knoop, 3 years ago

OpenACC port for SPEC

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