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

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

Bugfix in calculation of turbulent fluxes

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