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

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

various changes to avoid compiler warnings (mainly removal of unused variables)

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