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

Last change on this file since 1053 was 1053, checked in by hoffmann, 9 years ago

two-moment cloud physics implemented

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