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

Last change on this file since 1222 was 1222, checked in by raasch, 8 years ago

last commit documented

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