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

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

last commit documented

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