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

Last change on this file since 3298 was 3298, checked in by kanani, 3 years ago

Merge chemistry branch at r3297 to trunk

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