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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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