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

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

last commit documented

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