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

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

filter topography by filling holes resolved by only one grid point ;initialization of wall_flags_0 and wall_flags_00 moved to advec_ws

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