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

Last change on this file since 1310 was 1310, checked in by raasch, 7 years ago

update of GPL copyright

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