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

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

last commit documented

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