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

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

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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