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

Last change on this file since 2731 was 2731, checked in by suehring, 4 years ago

Enable vectorization by splitting the k-loop

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