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

Last change on this file since 2037 was 2037, checked in by knoop, 5 years ago

Anelastic approximation implemented

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