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

Last change on this file since 862 was 862, checked in by suehring, 11 years ago

Formatting adjustments.

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