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

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

last commit documented

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 191.4 KB
Line 
1 MODULE advec_ws
2
3!-----------------------------------------------------------------------------!
4! Current revisions:
5! ------------------
6!
7! Former revisions:
8! -----------------
9! $Id: advec_ws.f90 863 2012-03-26 14:28:17Z suehring $
10
11! 862 2012-03-26 14:21:38Z suehring
12! ws-scheme also work with topography in combination with vector version.
13! ws-scheme also work with outflow boundaries in combination with
14! vector version.
15! Degradation of the applied order of scheme is now steered by multiplying with
16! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each
17! grid point and mulitplied with the appropriate flag.
18! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
19! term derived according to the 4th and 6th order terms is applied. It turns
20! out that diss_2nd does not provide sufficient dissipation near walls.
21! Therefore, the function diss_2nd is removed.
22! Near walls a divergence correction is necessary to overcome numerical
23! instabilities due to too less divergence reduction of the velocity field.
24! boundary_flags and logicals steering the degradation are removed.
25! Empty SUBROUTINE local_diss removed.
26! Further formatting adjustments.
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, k_mm, k_pp, k_ppp,  tn
277       REAL    ::  diss_d, div, flux_d, u_comp, v_comp
278       REAL, DIMENSION(:,:,:), POINTER    :: sk
279       REAL, DIMENSION(nzb:nzt+1)         :: diss_n, diss_r, diss_t, flux_n,  &
280                                             flux_r, flux_t
281       REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local,  &
282                                                          swap_flux_y_local
283       REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::             &
284                                                          swap_diss_x_local,  &
285                                                          swap_flux_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, k_mm, k_pp, k_ppp, tn
674       REAL    ::  diss_d, div, flux_d, gu, gv, u_comp_l, v_comp, w_comp
675       REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r,  &
676                                     flux_t, 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, k_mm, k_pp, k_ppp, tn
1061       REAL     ::  diss_d, div, flux_d, gu, gv, u_comp, v_comp_l, w_comp
1062       REAL, DIMENSION(nzb:nzt+1)  :: diss_n, diss_r, diss_t, flux_n, flux_r,  &
1063                                      flux_t, 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, k_mm, k_pp, k_ppp, tn
1454       REAL    ::  diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp
1455       REAL, DIMENSION(nzb:nzt+1)  :: diss_n, diss_r, diss_t, flux_n, flux_r, &
1456                                      flux_t
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!--    The lower flux has to be calculated explicetely for the tendency at
1571!--    the first w-level. For topography wall this is done implicitely by
1572!--    wall_flags_0.
1573       k         = nzb + 1
1574       w_comp    = w(k,j,i) + w(k-1,j,i)
1575       flux_t(0) = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
1576       diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
1577       flux_d    = flux_t(0)
1578       diss_d    = diss_t(0)
1579!
1580!--    Now compute the fluxes and tendency terms for the horizontal
1581!--    and vertical parts.
1582       DO  k = nzb+1, nzb_max
1583
1584          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1585          flux_r(k) = u_comp * (                                             &
1586                     ( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
1587                  +     7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
1588                  +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1    &
1589                     ) *                                                     &
1590                                 ( w(k,j,i+1) + w(k,j,i)   )                 &
1591              -      (  8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
1592                  +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
1593                     ) *                                                     &
1594                                 ( w(k,j,i+2) + w(k,j,i-1) )                 &
1595              +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
1596                     ) *                                                     &
1597                                 ( w(k,j,i+3) + w(k,j,i-2) )                 &
1598                               )
1599
1600          diss_r(k) = - ABS( u_comp ) * (                                    &
1601                     ( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
1602                  +     3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
1603                  +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1    &
1604                     ) *                                                     &
1605                                 ( w(k,j,i+1) - w(k,j,i)  )                  &
1606              -      (  5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
1607                  +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
1608                     ) *                                                     &
1609                                 ( w(k,j,i+2) - w(k,j,i-1) )                 &
1610              +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
1611                     ) *                                                     &
1612                                 ( w(k,j,i+3) - w(k,j,i-2) )                 &
1613                                        )
1614
1615          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1616          flux_n(k) = v_comp * (                                             &
1617                     ( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
1618                  +     7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
1619                  +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1    &
1620                     ) *                                                     &
1621                                 ( w(k,j+1,i) + w(k,j,i)   )                 &
1622              -      (  8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
1623                  +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
1624                     ) *                                                     &
1625                                 ( w(k,j+2,i) + w(k,j-1,i) )                 &
1626              +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
1627                     ) *                                                     &
1628                                 ( w(k,j+3,i) + w(k,j-2,i) )                 &
1629                               )
1630
1631          diss_n(k) = - ABS( v_comp ) * (                                    &
1632                     ( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
1633                  +     3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
1634                  +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1    &
1635                     ) *                                                     &
1636                                 ( w(k,j+1,i) - w(k,j,i)  )                  &
1637              -      (  5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
1638                  +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
1639                     ) *                                                     &
1640                                 ( w(k,j+2,i) - w(k,j-1,i) )                 &
1641              +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
1642                     ) *                                                     &
1643                                 ( w(k,j+3,i) - w(k,j-2,i) )                 &
1644                                        )
1645!
1646!--       k index has to be modified near bottom and top, else array
1647!--       subscripts will be exceeded.
1648          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1)
1649          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1)  )
1650          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),35,1)
1651
1652          w_comp    = w(k+1,j,i) + w(k,j,i)
1653          flux_t(k) = w_comp  * (                                            &
1654                     ( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1655                  +     7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1656                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
1657                     ) *                                                     &
1658                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
1659              -      (  8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1660                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1661                     ) *                                                     &
1662                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
1663              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1664                     ) *                                                     &
1665                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
1666                                 )
1667
1668          diss_t(k) = - ABS( w_comp ) * (                                    &
1669                     ( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1670                  +     3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1671                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
1672                     ) *                                                     &
1673                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
1674              -      (  5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1675                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1676                     ) *                                                     &
1677                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
1678              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1679                     ) *                                                     &
1680                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
1681                                         )
1682
1683!
1684!--       Calculate the divergence of the velocity field. A respective
1685!--       correction is needed to overcome numerical instabilities introduced
1686!--       by a not sufficient reduction of divergences near topography.
1687          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
1688              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
1689              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
1690                ) * 0.5
1691
1692          tend(k,j,i) = tend(k,j,i) - (                                       &
1693                      ( flux_r(k) + diss_r(k)                                 &
1694                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
1695                    + ( flux_n(k) + diss_n(k)                                 &
1696                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
1697                    + ( flux_t(k) + diss_t(k)                                 &
1698                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
1699                                      ) + div * w(k,j,i)
1700
1701          flux_l_w(k,j,tn) = flux_r(k)
1702          diss_l_w(k,j,tn) = diss_r(k)
1703          flux_s_w(k,tn)   = flux_n(k)
1704          diss_s_w(k,tn)   = diss_n(k)
1705          flux_d           = flux_t(k)
1706          diss_d           = diss_t(k)
1707!
1708!--        Statistical Evaluation of w'w'.
1709          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                         &
1710                             + ( flux_t(k) + diss_t(k) )                     &
1711                             *   weight_substep(intermediate_timestep_count)
1712
1713       ENDDO
1714
1715       DO  k = nzb_max+1, nzt
1716
1717          u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1718          flux_r(k) = u_comp * (                                            &
1719                      37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
1720                    -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                    &
1721                    +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1722
1723          diss_r(k) = - ABS( u_comp ) * (                                   &
1724                      10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
1725                    -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                    &
1726                    +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1727
1728          v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1729          flux_n(k) = v_comp * (                                            &
1730                      37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
1731                    -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                    &
1732                    +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1733
1734          diss_n(k) = - ABS( v_comp ) * (                                   &
1735                      10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
1736                    -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
1737                    +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1738!
1739!--       k index has to be modified near bottom and top, else array
1740!--       subscripts will be exceeded.
1741          k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1)
1742          k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1)  )
1743          k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),35,1)
1744
1745          w_comp    = w(k+1,j,i) + w(k,j,i)
1746          flux_t(k) = w_comp  * (                                            &
1747                     ( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1748                  +     7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1749                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
1750                     ) *                                                     &
1751                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
1752              -      (  8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1753                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1754                     ) *                                                     &
1755                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
1756              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1757                     ) *                                                     &
1758                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
1759                                 )
1760
1761          diss_t(k) = - ABS( w_comp ) * (                                    &
1762                     ( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1763                  +     3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1764                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
1765                     ) *                                                     &
1766                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
1767              -      (  5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1768                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
1769                     ) *                                                     &
1770                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
1771              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
1772                     ) *                                                     &
1773                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
1774                                         )
1775!
1776!--       Calculate the divergence of the velocity field. A respective
1777!--       correction is needed to overcome numerical instabilities introduced
1778!--       by a not sufficient reduction of divergences near topography.
1779          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
1780              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
1781              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
1782                ) * 0.5
1783
1784          tend(k,j,i) = tend(k,j,i) - (                                       &
1785                      ( flux_r(k) + diss_r(k)                                 &
1786                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx         &
1787                    + ( flux_n(k) + diss_n(k)                                 &
1788                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy         &
1789                    + ( flux_t(k) + diss_t(k)                                 &
1790                    -   flux_d    - diss_d                    ) * ddzu(k+1)   &
1791                                      ) + div * w(k,j,i)
1792
1793          flux_l_w(k,j,tn) = flux_r(k)
1794          diss_l_w(k,j,tn) = diss_r(k)
1795          flux_s_w(k,tn)   = flux_n(k)
1796          diss_s_w(k,tn)   = diss_n(k)
1797          flux_d           = flux_t(k)
1798          diss_d           = diss_t(k)
1799!
1800!--        Statistical Evaluation of w'w'.
1801          sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                         &
1802                             + ( flux_t(k) + diss_t(k) )                     &
1803                             *   weight_substep(intermediate_timestep_count)
1804
1805       ENDDO
1806
1807
1808    END SUBROUTINE advec_w_ws_ij
1809   
1810
1811!------------------------------------------------------------------------------!
1812! Scalar advection - Call for all grid points
1813!------------------------------------------------------------------------------!
1814    SUBROUTINE advec_s_ws( sk, sk_char )
1815
1816       USE arrays_3d
1817       USE constants
1818       USE control_parameters
1819       USE grid_variables
1820       USE indices
1821       USE statistics
1822
1823       IMPLICIT NONE
1824
1825       INTEGER ::  i, j, k, k_mm, k_pp, k_ppp, tn = 0
1826       REAL, DIMENSION(:,:,:), POINTER ::  sk
1827       REAL    :: diss_d, div, flux_d, u_comp, v_comp
1828       REAL, DIMENSION(nzb:nzt)   :: diss_n, diss_r, diss_t, flux_n, flux_r,  &
1829                                     flux_t
1830       REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local, swap_flux_y_local
1831       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local,               &
1832                                             swap_flux_x_local
1833       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
1834
1835!
1836!--    Compute the fluxes for the whole left boundary of the processor domain.
1837       i = nxl
1838       DO  j = nys, nyn
1839
1840          DO  k = nzb+1, nzb_max
1841
1842             u_comp                 = u(k,j,i) - u_gtrans
1843             swap_flux_x_local(k,j) = u_comp * (                              &
1844                         ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
1845                      +     7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
1846                      +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1  &
1847                         ) *                                                  &
1848                                     ( sk(k,j,i)   + sk(k,j,i-1)    )         &
1849                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
1850                      +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
1851                         ) *                                                  &
1852                                     ( sk(k,j,i+1) + sk(k,j,i-2)    )         &
1853                  +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
1854                         ) *         ( sk(k,j,i+2) + sk(k,j,i-3)    )         &
1855                                                  )
1856
1857              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
1858                         ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
1859                      +     3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
1860                      +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1  &
1861                         ) *                                                  &
1862                                     ( sk(k,j,i)   - sk(k,j,i-1)    )         &
1863                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
1864                      +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3  &
1865                         ) *                                                  &
1866                                     ( sk(k,j,i+1) - sk(k,j,i-2)  )           &
1867                  +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5  &
1868                         ) *                                                  &
1869                                     ( sk(k,j,i+2) - sk(k,j,i-3) )            &
1870                                                           )
1871
1872          ENDDO
1873
1874          DO  k = nzb_max+1, nzt
1875
1876             u_comp                 = u(k,j,i) - u_gtrans
1877             swap_flux_x_local(k,j) = u_comp * (                             &
1878                                      37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
1879                                    -  8.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             swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
1884                                      10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
1885                                    -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
1886                                    +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
1887                                                          ) * adv_sca_5
1888
1889          ENDDO
1890
1891       ENDDO
1892
1893       DO  i = nxl, nxr
1894
1895          j = nys
1896          DO  k = nzb+1, nzb_max
1897
1898             v_comp                  = v(k,j,i) - v_gtrans
1899             swap_flux_y_local(k) = v_comp * (                                &
1900                         ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
1901                      +     7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
1902                      +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1  &
1903                         ) *                                                  &
1904                                     ( sk(k,j,i)  + sk(k,j-1,i)     )         &
1905                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
1906                      +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
1907                         ) *                                                  &
1908                                     ( sk(k,j+1,i) + sk(k,j-2,i)    )         &
1909                  +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
1910                         ) *                                                  &
1911                                     ( sk(k,j+2,i) + sk(k,j-3,i)    )         &
1912                                             )
1913
1914             swap_diss_y_local(k) = -ABS( v_comp ) * (                        &
1915                         ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
1916                      +     3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
1917                      +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1  &
1918                         ) *                                                  &
1919                                     ( sk(k,j,i)   - sk(k,j-1,i)    )         &
1920                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
1921                      +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3  &
1922                         ) *                                                  &
1923                                     ( sk(k,j+1,i) - sk(k,j-2,i)  )           &
1924                  +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5  &
1925                         ) *                                                  &
1926                                     ( sk(k,j+2,i) - sk(k,j-3,i) )            &
1927                                                     )
1928
1929          ENDDO
1930!
1931!--       Above to the top of the highest topography. No degradation necessary.
1932          DO  k = nzb_max+1, nzt
1933
1934             v_comp               = v(k,j,i) - v_gtrans
1935             swap_flux_y_local(k) = v_comp * (                               &
1936                                    37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
1937                                  -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
1938                                  +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
1939                                             ) * adv_sca_5
1940              swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
1941                                    10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
1942                                  -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
1943                                  +          sk(k,j+2,i) - sk(k,j-3,i)       &
1944                                                      ) * adv_sca_5
1945
1946          ENDDO
1947
1948          DO  j = nys, nyn
1949
1950             flux_t(0) = 0.0
1951             diss_t(0) = 0.0
1952             flux_d    = 0.0
1953             diss_d    = 0.0
1954
1955             DO  k = nzb+1, nzb_max
1956
1957                u_comp    = u(k,j,i+1) - u_gtrans
1958                flux_r(k) = u_comp * (                                      &
1959                     ( 37.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
1960                  +     7.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
1961                  +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1    &
1962                     ) *                                                    &
1963                             ( sk(k,j,i+1) + sk(k,j,i)   )                  &
1964              -      (  8.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
1965                  +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
1966                     ) *                                                    &
1967                             ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
1968              +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
1969                     ) *     ( sk(k,j,i+3) + sk(k,j,i-2) )                  &
1970                                     )
1971
1972                diss_r(k) = -ABS( u_comp ) * (                              &
1973                     ( 10.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
1974                  +     3.0 * IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
1975                  +           IBITS(wall_flags_0(k,j,i),0,1) * adv_sca_1    &
1976                     ) *                                                    &
1977                             ( sk(k,j,i+1) - sk(k,j,i)  )                   &
1978              -      (  5.0 * IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
1979                  +           IBITS(wall_flags_0(k,j,i),1,1) * adv_sca_3    &
1980                     ) *                                                    &
1981                             ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
1982              +      (        IBITS(wall_flags_0(k,j,i),2,1) * adv_sca_5    &
1983                     ) *                                                    &
1984                             ( sk(k,j,i+3) - sk(k,j,i-2) )                  &
1985                                             )
1986
1987                v_comp    = v(k,j+1,i) - v_gtrans
1988                flux_n(k) = v_comp * (                                      &
1989                     ( 37.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
1990                  +     7.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
1991                  +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1    &
1992                     ) *                                                    &
1993                             ( sk(k,j+1,i) + sk(k,j,i)   )                  &
1994              -      (  8.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
1995                  +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
1996                     ) *                                                    &
1997                             ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
1998              +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
1999                     ) *                                                    &
2000                             ( sk(k,j+3,i) + sk(k,j-2,i) )                  &
2001                                     )
2002
2003                diss_n(k) = -ABS( v_comp ) * (                              &
2004                     ( 10.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
2005                  +     3.0 * IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
2006                  +           IBITS(wall_flags_0(k,j,i),3,1) * adv_sca_1    &
2007                     ) *                                                    &
2008                             ( sk(k,j+1,i) - sk(k,j,i)    )                 &
2009              -      (  5.0 * IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
2010                  +           IBITS(wall_flags_0(k,j,i),4,1) * adv_sca_3    &
2011                     ) *                                                    &
2012                             ( sk(k,j+2,i) - sk(k,j-1,i)  )                 &
2013              +      (        IBITS(wall_flags_0(k,j,i),5,1) * adv_sca_5    &
2014                     ) *                                                    &
2015                             ( sk(k,j+3,i) - sk(k,j-2,i) )                  &
2016                                             )
2017!
2018!--             k index has to be modified near bottom and top, else array
2019!--             subscripts will be exceeded.
2020                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
2021                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1)  )
2022                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
2023
2024
2025                flux_t(k) = w(k,j,i) * (                                    &
2026                     ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2027                  +     7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2028                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
2029                     ) *                                                    &
2030                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
2031              -      (  8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2032                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2033                     ) *                                                    &
2034                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
2035              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2036                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
2037                                       )
2038
2039                diss_t(k) = -ABS( w(k,j,i) ) * (                            &
2040                     ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2041                  +     3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2042                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
2043                     ) *                                                    &
2044                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
2045              -      (  5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2046                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2047                     ) *                                                    &
2048                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
2049              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2050                     ) *                                                    &
2051                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
2052                                               )
2053!
2054!--             Calculate the divergence of the velocity field. A respective
2055!--             correction is needed to overcome numerical instabilities caused
2056!--             by a not sufficient reduction of divergences near topography.
2057                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
2058                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
2059                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
2060
2061                tend(k,j,i) = tend(k,j,i) - (                                 &
2062                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
2063                          swap_diss_x_local(k,j)            ) * ddx           &
2064                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
2065                          swap_diss_y_local(k)              ) * ddy           &
2066                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
2067                                                               ) * ddzw(k)    &
2068                                            ) + sk(k,j,i) * div
2069
2070                swap_flux_y_local(k)   = flux_n(k)
2071                swap_diss_y_local(k)   = diss_n(k)
2072                swap_flux_x_local(k,j) = flux_r(k)
2073                swap_diss_x_local(k,j) = diss_r(k)
2074                flux_d                 = flux_t(k)
2075                diss_d                 = diss_t(k)
2076
2077             ENDDO
2078
2079             DO  k = nzb_max+1, nzt
2080
2081                u_comp    = u(k,j,i+1) - u_gtrans
2082                flux_r(k) = u_comp * (                                      &
2083                      37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                  &
2084                    -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                  &
2085                    +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
2086                diss_r(k) = -ABS( u_comp ) * (                              &
2087                      10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                  &
2088                    -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                  &
2089                    +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
2090
2091                v_comp    = v(k,j+1,i) - v_gtrans
2092                flux_n(k) = v_comp * (                                      &
2093                      37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                  &
2094                    -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                  &
2095                    +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
2096                diss_n(k) = -ABS( v_comp ) * (                              &
2097                      10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                  &
2098                    -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                  &
2099                    +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
2100!
2101!--             k index has to be modified near bottom and top, else array
2102!--             subscripts will be exceeded.
2103                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),8,1)
2104                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),6,1)  )
2105                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),8,1)
2106
2107                flux_t(k) = w(k,j,i) * (                                    &
2108                     ( 37.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2109                  +     7.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2110                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
2111                     ) *                                                    &
2112                             ( sk(k+1,j,i)  + sk(k,j,i)   )                 &
2113              -      (  8.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2114                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2115                     ) *                                                    &
2116                             ( sk(k_pp,j,i) + sk(k-1,j,i) )                 &
2117              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2118                     ) *     ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )                &
2119                                       )
2120
2121                diss_t(k) = -ABS( w(k,j,i) ) * (                            &
2122                     ( 10.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2123                  +     3.0 * IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2124                  +           IBITS(wall_flags_0(k,j,i),6,1) * adv_sca_1    &
2125                     ) *                                                    &
2126                             ( sk(k+1,j,i)   - sk(k,j,i)    )               &
2127              -      (  5.0 * IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2128                  +           IBITS(wall_flags_0(k,j,i),7,1) * adv_sca_3    &
2129                     ) *                                                    &
2130                             ( sk(k_pp,j,i)  - sk(k-1,j,i)  )               &
2131              +      (        IBITS(wall_flags_0(k,j,i),8,1) * adv_sca_5    &
2132                     ) *                                                    &
2133                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )               &
2134                                               )
2135!
2136!--             Calculate the divergence of the velocity field. A respective
2137!--             correction is needed to overcome numerical instabilities introduced
2138!--             by a not sufficient reduction of divergences near topography.
2139                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
2140                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
2141                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
2142
2143                tend(k,j,i) = tend(k,j,i) - (                                 &
2144                        ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) -    &
2145                          swap_diss_x_local(k,j)            ) * ddx           &
2146                      + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   -    &
2147                          swap_diss_y_local(k)              ) * ddy           &
2148                      + ( flux_t(k) + diss_t(k) - flux_d - diss_d             &
2149                                                               ) * ddzw(k)    &
2150                                            ) + sk(k,j,i) * div
2151
2152                swap_flux_y_local(k)   = flux_n(k)
2153                swap_diss_y_local(k)   = diss_n(k)
2154                swap_flux_x_local(k,j) = flux_r(k)
2155                swap_diss_x_local(k,j) = diss_r(k)
2156                flux_d                 = flux_t(k)
2157                diss_d                 = diss_t(k)
2158
2159             ENDDO
2160!
2161!--          evaluation of statistics
2162             SELECT CASE ( sk_char )
2163
2164                 CASE ( 'pt' )
2165                    DO  k = nzb, nzt
2166                       sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)         &
2167                        + ( flux_t(k) + diss_t(k) )                          &
2168                        *   weight_substep(intermediate_timestep_count)
2169                    ENDDO
2170                 CASE ( 'sa' )
2171                    DO  k = nzb, nzt
2172                       sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)         &
2173                        + ( flux_t(k) + diss_t(k) )                          &
2174                        *   weight_substep(intermediate_timestep_count)
2175                    ENDDO
2176                 CASE ( 'q' )
2177                    DO  k = nzb, nzt
2178                       sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)           &
2179                        + ( flux_t(k) + diss_t(k) )                          &
2180                        *   weight_substep(intermediate_timestep_count)
2181                    ENDDO
2182
2183              END SELECT
2184
2185         ENDDO
2186      ENDDO
2187
2188    END SUBROUTINE advec_s_ws
2189
2190
2191!------------------------------------------------------------------------------!
2192! Advection of u - Call for all grid points
2193!------------------------------------------------------------------------------!
2194    SUBROUTINE advec_u_ws
2195
2196       USE arrays_3d
2197       USE constants
2198       USE control_parameters
2199       USE grid_variables
2200       USE indices
2201       USE statistics
2202
2203       IMPLICIT NONE
2204
2205       INTEGER ::  i, j, k, k_mm, k_pp, k_ppp, tn = 0
2206       REAL    ::  diss_d, div, flux_d, gu, gv, v_comp, w_comp
2207       REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u, swap_flux_y_local_u
2208       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u,              &
2209                                             swap_flux_x_local_u
2210       REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r,     &
2211                                   flux_t, u_comp
2212 
2213       gu = 2.0 * u_gtrans
2214       gv = 2.0 * v_gtrans
2215
2216!
2217!--    Compute the fluxes for the whole left boundary of the processor domain.
2218       i = nxlu
2219       DO  j = nys, nyn
2220          DO  k = nzb+1, nzb_max
2221
2222             u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
2223             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
2224                          ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2225                       +     7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2226                       +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1 &
2227                          ) *                                                  &
2228                                     ( u(k,j,i)   + u(k,j,i-1) )               &
2229                   -      (  8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2230                       +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2231                          ) *                                                  &
2232                                     ( u(k,j,i+1) + u(k,j,i-2) )               &
2233                   +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2234                          ) *                                                  &
2235                                     ( u(k,j,i+2) + u(k,j,i-3) )               &
2236                                                   )
2237
2238              swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * (                &
2239                          ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2240                       +     3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2241                       +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1 &
2242                          ) *                                                  &
2243                                     ( u(k,j,i)   - u(k,j,i-1) )               &
2244                   -      (  5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2245                       +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3 &
2246                          ) *                                                  &
2247                                     ( u(k,j,i+1) - u(k,j,i-2) )               &
2248                   +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5 &
2249                          ) *                                                  &
2250                                     ( u(k,j,i+2) - u(k,j,i-3) )               &
2251                                                             )
2252
2253          ENDDO
2254
2255          DO  k = nzb_max+1, nzt
2256
2257             u_comp(k)         = u(k,j,i) + u(k,j,i-1) - gu
2258             swap_flux_x_local_u(k,j) = u_comp(k) * (                          &
2259                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )                &
2260                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )                &
2261                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
2262             swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                   &
2263                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )                &
2264                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )                &
2265                           +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
2266
2267          ENDDO
2268       ENDDO
2269
2270       DO i = nxlu, nxr
2271!       
2272!--       The following loop computes the fluxes for the south boundary points
2273          j = nys
2274          DO  k = nzb+1, nzb_max
2275
2276             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2277             swap_flux_y_local_u(k) = v_comp * (                              &
2278                         ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2279                      +     7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2280                      +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
2281                         ) *                                                  &
2282                                     ( u(k,j,i)   + u(k,j-1,i) )              &
2283                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2284                      +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2285                         ) *                                                  &
2286                                     ( u(k,j+1,i) + u(k,j-2,i) )              &
2287                  +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2288                         ) *                                                  &
2289                                     ( u(k,j+2,i) + u(k,j-3,i) )              &
2290                                               )
2291
2292             swap_diss_y_local_u(k) = - ABS ( v_comp ) * (                    &
2293                         ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2294                      +     3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2295                      +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1 &
2296                         ) *                                                  &
2297                                     ( u(k,j,i)   - u(k,j-1,i) )              &
2298                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2299                      +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3 &
2300                         ) *                                                  &
2301                                     ( u(k,j+1,i) - u(k,j-2,i) )              &
2302                  +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5 &
2303                         ) *                                                  &
2304                                     ( u(k,j+2,i) - u(k,j-3,i) )              &
2305                                                         )
2306
2307          ENDDO
2308
2309          DO  k = nzb_max+1, nzt
2310
2311             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2312             swap_flux_y_local_u(k) = v_comp * (                              &
2313                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
2314                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
2315                         +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
2316             swap_diss_y_local_u(k) = - ABS(v_comp) * (                       &
2317                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
2318                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
2319                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
2320
2321          ENDDO
2322!
2323!--       Computation of interior fluxes and tendency terms
2324          DO  j = nys, nyn
2325
2326             flux_t(0) = 0.0
2327             diss_t(0) = 0.0
2328             flux_d    = 0.0
2329             diss_d    = 0.0
2330
2331             DO  k = nzb+1, nzb_max
2332
2333                u_comp(k) = u(k,j,i+1) + u(k,j,i)
2334                flux_r(k) = ( u_comp(k) - gu ) * (                           &
2335                     ( 37.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
2336                  +     7.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
2337                  +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1    &
2338                     ) *                                                     &
2339                                 ( u(k,j,i+1) + u(k,j,i)   )                 &
2340              -      (  8.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
2341                  +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
2342                     ) *                                                     &
2343                                 ( u(k,j,i+2) + u(k,j,i-1) )                 &
2344              +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
2345                     ) *                                                     &
2346                                 ( u(k,j,i+3) + u(k,j,i-2) )                 &
2347                                                 )
2348
2349                diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
2350                     ( 10.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
2351                  +     3.0 * IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
2352                  +           IBITS(wall_flags_0(k,j,i),9,1)  * adv_mom_1    &
2353                     ) *                                                     &
2354                                 ( u(k,j,i+1) - u(k,j,i)  )                  &
2355              -      (  5.0 * IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
2356                  +           IBITS(wall_flags_0(k,j,i),10,1) * adv_mom_3    &
2357                     ) *                                                     &
2358                                 ( u(k,j,i+2) - u(k,j,i-1) )                 &
2359              +      (        IBITS(wall_flags_0(k,j,i),11,1) * adv_mom_5    &
2360                     ) *                                                     &
2361                                 ( u(k,j,i+3) - u(k,j,i-2) )                 &
2362                                                     )
2363
2364                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2365                flux_n(k) = v_comp * (                                       &
2366                     ( 37.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
2367                  +     7.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
2368                  +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1    &
2369                     ) *                                                     &
2370                                 ( u(k,j+1,i) + u(k,j,i)   )                 &
2371              -      (  8.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
2372                  +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
2373                     ) *                                                     &
2374                                 ( u(k,j+2,i) + u(k,j-1,i) )                 &
2375              +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
2376                     ) *                                                     &
2377                                 ( u(k,j+3,i) + u(k,j-2,i) )                 &
2378                                                 )
2379
2380                diss_n(k) = - ABS ( v_comp ) * (                             &
2381                     ( 10.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
2382                  +     3.0 * IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
2383                  +           IBITS(wall_flags_0(k,j,i),12,1) * adv_mom_1    &
2384                     ) *                                                     &
2385                                 ( u(k,j+1,i) - u(k,j,i)  )                  &
2386              -      (  5.0 * IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
2387                  +           IBITS(wall_flags_0(k,j,i),13,1) * adv_mom_3    &
2388                     ) *                                                     &
2389                                 ( u(k,j+2,i) - u(k,j-1,i) )                 &
2390              +      (        IBITS(wall_flags_0(k,j,i),14,1) * adv_mom_5    &
2391                     ) *                                                     &
2392                                 ( u(k,j+3,i) - u(k,j-2,i) )                 &
2393                                                      )
2394!
2395!--             k index has to be modified near bottom and top, else array
2396!--             subscripts will be exceeded.
2397                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
2398                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1)  )
2399                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
2400
2401                w_comp    = w(k,j,i) + w(k,j,i-1)
2402                flux_t(k) = w_comp  * (                                      &
2403                     ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2404                  +     7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2405                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
2406                     ) *                                                     &
2407                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
2408              -      (  8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2409                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2410                     ) *                                                     &
2411                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
2412              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2413                     ) *                                                     &
2414                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
2415                                      )
2416
2417                diss_t(k) = - ABS( w_comp ) * (                              &
2418                     ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2419                  +     3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2420                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
2421                     ) *                                                     &
2422                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
2423              -      (  5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2424                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2425                     ) *                                                     &
2426                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
2427              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2428                     ) *                                                     &
2429                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
2430                                              )
2431!
2432!--             Calculate the divergence of the velocity field. A respective
2433!--             correction is needed to overcome numerical instabilities caused
2434!--             by a not sufficient reduction of divergences near topography.
2435                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
2436                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
2437                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
2438                                                                    * ddzw(k)  &
2439                      ) * 0.5
2440
2441                tend(k,j,i) = tend(k,j,i) - (                                  &
2442                 ( flux_r(k) + diss_r(k)                                       &
2443               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
2444               + ( flux_n(k) + diss_n(k)                                       &
2445               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
2446               + ( flux_t(k) + diss_t(k)                                       &
2447               -   flux_d    - diss_d                                          &
2448                                                                  ) * ddzw(k)  &
2449                                           ) + div * u(k,j,i)
2450
2451                swap_flux_x_local_u(k,j) = flux_r(k)
2452                swap_diss_x_local_u(k,j) = diss_r(k)
2453                swap_flux_y_local_u(k)   = flux_n(k)
2454                swap_diss_y_local_u(k)   = diss_n(k)
2455                flux_d                   = flux_t(k)
2456                diss_d                   = diss_t(k)
2457!
2458!--             Statistical Evaluation of u'u'. The factor has to be applied
2459!--             for right evaluation when gallilei_trans = .T. .
2460                sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                     &
2461                              + ( flux_r(k) *                                 &
2462                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
2463                              / ( u_comp(k) - gu + 1.0E-20      )             &
2464                              +   diss_r(k) *                                 &
2465                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
2466                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
2467                              *   weight_substep(intermediate_timestep_count)
2468!
2469!--             Statistical Evaluation of w'u'.
2470                sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                   &
2471                              + ( flux_t(k) + diss_t(k) )                     &
2472                              *   weight_substep(intermediate_timestep_count)
2473             ENDDO
2474
2475             DO  k = nzb_max+1, nzt
2476
2477                u_comp(k) = u(k,j,i+1) + u(k,j,i)
2478                flux_r(k) = ( u_comp(k) - gu ) * (                            &
2479                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
2480                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
2481                       +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2482                diss_r(k) = - ABS( u_comp(k) - gu ) * (                       &
2483                         10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
2484                       -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
2485                       +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2486
2487                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2488                flux_n(k) = v_comp * (                                        &
2489                         37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
2490                       -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                   &
2491                       +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2492                diss_n(k) = - ABS( v_comp ) * (                               &
2493                         10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
2494                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
2495                       +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2496!
2497!--             k index has to be modified near bottom and top, else array
2498!--             subscripts will be exceeded.
2499                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),17,1)
2500                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),15,1)  )
2501                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),17,1)
2502
2503                w_comp    = w(k,j,i) + w(k,j,i-1)
2504                flux_t(k) = w_comp  * (                                      &
2505                     ( 37.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2506                  +     7.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2507                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
2508                     ) *                                                     &
2509                             ( u(k+1,j,i)  + u(k,j,i)     )                  &
2510              -      (  8.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2511                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2512                     ) *                                                     &
2513                             ( u(k_pp,j,i) + u(k-1,j,i)   )                  &
2514              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2515                     ) *                                                     &
2516                             ( u(k_ppp,j,i) + u(k_mm,j,i) )                  &
2517                                      )
2518
2519                 diss_t(k) = - ABS( w_comp ) * (                             &
2520                     ( 10.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2521                  +     3.0 * IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2522                  +           IBITS(wall_flags_0(k,j,i),15,1) * adv_mom_1    &
2523                     ) *                                                     &
2524                             ( u(k+1,j,i)   - u(k,j,i)    )                  &
2525              -      (  5.0 * IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2526                  +           IBITS(wall_flags_0(k,j,i),16,1) * adv_mom_3    &
2527                     ) *                                                     &
2528                             ( u(k_pp,j,i)  - u(k-1,j,i)  )                  &
2529              +      (        IBITS(wall_flags_0(k,j,i),17,1) * adv_mom_5    &
2530                     ) *                                                     &
2531                             ( u(k_ppp,j,i) - u(k_mm,j,i) )                  &
2532                                               )
2533!
2534!--             Calculate the divergence of the velocity field. A respective
2535!--             correction is needed to overcome numerical instabilities caused
2536!--             by a not sufficient reduction of divergences near topography.
2537                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
2538                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy &
2539                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )       &
2540                                                                    * ddzw(k) &
2541                      ) * 0.5
2542
2543                 tend(k,j,i) = tend(k,j,i) - (                                 &
2544                 ( flux_r(k) + diss_r(k)                                       &
2545               -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
2546               + ( flux_n(k) + diss_n(k)                                       &
2547               -   swap_flux_y_local_u(k)   - swap_diss_y_local_u(k)   ) * ddy &
2548               + ( flux_t(k) + diss_t(k)                                       &
2549               -   flux_d    - diss_d                                          &
2550                                                                  ) * ddzw(k)  &
2551                                           ) + div * u(k,j,i)
2552
2553                 swap_flux_x_local_u(k,j) = flux_r(k)
2554                 swap_diss_x_local_u(k,j) = diss_r(k)
2555                 swap_flux_y_local_u(k)   = flux_n(k)
2556                 swap_diss_y_local_u(k)   = diss_n(k)
2557                 flux_d                   = flux_t(k)
2558                 diss_d                   = diss_t(k)
2559!
2560!--              Statistical Evaluation of u'u'. The factor has to be applied
2561!--              for right evaluation when gallilei_trans = .T. .
2562                 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                    &
2563                              + ( flux_r(k) *                                 &
2564                                ( u_comp(k) - 2.0 * hom(k,1,1,0) )            &
2565                              / ( u_comp(k) - gu + 1.0E-20      )             &
2566                              +   diss_r(k) *                                 &
2567                                  ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )       &
2568                              / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )         &
2569                              *   weight_substep(intermediate_timestep_count)
2570!
2571!--              Statistical Evaluation of w'u'.
2572                 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                  &
2573                              + ( flux_t(k) + diss_t(k) )                     &
2574                              *   weight_substep(intermediate_timestep_count)
2575       ENDDO
2576          ENDDO
2577       ENDDO
2578       sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
2579
2580
2581    END SUBROUTINE advec_u_ws
2582   
2583   
2584!------------------------------------------------------------------------------!
2585! Advection of v - Call for all grid points
2586!------------------------------------------------------------------------------!
2587    SUBROUTINE advec_v_ws
2588
2589       USE arrays_3d
2590       USE constants
2591       USE control_parameters
2592       USE grid_variables
2593       USE indices
2594       USE statistics
2595
2596       IMPLICIT NONE
2597
2598
2599       INTEGER ::  i, j, k, k_mm, k_pp, k_ppp, tn = 0
2600       REAL    ::  diss_d, div, flux_d, gu, gv, u_comp, w_comp
2601       REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v, swap_flux_y_local_v
2602       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v,             &
2603                                             swap_flux_x_local_v
2604       REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r,    &
2605                                   flux_t, v_comp
2606
2607       gu = 2.0 * u_gtrans
2608       gv = 2.0 * v_gtrans
2609!
2610!--    First compute the whole left boundary of the processor domain
2611       i = nxl
2612       DO  j = nysv, nyn
2613          DO  k = nzb+1, nzb_max
2614
2615             u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
2616             swap_flux_x_local_v(k,j) = u_comp * (                             &
2617                          ( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2618                       +     7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2619                       +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 &
2620                          ) *                                                  &
2621                                     ( v(k,j,i)   + v(k,j,i-1) )               &
2622                   -      (  8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2623                       +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2624                          ) *                                                  &
2625                                     ( v(k,j,i+1) + v(k,j,i-2) )               &
2626                   +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2627                          ) *                                                  &
2628                                     ( v(k,j,i+2) + v(k,j,i-3) )               &
2629                                                 )
2630
2631              swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                   &
2632                          ( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2633                       +     3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2634                       +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1 &
2635                          ) *                                                  &
2636                                     ( v(k,j,i)   - v(k,j,i-1) )               &
2637                   -      (  5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2638                       +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3 &
2639                          ) *                                                  &
2640                                     ( v(k,j,i+1) - v(k,j,i-2) )               &
2641                   +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5 &
2642                          ) *                                                  &
2643                                     ( v(k,j,i+2) - v(k,j,i-3) )               &
2644                                                           )
2645
2646          ENDDO
2647
2648          DO  k = nzb_max+1, nzt
2649
2650             u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
2651             swap_flux_x_local_v(k,j) = u_comp * (                            &
2652                             37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
2653                           -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
2654                           +        ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
2655             swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                   &
2656                             10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
2657                           -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
2658                           +        ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
2659
2660          ENDDO
2661
2662       ENDDO
2663
2664       DO i = nxl, nxr
2665
2666          j = nysv
2667          DO  k = nzb+1, nzb_max
2668
2669             v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
2670             swap_flux_y_local_v(k) = v_comp(k) * (                           &
2671                         ( 37.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
2672                      +     7.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
2673                      +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 &
2674                         ) *                                                  &
2675                                     ( v(k,j,i)   + v(k,j-1,i) )              &
2676                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
2677                      +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
2678                         ) *                                                  &
2679                                     ( v(k,j+1,i) + v(k,j-2,i) )              &
2680                  +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
2681                         ) *                                                  &
2682                                     ( v(k,j+2,i) + v(k,j-3,i) )              &
2683                                                 )
2684
2685             swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * (                  &
2686                         ( 10.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
2687                      +     3.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
2688                      +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1 &
2689                         ) *                                                  &
2690                                     ( v(k,j,i)   - v(k,j-1,i) )              &
2691                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
2692                      +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3 &
2693                         ) *                                                  &
2694                                     ( v(k,j+1,i) - v(k,j-2,i) )              &
2695                  +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5 &
2696                         ) *                                                  &
2697                                     ( v(k,j+2,i) - v(k,j-3,i) )              &
2698                                                          )
2699
2700          ENDDO
2701
2702          DO  k = nzb_max+1, nzt
2703
2704             v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
2705             swap_flux_y_local_v(k) = v_comp(k) * (                           &
2706                           37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
2707                         -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
2708                         +        ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
2709             swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * (                  &
2710                           10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
2711                         -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
2712                         +        ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
2713
2714          ENDDO
2715
2716          DO  j = nysv, nyn
2717
2718             flux_t(0) = 0.0
2719             diss_t(0) = 0.0
2720             flux_d    = 0.0
2721             diss_d    = 0.0
2722
2723             DO  k = nzb+1, nzb_max
2724
2725                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2726                flux_r(k) = u_comp * (                                       &
2727                     ( 37.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
2728                  +     7.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
2729                  +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1    &
2730                     ) *                                                     &
2731                                 ( v(k,j,i+1) + v(k,j,i)   )                 &
2732              -      (  8.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
2733                  +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
2734                     ) *                                                     &
2735                                 ( v(k,j,i+2) + v(k,j,i-1) )                 &
2736              +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
2737                     ) *                                                     &
2738                                 ( v(k,j,i+3) + v(k,j,i-2) )                 &
2739                                     )
2740
2741                diss_r(k) = - ABS( u_comp ) * (                              &
2742                     ( 10.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
2743                  +     3.0 * IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
2744                  +           IBITS(wall_flags_0(k,j,i),18,1) * adv_mom_1    &
2745                     ) *                                                     &
2746                                 ( v(k,j,i+1) - v(k,j,i)  )                  &
2747              -      (  5.0 * IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
2748                  +           IBITS(wall_flags_0(k,j,i),19,1) * adv_mom_3    &
2749                     ) *                                                     &
2750                                 ( v(k,j,i+2) - v(k,j,i-1) )                 &
2751              +      (        IBITS(wall_flags_0(k,j,i),20,1) * adv_mom_5    &
2752                     ) *                                                     &
2753                                 ( v(k,j,i+3) - v(k,j,i-2) )                 &
2754                                              )
2755
2756                v_comp(k) = v(k,j+1,i) + v(k,j,i)
2757                flux_n(k) = ( v_comp(k) - gv ) * (                           &
2758                     ( 37.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
2759                  +     7.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
2760                  +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1    &
2761                     ) *                                                     &
2762                                 ( v(k,j+1,i) + v(k,j,i)   )                 &
2763              -      (  8.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
2764                  +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
2765                     ) *                                                     &
2766                                 ( v(k,j+2,i) + v(k,j-1,i) )                 &
2767              +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
2768                     ) *                                                     &
2769                                 ( v(k,j+3,i) + v(k,j-2,i) )                 &
2770                                     )
2771
2772                diss_n(k) = - ABS( v_comp(k) - gv ) * (                      &
2773                     ( 10.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
2774                  +     3.0 * IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
2775                  +           IBITS(wall_flags_0(k,j,i),21,1) * adv_mom_1    &
2776                     ) *                                                     &
2777                                 ( v(k,j+1,i) - v(k,j,i)  )                  &
2778              -      (  5.0 * IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
2779                  +           IBITS(wall_flags_0(k,j,i),22,1) * adv_mom_3    &
2780                     ) *                                                     &
2781                                 ( v(k,j+2,i) - v(k,j-1,i) )                 &
2782              +      (        IBITS(wall_flags_0(k,j,i),23,1) * adv_mom_5    &
2783                     ) *                                                     &
2784                                 ( v(k,j+3,i) - v(k,j-2,i) )                 &
2785                                                      )
2786!
2787!--             k index has to be modified near bottom and top, else array
2788!--             subscripts will be exceeded.
2789                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),26,1)
2790                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),24,1)  )
2791                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),26,1)
2792
2793                w_comp    = w(k,j-1,i) + w(k,j,i)
2794                flux_t(k) = w_comp  * (                                      &
2795                     ( 37.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2796                  +     7.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2797                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
2798                     ) *                                                     &
2799                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
2800              -      (  8.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2801                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2802                     ) *                                                     &
2803                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
2804              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2805                     ) *                                                     &
2806                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
2807                                      )
2808
2809                diss_t(k) = - ABS( w_comp ) * (                              &
2810                     ( 10.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2811                  +     3.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2812                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
2813                     ) *                                                     &
2814                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
2815              -      (  5.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2816                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2817                     ) *                                                     &
2818                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
2819              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2820                     ) *                                                     &
2821                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
2822                                               )
2823!
2824!--             Calculate the divergence of the velocity field. A respective
2825!--             correction is needed to overcome numerical instabilities caused
2826!--             by a not sufficient reduction of divergences near topography.
2827                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
2828                     +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
2829                     +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) )         &
2830                                                                  ) * ddzw(k) &
2831                      ) * 0.5
2832
2833                tend(k,j,i) = tend(k,j,i) - (                                 &
2834                       ( flux_r(k) + diss_r(k)                                &
2835                     -   swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j)  &
2836                       ) * ddx                                                &
2837                     + ( flux_n(k) + diss_n(k)                                &
2838                     -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
2839                       ) * ddy                                                &
2840                     + ( flux_t(k) + diss_t(k)                                &
2841                     -   flux_d    - diss_d                                   &
2842                       ) * ddzw(k)                                            &
2843                                            )  + v(k,j,i) * div
2844
2845                swap_flux_x_local_v(k,j) = flux_r(k)
2846                swap_diss_x_local_v(k,j) = diss_r(k)
2847                swap_flux_y_local_v(k)   = flux_n(k)
2848                swap_diss_y_local_v(k)   = diss_n(k)
2849                flux_d                   = flux_t(k)
2850                diss_d                   = diss_t(k)
2851
2852!
2853!--             Statistical Evaluation of v'v'. The factor has to be applied
2854!--             for right evaluation when gallilei_trans = .T. .
2855                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                     &
2856                      + ( flux_n(k)                                           &
2857                      * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                    &
2858                      / ( v_comp(k) - gv + 1.0E-20 )                          &
2859                      +   diss_n(k)                                           &
2860                      *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )               &
2861                      / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                  &
2862                      *   weight_substep(intermediate_timestep_count)
2863!
2864!--              Statistical Evaluation of w'v'.
2865                 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                  &
2866                              + ( flux_t(k) + diss_t(k) )                     &
2867                              *   weight_substep(intermediate_timestep_count)
2868
2869             ENDDO
2870
2871             DO  k = nzb_max+1, nzt
2872
2873                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2874                flux_r(k) = u_comp * (                                        &
2875                      37.0 * ( v(k,j,i+1) + v(k,j,i)   )                      &
2876                    -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                      &
2877                    +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
2878
2879                diss_r(k) = - ABS( u_comp ) * (                               &
2880                      10.0 * ( v(k,j,i+1) - v(k,j,i) )                        &
2881                    -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                      &
2882                    +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
2883
2884
2885                v_comp(k) = v(k,j+1,i) + v(k,j,i)
2886                flux_n(k) = ( v_comp(k) - gv ) * (                            &
2887                      37.0 * ( v(k,j+1,i) + v(k,j,i)   )                      &
2888                    -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                      &
2889                      +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
2890
2891                diss_n(k) = - ABS( v_comp(k) - gv ) * (                       &
2892                      10.0 * ( v(k,j+1,i) - v(k,j,i)   )                      &
2893                    -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                      &
2894                    +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
2895!
2896!--             k index has to be modified near bottom and top, else array
2897!--             subscripts will be exceeded.
2898                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),26,1)
2899                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),24,1)  )
2900                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),26,1)
2901
2902                w_comp    = w(k,j-1,i) + w(k,j,i)
2903                flux_t(k) = w_comp  * (                                      &
2904                     ( 37.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2905                  +     7.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2906                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
2907                     ) *                                                     &
2908                             ( v(k+1,j,i)   + v(k,j,i)    )                  &
2909              -      (  8.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2910                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2911                     ) *                                                     &
2912                             ( v(k_pp,j,i)  + v(k-1,j,i)  )                  &
2913              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2914                     ) *                                                     &
2915                             ( v(k_ppp,j,i) + v(k_mm,j,i) )                  &
2916                                 )
2917
2918                diss_t(k) = - ABS( w_comp ) * (                              &
2919                     ( 10.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2920                  +     3.0 * IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2921                  +           IBITS(wall_flags_0(k,j,i),24,1) * adv_mom_1    &
2922                     ) *                                                     &
2923                             ( v(k+1,j,i)   - v(k,j,i)    )                  &
2924              -      (  5.0 * IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2925                  +           IBITS(wall_flags_0(k,j,i),25,1) * adv_mom_3    &
2926                     ) *                                                     &
2927                             ( v(k_pp,j,i)  - v(k-1,j,i)  )                  &
2928              +      (        IBITS(wall_flags_0(k,j,i),26,1) * adv_mom_5    &
2929                     ) *                                                     &
2930                             ( v(k_ppp,j,i) - v(k_mm,j,i) )                  &
2931                                         )
2932!
2933!--             Calculate the divergence of the velocity field. A respective
2934!--             correction is needed to overcome numerical instabilities caused
2935!--             by a not sufficient reduction of divergences near topography.
2936                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
2937                     +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
2938                     +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) )       &
2939                                                                    * ddzw(k) &
2940                      ) * 0.5
2941 
2942                tend(k,j,i) = tend(k,j,i) - (                                 &
2943                       ( flux_r(k) + diss_r(k)                                &
2944                     -   swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j)  &
2945                       ) * ddx                                                &
2946                     + ( flux_n(k) + diss_n(k)                                &
2947                     -   swap_flux_y_local_v(k) - swap_diss_y_local_v(k)      &
2948                       ) * ddy                                                &
2949                     + ( flux_t(k) + diss_t(k)                                &
2950                     -   flux_d    - diss_d                                   &
2951                       ) * ddzw(k)                                            &
2952                                            )  + v(k,j,i) * div
2953
2954                swap_flux_x_local_v(k,j) = flux_r(k)
2955                swap_diss_x_local_v(k,j) = diss_r(k)
2956                swap_flux_y_local_v(k)   = flux_n(k)
2957                swap_diss_y_local_v(k)   = diss_n(k)
2958                flux_d                   = flux_t(k)
2959                diss_d                   = diss_t(k)
2960
2961!
2962!--             Statistical Evaluation of v'v'. The factor has to be applied
2963!--             for right evaluation when gallilei_trans = .T. .
2964                sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                     &
2965                         + ( flux_n(k)                                        &
2966                         * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                 &
2967                         / ( v_comp(k) - gv + 1.0E-20 )                       &
2968                         +   diss_n(k)                                        &
2969                         *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )            &
2970                         / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )               &
2971                         *   weight_substep(intermediate_timestep_count)
2972!
2973!--             Statistical Evaluation of w'v'.
2974                sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                    &
2975                              + ( flux_t(k) + diss_t(k) )                      &
2976                              *   weight_substep(intermediate_timestep_count)
2977
2978             ENDDO
2979          ENDDO
2980       ENDDO
2981       sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
2982
2983
2984    END SUBROUTINE advec_v_ws
2985   
2986   
2987!------------------------------------------------------------------------------!
2988! Advection of w - Call for all grid points
2989!------------------------------------------------------------------------------!
2990    SUBROUTINE advec_w_ws
2991
2992       USE arrays_3d
2993       USE constants
2994       USE control_parameters
2995       USE grid_variables
2996       USE indices
2997       USE statistics
2998
2999       IMPLICIT NONE
3000
3001       INTEGER ::  i, j, k, k_mm, k_pp, k_ppp, tn = 0
3002       REAL    ::  diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp
3003       REAL, DIMENSION(nzb:nzt)    ::  diss_t, flux_t
3004       REAL, DIMENSION(nzb+1:nzt)  ::  diss_n, diss_r, flux_n, flux_r,        &
3005                                       swap_diss_y_local_w,                   &
3006                                       swap_flux_y_local_w
3007       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_w,             &
3008                                             swap_flux_x_local_w
3009 
3010       gu = 2.0 * u_gtrans
3011       gv = 2.0 * v_gtrans
3012!
3013!--   compute the whole left boundary of the processor domain
3014       i = nxl
3015       DO  j = nys, nyn
3016          DO  k = nzb+1, nzb_max
3017
3018             u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
3019             swap_flux_x_local_w(k,j) = u_comp * (                             &
3020                          ( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
3021                       +     7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
3022                       +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 &
3023                          ) *                                                  &
3024                                     ( w(k,j,i)   + w(k,j,i-1) )               &
3025                   -      (  8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
3026                       +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
3027                          ) *                                                  &
3028                                     ( w(k,j,i+1) + w(k,j,i-2) )               &
3029                   +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
3030                          ) *                                                  &
3031                                     ( w(k,j,i+2) + w(k,j,i-3) )               &
3032                                                 )
3033
3034               swap_diss_x_local_w(k,j) = - ABS( u_comp ) * (                  &
3035                          ( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
3036                       +     3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
3037                       +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1 &
3038                          ) *                                                  &
3039                                     ( w(k,j,i)   - w(k,j,i-1) )               &
3040                   -      (  5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
3041                       +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3 &
3042                          ) *                                                  &
3043                                     ( w(k,j,i+1) - w(k,j,i-2) )               &
3044                   +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5 &
3045                          ) *                                                  &
3046                                     ( w(k,j,i+2) - w(k,j,i-3) )               &
3047                                                            )
3048
3049          ENDDO
3050
3051          DO  k = nzb_max+1, nzt
3052
3053             u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
3054             swap_flux_x_local_w(k,j) = u_comp * (                             &
3055                            37.0 * ( w(k,j,i) + w(k,j,i-1)   )                 &
3056                          -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                 &
3057                          +        ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
3058             swap_diss_x_local_w(k,j) = - ABS( u_comp ) * (                    &
3059                            10.0 * ( w(k,j,i) - w(k,j,i-1)   )                 &
3060                          -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                 &
3061                          +        ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
3062
3063          ENDDO
3064
3065       ENDDO
3066
3067       DO i = nxl, nxr
3068
3069          j = nys
3070          DO  k = nzb+1, nzb_max
3071
3072             v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
3073             swap_flux_y_local_w(k) = v_comp * (                              &
3074                         ( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
3075                      +     7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
3076                      +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 &
3077                         ) *                                                  &
3078                                     ( w(k,j,i)   + w(k,j-1,i) )              &
3079                  -      (  8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
3080                      +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
3081                         ) *                                                  &
3082                                     ( w(k,j+1,i) + w(k,j-2,i) )              &
3083                  +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
3084                         ) *                                                  &
3085                                     ( w(k,j+2,i) + w(k,j-3,i) )              &
3086                                               )
3087
3088             swap_diss_y_local_w(k) = - ABS( v_comp ) * (                     &
3089                         ( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
3090                      +     3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
3091                      +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1 &
3092                         ) *                                                  &
3093                                     ( w(k,j,i)   - w(k,j-1,i) )              &
3094                  -      (  5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
3095                      +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3 &
3096                         ) *                                                  &
3097                                     ( w(k,j+1,i) - w(k,j-2,i) )              &
3098                  +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5 &
3099                         ) *                                                  &
3100                                     ( w(k,j+2,i) - w(k,j-3,i) )              &
3101                                                        )
3102
3103          ENDDO
3104
3105          DO  k = nzb_max+1, nzt
3106
3107             v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
3108             swap_flux_y_local_w(k) = v_comp * (                              &
3109                           37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
3110                         -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
3111                         +        ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
3112             swap_diss_y_local_w(k) = - ABS( v_comp ) * (                     &
3113                           10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
3114                         -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
3115                         +        ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
3116
3117          ENDDO
3118
3119          DO  j = nys, nyn
3120
3121!
3122!--          The lower flux has to be calculated explicetely for the tendency
3123!--          at the first w-level. For topography wall this is done implicitely
3124!--          by wall_flags_0.
3125             k         = nzb + 1
3126             w_comp    = w(k,j,i) + w(k-1,j,i)
3127             flux_t(0) = w_comp       * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
3128             diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
3129             flux_d    = flux_t(0)
3130             diss_d    = diss_t(0)
3131
3132             DO  k = nzb+1, nzb_max
3133
3134                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
3135                flux_r(k) = u_comp * (                                       &
3136                     ( 37.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
3137                  +     7.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
3138                  +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1    &
3139                     ) *                                                     &
3140                                 ( w(k,j,i+1) + w(k,j,i)   )                 &
3141              -      (  8.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
3142                  +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
3143                     ) *                                                     &
3144                                 ( w(k,j,i+2) + w(k,j,i-1) )                 &
3145              +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
3146                     ) *                                                     &
3147                                 ( w(k,j,i+3) + w(k,j,i-2) )                 &
3148                                     )
3149
3150                diss_r(k) = - ABS( u_comp ) * (                              &
3151                     ( 10.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
3152                  +     3.0 * IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
3153                  +           IBITS(wall_flags_0(k,j,i),27,1) * adv_mom_1    &
3154                     ) *                                                     &
3155                                 ( w(k,j,i+1) - w(k,j,i)  )                  &
3156              -      (  5.0 * IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
3157                  +           IBITS(wall_flags_0(k,j,i),28,1) * adv_mom_3    &
3158                     ) *                                                     &
3159                                 ( w(k,j,i+2) - w(k,j,i-1) )                 &
3160              +      (        IBITS(wall_flags_0(k,j,i),29,1) * adv_mom_5    &
3161                     ) *                                                     &
3162                                 ( w(k,j,i+3) - w(k,j,i-2) )                 &
3163                                              )
3164
3165                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
3166                flux_n(k) = v_comp * (                                       &
3167                     ( 37.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
3168                  +     7.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
3169                  +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1    &
3170                     ) *                                                     &
3171                                 ( w(k,j+1,i) + w(k,j,i)   )                 &
3172              -      (  8.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
3173                  +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
3174                     ) *                                                     &
3175                                 ( w(k,j+2,i) + w(k,j-1,i) )                 &
3176              +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
3177                     ) *                                                     &
3178                                 ( w(k,j+3,i) + w(k,j-2,i) )                 &
3179                                     )
3180
3181                diss_n(k) = - ABS( v_comp ) * (                              &
3182                     ( 10.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
3183                  +     3.0 * IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
3184                  +           IBITS(wall_flags_0(k,j,i),30,1) * adv_mom_1    &
3185                     ) *                                                     &
3186                                 ( w(k,j+1,i) - w(k,j,i)  )                  &
3187              -      (  5.0 * IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
3188                  +           IBITS(wall_flags_0(k,j,i),31,1) * adv_mom_3    &
3189                     ) *                                                     &
3190                                 ( w(k,j+2,i) - w(k,j-1,i) )                 &
3191              +      (        IBITS(wall_flags_0(k,j,i),32,1) * adv_mom_5    &
3192                     ) *                                                     &
3193                                 ( w(k,j+3,i) - w(k,j-2,i) )                 &
3194                                              )
3195!
3196!--             k index has to be modified near bottom and top, else array
3197!--             subscripts will be exceeded.
3198                k_ppp = k + 3 * IBITS(wall_flags_0(k,j,i),35,1)
3199                k_pp  = k + 2 * ( 1 - IBITS(wall_flags_0(k,j,i),33,1)  )
3200                k_mm  = k - 2 * IBITS(wall_flags_0(k,j,i),35,1)
3201
3202                w_comp    = w(k+1,j,i) + w(k,j,i)
3203                flux_t(k) = w_comp  * (                                      &
3204                     ( 37.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
3205                  +     7.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
3206                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
3207                     ) *                                                     &
3208                             ( w(k+1,j,i)  + w(k,j,i)     )                  &
3209              -      (  8.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
3210                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
3211                     ) *                                                     &
3212                             ( w(k_pp,j,i)  + w(k-1,j,i)  )                  &
3213              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
3214                     ) *                                                     &
3215                             ( w(k_ppp,j,i) + w(k_mm,j,i) )                  &
3216                                       )
3217
3218                diss_t(k) = - ABS( w_comp ) * (                              &
3219                     ( 10.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
3220                  +     3.0 * IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
3221                  +           IBITS(wall_flags_0(k,j,i),33,1) * adv_mom_1    &
3222                     ) *                                                     &
3223                             ( w(k+1,j,i)   - w(k,j,i)    )                  &
3224              -      (  5.0 * IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
3225                  +           IBITS(wall_flags_0(k,j,i),34,1) * adv_mom_3    &
3226                     ) *                                                     &
3227                             ( w(k_pp,j,i)  - w(k-1,j,i)  )                  &
3228              +      (        IBITS(wall_flags_0(k,j,i),35,1) * adv_mom_5    &
3229                     ) *                                                     &
3230                             ( w(k_ppp,j,i) - w(k_mm,j,i) )                  &
3231                                               )
3232!
3233!--             Calculate the divergence of the velocity field. A respective
3234!--             correction is needed to overcome numerical instabilities caused
3235!--             by a not sufficient reduction of divergences near topography.
3236                div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
3237                    +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
3238                    +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
3239                                                                 * ddzu(k+1) &
3240                      ) * 0.5
3241
3242                tend(k,j,i) = tend(k,j,i) - (                                 &
3243                      ( flux_r(k) + diss_r(k)                                 &
3244                    -   swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j)   &
3245                      ) * ddx                                                 &
3246                    + ( flux_n(k) + diss_n(k)                                 &
3247                    -   swap_flux_y_local_w(k)   - swap_diss_y_local_w(k)     &
3248                      ) * ddy                                                 &
3249                    + ( flux_t(k) + diss_t(k)                                 &
3250                    -   flux_d    - diss_d                                    &
3251                      ) * ddzu(k+1)                                           &
3252                                            )  + div * w(k,j,i)
3253
3254                swap_flux_x_local_w(k,j) = flux_r(k)
3255                swap_diss_x_local_w(k,j) = diss_r(k)
3256                swap_flux_y_local_w(k)   = flux_n(k)
3257                swap_diss_y_local_w(k)   = diss_n(k)
3258                flux_d                   = flux_t(k)
3259                diss_d                   = diss_t(k)
3260
3261                sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                    &
3262                               + ( flux_t(k) + diss_t(k) )                    &
3263                               *   weight_substep(intermediate_timestep_count)
3264
3265             ENDDO
3266
3267             DO  k = nzb_max+1, nzt
3268
3269                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
3270                flux_r(k) = u_comp * (                                      &
3271                      37.0 * ( w(k,j,i+1) + w(k,j,i)   )