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

Last change on this file since 1010 was 1010, checked in by raasch, 9 years ago

pointer free version can be generated with cpp switch nopointer

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