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

Last change on this file since 1027 was 1027, checked in by suehring, 12 years ago

Bugfix in calculation indices k_mm, k_pp in accelerator version

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 268.6 KB
<
Line 
1 MODULE advec_ws
2
3!-----------------------------------------------------------------------------!
4! Current revisions:
5! ------------------
6! Bugfix in calculation indices k_mm, k_pp in accelerator version
7!
8! Former revisions:
9! -----------------
10! $Id: advec_ws.f90 1027 2012-10-15 17:18:39Z suehring $
11!
12! 1019 2012-09-28 06:46:45Z raasch
13! small change in comment lines
14!
15! 1015 2012-09-27 09:23:24Z raasch
16! accelerator versions (*_acc) added
17!
18! 1010 2012-09-20 07:59:54Z raasch
19! cpp switch __nopointer added for pointer free version
20!
21! 888 2012-04-20 15:03:46Z suehring
22! Number of IBITS() calls with identical arguments is reduced.
23!
24! 862 2012-03-26 14:21:38Z suehring
25! ws-scheme also work with topography in combination with vector version.
26! ws-scheme also work with outflow boundaries in combination with
27! vector version.
28! Degradation of the applied order of scheme is now steered by multiplying with
29! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each
30! grid point and mulitplied with the appropriate flag.
31! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
32! term derived according to the 4th and 6th order terms is applied. It turns
33! out that diss_2nd does not provide sufficient dissipation near walls.
34! Therefore, the function diss_2nd is removed.
35! Near walls a divergence correction is necessary to overcome numerical
36! instabilities due to too less divergence reduction of the velocity field.
37! boundary_flags and logicals steering the degradation are removed.
38! Empty SUBROUTINE local_diss removed.
39! Further formatting adjustments.
40!
41! 801 2012-01-10 17:30:36Z suehring
42! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
43! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
44! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
45! dimension.
46!
47! 743 2011-08-18 16:10:16Z suehring
48! Evaluation of turbulent fluxes with WS-scheme only for the whole model
49! domain. Therefor dimension of arrays needed for statistical evaluation
50! decreased.
51!
52! 736 2011-08-17 14:13:26Z suehring
53! Bugfix concerning OpenMP parallelization. i_omp introduced, because first
54! index where fluxes on left side have to be calculated explicitly is
55! different on each thread. Furthermore the swapping of fluxes is now
56! thread-safe by adding an additional dimension.
57!
58! 713 2011-03-30 14:21:21Z suehring
59! File reformatted.
60! Bugfix in vertical advection of w concerning the optimized version for   
61! vector architecture.
62! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as
63! broken numbers.
64!
65! 709 2011-03-30 09:31:40Z raasch
66! formatting adjustments
67!
68! 705 2011-03-25 11:21:43 Z suehring $
69! Bugfix in declaration of logicals concerning outflow boundaries.
70!
71! 411 2009-12-11 12:31:43 Z suehring
72! Allocation of weight_substep moved to init_3d_model.
73! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters.
74! Setting bc for the horizontal velocity variances added (moved from
75! flow_statistics).
76!
77! Initial revision
78!
79! 411 2009-12-11 12:31:43 Z suehring
80!
81! Description:
82! ------------
83! Advection scheme for scalars and momentum using the flux formulation of
84! Wicker and Skamarock 5th order. Additionally the module contains of a
85! routine using for initialisation and steering of the statical evaluation.
86! The computation of turbulent fluxes takes place inside the advection
87! routines.
88! Near non-cyclic boundaries the order of the applied advection scheme is
89! degraded.
90! A divergence correction is applied. It is necessary for topography, since
91! the divergence is not sufficiently reduced, resulting in erroneous fluxes and
92! partly numerical instabilities.
93!-----------------------------------------------------------------------------!
94
95    PRIVATE
96    PUBLIC   advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, &
97             advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc, &
98             ws_init, ws_statistics
99
100    INTERFACE ws_init
101       MODULE PROCEDURE ws_init
102    END INTERFACE ws_init
103
104    INTERFACE ws_statistics
105       MODULE PROCEDURE ws_statistics
106    END INTERFACE ws_statistics
107
108    INTERFACE advec_s_ws
109       MODULE PROCEDURE advec_s_ws
110       MODULE PROCEDURE advec_s_ws_ij
111    END INTERFACE advec_s_ws
112
113    INTERFACE advec_u_ws
114       MODULE PROCEDURE advec_u_ws
115       MODULE PROCEDURE advec_u_ws_ij
116    END INTERFACE advec_u_ws
117
118    INTERFACE advec_u_ws_acc
119       MODULE PROCEDURE advec_u_ws_acc
120    END INTERFACE advec_u_ws_acc
121
122    INTERFACE advec_v_ws
123       MODULE PROCEDURE advec_v_ws
124       MODULE PROCEDURE advec_v_ws_ij
125    END INTERFACE advec_v_ws
126
127    INTERFACE advec_v_ws_acc
128       MODULE PROCEDURE advec_v_ws_acc
129    END INTERFACE advec_v_ws_acc
130
131    INTERFACE advec_w_ws
132       MODULE PROCEDURE advec_w_ws
133       MODULE PROCEDURE advec_w_ws_ij
134    END INTERFACE advec_w_ws
135
136    INTERFACE advec_w_ws_acc
137       MODULE PROCEDURE advec_w_ws_acc
138    END INTERFACE advec_w_ws_acc
139
140 CONTAINS
141
142
143!------------------------------------------------------------------------------!
144! Initialization of WS-scheme
145!------------------------------------------------------------------------------!
146    SUBROUTINE ws_init
147
148       USE arrays_3d
149       USE constants
150       USE control_parameters
151       USE indices
152       USE pegrid
153       USE statistics
154
155!
156!--    Set the appropriate factors for scalar and momentum advection.
157       adv_sca_5 = 1./60.
158       adv_sca_3 = 1./12.
159       adv_sca_1 = 1./2.
160       adv_mom_5 = 1./120.
161       adv_mom_3 = 1./24.
162       adv_mom_1 = 1./4.
163!         
164!--    Arrays needed for statical evaluation of fluxes.
165       IF ( ws_scheme_mom )  THEN
166
167          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1),  &
168                    sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1),  &
169                    sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1),   &
170                    sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1),   &
171                    sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
172
173          sums_wsus_ws_l = 0.0
174          sums_wsvs_ws_l = 0.0
175          sums_us2_ws_l  = 0.0
176          sums_vs2_ws_l  = 0.0
177          sums_ws2_ws_l  = 0.0
178
179       ENDIF
180
181       IF ( ws_scheme_sca )  THEN
182
183          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
184          sums_wspts_ws_l = 0.0
185
186          IF ( humidity .OR. passive_scalar )  THEN
187             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
188             sums_wsqs_ws_l = 0.0
189          ENDIF
190
191          IF ( ocean )  THEN
192             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
193             sums_wssas_ws_l = 0.0
194          ENDIF
195
196       ENDIF
197
198!
199!--    Arrays needed for reasons of speed optimization for cache version.
200!--    For the vector version the buffer arrays are not necessary,
201!--    because the the fluxes can swapped directly inside the loops of the
202!--    advection routines.
203       IF ( loop_optimization /= 'vector' )  THEN
204
205          IF ( ws_scheme_mom )  THEN
206
207             ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),            &
208                       flux_s_v(nzb+1:nzt,0:threads_per_task-1),            &
209                       flux_s_w(nzb+1:nzt,0:threads_per_task-1),            &
210                       diss_s_u(nzb+1:nzt,0:threads_per_task-1),            &
211                       diss_s_v(nzb+1:nzt,0:threads_per_task-1),            &
212                       diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
213             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
214                       flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
215                       flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
216                       diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
217                       diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
218                       diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
219
220          ENDIF
221
222          IF ( ws_scheme_sca )  THEN
223
224             ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
225                       flux_s_e(nzb+1:nzt,0:threads_per_task-1),            &
226                       diss_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
227                       diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) 
228             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
229                       flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
230                       diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
231                       diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
232
233             IF ( humidity .OR. passive_scalar )  THEN
234                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),         &
235                          diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
236                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
237                          diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
238             ENDIF
239
240             IF ( ocean )  THEN
241                ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1),         &
242                          diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
243                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
244                          diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
245             ENDIF
246
247          ENDIF
248
249       ENDIF
250
251    END SUBROUTINE ws_init
252
253
254!------------------------------------------------------------------------------!
255! Initialize variables used for storing statistic qauntities (fluxes, variances)
256!------------------------------------------------------------------------------!
257    SUBROUTINE ws_statistics
258   
259       USE control_parameters
260       USE statistics
261
262       IMPLICIT NONE
263
264!       
265!--    The arrays needed for statistical evaluation are set to to 0 at the
266!--    beginning of prognostic_equations.
267       IF ( ws_scheme_mom )  THEN
268          sums_wsus_ws_l = 0.0
269          sums_wsvs_ws_l = 0.0
270          sums_us2_ws_l  = 0.0
271          sums_vs2_ws_l  = 0.0
272          sums_ws2_ws_l  = 0.0
273       ENDIF
274
275       IF ( ws_scheme_sca )  THEN
276          sums_wspts_ws_l = 0.0
277          IF ( humidity .OR. passive_scalar )  sums_wsqs_ws_l = 0.0
278          IF ( ocean )  sums_wssas_ws_l = 0.0
279
280       ENDIF
281
282    END SUBROUTINE ws_statistics
283
284
285!------------------------------------------------------------------------------!
286! Scalar advection - Call for grid point i,j
287!------------------------------------------------------------------------------!
288    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
289                              swap_diss_y_local, swap_flux_x_local, &
290                              swap_diss_x_local, i_omp, tn )
291
292       USE arrays_3d
293       USE constants
294       USE control_parameters
295       USE grid_variables
296       USE indices
297       USE pegrid
298       USE statistics
299
300       IMPLICIT NONE
301
302       INTEGER ::  i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6,        &
303                   ibit7, ibit8, i_omp, j, k, k_mm, k_pp, k_ppp,  tn
304       REAL    ::  diss_d, div, flux_d, u_comp, v_comp
305#if defined( __nopointer )
306       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk
307#else
308       REAL, DIMENSION(:,:,:), POINTER    ::  sk
309#endif
310       REAL, DIMENSION(nzb:nzt+1)         ::  diss_n, diss_r, diss_t, flux_n,  &
311                                              flux_r, flux_t
312       REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local,  &
313                                                           swap_flux_y_local
314       REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::              &
315                                                          swap_diss_x_local,   &
316                                                          swap_flux_x_local
317       CHARACTER (LEN = *), INTENT(IN)    ::  sk_char
318
319!
320!--    Compute southside fluxes of the respective PE bounds.
321       IF ( j == nys )  THEN
322!
323!--       Up to the top of the highest topography.
324          DO  k = nzb+1, nzb_max
325
326             ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
327             ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
328             ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
329
330             v_comp                  = v(k,j,i) - v_gtrans
331             swap_flux_y_local(k,tn) = v_comp * (                             &
332                                                  ( 37.0 * ibit5 * adv_sca_5  &
333                                               +     7.0 * ibit4 * adv_sca_3  &
334                                               +           ibit3 * adv_sca_1  &
335                                                  ) *                         &
336                                           ( sk(k,j,i)  + sk(k,j-1,i)     )   &
337                                            -     (  8.0 * ibit5 * adv_sca_5  &
338                                               +           ibit4 * adv_sca_3  &
339                                                   ) *                        &
340                                           ( sk(k,j+1,i) + sk(k,j-2,i)    )   &
341                                           +      (        ibit5 * adv_sca_5  &
342                                                  ) *                         &
343                                           ( sk(k,j+2,i) + sk(k,j-3,i)    )   &
344                                                )
345
346             swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                     &
347                                                  ( 10.0 * ibit5 * adv_sca_5  &
348                                               +     3.0 * ibit4 * adv_sca_3  &
349                                               +           ibit3 * adv_sca_1  &
350                                                  ) *                         &
351                                            ( sk(k,j,i)   - sk(k,j-1,i)    )  &
352                                           -      (  5.0 * ibit5 * adv_sca_5  &
353                                               +           ibit4 * adv_sca_3  &
354                                            ) *                               &
355                                            ( sk(k,j+1,i) - sk(k,j-2,i)  )    &
356                                           +      (        ibit5 * adv_sca_5  &
357                                                  ) *                         &
358                                            ( sk(k,j+2,i) - sk(k,j-3,i) )     &
359                                                        )
360
361          ENDDO
362!
363!--       Above to the top of the highest topography. No degradation necessary.
364          DO  k = nzb_max+1, nzt
365
366             v_comp                  = v(k,j,i) - v_gtrans
367             swap_flux_y_local(k,tn) = v_comp * (                            &
368                                    37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
369                                  -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
370                                  +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
371                                             ) * adv_sca_5
372              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                   &
373                                    10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
374                                  -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
375                                  +          sk(k,j+2,i) - sk(k,j-3,i)       &
376                                                        ) * adv_sca_5
377
378          ENDDO
379
380       ENDIF
381!
382!--    Compute leftside fluxes of the respective PE bounds.
383       IF ( i == i_omp )  THEN
384       
385          DO  k = nzb+1, nzb_max
386
387             ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
388             ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
389             ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
390
391             u_comp                     = u(k,j,i) - u_gtrans
392             swap_flux_x_local(k,j,tn) = u_comp * (                           &
393                                                  ( 37.0 * ibit2 * adv_sca_5  &
394                                               +     7.0 * ibit1 * adv_sca_3  &
395                                               +           ibit0 * adv_sca_1  &
396                                                  ) *                         &
397                                            ( sk(k,j,i)   + sk(k,j,i-1)    )  &
398                                           -      (  8.0 * ibit2 * adv_sca_5  &
399                                               +           ibit1 * adv_sca_3  &
400                                                  ) *                         &
401                                            ( sk(k,j,i+1) + sk(k,j,i-2)    )  &
402                                           +      (        ibit2 * adv_sca_5  &
403                                                  ) *                         &
404                                            ( sk(k,j,i+2) + sk(k,j,i-3)    )  &
405                                                  )
406
407              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
408                                                  ( 10.0 * ibit2 * adv_sca_5  &
409                                               +     3.0 * ibit1 * adv_sca_3  &
410                                               +           ibit0 * adv_sca_1  &
411                                                  ) *                         &
412                                            ( sk(k,j,i)   - sk(k,j,i-1)    )  &
413                                           -      (  5.0 * ibit2 * adv_sca_5  &
414                                               +           ibit1 * adv_sca_3  &
415                                                  ) *                         &
416                                            ( sk(k,j,i+1) - sk(k,j,i-2)  )    &
417                                           +      (        ibit2 * adv_sca_5  &
418                                                  ) *                         &
419                                            ( sk(k,j,i+2) - sk(k,j,i-3) )     &
420                                                           )
421
422          ENDDO
423
424          DO  k = nzb_max+1, nzt
425
426             u_comp                 = u(k,j,i) - u_gtrans
427             swap_flux_x_local(k,j,tn) = u_comp * (                          &
428                                      37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
429                                    -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
430                                    +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
431                                                  ) * adv_sca_5
432
433             swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
434                                      10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
435                                    -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
436                                    +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
437                                                          ) * adv_sca_5
438
439          ENDDO
440           
441       ENDIF
442
443       flux_t(0) = 0.0
444       diss_t(0) = 0.0
445       flux_d    = 0.0
446       diss_d    = 0.0
447!       
448!--    Now compute the fluxes and tendency terms for the horizontal and
449!--    vertical parts up to the top of the highest topography.
450       DO  k = nzb+1, nzb_max
451!
452!--       Note: It is faster to conduct all multiplications explicitly, e.g.
453!--       * adv_sca_5 ... than to determine a factor and multiplicate the
454!--       flux at the end.
455
456          ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
457          ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
458          ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
459
460          u_comp    = u(k,j,i+1) - u_gtrans
461          flux_r(k) = u_comp * (                                            &
462                     ( 37.0 * ibit2 * adv_sca_5                             &
463                  +     7.0 * ibit1 * adv_sca_3                             &
464                  +           ibit0 * adv_sca_1                             &
465                     ) *                                                    &
466                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
467              -      (  8.0 * ibit2 * adv_sca_5                             &
468                  +           ibit1 * adv_sca_3                             &
469                     ) *                                                    &
470                             ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
471              +      (        ibit2 * adv_sca_5                             &
472                     ) *                                                    &
473                             ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
474                               )
475
476          diss_r(k) = -ABS( u_comp ) * (                                    &
477                     ( 10.0 * ibit2 * adv_sca_5                             &
478                  +     3.0 * ibit1 * adv_sca_3                             &
479                  +           ibit0 * adv_sca_1                             &
480                     ) *                                                    &
481                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
482              -      (  5.0 * ibit2 * adv_sca_5                             &
483                  +           ibit1 * adv_sca_3                             &
484                     ) *                                                    &
485                             ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
486              +      (        ibit2 * adv_sca_5                             &
487                     ) *                                                    &
488                             ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
489                                       )
490
491          ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
492          ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
493          ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
494
495          v_comp    = v(k,j+1,i) - v_gtrans
496          flux_n(k) = v_comp * (                                            &
497                     ( 37.0 * ibit5 * adv_sca_5                             &
498                  +     7.0 * ibit4 * adv_sca_3                             &
499                  +           ibit3 * adv_sca_1                             &
500                     ) *                                                    &
501                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
502              -      (  8.0 * ibit5 * adv_sca_5                             &
503                  +           ibit4 * adv_sca_3                             &
504                     ) *                                                    &
505                             ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
506              +      (        ibit5 * adv_sca_5                             &
507                     ) *                                                    &
508                             ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
509                               )
510
511          diss_n(k) = -ABS( v_comp ) * (                                    &
512                     ( 10.0 * ibit5 * adv_sca_5                             &
513                  +     3.0 * ibit4 * adv_sca_3                             &
514                  +           ibit3 * adv_sca_1                             &
515                     ) *                                                    &
516                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
517              -      (  5.0 * ibit5 * adv_sca_5                             &
518                  +           ibit4 * adv_sca_3                             &
519                     ) *                                                    &
520                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
521              +      (        ibit5 * adv_sca_5                             &
522                     ) *                                                    &
523                             ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
524                                       )
525!
526!--       k index has to be modified near bottom and top, else array
527!--       subscripts will be exceeded.
528          ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
529          ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
530          ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
531
532          k_ppp = k + 3 * ibit8
533          k_pp  = k + 2 * ( 1 - ibit6  )
534          k_mm  = k - 2 * ibit8
535
536
537          flux_t(k) = w(k,j,i) * (                                          &
538                     ( 37.0 * ibit8 * adv_sca_5                             &
539                  +     7.0 * ibit7 * adv_sca_3                             &
540                  +           ibit6 * adv_sca_1                             &
541                     ) *                                                    &
542                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
543              -      (  8.0 * ibit8 * adv_sca_5                             &
544                  +           ibit7 * adv_sca_3                             &
545                     ) *                                                    &
546                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
547              +      (        ibit8 * adv_sca_5                             &
548                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
549                                 )
550
551          diss_t(k) = -ABS( w(k,j,i) ) * (                                  &
552                     ( 10.0 * ibit8 * adv_sca_5                             &
553                  +     3.0 * ibit7 * adv_sca_3                             &
554                  +           ibit6 * adv_sca_1                             &
555                     ) *                                                    &
556                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
557              -      (  5.0 * ibit8 * adv_sca_5                             &
558                  +           ibit7 * adv_sca_3                             &
559                     ) *                                                    &
560                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
561              +      (        ibit8 * adv_sca_5                             &
562                     ) *                                                    &
563                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
564                                         )
565!
566!--       Calculate the divergence of the velocity field. A respective
567!--       correction is needed to overcome numerical instabilities caused
568!--       by a not sufficient reduction of divergences near topography.
569          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
570                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
571                        + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
572
573          tend(k,j,i) = tend(k,j,i) - (                                       &
574                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
575                          swap_diss_x_local(k,j,tn)            ) * ddx        &
576                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
577                          swap_diss_y_local(k,tn)              ) * ddy        &
578                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
579                                                               ) * ddzw(k)    &
580                                      ) + sk(k,j,i) * div
581
582          swap_flux_y_local(k,tn)   = flux_n(k)
583          swap_diss_y_local(k,tn)   = diss_n(k)
584          swap_flux_x_local(k,j,tn) = flux_r(k)
585          swap_diss_x_local(k,j,tn) = diss_r(k)
586          flux_d                    = flux_t(k)
587          diss_d                    = diss_t(k)
588
589       ENDDO
590!
591!--    Now compute the fluxes and tendency terms for the horizontal and
592!--    vertical parts above the top of the highest topography. No degradation
593!--    for the horizontal parts, but for the vertical it is stell needed.
594       DO  k = nzb_max+1, nzt
595
596          u_comp    = u(k,j,i+1) - u_gtrans
597          flux_r(k) = u_comp * (                                            &
598                      37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                  &
599                    -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
600                    +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
601          diss_r(k) = -ABS( u_comp ) * (                                    &
602                      10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                  &
603                    -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
604                    +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
605
606          v_comp    = v(k,j+1,i) - v_gtrans
607          flux_n(k) = v_comp * (                                            &
608                      37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                  &
609                    -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
610                    +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
611          diss_n(k) = -ABS( v_comp ) * (                                    &
612                      10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                  &
613                    -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                  &
614                    +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
615!
616!--       k index has to be modified near bottom and top, else array
617!--       subscripts will be exceeded.
618          ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
619          ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
620          ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
621
622          k_ppp = k + 3 * ibit8
623          k_pp  = k + 2 * ( 1 - ibit6  )
624          k_mm  = k - 2 * ibit8
625
626
627          flux_t(k) = w(k,j,i) * (                                          &
628                     ( 37.0 * ibit8 * adv_sca_5                             &
629                  +     7.0 * ibit7 * adv_sca_3                             &
630                  +           ibit6 * adv_sca_1                             &
631                     ) *                                                    &
632                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
633              -      (  8.0 * ibit8 * adv_sca_5                             &
634                  +           ibit7 * adv_sca_3                             &
635                     ) *                                                    &
636                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
637              +      (        ibit8 * adv_sca_5                             &
638                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
639                                 )
640
641          diss_t(k) = -ABS( w(k,j,i) ) * (                                  &
642                     ( 10.0 * ibit8 * adv_sca_5                             &
643                  +     3.0 * ibit7 * adv_sca_3                             &
644                  +           ibit6 * adv_sca_1                             &
645                     ) *                                                    &
646                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
647              -      (  5.0 * ibit8 * adv_sca_5                             &
648                  +           ibit7 * adv_sca_3                             &
649                     ) *                                                    &
650                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
651              +      (        ibit8 * adv_sca_5                             &
652                     ) *                                                    &
653                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
654                                         )
655!
656!--       Calculate the divergence of the velocity field. A respective
657!--       correction is needed to overcome numerical instabilities introduced
658!--       by a not sufficient reduction of divergences near topography.
659          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
660                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
661                        + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
662
663          tend(k,j,i) = tend(k,j,i) - (                                       &
664                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
665                          swap_diss_x_local(k,j,tn)            ) * ddx        &
666                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
667                          swap_diss_y_local(k,tn)              ) * ddy        &
668                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
669                                                               ) * ddzw(k)    &
670                                      ) + sk(k,j,i) * div
671
672          swap_flux_y_local(k,tn)   = flux_n(k)
673          swap_diss_y_local(k,tn)   = diss_n(k)
674          swap_flux_x_local(k,j,tn) = flux_r(k)
675          swap_diss_x_local(k,j,tn) = diss_r(k)
676          flux_d                    = flux_t(k)
677          diss_d                    = diss_t(k)
678
679       ENDDO
680
681!
682!--    Evaluation of statistics
683       SELECT CASE ( sk_char )
684
685          CASE ( 'pt' )
686
687             DO  k = nzb, nzt
688                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                &
689                                       ( flux_t(k) + diss_t(k) )               &
690                                 * weight_substep(intermediate_timestep_count)
691             ENDDO
692             
693          CASE ( 'sa' )
694
695             DO  k = nzb, nzt
696                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                &
697                                       ( flux_t(k) + diss_t(k) )               &
698                                 * weight_substep(intermediate_timestep_count)
699             ENDDO
700             
701          CASE ( 'q' )
702
703             DO  k = nzb, nzt
704                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                 &
705                                      ( flux_t(k) + diss_t(k) )                &
706                                 * weight_substep(intermediate_timestep_count)
707             ENDDO
708
709         END SELECT
710
711    END SUBROUTINE advec_s_ws_ij
712
713
714
715
716!------------------------------------------------------------------------------!
717! Advection of u-component - Call for grid point i,j
718!------------------------------------------------------------------------------!
719    SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
720
721       USE arrays_3d
722       USE constants
723       USE control_parameters
724       USE grid_variables
725       USE indices
726       USE statistics
727
728       IMPLICIT NONE
729
730       INTEGER ::  i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15,  &
731                   ibit16, ibit17, i_omp, j, k, k_mm, k_pp, k_ppp, tn
732       REAL    ::  diss_d, div, flux_d, gu, gv, u_comp_l, v_comp, w_comp
733       REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r,  &
734                                     flux_t, u_comp
735
736       gu = 2.0 * u_gtrans
737       gv = 2.0 * v_gtrans
738!
739!--    Compute southside fluxes for the respective boundary of PE
740       IF ( j == nys  )  THEN
741       
742          DO  k = nzb+1, nzb_max
743
744             ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
745             ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
746             ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
747
748             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
749             flux_s_u(k,tn) = v_comp * (                                      &
750                            ( 37.0 * ibit14 * adv_mom_5                       &
751                         +     7.0 * ibit13 * adv_mom_3                       &
752                         +           ibit12 * adv_mom_1                       &
753                            ) *                                               &
754                                     ( u(k,j,i)   + u(k,j-1,i) )              &
755                     -      (  8.0 * ibit14 * adv_mom_5                       &
756                         +           ibit13 * adv_mom_3                       &
757                            ) *                                               &
758                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
759                     +      (        ibit14 * adv_mom_5                       &
760                            ) *                                               &
761                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
762                                       )
763
764             diss_s_u(k,tn) = - ABS ( v_comp ) * (                            &
765                            ( 10.0 * ibit14 * adv_mom_5                       &
766                         +     3.0 * ibit13 * adv_mom_3                       &
767                         +           ibit12 * adv_mom_1                       &
768                            ) *                                               &
769                                     ( u(k,j,i)   - u(k,j-1,i) )              &
770                     -      (  5.0 * ibit14 * adv_mom_5                       &
771                         +           ibit13 * adv_mom_3                       &
772                            ) *                                               &
773                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
774                     +      (        ibit14 * adv_mom_5                       &
775                            ) *                                               &
776                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
777                                                 )
778
779          ENDDO
780
781          DO  k = nzb_max+1, nzt
782
783             v_comp         = v(k,j,i) + v(k,j,i-1) - gv
784             flux_s_u(k,tn) = v_comp * (                                      &
785                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
786                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
787                         +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
788             diss_s_u(k,tn) = - ABS(v_comp) * (                               &
789                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
790                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
791                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
792
793          ENDDO
794         
795       ENDIF
796!
797!--    Compute leftside fluxes for the respective boundary of PE
798       IF ( i == i_omp )  THEN
799       
800          DO  k = nzb+1, nzb_max
801
802             ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
803             ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
804             ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
805
806             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
807             flux_l_u(k,j,tn) = u_comp_l * (                                   &
808                              ( 37.0 * ibit11 * adv_mom_5                      &
809                           +     7.0 * ibit10 * adv_mom_3                      &
810                           +           ibit9  * adv_mom_1                      &
811                              ) *                                              &
812                                     ( u(k,j,i)   + u(k,j,i-1) )               &
813                       -      (  8.0 * ibit11 * adv_mom_5                      &
814                           +           ibit10 * adv_mom_3                      &
815                              ) *                                              &
816                                     ( u(k,j,i+1) + u(k,j,i-2) )               &
817                       +      (        ibit11 * adv_mom_5                      &
818                              ) *                                              &
819                                     ( u(k,j,i+2) + u(k,j,i-3) )               &
820                                           )
821
822              diss_l_u(k,j,tn) = - ABS( u_comp_l ) * (                         &
823                              ( 10.0 * ibit11 * adv_mom_5                      &
824                           +     3.0 * ibit10 * adv_mom_3                      &
825                           +           ibit9  * adv_mom_1                      &
826                              ) *                                              &
827                                     ( u(k,j,i)   - u(k,j,i-1) )               &
828                       -      (  5.0 * ibit11 * adv_mom_5                      &
829                           +           ibit10 * adv_mom_3                      &
830                              ) *                                              &
831                                     ( u(k,j,i+1) - u(k,j,i-2) )               &
832                       +      (        ibit11 * adv_mom_5                      &
833                              ) *                                              &
834                                     ( u(k,j,i+2) - u(k,j,i-3) )               &
835                                                     )
836
837          ENDDO
838
839          DO  k = nzb_max+1, nzt
840
841             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
842             flux_l_u(k,j,tn) = u_comp_l * (                                   &
843                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )                &
844                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )                &
845                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
846             diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                            &
847                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )                &
848                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )                &
849                           +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
850
851          ENDDO
852         
853       ENDIF
854
855       flux_t(0) = 0.0
856       diss_t(0) = 0.0
857       flux_d    = 0.0
858       diss_d    = 0.0
859!
860!--    Now compute the fluxes tendency terms for the horizontal and
861!--    vertical parts.
862       DO  k = nzb+1, nzb_max
863
864          ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
865          ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
866          ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
867
868          u_comp(k) = u(k,j,i+1) + u(k,j,i)
869          flux_r(k) = ( u_comp(k) - gu ) * (                                 &
870                     ( 37.0 * ibit11 * adv_mom_5                             &
871                  +     7.0 * ibit10 * adv_mom_3                             &
872                  +           ibit9  * adv_mom_1                             &
873                     ) *                                                     &
874                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
875              -      (  8.0 * ibit11 * adv_mom_5                             &
876                  +           ibit10 * adv_mom_3                             &
877                     ) *                                                     &
878                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
879              +      (        ibit11 * adv_mom_5                             &
880                     ) *                                                     &
881                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
882                                           )
883
884          diss_r(k) = - ABS( u_comp(k) - gu ) * (                            &
885                     ( 10.0 * ibit11 * adv_mom_5                             &
886                  +     3.0 * ibit10 * adv_mom_3                             &
887                  +           ibit9  * adv_mom_1                             &
888                     ) *                                                     &
889                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
890              -      (  5.0 * ibit11 * adv_mom_5                             &
891                  +           ibit10 * adv_mom_3                             &
892                     ) *                                                     &
893                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
894              +      (        ibit11 * adv_mom_5                             &
895                     ) *                                                     &
896                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
897                                                )
898
899          ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
900          ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
901          ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
902
903          v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
904          flux_n(k) = v_comp * (                                             &
905                     ( 37.0 * ibit14 * adv_mom_5                             &
906                  +     7.0 * ibit13 * adv_mom_3                             &
907                  +           ibit12 * adv_mom_1                             &
908                     ) *                                                     &
909                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
910              -      (  8.0 * ibit14 * adv_mom_5                             &
911                  +           ibit13 * adv_mom_3                             &
912                     ) *                                                     &
913                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
914              +      (        ibit14 * adv_mom_5                             &
915                     ) *                                                     &
916                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
917                               )
918
919          diss_n(k) = - ABS ( v_comp ) * (                                   &
920                     ( 10.0 * ibit14 * adv_mom_5                             &
921                  +     3.0 * ibit13 * adv_mom_3                             &
922                  +           ibit12 * adv_mom_1                             &
923                     ) *                                                     &
924                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
925              -      (  5.0 * ibit14 * adv_mom_5                             &
926                  +           ibit13 * adv_mom_3                             &
927                     ) *                                                     &
928                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
929              +      (        ibit14 * adv_mom_5                             &
930                     ) *                                                     &
931                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
932                                         )
933!
934!--       k index has to be modified near bottom and top, else array
935!--       subscripts will be exceeded.
936          ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
937          ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
938          ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
939
940          k_ppp = k + 3 * ibit17
941          k_pp  = k + 2 * ( 1 - ibit15 )
942          k_mm  = k - 2 * ibit17
943
944          w_comp    = w(k,j,i) + w(k,j,i-1)
945          flux_t(k) = w_comp  * (                                            &
946                     ( 37.0 * ibit17 * adv_mom_5                             &
947                  +     7.0 * ibit16 * adv_mom_3                             &
948                  +           ibit15 * adv_mom_1                             &
949                     ) *                                                     &
950                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
951              -      (  8.0 * ibit17 * adv_mom_5                             &
952                  +           ibit16 * adv_mom_3                             &
953                     ) *                                                     &
954                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
955              +      (        ibit17 * adv_mom_5                             &
956                     ) *                                                     &
957                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
958                                 )
959
960          diss_t(k) = - ABS( w_comp ) * (                                    &
961                     ( 10.0 * ibit17 * adv_mom_5                             &
962                  +     3.0 * ibit16 * adv_mom_3                             &
963                  +           ibit15 * adv_mom_1                             &
964                     ) *                                                     &
965                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
966              -      (  5.0 * ibit17 * adv_mom_5                             &
967                  +           ibit16 * adv_mom_3                             &
968                     ) *                                                     &
969                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
970              +      (        ibit17 * adv_mom_5                             &
971                     ) *                                                     &
972                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
973                                         )
974!
975!--       Calculate the divergence of the velocity field. A respective
976!--       correction is needed to overcome numerical instabilities introduced
977!--       by a not sufficient reduction of divergences near topography.
978          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx      &
979               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy      &
980               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)  &
981                ) * 0.5
982
983          tend(k,j,i) = tend(k,j,i) - (                                       &
984                            ( flux_r(k) + diss_r(k)                           &
985                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
986                          + ( flux_n(k) + diss_n(k)                           &
987                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
988                          + ( flux_t(k) + diss_t(k)                           &
989                          -   flux_d    - diss_d                  ) * ddzw(k) &
990                                       ) + div * u(k,j,i)
991
992           flux_l_u(k,j,tn) = flux_r(k)
993           diss_l_u(k,j,tn) = diss_r(k)
994           flux_s_u(k,tn)   = flux_n(k)
995           diss_s_u(k,tn)   = diss_n(k)
996           flux_d           = flux_t(k)
997           diss_d           = diss_t(k)
998!
999!--        Statistical Evaluation of u'u'. The factor has to be applied for
1000!--        right evaluation when gallilei_trans = .T. .
1001           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
1002                              + ( flux_r(k) *                                 &
1003                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
1004                              / ( u_comp(k) - gu + 1.0E-20      )             &
1005                              +   diss_r(k) *                                 &
1006                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
1007                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
1008                              *   weight_substep(intermediate_timestep_count)
1009!
1010!--        Statistical Evaluation of w'u'.
1011           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                        &
1012                              + ( flux_t(k) + diss_t(k) )                     &
1013                              *   weight_substep(intermediate_timestep_count)
1014       ENDDO
1015
1016       DO  k = nzb_max+1, nzt
1017
1018          u_comp(k) = u(k,j,i+1) + u(k,j,i)
1019          flux_r(k) = ( u_comp(k) - gu ) * (                                  &
1020                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
1021                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
1022                       +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
1023             diss_r(k) = - ABS( u_comp(k) - gu ) * (                          &
1024                         10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
1025                       -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
1026                       +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
1027
1028             v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
1029             flux_n(k) = v_comp * (                                           &
1030                         37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
1031                       -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                   &
1032                       +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
1033             diss_n(k) = - ABS( v_comp ) * (                                  &
1034                         10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
1035                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
1036                       +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
1037!
1038!--       k index has to be modified near bottom and top, else array
1039!--       subscripts will be exceeded.
1040          ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
1041          ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
1042          ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
1043
1044          k_ppp = k + 3 * ibit17
1045          k_pp  = k + 2 * ( 1 - ibit15 )
1046          k_mm  = k - 2 * ibit17
1047
1048          w_comp    = w(k,j,i) + w(k,j,i-1)
1049          flux_t(k) = w_comp  * (                                            &
1050                     ( 37.0 * ibit17 * adv_mom_5                             &
1051                  +     7.0 * ibit16 * adv_mom_3                             &
1052                  +           ibit15 * adv_mom_1                             &
1053                     ) *                                                     &
1054                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
1055              -      (  8.0 * ibit17 * adv_mom_5                             &
1056                  +           ibit16 * adv_mom_3                             &
1057                     ) *                                                     &
1058                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
1059              +      (        ibit17 * adv_mom_5                             &
1060                     ) *                                                     &
1061                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
1062                                 )
1063
1064          diss_t(k) = - ABS( w_comp ) * (                                    &
1065                     ( 10.0 * ibit17 * adv_mom_5                             &
1066                  +     3.0 * ibit16 * adv_mom_3                             &
1067                  +           ibit15 * adv_mom_1                             &
1068                     ) *                                                     &
1069                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
1070              -      (  5.0 * ibit17 * adv_mom_5                             &
1071                  +           ibit16 * adv_mom_3                             &
1072                     ) *                                                     &
1073                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
1074              +      (        ibit17 * adv_mom_5                             &
1075                     ) *                                                     &
1076                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
1077                                         )
1078!
1079!--       Calculate the divergence of the velocity field. A respective
1080!--       correction is needed to overcome numerical instabilities introduced
1081!--       by a not sufficient reduction of divergences near topography.
1082          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx      &
1083               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy      &
1084               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)  &
1085                ) * 0.5
1086
1087          tend(k,j,i) = tend(k,j,i) - (                                       &
1088                            ( flux_r(k) + diss_r(k)                           &
1089                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx     &
1090                          + ( flux_n(k) + diss_n(k)                           &
1091                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy     &
1092                          + ( flux_t(k) + diss_t(k)                           &
1093                          -   flux_d    - diss_d                  ) * ddzw(k) &
1094                                       ) + div * u(k,j,i)
1095
1096           flux_l_u(k,j,tn) = flux_r(k)
1097           diss_l_u(k,j,tn) = diss_r(k)
1098           flux_s_u(k,tn)   = flux_n(k)
1099           diss_s_u(k,tn)   = diss_n(k)
1100           flux_d           = flux_t(k)
1101           diss_d           = diss_t(k)
1102!
1103!--        Statistical Evaluation of u'u'. The factor has to be applied for
1104!--        right evaluation when gallilei_trans = .T. .
1105           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                          &
1106                              + ( flux_r(k) *                                 &
1107                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
1108                              / ( u_comp(k) - gu + 1.0E-20      )             &
1109                              +   diss_r(k) *                                 &
1110                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
1111                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
1112                              *   weight_substep(intermediate_timestep_count)
1113!
1114!--        Statistical Evaluation of w'u'.
1115           sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                        &
1116                              + ( flux_t(k) + diss_t(k) )                     &
1117                              *   weight_substep(intermediate_timestep_count)
1118       ENDDO
1119
1120       sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
1121
1122
1123
1124    END SUBROUTINE advec_u_ws_ij
1125
1126
1127
1128!-----------------------------------------------------------------------------!
1129! Advection of v-component - Call for grid point i,j
1130!-----------------------------------------------------------------------------!
1131   SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
1132
1133       USE arrays_3d
1134       USE constants
1135       USE control_parameters
1136       USE grid_variables
1137       USE indices
1138       USE statistics
1139
1140       IMPLICIT NONE
1141
1142       INTEGER  ::  i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, &
1143                    ibit25, ibit26, i_omp, j, k, k_mm, k_pp, k_ppp, tn
1144       REAL     ::  diss_d, div, flux_d, gu, gv, u_comp, v_comp_l, w_comp
1145       REAL, DIMENSION(nzb:nzt+1)  :: diss_n, diss_r, diss_t, flux_n, flux_r,  &
1146                                      flux_t, v_comp
1147
1148       gu = 2.0 * u_gtrans
1149       gv = 2.0 * v_gtrans
1150
1151!       
1152!--    Compute leftside fluxes for the respective boundary.
1153       IF ( i == i_omp )  THEN
1154
1155          DO  k = nzb+1, nzb_max
1156
1157             ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
1158             ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
1159             ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
1160
1161             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
1162             flux_l_v(k,j,tn) = u_comp * (                                     &
1163                              ( 37.0 * ibit20 * adv_mom_5                      &
1164                           +     7.0 * ibit19 * adv_mom_3                      &
1165                           +           ibit18 * adv_mom_1                      &
1166                              ) *                                              &
1167                                     ( v(k,j,i)   + v(k,j,i-1) )               &
1168                       -      (  8.0 * ibit20 * adv_mom_5                      &
1169                           +           ibit19 * adv_mom_3                      &
1170                              ) *                                              &
1171                                     ( v(k,j,i+1) + v(k,j,i-2) )               &
1172                       +      (        ibit20 * adv_mom_5                      &
1173                              ) *                                              &
1174                                     ( v(k,j,i+2) + v(k,j,i-3) )               &
1175                                         )
1176
1177              diss_l_v(k,j,tn) = - ABS( u_comp ) * (                           &
1178                              ( 10.0 * ibit20 * adv_mom_5                      &
1179                           +     3.0 * ibit19 * adv_mom_3                      &
1180                           +           ibit18 * adv_mom_1                      &
1181                              ) *                                              &
1182                                     ( v(k,j,i)   - v(k,j,i-1) )               &
1183                       -      (  5.0 * ibit20 * adv_mom_5                      &
1184                           +           ibit19 * adv_mom_3                      &
1185                              ) *                                              &
1186                                     ( v(k,j,i+1) - v(k,j,i-2) )               &
1187                       +      (        ibit20 * adv_mom_5                      &
1188                              ) *                                              &
1189                                     ( v(k,j,i+2) - v(k,j,i-3) )               &
1190                                                   )
1191
1192          ENDDO
1193
1194          DO  k = nzb_max+1, nzt
1195
1196             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
1197             flux_l_v(k,j,tn) = u_comp * (                                    &
1198                             37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
1199                           -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
1200                           +        ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
1201             diss_l_v(k,j,tn) = - ABS( u_comp ) * (                           &
1202                             10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
1203                           -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
1204                           +        ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
1205
1206          ENDDO
1207         
1208       ENDIF
1209!
1210!--    Compute southside fluxes for the respective boundary.
1211       IF ( j == nysv )  THEN
1212       
1213          DO  k = nzb+1, nzb_max
1214
1215             ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
1216             ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
1217             ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
1218
1219             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
1220             flux_s_v(k,tn) = v_comp_l * (                                    &
1221                            ( 37.0 * ibit23 * adv_mom_5                       &
1222                         +     7.0 * ibit22 * adv_mom_3                       &
1223                         +           ibit21 * adv_mom_1                       &
1224                            ) *                                               &
1225                                     ( v(k,j,i)   + v(k,j-1,i) )              &
1226                     -      (  8.0 * ibit23 * adv_mom_5                       &
1227                         +           ibit22 * adv_mom_3                       &
1228                            ) *                                               &
1229                                     ( v(k,j+1,i) + v(k,j-2,i) )              &
1230                     +      (        ibit23 * adv_mom_5                       &
1231                            ) *                                               &
1232                                     ( v(k,j+2,i) + v(k,j-3,i) )              &
1233                                         )
1234
1235             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
1236                            ( 10.0 * ibit23 * adv_mom_5                       &
1237                         +     3.0 * ibit22 * adv_mom_3                       &
1238                         +           ibit21 * adv_mom_1                       &
1239                            ) *                                               &
1240                                     ( v(k,j,i)   - v(k,j-1,i) )              &
1241                     -      (  5.0 * ibit23 * adv_mom_5                       &
1242                         +           ibit22 * adv_mom_3                       &
1243                            ) *                                               &
1244                                     ( v(k,j+1,i) - v(k,j-2,i) )              &
1245                     +      (        ibit23 * adv_mom_5                       &
1246                            ) *                                               &
1247                                     ( v(k,j+2,i) - v(k,j-3,i) )              &
1248                                                 )
1249
1250          ENDDO
1251
1252          DO  k = nzb_max+1, nzt
1253
1254             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
1255             flux_s_v(k,tn) = v_comp_l * (                                    &
1256                           37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
1257                         -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
1258                         +        ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
1259             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                           &
1260                           10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
1261                         -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
1262                         +        ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
1263
1264          ENDDO
1265         
1266       ENDIF
1267
1268       flux_t(0) = 0.0
1269       diss_t(0) = 0.0
1270       flux_d    = 0.0
1271       diss_d    = 0.0
1272!
1273!--    Now compute the fluxes and tendency terms for the horizontal and
1274!--    verical parts.
1275       DO  k = nzb+1, nzb_max
1276
1277          ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
1278          ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
1279          ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
1280 
1281          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1282          flux_r(k) = u_comp * (                                             &
1283                     ( 37.0 * ibit20 * adv_mom_5                             &
1284                  +     7.0 * ibit19 * adv_mom_3                             &
1285                  +           ibit18 * adv_mom_1                             &
1286                     ) *                                                     &
1287                                 ( v(k,j,i+1) + v(k,j,i)   )                 &
1288              -      (  8.0 * ibit20 * adv_mom_5                             &
1289                  +           ibit19 * adv_mom_3                             &
1290                     ) *                                                     &
1291                                 ( v(k,j,i+2) + v(k,j,i-1) )                 &
1292              +      (        ibit20 * adv_mom_5                             &
1293                     ) *                                                     &
1294                                 ( v(k,j,i+3) + v(k,j,i-2) )                 &
1295                               )
1296
1297          diss_r(k) = - ABS( u_comp ) * (                                    &
1298                     ( 10.0 * ibit20 * adv_mom_5                             &
1299                  +     3.0 * ibit19 * adv_mom_3                             &
1300                  +           ibit18 * adv_mom_1                             &
1301                     ) *                                                     &
1302                                 ( v(k,j,i+1) - v(k,j,i)  )                  &
1303              -      (  5.0 * ibit20 * adv_mom_5                             &
1304                  +           ibit19 * adv_mom_3                             &
1305                     ) *                                                     &
1306                                 ( v(k,j,i+2) - v(k,j,i-1) )                 &
1307              +      (        ibit20 * adv_mom_5                             &
1308                     ) *                                                     &
1309                                 ( v(k,j,i+3) - v(k,j,i-2) )                 &
1310                                        )
1311
1312          ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
1313          ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
1314          ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
1315
1316
1317          v_comp(k) = v(k,j+1,i) + v(k,j,i)
1318          flux_n(k) = ( v_comp(k) - gv ) * (                                 &
1319                     ( 37.0 * ibit23 * adv_mom_5                             &
1320                  +     7.0 * ibit22 * adv_mom_3                             &
1321                  +           ibit21 * adv_mom_1                             &
1322                     ) *                                                     &
1323                                 ( v(k,j+1,i) + v(k,j,i)   )                 &
1324              -      (  8.0 * ibit23 * adv_mom_5                             &
1325                  +           ibit22 * adv_mom_3                             &
1326                     ) *                                                     &
1327                                 ( v(k,j+2,i) + v(k,j-1,i) )                 &
1328              +      (        ibit23 * adv_mom_5                             &
1329                     ) *                                                     &
1330                                 ( v(k,j+3,i) + v(k,j-2,i) )                 &
1331                               )
1332
1333          diss_n(k) = - ABS( v_comp(k) - gv ) * (                            &
1334                     ( 10.0 * ibit23 * adv_mom_5                             &
1335                  +     3.0 * ibit22 * adv_mom_3                             &
1336                  +           ibit21 * adv_mom_1                             &
1337                     ) *                                                     &
1338                                 ( v(k,j+1,i) - v(k,j,i)  )                  &
1339              -      (  5.0 * ibit23 * adv_mom_5                             &
1340                  +           ibit22 * adv_mom_3                             &
1341                     ) *                                                     &
1342                                 ( v(k,j+2,i) - v(k,j-1,i) )                 &
1343              +      (        ibit23 * adv_mom_5                             &
1344                     ) *                                                     &
1345                                 ( v(k,j+3,i) - v(k,j-2,i) )                 &
1346                                        )
1347!
1348!--       k index has to be modified near bottom and top, else array
1349!--       subscripts will be exceeded.
1350          ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
1351          ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
1352          ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
1353
1354          k_ppp = k + 3 * ibit26
1355          k_pp  = k + 2 * ( 1 - ibit24  )
1356          k_mm  = k - 2 * ibit26
1357
1358          w_comp    = w(k,j-1,i) + w(k,j,i)
1359          flux_t(k) = w_comp  * (                                            &
1360                     ( 37.0 * ibit26 * adv_mom_5                             &
1361                  +     7.0 * ibit25 * adv_mom_3                             &
1362                  +           ibit24 * adv_mom_1                             &
1363                     ) *                                                     &
1364                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
1365              -      (  8.0 * ibit26 * adv_mom_5                             &
1366                  +           ibit25 * adv_mom_3                             &
1367                     ) *                                                     &
1368                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
1369              +      (        ibit26 * adv_mom_5                             &
1370                     ) *                                                     &
1371                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
1372                                 )
1373
1374          diss_t(k) = - ABS( w_comp ) * (                                    &
1375                     ( 10.0 * ibit26 * adv_mom_5                             &
1376                  +     3.0 * ibit25 * adv_mom_3                             &
1377                  +           ibit24 * adv_mom_1                             &
1378                     ) *                                                     &
1379                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
1380              -      (  5.0 * ibit26 * adv_mom_5                             &
1381                  +           ibit25 * adv_mom_3                             &
1382                     ) *                                                     &
1383                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
1384              +      (        ibit26 * adv_mom_5                             &
1385                     ) *                                                     &
1386                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
1387                                         )
1388!
1389!--       Calculate the divergence of the velocity field. A respective
1390!--       correction is needed to overcome numerical instabilities introduced
1391!--       by a not sufficient reduction of divergences near topography.
1392          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx      &
1393               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy      &
1394               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)  &
1395                ) * 0.5
1396
1397          tend(k,j,i) = tend(k,j,i) - (                                       &
1398                         ( flux_r(k) + diss_r(k)                              &
1399                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx      &
1400                       + ( flux_n(k) + diss_n(k)                              &
1401                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
1402                       + ( flux_t(k) + diss_t(k)                              &
1403                       -   flux_d    - diss_d                    ) * ddzw(k)  &
1404                                      ) + v(k,j,i) * div
1405
1406           flux_l_v(k,j,tn) = flux_r(k)
1407           diss_l_v(k,j,tn) = diss_r(k)
1408           flux_s_v(k,tn)   = flux_n(k)
1409           diss_s_v(k,tn)   = diss_n(k)
1410           flux_d           = flux_t(k)
1411           diss_d           = diss_t(k)
1412
1413!
1414!--        Statistical Evaluation of v'v'. The factor has to be applied for
1415!--        right evaluation when gallilei_trans = .T. .
1416           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
1417             + ( flux_n(k)                                                    &
1418             * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
1419             / ( v_comp(k) - gv + 1.0E-20 )                                   &
1420             +   diss_n(k)                                                    &
1421             *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
1422             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
1423             *   weight_substep(intermediate_timestep_count)
1424!
1425!--        Statistical Evaluation of w'v'.
1426           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                        &
1427                              + ( flux_t(k) + diss_t(k) )                     &
1428                              *   weight_substep(intermediate_timestep_count)
1429
1430       ENDDO
1431
1432       DO  k = nzb_max+1, nzt
1433
1434          u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1435          flux_r(k) = u_comp * (                                           &
1436                      37.0 * ( v(k,j,i+1) + v(k,j,i)   )                   &
1437                    -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                   &
1438                    +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1439
1440          diss_r(k) = - ABS( u_comp ) * (                                  &
1441                      10.0 * ( v(k,j,i+1) - v(k,j,i) )                     &
1442                    -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                   &
1443                    +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1444
1445
1446          v_comp(k) = v(k,j+1,i) + v(k,j,i)
1447          flux_n(k) = ( v_comp(k) - gv ) * (                               &
1448                      37.0 * ( v(k,j+1,i) + v(k,j,i)   )                   &
1449                    -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                   &
1450                      +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1451
1452          diss_n(k) = - ABS( v_comp(k) - gv ) * (                          &
1453                      10.0 * ( v(k,j+1,i) - v(k,j,i)   )                   &
1454                    -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                   &
1455                    +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
1456!
1457!--       k index has to be modified near bottom and top, else array
1458!--       subscripts will be exceeded.
1459          ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
1460          ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
1461          ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
1462
1463          k_ppp = k + 3 * ibit26
1464          k_pp  = k + 2 * ( 1 - ibit24  )
1465          k_mm  = k - 2 * ibit26
1466
1467          w_comp    = w(k,j-1,i) + w(k,j,i)
1468          flux_t(k) = w_comp  * (                                            &
1469                     ( 37.0 * ibit26 * adv_mom_5                             &
1470                  +     7.0 * ibit25 * adv_mom_3                             &
1471                  +           ibit24 * adv_mom_1                             &
1472                     ) *                                                     &
1473                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
1474              -      (  8.0 * ibit26 * adv_mom_5                             &
1475                  +           ibit25 * adv_mom_3                             &
1476                     ) *                                                     &
1477                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
1478              +      (        ibit26 * adv_mom_5                             &
1479                     ) *                                                     &
1480                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
1481                                 )
1482
1483          diss_t(k) = - ABS( w_comp ) * (                                    &
1484                     ( 10.0 * ibit26 * adv_mom_5                             &
1485                  +     3.0 * ibit25 * adv_mom_3                             &
1486                  +           ibit24 * adv_mom_1                             &
1487                     ) *                                                     &
1488                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
1489              -      (  5.0 * ibit26 * adv_mom_5                             &
1490                  +           ibit25 * adv_mom_3                             &
1491                     ) *                                                     &
1492                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
1493              +      (        ibit26 * adv_mom_5                             &
1494                     ) *                                                     &
1495                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
1496                                         )
1497!
1498!--       Calculate the divergence of the velocity field. A respective
1499!--       correction is needed to overcome numerical instabilities introduced
1500!--       by a not sufficient reduction of divergences near topography.
1501          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx      &
1502               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy      &
1503               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)  &
1504                ) * 0.5
1505
1506          tend(k,j,i) = tend(k,j,i) - (                                       &
1507                         ( flux_r(k) + diss_r(k)                              &
1508                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx      &
1509                       + ( flux_n(k) + diss_n(k)                              &
1510                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy      &
1511                       + ( flux_t(k) + diss_t(k)                              &
1512                       -   flux_d    - diss_d                    ) * ddzw(k)  &
1513                                      ) + v(k,j,i) * div
1514
1515           flux_l_v(k,j,tn) = flux_r(k)
1516           diss_l_v(k,j,tn) = diss_r(k)
1517           flux_s_v(k,tn)   = flux_n(k)
1518           diss_s_v(k,tn)   = diss_n(k)
1519           flux_d           = flux_t(k)
1520           diss_d           = diss_t(k)
1521
1522!
1523!--        Statistical Evaluation of v'v'. The factor has to be applied for
1524!--        right evaluation when gallilei_trans = .T. .
1525           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                          &
1526             + ( flux_n(k)                                                    &
1527             * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
1528             / ( v_comp(k) - gv + 1.0E-20 )                                   &
1529             +   diss_n(k)                                                    &
1530             *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
1531             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
1532             *   weight_substep(intermediate_timestep_count)
1533!
1534!--        Statistical Evaluation of w'v'.
1535           sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                         &
1536                              + ( flux_t(k) + diss_t(k) )                      &
1537                              *   weight_substep(intermediate_timestep_count)
1538
1539       ENDDO
1540       sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
1541
1542
1543    END SUBROUTINE advec_v_ws_ij
1544
1545
1546
1547!------------------------------------------------------------------------------!
1548! Advection of w-component - Call for grid point i,j
1549!------------------------------------------------------------------------------!
1550    SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
1551
1552       USE arrays_3d
1553       USE constants
1554       USE control_parameters
1555       USE grid_variables
1556       USE indices
1557       USE statistics
1558
1559       IMPLICIT NONE
1560
1561       INTEGER ::  i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, &
1562                   ibit34, ibit35, i_omp, j, k, k_mm, k_pp, k_ppp, tn
1563       REAL    ::  diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp
1564       REAL, DIMENSION(nzb:nzt+1)  :: diss_n, diss_r, diss_t, flux_n, flux_r, &
1565                                      flux_t
1566
1567       gu = 2.0 * u_gtrans
1568       gv = 2.0 * v_gtrans
1569
1570!
1571!--    Compute southside fluxes for the respective boundary.
1572       IF ( j == nys )  THEN
1573
1574          DO  k = nzb+1, nzb_max
1575             ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
1576             ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
1577             ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
1578
1579             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
1580             flux_s_w(k,tn) = v_comp * (                                      &
1581                            ( 37.0 * ibit32 * adv_mom_5                       &
1582                         +     7.0 * ibit31 * adv_mom_3                       &
1583                         +           ibit30 * adv_mom_1                       &
1584                            ) *                                               &
1585                                     ( w(k,j,i)   + w(k,j-1,i) )              &
1586                     -      (  8.0 * ibit32 * adv_mom_5                       &
1587                         +           ibit31 * adv_mom_3                       &
1588                            ) *                                               &
1589                                     ( w(k,j+1,i) + w(k,j-2,i) )              &
1590                     +      (        ibit32 * adv_mom_5                       &
1591                            ) *                                               &
1592                                     ( w(k,j+2,i) + w(k,j-3,i) )              &
1593                                         )
1594
1595             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
1596                            ( 10.0 * ibit32 * adv_mom_5                       &
1597                         +     3.0 * ibit31 * adv_mom_3                       &
1598                         +           ibit30 * adv_mom_1                       &
1599                            ) *                                               &
1600                                     ( w(k,j,i)   - w(k,j-1,i) )              &
1601                     -      (  5.0 * ibit32 * adv_mom_5                       &
1602                         +           ibit31 * adv_mom_3                       &
1603                            ) *                                               &
1604                                     ( w(k,j+1,i) - w(k,j-2,i) )              &
1605                     +      (        ibit32 * adv_mom_5                       &
1606                            ) *                                               &
1607                                     ( w(k,j+2,i) - w(k,j-3,i) )              &
1608                                                 )
1609
1610          ENDDO
1611
1612          DO  k = nzb_max+1, nzt
1613
1614             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
1615             flux_s_w(k,tn) = v_comp * (                                      &
1616                           37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
1617                         -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
1618                         +        ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
1619             diss_s_w(k,tn) = - ABS( v_comp ) * (                             &
1620                           10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
1621                         -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
1622                         +        ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
1623
1624          ENDDO
1625
1626       ENDIF
1627!
1628!--    Compute leftside fluxes for the respective boundary.
1629       IF ( i == i_omp ) THEN
1630
1631          DO  k = nzb+1, nzb_max
1632
1633             ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
1634             ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
1635             ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
1636
1637             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
1638             flux_l_w(k,j,tn) = u_comp * (                                     &
1639                             ( 37.0 * ibit29 * adv_mom_5                       &
1640                          +     7.0 * ibit28 * adv_mom_3                       &
1641                          +           ibit27 * adv_mom_1                       &
1642                             ) *                                               &
1643                                     ( w(k,j,i)   + w(k,j,i-1) )               &
1644                      -      (  8.0 * ibit29 * adv_mom_5                       &
1645                          +           ibit28 * adv_mom_3                       &
1646                             ) *                                               &
1647                                     ( w(k,j,i+1) + w(k,j,i-2) )               &
1648                      +      (        ibit29 * adv_mom_5                       &
1649                             ) *                                               &
1650                                     ( w(k,j,i+2) + w(k,j,i-3) )               &
1651                                         )
1652
1653               diss_l_w(k,j,tn) = - ABS( u_comp ) * (                          &
1654                             ( 10.0 * ibit29 * adv_mom_5                       &
1655                          +     3.0 * ibit28 * adv_mom_3                       &
1656                          +           ibit27 * adv_mom_1                       &
1657                             ) *                                               &
1658                                     ( w(k,j,i)   - w(k,j,i-1) )               &
1659                      -      (  5.0 * ibit29 * adv_mom_5                       &
1660                          +           ibit28 * adv_mom_3                       &
1661                             ) *                                               &
1662                                     ( w(k,j,i+1) - w(k,j,i-2) )               &
1663                      +      (        ibit29 * adv_mom_5                       &
1664                             ) *                                               &
1665                                     ( w(k,j,i+2) - w(k,j,i-3) )               &
1666                                                    )
1667
1668          ENDDO
1669
1670          DO  k = nzb_max+1, nzt
1671
1672             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
1673             flux_l_w(k,j,tn) = u_comp * (                                    &
1674                            37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
1675                          -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                &
1676                          +        ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
1677             diss_l_w(k,j,tn) = - ABS( u_comp ) * (                           &
1678                            10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
1679                          -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
1680                          +        ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 
1681
1682          ENDDO
1683
1684       ENDIF
1685!
1686!--    The lower flux has to be calculated explicetely for the tendency at
1687!--    the first w-level. For topography wall this is done implicitely by
1688!--    wall_flags_0.
1689       k         = nzb + 1
1690       w_comp    = w(k,j,i) + w(k-1,j,i)
1691       flux_t(0) = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
1692       diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
1693       flux_d    = flux_t(0)
1694       diss_d    = diss_t(0)
1695!
1696!--    Now compute the fluxes and tendency terms for the horizontal
1697!--    and vertical parts.
1698       DO  k = nzb+1, nzb_max
1699
1700          ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
1701          ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
1702          ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
1703
1704          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1705          flux_r(k) = u_comp * (                                             &
1706                     ( 37.0 * ibit29 * adv_mom_5                             &
1707                  +     7.0 * ibit28 * adv_mom_3                             &
1708                  +           ibit27 * adv_mom_1                             &
1709                     ) *                                                     &
1710                                 ( w(k,j,i+1) + w(k,j,i)   )                 &
1711              -      (  8.0 * ibit29 * adv_mom_5                             &
1712                  +           ibit28 * adv_mom_3                             &
1713                     ) *                                                     &
1714                                 ( w(k,j,i+2) + w(k,j,i-1) )                 &
1715              +      (        ibit29 * adv_mom_5                             &
1716                     ) *                                                     &
1717                                 ( w(k,j,i+3) + w(k,j,i-2) )                 &
1718                                )
1719
1720          diss_r(k) = - ABS( u_comp ) * (                                    &
1721                     ( 10.0 * ibit29 * adv_mom_5                             &
1722                  +     3.0 * ibit28 * adv_mom_3                             &
1723                  +           ibit27 * adv_mom_1                             &
1724                     ) *                                                     &
1725                                 ( w(k,j,i+1) - w(k,j,i)  )                  &
1726              -      (  5.0 * ibit29 * adv_mom_5                             &
1727                  +           ibit28 * adv_mom_3                             &
1728                     ) *                                                     &
1729                                 ( w(k,j,i+2) - w(k,j,i-1) )                 &
1730              +      (        ibit29 * adv_mom_5                             &
1731                     ) *                                                     &
1732                                 ( w(k,j,i+3) - w(k,j,i-2) )                 &
1733                                        )
1734
1735          ibit32 = IBITS(wall_flags_0(k,j,i),32,1)
1736          ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
1737          ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
1738
1739          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1740          flux_n(k) = v_comp * (                                             &
1741                     ( 37.0 * ibit32 * adv_mom_5                             &
1742                  +     7.0 * ibit31 * adv_mom_3                             &
1743                  +           ibit30 * adv_mom_1                             &
1744                     ) *                                                     &
1745                                 ( w(k,j+1,i) + w(k,j,i)   )                 &
1746              -      (  8.0 * ibit32 * adv_mom_5                             &
1747                  +           ibit31 * adv_mom_3                             &
1748                     ) *                                                     &
1749                                 ( w(k,j+2,i) + w(k,j-1,i) )                 &
1750              +      (        ibit32 * adv_mom_5                             &
1751                     ) *                                                     &
1752                                 ( w(k,j+3,i) + w(k,j-2,i) )                 &
1753                               )
1754
1755          diss_n(k) = - ABS( v_comp ) * (                                    &
1756                     ( 10.0 * ibit32 * adv_mom_5                             &
1757                  +     3.0 * ibit31 * adv_mom_3                             &
1758                  +           ibit30 * adv_mom_1                             &
1759                     ) *                                                     &
1760                                 ( w(k,j+1,i) - w(k,j,i)  )                  &
1761              -      (  5.0 * ibit32 * adv_mom_5                             &
1762                  +           ibit31 * adv_mom_3                             &
1763                     ) *                                                     &
1764                                 ( w(k,j+2,i) - w(k,j-1,i) )                 &
1765              +      (        ibit32 * adv_mom_5                             &
1766                     ) *                                                     &
1767                                 ( w(k,j+3,i) - w(k,j-2,i) )                 &
1768                                        )
1769!
1770!--       k index has to be modified near bottom and top, else array
1771!--       subscripts will be exceeded.
1772          ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
1773          ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
1774          ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
1775
1776          k_ppp = k + 3 * ibit35
1777          k_pp  = k + 2 * ( 1 - ibit33  )
1778          k_mm  = k - 2 * ibit35
1779
1780          w_comp    = w(k+1,j,i) + w(k,j,i)
1781          flux_t(k) = w_comp  * (                                            &
1782                     ( 37.0 * ibit35 * adv_mom_5                             &
1783                  +     7.0 * ibit34 * adv_mom_3                             &
1784                  +           ibit33 * adv_mom_1                             &
1785                     ) *                                                     &
1786                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
1787              -      (  8.0 * ibit35 * adv_mom_5                             &
1788                  +           ibit34 * adv_mom_3                             &
1789                     ) *                                                     &
1790                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
1791              +      (        ibit35 * adv_mom_5                             &
1792                     ) *                                                     &
1793                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
1794                                 )
1795
1796          diss_t(k) = - ABS( w_comp ) * (                                    &
1797                     ( 10.0 * ibit35 * adv_mom_5                             &
1798                  +     3.0 * ibit34 * adv_mom_3                             &
1799                  +           ibit33 * adv_mom_1                             &
1800                     ) *                                                     &
1801                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
1802              -      (  5.0 * ibit35 * adv_mom_5                             &
1803                  +           ibit34 * adv_mom_3                             &
1804                     ) *                                                     &
1805                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
1806              +      (        ibit35 * adv_mom_5                             &
1807                     ) *                                                     &
1808                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
1809                                         )
1810
1811!
1812!--       Calculate the divergence of the velocity field. A respective
1813!--       correction is needed to overcome numerical instabilities introduced
1814!--       by a not sufficient reduction of divergences near topography.
1815          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
1816              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
1817              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
1818                ) * 0.5
1819
1820          tend(k,j,i) = tend(k,j,i) - (                                       &
1821                      ( flux_r(k) + diss_r(k)                                 &
1822                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
1823                    + ( flux_n(k) + diss_n(k)                                 &
1824                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
1825                    + ( flux_t(k) + diss_t(k)                                 &
1826                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
1827                                      ) + div * w(k,j,i)
1828
1829          flux_l_w(k,j,tn) = flux_r(k)
1830          diss_l_w(k,j,tn) = diss_r(k)
1831          flux_s_w(k,tn)   = flux_n(k)
1832          diss_s_w(k,tn)   = diss_n(k)
1833          flux_d           = flux_t(k)
1834          diss_d           = diss_t(k)
1835!
1836!--        Statistical Evaluation of w'w'.
1837          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                         &
1838                             + ( flux_t(k) + diss_t(k) )                     &
1839                             *   weight_substep(intermediate_timestep_count)
1840
1841       ENDDO
1842
1843       DO  k = nzb_max+1, nzt
1844
1845          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1846          flux_r(k) = u_comp * (                                            &
1847                      37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
1848                    -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                    &
1849                    +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1850
1851          diss_r(k) = - ABS( u_comp ) * (                                   &
1852                      10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
1853                    -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                    &
1854                    +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1855
1856          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1857          flux_n(k) = v_comp * (                                            &
1858                      37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
1859                    -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                    &
1860                    +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1861
1862          diss_n(k) = - ABS( v_comp ) * (                                   &
1863                      10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
1864                    -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
1865                    +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1866!
1867!--       k index has to be modified near bottom and top, else array
1868!--       subscripts will be exceeded.
1869          ibit35 = IBITS(wall_flags_0(k,j,i),35,1)
1870          ibit34 = IBITS(wall_flags_0(k,j,i),34,1)
1871          ibit33 = IBITS(wall_flags_0(k,j,i),33,1)
1872
1873          k_ppp = k + 3 * ibit35
1874          k_pp  = k + 2 * ( 1 - ibit33  )
1875          k_mm  = k - 2 * ibit35
1876
1877          w_comp    = w(k+1,j,i) + w(k,j,i)
1878          flux_t(k) = w_comp  * (                                            &
1879                     ( 37.0 * ibit35 * adv_mom_5                             &
1880                  +     7.0 * ibit34 * adv_mom_3                             &
1881                  +           ibit33 * adv_mom_1                             &
1882                     ) *                                                     &
1883                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
1884              -      (  8.0 * ibit35 * adv_mom_5                             &
1885                  +           ibit34 * adv_mom_3                             &
1886                     ) *                                                     &
1887                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
1888              +      (        ibit35 * adv_mom_5                             &
1889                     ) *                                                     &
1890                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
1891                                 )
1892
1893          diss_t(k) = - ABS( w_comp ) * (                                    &
1894                     ( 10.0 * ibit35 * adv_mom_5                             &
1895                  +     3.0 * ibit34 * adv_mom_3                             &
1896                  +           ibit33 * adv_mom_1                             &
1897                     ) *                                                     &
1898                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
1899              -      (  5.0 * ibit35 * adv_mom_5                             &
1900                  +           ibit34 * adv_mom_3                             &
1901                     ) *                                                     &
1902                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
1903              +      (        ibit35 * adv_mom_5                             &
1904                     ) *                                                     &
1905                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
1906                                         )
1907!
1908!--       Calculate the divergence of the velocity field. A respective
1909!--       correction is needed to overcome numerical instabilities introduced
1910!--       by a not sufficient reduction of divergences near topography.
1911          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
1912              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
1913              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
1914                ) * 0.5
1915
1916          tend(k,j,i) = tend(k,j,i) - (                                       &
1917                      ( flux_r(k) + diss_r(k)                                 &
1918                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
1919                    + ( flux_n(k) + diss_n(k)                                 &
1920                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
1921                    + ( flux_t(k) + diss_t(k)                                 &
1922                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
1923                                      ) + div * w(k,j,i)
1924
1925          flux_l_w(k,j,tn) = flux_r(k)
1926          diss_l_w(k,j,tn) = diss_r(k)
1927          flux_s_w(k,tn)   = flux_n(k)
1928          diss_s_w(k,tn)   = diss_n(k)
1929          flux_d           = flux_t(k)
1930          diss_d           = diss_t(k)
1931!
1932!--        Statistical Evaluation of w'w'.
1933          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                         &
1934                             + ( flux_t(k) + diss_t(k) )                     &
1935                             *   weight_substep(intermediate_timestep_count)
1936
1937       ENDDO
1938
1939
1940    END SUBROUTINE advec_w_ws_ij
1941   
1942
1943!------------------------------------------------------------------------------!
1944! Scalar advection - Call for all grid points
1945!------------------------------------------------------------------------------!
1946    SUBROUTINE advec_s_ws( sk, sk_char )
1947
1948       USE arrays_3d
1949       USE constants
1950       USE control_parameters
1951       USE grid_variables
1952       USE indices
1953       USE statistics
1954
1955       IMPLICIT NONE
1956
1957       INTEGER ::  i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6,        &
1958                   ibit7, ibit8, j, k, k_mm, k_pp, k_ppp, tn = 0
1959#if defined( __nopointer )
1960       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk
1961#else
1962       REAL, DIMENSION(:,:,:), POINTER ::  sk
1963#endif
1964       REAL ::  diss_d, div, flux_d, u_comp, v_comp
1965       REAL, DIMENSION(nzb:nzt)   ::  diss_n, diss_r, diss_t, flux_n, flux_r,  &
1966                                      flux_t
1967       REAL, DIMENSION(nzb+1:nzt) ::  swap_diss_y_local, swap_flux_y_local
1968       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  swap_diss_x_local,               &
1969                                              swap_flux_x_local
1970       CHARACTER (LEN = *), INTENT(IN)    ::  sk_char
1971
1972!
1973!--    Compute the fluxes for the whole left boundary of the processor domain.
1974       i = nxl
1975       DO  j = nys, nyn
1976
1977          DO  k = nzb+1, nzb_max
1978
1979             ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
1980             ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
1981             ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
1982
1983             u_comp                 = u(k,j,i) - u_gtrans
1984             swap_flux_x_local(k,j) = u_comp * (                              &
1985                                                  ( 37.0 * ibit2 * adv_sca_5  &
1986                                               +     7.0 * ibit1 * adv_sca_3  &
1987                                               +           ibit0 * adv_sca_1  &
1988                                                  ) *                         &
1989                                            ( sk(k,j,i)   + sk(k,j,i-1)    )  &
1990                                           -      (  8.0 * ibit2 * adv_sca_5  &
1991                                               +           ibit1 * adv_sca_3  &
1992                                                  ) *                         &
1993                                            ( sk(k,j,i+1) + sk(k,j,i-2)    )  &
1994                                           +      (        ibit2 * adv_sca_5  &
1995                                                  ) *                         &
1996                                            ( sk(k,j,i+2) + sk(k,j,i-3)    )  &
1997                                               )
1998
1999              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
2000                                                  ( 10.0 * ibit2 * adv_sca_5  &
2001                                               +     3.0 * ibit1 * adv_sca_3  &
2002                                               +           ibit0 * adv_sca_1  &
2003                                                  ) *                         &
2004                                            ( sk(k,j,i)   - sk(k,j,i-1)    )  &
2005                                           -      (  5.0 * ibit2 * adv_sca_5  &
2006                                               +           ibit1 * adv_sca_3  &
2007                                                  ) *                         &
2008                                            ( sk(k,j,i+1) - sk(k,j,i-2)  )    &
2009                                           +      (        ibit2 * adv_sca_5  &
2010                                                  ) *                         &
2011                                            ( sk(k,j,i+2) - sk(k,j,i-3) )     &
2012                                                        )
2013
2014          ENDDO
2015
2016          DO  k = nzb_max+1, nzt
2017
2018             u_comp                 = u(k,j,i) - u_gtrans
2019             swap_flux_x_local(k,j) = u_comp * (                             &
2020                                      37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
2021                                    -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
2022                                    +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
2023                                               ) * adv_sca_5
2024
2025             swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
2026                                      10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
2027                                    -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
2028                                    +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
2029                                                       ) * adv_sca_5
2030
2031          ENDDO
2032
2033       ENDDO
2034
2035       DO  i = nxl, nxr
2036
2037          j = nys
2038          DO  k = nzb+1, nzb_max
2039
2040             ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
2041             ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
2042             ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
2043
2044             v_comp               = v(k,j,i) - v_gtrans
2045             swap_flux_y_local(k) = v_comp * (                                &
2046                                                  ( 37.0 * ibit5 * adv_sca_5  &
2047                                               +     7.0 * ibit4 * adv_sca_3  &
2048                                               +           ibit3 * adv_sca_1  &
2049                                                  ) *                         &
2050                                           ( sk(k,j,i)  + sk(k,j-1,i)     )   &
2051                                            -     (  8.0 * ibit5 * adv_sca_5  &
2052                                               +           ibit4 * adv_sca_3  &
2053                                                   ) *                        &
2054                                           ( sk(k,j+1,i) + sk(k,j-2,i)    )   &
2055                                           +      (        ibit5 * adv_sca_5  &
2056                                                  ) *                         &
2057                                           ( sk(k,j+2,i) + sk(k,j-3,i)    )   &
2058                                             )
2059
2060             swap_diss_y_local(k) = -ABS( v_comp ) * (                        &
2061                                                  ( 10.0 * ibit5 * adv_sca_5  &
2062                                               +     3.0 * ibit4 * adv_sca_3  &
2063                                               +           ibit3 * adv_sca_1  &
2064                                                  ) *                         &
2065                                            ( sk(k,j,i)   - sk(k,j-1,i)    )  &
2066                                           -      (  5.0 * ibit5 * adv_sca_5  &
2067                                               +           ibit4 * adv_sca_3  &
2068                                            ) *                               &
2069                                            ( sk(k,j+1,i) - sk(k,j-2,i)  )    &
2070                                           +      (        ibit5 * adv_sca_5  &
2071                                                  ) *                         &
2072                                            ( sk(k,j+2,i) - sk(k,j-3,i) )     &
2073                                                     )
2074
2075          ENDDO
2076!
2077!--       Above to the top of the highest topography. No degradation necessary.
2078          DO  k = nzb_max+1, nzt
2079
2080             v_comp               = v(k,j,i) - v_gtrans
2081             swap_flux_y_local(k) = v_comp * (                               &
2082                                    37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
2083                                  -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
2084                                  +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
2085                                             ) * adv_sca_5
2086              swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
2087                                    10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
2088                                  -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
2089                                  +          sk(k,j+2,i) - sk(k,j-3,i)       &
2090                                                      ) * adv_sca_5
2091
2092          ENDDO
2093
2094          DO  j = nys, nyn
2095
2096             flux_t(0) = 0.0
2097             diss_t(0) = 0.0
2098             flux_d    = 0.0
2099             diss_d    = 0.0
2100
2101             DO  k = nzb+1, nzb_max
2102
2103                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
2104                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
2105                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
2106
2107                u_comp    = u(k,j,i+1) - u_gtrans
2108                flux_r(k) = u_comp * (                                      &
2109                          ( 37.0 * ibit2 * adv_sca_5                        &
2110                      +      7.0 * ibit1 * adv_sca_3                        &
2111                      +           ibit0 * adv_sca_1                         &
2112                          ) *                                               &
2113                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
2114                   -      (  8.0 * ibit2 * adv_sca_5                        &
2115                       +           ibit1 * adv_sca_3                        &
2116                          ) *                                               &
2117                             ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
2118                   +      (        ibit2 * adv_sca_5                        &
2119                          ) *                                               &
2120                             ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
2121                                     )
2122
2123                diss_r(k) = -ABS( u_comp ) * (                              &
2124                          ( 10.0 * ibit2 * adv_sca_5                        &
2125                       +     3.0 * ibit1 * adv_sca_3                        &
2126                       +           ibit0 * adv_sca_1                        &
2127                          ) *                                               &
2128                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
2129                   -      (  5.0 * ibit2 * adv_sca_5                        &
2130                       +           ibit1 * adv_sca_3                        &
2131                          ) *                                               &
2132                             ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
2133                   +      (        ibit2 * adv_sca_5                        &
2134                          ) *                                               &
2135                             ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
2136                                             )
2137
2138                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
2139                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
2140                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
2141
2142                v_comp    = v(k,j+1,i) - v_gtrans
2143                flux_n(k) = v_comp * (                                      &
2144                          ( 37.0 * ibit5 * adv_sca_5                        &
2145                       +     7.0 * ibit4 * adv_sca_3                        &
2146                       +           ibit3 * adv_sca_1                        &
2147                          ) *                                               &
2148                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
2149                   -      (  8.0 * ibit5 * adv_sca_5                        &
2150                       +           ibit4 * adv_sca_3                        &
2151                          ) *                                               &
2152                             ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
2153                   +      (        ibit5 * adv_sca_5                        &
2154                          ) *                                               &
2155                             ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
2156                                     )
2157
2158                diss_n(k) = -ABS( v_comp ) * (                              &
2159                          ( 10.0 * ibit5 * adv_sca_5                        &
2160                       +     3.0 * ibit4 * adv_sca_3                        &
2161                       +           ibit3 * adv_sca_1                        &
2162                          ) *                                               &
2163                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
2164                   -      (  5.0 * ibit5 * adv_sca_5                        &
2165                       +           ibit4 * adv_sca_3                        &
2166                          ) *                                               &
2167                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
2168                   +      (        ibit5 * adv_sca_5                        &
2169                          ) *                                               &
2170                             ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
2171                                             )
2172!
2173!--             k index has to be modified near bottom and top, else array
2174!--             subscripts will be exceeded.
2175                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
2176                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
2177                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
2178
2179                k_ppp = k + 3 * ibit8
2180                k_pp  = k + 2 * ( 1 - ibit6  )
2181                k_mm  = k - 2 * ibit8
2182
2183
2184                flux_t(k) = w(k,j,i) * (                                      &
2185                           ( 37.0 * ibit8 * adv_sca_5                         &
2186                        +     7.0 * ibit7 * adv_sca_3                         &
2187                        +           ibit6 * adv_sca_1                         &
2188                           ) *                                                &
2189                                   ( sk(k+1,j,i)  + sk(k,j,i)   )             &
2190                          -      (  8.0 * ibit8 * adv_sca_5                   &
2191                        +           ibit7 * adv_sca_3                         &
2192                           ) *                                                &
2193                                   ( sk(k_pp,j,i) + sk(k-1,j,i) )             &
2194                    +      (        ibit8 * adv_sca_5                         &
2195                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
2196                                       )
2197
2198                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
2199                           ( 10.0 * ibit8 * adv_sca_5                         &
2200                        +     3.0 * ibit7 * adv_sca_3                         &
2201                        +           ibit6 * adv_sca_1                         &
2202                           ) *                                                &
2203                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
2204                    -      (  5.0 * ibit8 * adv_sca_5                         &
2205                        +           ibit7 * adv_sca_3                         &
2206                           ) *                                                &
2207                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
2208                    +      (        ibit8 * adv_sca_5                         &
2209                           ) *                                                &
2210                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
2211                                         )
2212!
2213!--             Calculate the divergence of the velocity field. A respective
2214!--             correction is needed to overcome numerical instabilities caused
2215!--             by a not sufficient reduction of divergences near topography.
2216                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
2217                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
2218                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
2219
2220                tend(k,j,i) = tend(k,j,i) - (                                 &
2221                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
2222                          swap_diss_x_local(k,j)            ) * ddx           &
2223                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
2224                          swap_diss_y_local(k)              ) * ddy           &
2225                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
2226                                                               ) * ddzw(k)    &
2227                                            ) + sk(k,j,i) * div
2228
2229                swap_flux_y_local(k)   = flux_n(k)
2230                swap_diss_y_local(k)   = diss_n(k)
2231                swap_flux_x_local(k,j) = flux_r(k)
2232                swap_diss_x_local(k,j) = diss_r(k)
2233                flux_d                 = flux_t(k)
2234                diss_d                 = diss_t(k)
2235
2236             ENDDO
2237
2238             DO  k = nzb_max+1, nzt
2239
2240                u_comp    = u(k,j,i+1) - u_gtrans
2241                flux_r(k) = u_comp * (                                      &
2242                      37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                  &
2243                    -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
2244                    +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
2245                diss_r(k) = -ABS( u_comp ) * (                              &
2246                      10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                  &
2247                    -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
2248                    +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
2249
2250                v_comp    = v(k,j+1,i) - v_gtrans
2251                flux_n(k) = v_comp * (                                      &
2252                      37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                  &
2253                    -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
2254                    +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
2255                diss_n(k) = -ABS( v_comp ) * (                              &
2256                      10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                  &
2257                    -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                  &
2258                    +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
2259!
2260!--             k index has to be modified near bottom and top, else array
2261!--             subscripts will be exceeded.
2262                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
2263                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
2264                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
2265
2266                k_ppp = k + 3 * ibit8
2267                k_pp  = k + 2 * ( 1 - ibit6  )
2268                k_mm  = k - 2 * ibit8
2269
2270
2271                flux_t(k) = w(k,j,i) * (                                      &
2272                           ( 37.0 * ibit8 * adv_sca_5                         &
2273                        +     7.0 * ibit7 * adv_sca_3                         &
2274                        +           ibit6 * adv_sca_1                         &
2275                           ) *                                                &
2276                                   ( sk(k+1,j,i)  + sk(k,j,i)   )             &
2277                          -      (  8.0 * ibit8 * adv_sca_5                   &
2278                        +           ibit7 * adv_sca_3                         &
2279                           ) *                                                &
2280                                   ( sk(k_pp,j,i) + sk(k-1,j,i) )             &
2281                    +      (        ibit8 * adv_sca_5                         &
2282                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
2283                                       )
2284
2285                diss_t(k) = -ABS( w(k,j,i) ) * (                              &
2286                           ( 10.0 * ibit8 * adv_sca_5                         &
2287                        +     3.0 * ibit7 * adv_sca_3                         &
2288                        +           ibit6 * adv_sca_1                         &
2289                           ) *                                                &
2290                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
2291                    -      (  5.0 * ibit8 * adv_sca_5                         &
2292                        +           ibit7 * adv_sca_3                         &
2293                           ) *                                                &
2294                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
2295                    +      (        ibit8 * adv_sca_5                         &
2296                           ) *                                                &
2297                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
2298                                         )
2299!
2300!--             Calculate the divergence of the velocity field. A respective
2301!--             correction is needed to overcome numerical instabilities introduced
2302!--             by a not sufficient reduction of divergences near topography.
2303                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
2304                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
2305                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
2306
2307                tend(k,j,i) = tend(k,j,i) - (                                 &
2308                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
2309                          swap_diss_x_local(k,j)            ) * ddx           &
2310                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
2311                          swap_diss_y_local(k)              ) * ddy           &
2312                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
2313                                                               ) * ddzw(k)    &
2314                                            ) + sk(k,j,i) * div
2315
2316                swap_flux_y_local(k)   = flux_n(k)
2317                swap_diss_y_local(k)   = diss_n(k)
2318                swap_flux_x_local(k,j) = flux_r(k)
2319                swap_diss_x_local(k,j) = diss_r(k)
2320                flux_d                 = flux_t(k)
2321                diss_d                 = diss_t(k)
2322
2323             ENDDO
2324!
2325!--          evaluation of statistics
2326             SELECT CASE ( sk_char )
2327
2328                 CASE ( 'pt' )
2329                    DO  k = nzb, nzt
2330                       sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
2331                        + ( flux_t(k) + diss_t(k) )                          &
2332                        *   weight_substep(intermediate_timestep_count)
2333                    ENDDO
2334                 CASE ( 'sa' )
2335                    DO  k = nzb, nzt
2336                       sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
2337                        + ( flux_t(k) + diss_t(k) )                          &
2338                        *   weight_substep(intermediate_timestep_count)
2339                    ENDDO
2340                 CASE ( 'q' )
2341                    DO  k = nzb, nzt
2342                       sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
2343                        + ( flux_t(k) + diss_t(k) )                          &
2344                        *   weight_substep(intermediate_timestep_count)
2345                    ENDDO
2346
2347              END SELECT
2348
2349         ENDDO
2350      ENDDO
2351
2352    END SUBROUTINE advec_s_ws
2353
2354
2355!------------------------------------------------------------------------------!
2356! Scalar advection - Call for all grid points - accelerator version
2357!------------------------------------------------------------------------------!
2358    SUBROUTINE advec_s_ws_acc ( sk, sk_char )
2359
2360       USE arrays_3d
2361       USE constants
2362       USE control_parameters
2363       USE grid_variables
2364       USE indices
2365       USE statistics
2366
2367       IMPLICIT NONE
2368
2369       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
2370
2371       INTEGER ::  i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6,        &
2372                   ibit7, ibit8, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
2373
2374       REAL    :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, flux_d, &
2375                  flux_l, flux_n, flux_r, flux_s, flux_t, u_comp, v_comp
2376
2377       REAL, INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  ::  sk
2378
2379!
2380!--    Computation of fluxes and tendency terms
2381       !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 )
2382       !$acc loop
2383       DO  i = nxl, nxr
2384          DO  j = nys, nyn
2385             !$acc loop vector( 32 )
2386             DO  k = nzb+1, nzt
2387
2388                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
2389                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
2390                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
2391
2392                u_comp              = u(k,j,i) - u_gtrans
2393                flux_l              = u_comp * (                           &
2394                                               ( 37.0 * ibit2 * adv_sca_5  &
2395                                            +     7.0 * ibit1 * adv_sca_3  &
2396                                            +           ibit0 * adv_sca_1  &
2397                                               ) *                         &
2398                                         ( sk(k,j,i)   + sk(k,j,i-1)    )  &
2399                                        -      (  8.0 * ibit2 * adv_sca_5  &
2400                                            +           ibit1 * adv_sca_3  &
2401                                               ) *                         &
2402                                         ( sk(k,j,i+1) + sk(k,j,i-2)    )  &
2403                                        +      (        ibit2 * adv_sca_5  &
2404                                               ) *                         &
2405                                         ( sk(k,j,i+2) + sk(k,j,i-3)    )  &
2406                                            )
2407
2408                diss_l              = -ABS( u_comp ) * (                   &
2409                                               ( 10.0 * ibit2 * adv_sca_5  &
2410                                            +     3.0 * ibit1 * adv_sca_3  &
2411                                            +           ibit0 * adv_sca_1  &
2412                                               ) *                         &
2413                                         ( sk(k,j,i)   - sk(k,j,i-1)    )  &
2414                                        -      (  5.0 * ibit2 * adv_sca_5  &
2415                                            +           ibit1 * adv_sca_3  &
2416                                               ) *                         &
2417                                         ( sk(k,j,i+1) - sk(k,j,i-2)  )    &
2418                                        +      (        ibit2 * adv_sca_5  &
2419                                               ) *                         &
2420                                         ( sk(k,j,i+2) - sk(k,j,i-3) )     &
2421                                                    )
2422
2423                u_comp    = u(k,j,i+1) - u_gtrans
2424                flux_r    = u_comp * (                                      &
2425                          ( 37.0 * ibit2 * adv_sca_5                        &
2426                      +      7.0 * ibit1 * adv_sca_3                        &
2427                      +            ibit0 * adv_sca_1                        &
2428                          ) *                                               &
2429                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
2430                   -      (  8.0 * ibit2 * adv_sca_5                        &
2431                       +           ibit1 * adv_sca_3                        &
2432                          ) *                                               &
2433                             ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
2434                   +      (        ibit2 * adv_sca_5                        &
2435                          ) *                                               &
2436                             ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
2437                                     )
2438
2439                diss_r    = -ABS( u_comp ) * (                              &
2440                          ( 10.0 * ibit2 * adv_sca_5                        &
2441                       +     3.0 * ibit1 * adv_sca_3                        &
2442                       +           ibit0 * adv_sca_1                        &
2443                          ) *                                               &
2444                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
2445                   -      (  5.0 * ibit2 * adv_sca_5                        &
2446                       +           ibit1 * adv_sca_3                        &
2447                          ) *                                               &
2448                             ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
2449                   +      (        ibit2 * adv_sca_5                        &
2450                          ) *                                               &
2451                             ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
2452                                             )
2453
2454                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
2455                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
2456                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
2457
2458                v_comp               = v(k,j,i) - v_gtrans
2459                flux_s               = v_comp * (                             &
2460                                                  ( 37.0 * ibit5 * adv_sca_5  &
2461                                               +     7.0 * ibit4 * adv_sca_3  &
2462                                               +           ibit3 * adv_sca_1  &
2463                                                  ) *                         &
2464                                           ( sk(k,j,i)  + sk(k,j-1,i)     )   &
2465                                            -     (  8.0 * ibit5 * adv_sca_5  &
2466                                               +           ibit4 * adv_sca_3  &
2467                                                   ) *                        &
2468                                           ( sk(k,j+1,i) + sk(k,j-2,i)    )   &
2469                                           +      (        ibit5 * adv_sca_5  &
2470                                                  ) *                         &
2471                                           ( sk(k,j+2,i) + sk(k,j-3,i)    )   &
2472                                             )
2473
2474                diss_s               = -ABS( v_comp ) * (                     &
2475                                                  ( 10.0 * ibit5 * adv_sca_5  &
2476                                               +     3.0 * ibit4 * adv_sca_3  &
2477                                               +           ibit3 * adv_sca_1  &
2478                                                  ) *                         &
2479                                            ( sk(k,j,i)   - sk(k,j-1,i)    )  &
2480                                           -      (  5.0 * ibit5 * adv_sca_5  &
2481                                               +           ibit4 * adv_sca_3  &
2482                                            ) *                               &
2483                                            ( sk(k,j+1,i) - sk(k,j-2,i)  )    &
2484                                           +      (        ibit5 * adv_sca_5  &
2485                                                  ) *                         &
2486                                            ( sk(k,j+2,i) - sk(k,j-3,i) )     &
2487                                                     )
2488
2489
2490                v_comp    = v(k,j+1,i) - v_gtrans
2491                flux_n    = v_comp * (                                      &
2492                          ( 37.0 * ibit5 * adv_sca_5                        &
2493                       +     7.0 * ibit4 * adv_sca_3                        &
2494                       +           ibit3 * adv_sca_1                        &
2495                          ) *                                               &
2496                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
2497                   -      (  8.0 * ibit5 * adv_sca_5                        &
2498                       +           ibit4 * adv_sca_3                        &
2499                          ) *                                               &
2500                             ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
2501                   +      (        ibit5 * adv_sca_5                        &
2502                          ) *                                               &
2503                             ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
2504                                     )
2505
2506                diss_n    = -ABS( v_comp ) * (                              &
2507                          ( 10.0 * ibit5 * adv_sca_5                        &
2508                       +     3.0 * ibit4 * adv_sca_3                        &
2509                       +           ibit3 * adv_sca_1                        &
2510                          ) *                                               &
2511                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
2512                   -      (  5.0 * ibit5 * adv_sca_5                        &
2513                       +           ibit4 * adv_sca_3                        &
2514                          ) *                                               &
2515                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
2516                   +      (        ibit5 * adv_sca_5                        &
2517                          ) *                                               &
2518                             ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
2519                                             )
2520
2521!
2522!--             indizes k_m, k_mm, ... should be known at these point
2523                ibit8 = IBITS(wall_flags_0(k-1,j,i),8,1)
2524                ibit7 = IBITS(wall_flags_0(k-1,j,i),7,1)
2525                ibit6 = IBITS(wall_flags_0(k-1,j,i),6,1)
2526
2527                k_pp  = k + 2 * ibit8
2528                k_mm  = k - 2 * ( ibit7 + ibit8 )
2529                k_mmm = k - 3 * ibit8
2530
2531                flux_d    = w(k-1,j,i) * (                                    &
2532                           ( 37.0 * ibit8 * adv_sca_5                         &
2533                        +     7.0 * ibit7 * adv_sca_3                         &
2534                        +           ibit6 * adv_sca_1                         &
2535                           ) *                                                &
2536                                   ( sk(k,j,i)    + sk(k-1,j,i) )             &
2537                          -      (  8.0 * ibit8 * adv_sca_5                   &
2538                          +               ibit7 * adv_sca_3                   &
2539                           ) *                                                &
2540                                   ( sk(k+1,j,i) + sk(k_mm,j,i) )             &
2541                    +      (        ibit8 * adv_sca_5                         &
2542                           ) *     ( sk(k_pp,j,i)+ sk(k_mmm,j,i) )            &
2543                                       )
2544
2545                diss_d    = -ABS( w(k-1,j,i) ) * (                            &
2546                           ( 10.0 * ibit8 * adv_sca_5                         &
2547                        +     3.0 * ibit7 * adv_sca_3                         &
2548                        +           ibit6 * adv_sca_1                         &
2549                           ) *                                                &
2550                                   ( sk(k,j,i)    - sk(k-1,j,i)   )           &
2551                    -      (  5.0 * ibit8 * adv_sca_5                         &
2552                        +           ibit7 * adv_sca_3                         &
2553                           ) *                                                &
2554                                   ( sk(k+1,j,i)  - sk(k_mm,j,i)  )           &
2555                    +      (        ibit8 * adv_sca_5                         &
2556                           ) *                                                &
2557                                   ( sk(k_pp,j,i) - sk(k_mmm,j,i) )           &
2558                                         )
2559
2560                ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
2561                ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
2562                ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
2563
2564                k_ppp = k + 3 * ibit8
2565                k_pp  = k + 2 * ( 1 - ibit6  )
2566                k_mm  = k - 2 * ibit8
2567
2568                flux_t    = w(k,j,i) * (                                      &
2569                           ( 37.0 * ibit8 * adv_sca_5                         &
2570                        +     7.0 * ibit7 * adv_sca_3                         &
2571                        +           ibit6 * adv_sca_1                         &
2572                           ) *                                                &
2573                                   ( sk(k+1,j,i)  + sk(k,j,i)   )             &
2574                          -      (  8.0 * ibit8 * adv_sca_5                   &
2575                        +                 ibit7 * adv_sca_3                   &
2576                           ) *                                                &
2577                                   ( sk(k_pp,j,i) + sk(k-1,j,i) )             &
2578                    +      (        ibit8 * adv_sca_5                         &
2579                           ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )            &
2580                                       )
2581
2582                diss_t    = -ABS( w(k,j,i) ) * (                              &
2583                           ( 10.0 * ibit8 * adv_sca_5                         &
2584                        +     3.0 * ibit7 * adv_sca_3                         &
2585                        +           ibit6 * adv_sca_1                         &
2586                           ) *                                                &
2587                                   ( sk(k+1,j,i)   - sk(k,j,i)    )           &
2588                    -      (  5.0 * ibit8 * adv_sca_5                         &
2589                        +           ibit7 * adv_sca_3                         &
2590                           ) *                                                &
2591                                   ( sk(k_pp,j,i)  - sk(k-1,j,i)  )           &
2592                    +      (        ibit8 * adv_sca_5                         &
2593                           ) *                                                &
2594                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
2595                                         )
2596!
2597!--             Calculate the divergence of the velocity field. A respective
2598!--             correction is needed to overcome numerical instabilities caused
2599!--             by a not sufficient reduction of divergences near topography.
2600                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
2601                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
2602                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
2603
2604                tend(k,j,i) = - (                                             &
2605                               ( flux_r + diss_r - flux_l - diss_l ) * ddx    &
2606                             + ( flux_n + diss_n - flux_s - diss_s ) * ddy    &
2607                             + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)&
2608                                ) + div * sk(k,j,i)
2609
2610!++
2611!--             Evaluation of statistics
2612!                SELECT CASE ( sk_char )
2613!
2614!                   CASE ( 'pt' )
2615!                      sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
2616!                       + ( flux_t + diss_t )                                &
2617!                       *   weight_substep(intermediate_timestep_count)
2618!                   CASE ( 'sa' )
2619!                      sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
2620!                       + ( flux_t + diss_t )                                &
2621!                       *   weight_substep(intermediate_timestep_count)
2622!                   CASE ( 'q' )
2623!                      sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
2624!                      + ( flux_t + diss_t )                                &
2625!                      *   weight_substep(intermediate_timestep_count)
2626!
2627!                END SELECT
2628
2629             ENDDO
2630         ENDDO
2631      ENDDO
2632      !$acc end kernels
2633
2634    END SUBROUTINE advec_s_ws_acc
2635
2636
2637!------------------------------------------------------------------------------!
2638! Advection of u - Call for all grid points
2639!------------------------------------------------------------------------------!
2640    SUBROUTINE advec_u_ws
2641
2642       USE arrays_3d
2643       USE constants
2644       USE control_parameters
2645       USE grid_variables
2646       USE indices
2647       USE statistics
2648
2649       IMPLICIT NONE
2650
2651       INTEGER ::  i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15,   &
2652                   ibit16, ibit17, j, k, k_mm, k_pp, k_ppp, tn = 0
2653       REAL    ::  diss_d, div, flux_d, gu, gv, v_comp, w_comp
2654       REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u, swap_flux_y_local_u
2655       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u,              &
2656                                             swap_flux_x_local_u
2657       REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r,     &
2658                                   flux_t, u_comp
2659 
2660       gu = 2.0 * u_gtrans
2661       gv = 2.0 * v_gtrans
2662
2663!
2664!--    Compute the fluxes for the whole left boundary of the processor domain.
2665       i = nxlu
2666       DO  j = nys, nyn
2667          DO  k = nzb+1, nzb_max
2668
2669             ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
2670             ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
2671             ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
2672
2673             u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
2674             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
2675                                       ( 37.0 * ibit11 * adv_mom_5             &
2676                                    +     7.0 * ibit10 * adv_mom_3             &
2677                                    +           ibit9  * adv_mom_1             &
2678                                       ) *                                     &
2679                                     ( u(k,j,i)   + u(k,j,i-1) )               &
2680                                -      (  8.0 * ibit11 * adv_mom_5             &
2681                                    +           ibit10 * adv_mom_3             &
2682                                       ) *                                     &
2683                                     ( u(k,j,i+1) + u(k,j,i-2) )               &
2684                                +      (        ibit11 * adv_mom_5             &
2685                                       ) *                                     &
2686                                     ( u(k,j,i+2) + u(k,j,i-3) )               &
2687                                                   )
2688
2689              swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * (                &
2690                                       ( 10.0 * ibit11 * adv_mom_5             &
2691                                    +     3.0 * ibit10 * adv_mom_3             &
2692                                    +           ibit9  * adv_mom_1             &
2693                                       ) *                                     &
2694                                     ( u(k,j,i)   - u(k,j,i-1) )               &
2695                                -      (  5.0 * ibit11 * adv_mom_5             &
2696                                    +           ibit10 * adv_mom_3             &
2697                                       ) *                                     &
2698                                     ( u(k,j,i+1) - u(k,j,i-2) )               &
2699                                +      (        ibit11 * adv_mom_5             &
2700                                       ) *                                     &
2701                                     ( u(k,j,i+2) - u(k,j,i-3) )               &
2702                                                             )
2703
2704          ENDDO
2705
2706          DO  k = nzb_max+1, nzt
2707
2708             u_comp(k)         = u(k,j,i) + u(k,j,i-1) - gu
2709             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
2710                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )                &
2711                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )                &
2712                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
2713             swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                   &
2714                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )                &
2715                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )                &
2716                           +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
2717
2718          ENDDO
2719       ENDDO
2720
2721       DO i = nxlu, nxr
2722!       
2723!--       The following loop computes the fluxes for the south boundary points
2724          j = nys
2725          DO  k = nzb+1, nzb_max
2726
2727             ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
2728             ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
2729             ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
2730
2731             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2732             swap_flux_y_local_u(k) = v_comp * (                              &
2733                                   ( 37.0 * ibit14 * adv_mom_5                &
2734                                +     7.0 * ibit13 * adv_mom_3                &
2735                                +           ibit12 * adv_mom_1                &
2736                                   ) *                                        &
2737                                     ( u(k,j,i)   + u(k,j-1,i) )              &
2738                            -      (  8.0 * ibit14 * adv_mom_5                &
2739                            +           ibit13 * adv_mom_3                    &
2740                                   ) *                                        &
2741                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
2742                        +      (        ibit14 * adv_mom_5                    &
2743                               ) *                                            &
2744                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
2745                                               )
2746
2747             swap_diss_y_local_u(k) = - ABS ( v_comp ) * (                    &
2748                                   ( 10.0 * ibit14 * adv_mom_5                &
2749                                +     3.0 * ibit13 * adv_mom_3                &
2750                                +           ibit12 * adv_mom_1                &
2751                                   ) *                                        &
2752                                     ( u(k,j,i)   - u(k,j-1,i) )              &
2753                            -      (  5.0 * ibit14 * adv_mom_5                &
2754                                +           ibit13 * adv_mom_3                &
2755                                   ) *                                        &
2756                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
2757                            +      (        ibit14 * adv_mom_5                &
2758                                   ) *                                        &
2759                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
2760                                                         )
2761
2762          ENDDO
2763
2764          DO  k = nzb_max+1, nzt
2765
2766             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2767             swap_flux_y_local_u(k) = v_comp * (                              &
2768                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
2769                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
2770                         +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
2771             swap_diss_y_local_u(k) = - ABS(v_comp) * (                       &
2772                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
2773                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
2774                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
2775
2776          ENDDO
2777!
2778!--       Computation of interior fluxes and tendency terms
2779          DO  j = nys, nyn
2780
2781             flux_t(0) = 0.0
2782             diss_t(0) = 0.0
2783             flux_d    = 0.0
2784             diss_d    = 0.0
2785
2786             DO  k = nzb+1, nzb_max
2787
2788                ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
2789                ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
2790                ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
2791
2792                u_comp(k) = u(k,j,i+1) + u(k,j,i)
2793                flux_r(k) = ( u_comp(k) - gu ) * (                           &
2794                          ( 37.0 * ibit11 * adv_mom_5                        &
2795                       +     7.0 * ibit10 * adv_mom_3                        &
2796                       +           ibit9  * adv_mom_1                        &
2797                          ) *                                                &
2798                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
2799                   -      (  8.0 * ibit11 * adv_mom_5                        &
2800                       +           ibit10 * adv_mom_3                        &
2801                          ) *                                                &
2802                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
2803                   +      (        ibit11 * adv_mom_5                        &
2804                          ) *                                                &
2805                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
2806                                                 )
2807
2808                diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
2809                          ( 10.0 * ibit11 * adv_mom_5                        &
2810                       +     3.0 * ibit10 * adv_mom_3                        & 
2811                       +           ibit9  * adv_mom_1                        &
2812                          ) *                                                &
2813                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
2814                   -      (  5.0 * ibit11 * adv_mom_5                        &
2815                       +           ibit10 * adv_mom_3                        &
2816                          ) *                                                &
2817                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
2818                   +      (        ibit11 * adv_mom_5                        &
2819                          ) *                                                &
2820                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
2821                                                     )
2822
2823                ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
2824                ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
2825                ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
2826
2827                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2828                flux_n(k) = v_comp * (                                       &
2829                          ( 37.0 * ibit14 * adv_mom_5                        &
2830                       +     7.0 * ibit13 * adv_mom_3                        &
2831                       +           ibit12 * adv_mom_1                        &
2832                          ) *                                                &
2833                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
2834                   -      (  8.0 * ibit14 * adv_mom_5                        &
2835                       +           ibit13 * adv_mom_3                        &
2836                          ) *                                                &
2837                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
2838                   +      (        ibit14 * adv_mom_5                        & 
2839                          ) *                                                &
2840                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
2841                                                 )
2842
2843                diss_n(k) = - ABS ( v_comp ) * (                             &
2844                          ( 10.0 * ibit14 * adv_mom_5                        &
2845                       +     3.0 * ibit13 * adv_mom_3                        &
2846                       +           ibit12 * adv_mom_1                        &
2847                          ) *                                                &
2848                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
2849                   -      (  5.0 * ibit14 * adv_mom_5                        &
2850                       +           ibit13 * adv_mom_3                        &
2851                          ) *                                                &
2852                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
2853                   +      (        ibit14 * adv_mom_5                        &
2854                          ) *                                                &
2855                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
2856                                                      )
2857!
2858!--             k index has to be modified near bottom and top, else array
2859!--             subscripts will be exceeded.
2860                ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
2861                ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
2862                ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
2863
2864                k_ppp = k + 3 * ibit17
2865                k_pp  = k + 2 * ( 1 - ibit15  )
2866                k_mm  = k - 2 * ibit17
2867
2868                w_comp    = w(k,j,i) + w(k,j,i-1)
2869                flux_t(k) = w_comp  * (                                      &
2870                          ( 37.0 * ibit17 * adv_mom_5                        &
2871                       +     7.0 * ibit16 * adv_mom_3                        &
2872                       +           ibit15 * adv_mom_1                        & 
2873                          ) *                                                &
2874                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
2875                   -      (  8.0 * ibit17 * adv_mom_5                        &
2876                       +           ibit16 * adv_mom_3                        &
2877                          ) *                                                &
2878                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
2879                   +      (        ibit17 * adv_mom_5                        &
2880                          ) *                                                &
2881                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
2882                                      )
2883
2884                diss_t(k) = - ABS( w_comp ) * (                              &
2885                          ( 10.0 * ibit17 * adv_mom_5                        &
2886                       +     3.0 * ibit16 * adv_mom_3                        &
2887                       +           ibit15 * adv_mom_1                        &
2888                          ) *                                                &
2889                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
2890                   -      (  5.0 * ibit17 * adv_mom_5                        &
2891                       +           ibit16 * adv_mom_3                        &
2892                          ) *                                                &
2893                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
2894                   +      (        ibit17 * adv_mom_5                        &
2895                           ) *                                               &
2896                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
2897                                              )
2898!
2899!--             Calculate the divergence of the velocity field. A respective
2900!--             correction is needed to overcome numerical instabilities caused
2901!--             by a not sufficient reduction of divergences near topography.
2902                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
2903                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
2904                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
2905                                                                    * ddzw(k)  &
2906                      ) * 0.5
2907
2908                tend(k,j,i) = tend(k,j,i) - (                                  &
2909                 ( flux_r(k) + diss_r(k)                                       &
2910               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
2911               + ( flux_n(k) + diss_n(k)                                       &
2912               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
2913               + ( flux_t(k) + diss_t(k)                                       &
2914               -   flux_d    - diss_d                                          &
2915                                                                  ) * ddzw(k)  &
2916                                           ) + div * u(k,j,i)
2917
2918                swap_flux_x_local_u(k,j) = flux_r(k)
2919                swap_diss_x_local_u(k,j) = diss_r(k)
2920                swap_flux_y_local_u(k)   = flux_n(k)
2921                swap_diss_y_local_u(k)   = diss_n(k)
2922                flux_d                   = flux_t(k)
2923                diss_d                   = diss_t(k)
2924!
2925!--             Statistical Evaluation of u'u'. The factor has to be applied
2926!--             for right evaluation when gallilei_trans = .T. .
2927                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                     &
2928                              + ( flux_r(k) *                                 &
2929                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
2930                              / ( u_comp(k) - gu + 1.0E-20      )             &
2931                              +   diss_r(k) *                                 &
2932                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
2933                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
2934                              *   weight_substep(intermediate_timestep_count)
2935!
2936!--             Statistical Evaluation of w'u'.
2937                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                   &
2938                              + ( flux_t(k) + diss_t(k) )                     &
2939                              *   weight_substep(intermediate_timestep_count)
2940             ENDDO
2941
2942             DO  k = nzb_max+1, nzt
2943
2944                u_comp(k) = u(k,j,i+1) + u(k,j,i)
2945                flux_r(k) = ( u_comp(k) - gu ) * (                            &
2946                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
2947                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
2948                       +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2949                diss_r(k) = - ABS( u_comp(k) - gu ) * (                       &
2950                         10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
2951                       -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
2952                       +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2953
2954                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2955                flux_n(k) = v_comp * (                                        &
2956                         37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
2957                       -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                   &
2958                       +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2959                diss_n(k) = - ABS( v_comp ) * (                               &
2960                         10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
2961                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
2962                       +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2963!
2964!--             k index has to be modified near bottom and top, else array
2965!--             subscripts will be exceeded.
2966                ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
2967                ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
2968                ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
2969
2970                k_ppp = k + 3 * ibit17
2971                k_pp  = k + 2 * ( 1 - ibit15  )
2972                k_mm  = k - 2 * ibit17
2973
2974                w_comp    = w(k,j,i) + w(k,j,i-1)
2975                flux_t(k) = w_comp  * (                                      &
2976                          ( 37.0 * ibit17 * adv_mom_5                        &
2977                       +     7.0 * ibit16 * adv_mom_3                        &
2978                       +           ibit15 * adv_mom_1                        &
2979                          ) *                                                &
2980                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
2981                   -      (  8.0 * ibit17 * adv_mom_5                        &
2982                       +           ibit16 * adv_mom_3                        &
2983                          ) *                                                &
2984                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
2985                   +      (        ibit17 * adv_mom_5                        &
2986                          ) *                                                &
2987                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
2988                                      )
2989
2990                diss_t(k) = - ABS( w_comp ) * (                              &
2991                          ( 10.0 * ibit17 * adv_mom_5                        &
2992                       +     3.0 * ibit16 * adv_mom_3                        &
2993                       +           ibit15 * adv_mom_1                        &
2994                          ) *                                                &
2995                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
2996                   -      (  5.0 * ibit17 * adv_mom_5                        &
2997                       +           ibit16 * adv_mom_3                        &
2998                          ) *                                                &
2999                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
3000                   +      (        ibit17 * adv_mom_5                        &
3001                           ) *                                               &
3002                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
3003                                              )
3004!
3005!--             Calculate the divergence of the velocity field. A respective
3006!--             correction is needed to overcome numerical instabilities caused
3007!--             by a not sufficient reduction of divergences near topography.
3008                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
3009                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy &
3010                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )       &
3011                                                                    * ddzw(k) &
3012                      ) * 0.5
3013
3014                 tend(k,j,i) = tend(k,j,i) - (                                 &
3015                 ( flux_r(k) + diss_r(k)                                       &
3016               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
3017               + ( flux_n(k) + diss_n(k)                                       &
3018               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
3019               + ( flux_t(k) + diss_t(k)                                       &
3020               -   flux_d    - diss_d                                          &
3021                                                                  ) * ddzw(k)  &
3022                                           ) + div * u(k,j,i)
3023
3024                 swap_flux_x_local_u(k,j) = flux_r(k)
3025                 swap_diss_x_local_u(k,j) = diss_r(k)
3026                 swap_flux_y_local_u(k)   = flux_n(k)
3027                 swap_diss_y_local_u(k)   = diss_n(k)
3028                 flux_d                   = flux_t(k)
3029                 diss_d                   = diss_t(k)
3030!
3031!--              Statistical Evaluation of u'u'. The factor has to be applied
3032!--              for right evaluation when gallilei_trans = .T. .
3033                 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                    &
3034                              + ( flux_r(k) *                                 &
3035                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
3036                              / ( u_comp(k) - gu + 1.0E-20      )             &
3037                              +   diss_r(k) *                                 &
3038                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
3039                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
3040                              *   weight_substep(intermediate_timestep_count)
3041!
3042!--              Statistical Evaluation of w'u'.
3043                 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                  &
3044                              + ( flux_t(k) + diss_t(k) )                     &
3045                              *   weight_substep(intermediate_timestep_count)
3046       ENDDO