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

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

last commit documented

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