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

Last change on this file since 1028 was 1028, checked in by suehring, 9 years ago

last commit documented

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