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

Last change on this file since 3421 was 3302, checked in by raasch, 3 years ago

Craik-Leibovich force and wave breaking effects added to ocean mode, header output for ocean mode

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