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

Last change on this file since 2118 was 2118, checked in by raasch, 5 years ago

all OpenACC directives and related parts removed from the code

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