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

Last change on this file since 1997 was 1997, checked in by suehring, 5 years ago

last commit documented

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