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

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

last commit documented

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