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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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