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

Last change on this file since 1116 was 1116, checked in by hoffmann, 12 years ago

last commit documented

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