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

Last change on this file since 3467 was 3467, checked in by suehring, 3 years ago

Branch salsa @3446 re-integrated into trunk

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