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

Last change on this file since 1221 was 1221, checked in by raasch, 11 years ago

New:


openACC porting of reduction operations
additional 3D-flag arrays for replacing the 2D-index arrays nzb_s_inner and nzb_diff_s_inner
(flow_statistics, init_grid, init_3d_model, modules, palm, pres, time_integration)

Changed:


for PGI/openACC performance reasons set default compile options for openACC to "-ta=nocache",
and set environment variable PGI_ACC_SYNCHRONOUS=1
(MAKE.inc.pgi.openacc, palm_simple_run)

wall_flags_0 changed to 32bit INTEGER, additional array wall_flags_00 introduced to hold
bits 32-63
(advec_ws, init_grid, modules, palm)

Errors:


dummy argument tri in 1d-routines replaced by tri_for_1d because of name
conflict with arry tri in module arrays_3d
(tridia_solver)

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