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

Last change on this file since 888 was 888, checked in by suehring, 9 years ago

IBITS() calls with identical arguments are replaced by a temporary variable

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