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

Last change on this file since 861 was 861, checked in by suehring, 10 years ago

WS5 is available in combination with topography. Version number changed from 3.8 to 3.8a. Modification in init_pt_anomaly.

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