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

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

Separate balance equations for humidity and passive_scalar

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