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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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