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

Last change on this file since 802 was 802, checked in by suehring, 11 years ago

Last commit documented.

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 161.1 KB
Line 
1 MODULE advec_ws
2
3!-----------------------------------------------------------------------------!
4! Current revisions:
5! ------------------
6!
7! Former revisions:
8! -----------------
9! $Id: advec_ws.f90 802 2012-01-10 17:34:03Z suehring $
10!
11! 801 2012-01-10 17:30:36Z suehring
12! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
13! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
14! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
15! dimension.
16!
17! 743 2011-08-18 16:10:16Z suehring
18! Evaluation of turbulent fluxes with WS-scheme only for the whole model
19! domain. Therefor dimension of arrays needed for statistical evaluation
20! decreased.
21!
22! 736 2011-08-17 14:13:26Z suehring
23! Bugfix concerning OpenMP parallelization. i_omp introduced, because first
24! index where fluxes on left side have to be calculated explicitly is
25! different on each thread. Furthermore the swapping of fluxes is now
26! thread-safe by adding an additional dimension.
27!
28! 713 2011-03-30 14:21:21Z suehring
29! File reformatted.
30! Bugfix in vertical advection of w concerning the optimized version for   
31! vector architecture.
32! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as
33! broken numbers.
34!
35! 709 2011-03-30 09:31:40Z raasch
36! formatting adjustments
37!
38! 705 2011-03-25 11:21:43 Z suehring $
39! Bugfix in declaration of logicals concerning outflow boundaries.
40!
41! 411 2009-12-11 12:31:43 Z suehring
42! Allocation of weight_substep moved to init_3d_model.
43! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters.
44! Setting bc for the horizontal velocity variances added (moved from
45! flow_statistics).
46!
47! Initial revision
48!
49! 411 2009-12-11 12:31:43 Z suehring
50!
51! Description:
52! ------------
53! Advection scheme for scalars and momentum using the flux formulation of
54! Wicker and Skamarock 5th order. Additionally the module contains of a
55! routine using for initialisation and steering of the statical evaluation.
56! The computation of turbulent fluxes takes place inside the advection
57! routines.
58! In case of vector architectures Dirichlet and Radiation boundary conditions
59! are outstanding and not available.
60! A further routine local_diss_ij is available (next weeks) and is used if a
61! control of dissipative fluxes is desired.
62!
63! OUTSTANDING: - Dirichlet and Radiation boundary conditions for
64!                vector architectures
65!              - dissipation control for cache architectures ( next weeks )
66!              - Topography ( next weeks )
67!-----------------------------------------------------------------------------!
68
69    PRIVATE
70    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
71             local_diss, ws_init, ws_statistics
72
73    INTERFACE ws_init
74       MODULE PROCEDURE ws_init
75    END INTERFACE ws_init
76
77    INTERFACE ws_statistics
78       MODULE PROCEDURE ws_statistics
79    END INTERFACE ws_statistics
80
81    INTERFACE advec_s_ws
82       MODULE PROCEDURE advec_s_ws
83       MODULE PROCEDURE advec_s_ws_ij
84    END INTERFACE advec_s_ws
85
86    INTERFACE advec_u_ws
87       MODULE PROCEDURE advec_u_ws
88       MODULE PROCEDURE advec_u_ws_ij
89    END INTERFACE advec_u_ws
90
91    INTERFACE advec_v_ws
92       MODULE PROCEDURE advec_v_ws
93       MODULE PROCEDURE advec_v_ws_ij
94    END INTERFACE advec_v_ws
95
96    INTERFACE advec_w_ws
97       MODULE PROCEDURE advec_w_ws
98       MODULE PROCEDURE advec_w_ws_ij
99    END INTERFACE advec_w_ws
100
101    INTERFACE local_diss
102       MODULE PROCEDURE local_diss
103       MODULE PROCEDURE local_diss_ij
104    END INTERFACE local_diss
105
106 CONTAINS
107
108
109!------------------------------------------------------------------------------!
110! Initialization of WS-scheme
111!------------------------------------------------------------------------------!
112    SUBROUTINE ws_init
113
114       USE arrays_3d
115       USE constants
116       USE control_parameters
117       USE indices
118       USE pegrid
119       USE statistics
120
121!
122!--    Allocate arrays needed for dissipation control.
123       IF ( dissipation_control )  THEN
124!          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),        &
125!                   var_y(nzb+1:nzt,nys:nyn,nxl:nxr),        &
126!                   var_z(nzb+1:nzt,nys:nyn,nxl:nxr),        &
127!                   gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
128!                   gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
129!                   gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
130       ENDIF
131
132!
133!--    Set the appropriate factors for scalar and momentum advection.
134       adv_sca_5 = 1./60.
135       adv_sca_3 = 1./12.
136       adv_mom_5 = 1./120.
137       adv_mom_3 = 1./24.
138
139!         
140!--    Arrays needed for statical evaluation of fluxes.
141       IF ( ws_scheme_mom )  THEN
142
143          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1),  &
144                    sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1),  &
145                    sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1),   &
146                    sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1),   &
147                    sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
148
149          sums_wsus_ws_l = 0.0
150          sums_wsvs_ws_l = 0.0
151          sums_us2_ws_l  = 0.0
152          sums_vs2_ws_l  = 0.0
153          sums_ws2_ws_l  = 0.0
154
155       ENDIF
156
157       IF ( ws_scheme_sca )  THEN
158
159          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
160          sums_wspts_ws_l = 0.0
161
162          IF ( humidity .OR. passive_scalar )  THEN
163             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
164             sums_wsqs_ws_l = 0.0
165          ENDIF
166
167          IF ( ocean )  THEN
168             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
169             sums_wssas_ws_l = 0.0
170          ENDIF
171
172       ENDIF
173
174!
175!--    Arrays needed for reasons of speed optimization for cache and noopt
176!--    version. For the vector version the buffer arrays are not necessary,
177!--    because the the fluxes can swapped directly inside the loops of the
178!--    advection routines.
179       IF ( loop_optimization /= 'vector' )  THEN
180
181          IF ( ws_scheme_mom )  THEN
182
183             ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),            &
184                       flux_s_v(nzb+1:nzt,0:threads_per_task-1),            &
185                       flux_s_w(nzb+1:nzt,0:threads_per_task-1),            &
186                       diss_s_u(nzb+1:nzt,0:threads_per_task-1),            &
187                       diss_s_v(nzb+1:nzt,0:threads_per_task-1),            &
188                       diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
189             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
190                       flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
191                       flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
192                       diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
193                       diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
194                       diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
195
196          ENDIF
197
198          IF ( ws_scheme_sca )  THEN
199
200             ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
201                       flux_s_e(nzb+1:nzt,0:threads_per_task-1),            &
202                       diss_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
203                       diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) 
204             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
205                       flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
206                       diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
207                       diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
208
209             IF ( humidity .OR. passive_scalar )  THEN
210                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),         &
211                          diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
212                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
213                          diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
214             ENDIF
215
216             IF ( ocean )  THEN
217                ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1),         &
218                          diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
219                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
220                          diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
221             ENDIF
222
223          ENDIF
224
225       ENDIF
226
227!
228!--    Determine the flags where the order of the scheme for horizontal
229!--    advection has to be degraded.
230       ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
231       DO  i = nxl, nxr
232          DO  j = nys, nyn
233
234             boundary_flags(j,i) = 0
235             IF ( outflow_l )  THEN
236                IF ( i == nxlu  )  boundary_flags(j,i) = 5
237                IF ( i == nxl   )  boundary_flags(j,i) = 6
238             ELSEIF ( outflow_r )  THEN
239                IF ( i == nxr-1 )  boundary_flags(j,i) = 1
240                IF ( i == nxr   )  boundary_flags(j,i) = 2
241             ELSEIF ( outflow_n )  THEN
242                IF ( j == nyn-1 )  boundary_flags(j,i) = 3
243                IF ( j == nyn   )  boundary_flags(j,i) = 4
244             ELSEIF ( outflow_s )  THEN
245                IF ( j == nysv  )  boundary_flags(j,i) = 7
246                IF ( j == nys   )  boundary_flags(j,i) = 8
247             ENDIF
248
249          ENDDO
250       ENDDO
251       
252    END SUBROUTINE ws_init
253
254
255!------------------------------------------------------------------------------!
256! Initialize variables used for storing statistic qauntities (fluxes, variances)
257!------------------------------------------------------------------------------!
258    SUBROUTINE ws_statistics
259   
260       USE control_parameters
261       USE statistics
262
263       IMPLICIT NONE
264
265!       
266!--    The arrays needed for statistical evaluation are set to to 0 at the
267!--    begin of prognostic_equations.
268       IF ( ws_scheme_mom )  THEN
269          sums_wsus_ws_l = 0.0
270          sums_wsvs_ws_l = 0.0
271          sums_us2_ws_l  = 0.0
272          sums_vs2_ws_l  = 0.0
273          sums_ws2_ws_l  = 0.0
274       ENDIF
275
276       IF ( ws_scheme_sca )  THEN
277          sums_wspts_ws_l = 0.0
278          IF ( humidity .OR. passive_scalar )  sums_wsqs_ws_l = 0.0
279          IF ( ocean )  sums_wssas_ws_l = 0.0
280
281       ENDIF
282
283    END SUBROUTINE ws_statistics
284
285
286!------------------------------------------------------------------------------!
287! Scalar advection - Call for grid point i,j
288!------------------------------------------------------------------------------!
289    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
290                              swap_diss_y_local, swap_flux_x_local, &
291                              swap_diss_x_local, i_omp, tn )
292
293       USE arrays_3d
294       USE constants
295       USE control_parameters
296       USE grid_variables
297       USE indices
298       USE pegrid
299       USE statistics
300
301       IMPLICIT NONE
302
303       INTEGER ::  i, i_omp, j, k, tn
304       LOGICAL ::  degraded_l, degraded_s
305       REAL    ::  flux_d, diss_d, u_comp, v_comp
306       REAL, DIMENSION(:,:,:), POINTER    :: sk
307       REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
308                                             flux_n, diss_n
309       REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local,  &
310                                                          swap_diss_y_local
311       REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::             &
312                                                          swap_flux_x_local,  &
313                                                          swap_diss_x_local
314       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
315
316
317       degraded_l = .FALSE.
318       degraded_s = .FALSE.
319
320       IF ( boundary_flags(j,i) /= 0 )  THEN
321!       
322!--       Degrade the order for Dirichlet bc. at the outflow boundary
323          SELECT CASE ( boundary_flags(j,i) )
324
325             CASE ( 1 )
326
327                DO  k = nzb_s_inner(j,i)+1, nzt
328                   u_comp    = u(k,j,i+1) - u_gtrans
329                   flux_r(k) = u_comp * (                                      &
330                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
331                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
332                   diss_r(k) = -ABS( u_comp ) * (                              &
333                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
334                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
335                ENDDO
336
337             CASE ( 2 )
338
339                DO  k = nzb_s_inner(j,i)+1, nzt
340                   u_comp    = u(k,j,i+1) - u_gtrans
341                   flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
342                   diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
343                                         sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
344                                         0.5, ddx )
345                ENDDO
346
347             CASE ( 3 )
348
349                DO  k = nzb_s_inner(j,i)+1, nzt
350                   v_comp = v(k,j+1,i) - v_gtrans
351                   flux_n(k) = v_comp * (                                      &
352                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
353                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
354                   diss_n(k) = -ABS( v_comp ) * (                              &
355                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)     )           &
356                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
357               ENDDO
358
359             CASE ( 4 )
360
361                DO  k = nzb_s_inner(j,i)+1, nzt
362                   v_comp    = v(k,j+1,i) - v_gtrans
363                   flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
364                   diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
365                                         sk(k,j-1,i), sk(k,j-2,i), v_comp,    &
366                                         0.5, ddy )     
367                ENDDO
368
369             CASE ( 5 )
370
371                DO  k = nzb_w_inner(j,i)+1, nzt
372                   u_comp    = u(k,j,i+1) - u_gtrans 
373                   flux_r(k) = u_comp  * (                                     &
374                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
375                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
376                   diss_r(k) = -ABS( u_comp ) * (                              &
377                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
378                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
379                ENDDO 
380
381             CASE ( 6 )
382
383                DO  k = nzb_s_inner(j,i)+1, nzt
384                   u_comp    = u(k,j,i+1) - u_gtrans
385                   flux_r(k) = u_comp * (                                      &
386                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
387                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
388                   diss_r(k) = -ABS( u_comp ) * (                              &
389                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
390                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
391!                             
392!--                Compute leftside fluxes for the left boundary of PE domain
393                   u_comp                 = u(k,j,i) - u_gtrans
394                   swap_flux_x_local(k,j,tn) = u_comp * (                         &
395                                            sk(k,j,i) + sk(k,j,i-1) ) * 0.5
396                   swap_diss_x_local(k,j,tn) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
397                                                      sk(k,j,i), sk(k,j,i-1),  &
398                                                      sk(k,j,i-1), u_comp,     &
399                                                      0.5, ddx ) 
400                ENDDO
401                degraded_l = .TRUE.
402
403             CASE ( 7 )
404
405                DO  k = nzb_s_inner(j,i)+1, nzt
406                   v_comp    = v(k,j+1,i)-v_gtrans
407                   flux_n(k) = v_comp * (                                      &
408                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
409                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
410                   diss_n(k) = -ABS( v_comp ) * (                              &
411                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
412                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
413                ENDDO
414
415             CASE ( 8 )
416
417                DO  k = nzb_s_inner(j,i)+1, nzt
418                   v_comp    = v(k,j+1,i) - v_gtrans 
419                   flux_n(k) = v_comp * (                                      &
420                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
421                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
422                   diss_n(k) = -ABS( v_comp ) * (                              &
423                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
424                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
425!                             
426!--                Compute southside fluxes for the south boundary of PE domain
427                   v_comp               = v(k,j,i) - v_gtrans
428                   swap_flux_y_local(k,tn) = v_comp *                             &
429                                          ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 
430                   swap_diss_y_local(k,tn) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
431                                                    sk(k,j,i), sk(k,j-1,i),    &
432                                                    sk(k,j-1,i), v_comp,       &
433                                                    0.5, ddy )
434                ENDDO
435                degraded_s = .TRUE.
436               
437             CASE DEFAULT
438           
439          END SELECT
440
441!         
442!--       Compute the crosswise 5th order fluxes at the outflow
443          IF ( boundary_flags(j,i) == 1  .OR.  boundary_flags(j,i) == 2  .OR. &
444               boundary_flags(j,i) == 5  .OR.  boundary_flags(j,i) == 6 )  THEN
445
446             DO  k = nzb_s_inner(j,i)+1, nzt
447                v_comp    = v(k,j+1,i) - v_gtrans
448                flux_n(k) = v_comp * (                                         &
449                            37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )               &
450                          -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
451                          +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
452                diss_n(k) = -ABS( v_comp ) * (                                 &
453                            10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )               &
454                          -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
455                          +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
456             ENDDO
457           
458          ELSE
459         
460             DO  k = nzb_s_inner(j,i)+1, nzt
461                u_comp    = u(k,j,i+1) - u_gtrans 
462                flux_r(k) = u_comp * (                                         &
463                            37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )               &
464                          -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
465                          +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
466                diss_r(k) = -ABS( u_comp ) * (                                 &
467                            10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )               &
468                          -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
469                          +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
470             ENDDO
471           
472          ENDIF
473
474       ELSE
475               
476!
477!--       Compute the fifth order fluxes for the interior of PE domain.
478          DO  k = nzb_u_inner(j,i)+1, nzt
479             u_comp    = u(k,j,i+1) - u_gtrans
480             flux_r(k) = u_comp * (                                           &
481                         37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
482                       -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
483                       +        ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
484             diss_r(k) = -ABS( u_comp ) * (                                   &
485                         10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
486                       -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
487                       +        ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
488
489             v_comp    = v(k,j+1,i) - v_gtrans
490             flux_n(k) = v_comp * (                                           &
491                         37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
492                       -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
493                       +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
494             diss_n(k) = -ABS( v_comp ) * (                                   &
495                         10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
496                       -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
497                       +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
498          ENDDO
499         
500       ENDIF
501!       
502!--    Compute left- and southside fluxes of the respective PE bounds.       
503       IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
504       
505           DO  k = nzb_s_inner(j,i)+1, nzt
506              v_comp               = v(k,j,i) - v_gtrans
507              swap_flux_y_local(k,tn) = v_comp * (                               &
508                                     37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
509                                   -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
510                                   +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
511                                              ) * adv_sca_5
512              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                       &
513                                     10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
514                                   -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
515                                   +          sk(k,j+2,i) - sk(k,j-3,i)       &
516                                                      ) * adv_sca_5
517           ENDDO           
518         
519        ENDIF
520       
521        IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
522       
523           DO  k = nzb_s_inner(j,i)+1, nzt
524              u_comp                 = u(k,j,i) - u_gtrans
525              swap_flux_x_local(k,j,tn) = u_comp * (                             &
526                                       37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
527                                     -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
528                                     +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
529                                                ) * adv_sca_5
530              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                     &
531                                       10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
532                                     -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
533                                     +        ( sk(k,j,i+2) - sk(k,j,i-3) )   &
534                                                        ) * adv_sca_5
535           ENDDO
536           
537        ENDIF
538
539!       
540!--    Now compute the tendency terms for the horizontal parts
541       DO  k = nzb_s_inner(j,i)+1, nzt
542         
543          tend(k,j,i) = tend(k,j,i) - (                                      &
544                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
545                            swap_diss_x_local(k,j,tn) ) * ddx                   &
546                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
547                            swap_diss_y_local(k,tn)   ) * ddy                   &
548                                      )
549
550          swap_flux_y_local(k,tn)   = flux_n(k)
551          swap_diss_y_local(k,tn)   = diss_n(k)
552          swap_flux_x_local(k,j,tn) = flux_r(k)
553          swap_diss_x_local(k,j,tn) = diss_r(k)
554         
555       ENDDO
556
557!
558!--    Vertical advection, degradation of order near bottom and top.
559!--    The fluxes flux_d and diss_d at the surface are 0. Due to later
560!--    calculation of statistics the top flux at the surface should be 0.
561       flux_t(nzb_s_inner(j,i)) = 0.0
562       diss_t(nzb_s_inner(j,i)) = 0.0
563
564!       
565!--    2nd-order scheme (bottom)
566       k         = nzb_s_inner(j,i)+1
567       flux_d    = flux_t(k-1)
568       diss_d    = diss_t(k-1)
569       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
570
571!       
572!--    sk(k,j,i) is referenced three times to avoid an access below surface
573       diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i),  &
574                             sk(k,j,i), w(k,j,i), 0.5, ddzw(k) )
575                   
576       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
577                                   * ddzw(k)
578!
579!--    WS3 as an intermediate step (bottom)
580       k         = nzb_s_inner(j,i) + 2
581       flux_d    = flux_t(k-1)
582       diss_d    = diss_t(k-1)
583       flux_t(k) = w(k,j,i) * (                                               &
584                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )           &
585                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )           &
586                              ) * adv_sca_3
587       diss_t(k) = -ABS( w(k,j,i) ) * (                                       &
588                                3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )           & 
589                                    - ( sk(k+2,j,i) - sk(k-1,j,i) )           &
590                                      ) * adv_sca_3
591
592       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
593                                   * ddzw(k)
594!
595!--    WS5
596       DO  k = nzb_s_inner(j,i)+3, nzt-2
597       
598          flux_d    = flux_t(k-1)
599          diss_d    = diss_t(k-1)         
600          flux_t(k) = w(k,j,i) * (                                            &
601                                   37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )       &
602                                 -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )       &
603                                 +        ( sk(k+3,j,i) + sk(k-2,j,i) )       &
604                                 ) * adv_sca_5
605          diss_t(k) = -ABS( w(k,j,i) ) * (                                     &
606                                           10.0 * ( sk(k+1,j,i) - sk(k,j,i)   )&
607                                         -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )&
608                                         +        ( sk(k+3,j,i) - sk(k-2,j,i) )&
609                                         ) * adv_sca_5
610
611          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
612                                    - ( flux_d + diss_d ) ) * ddzw(k)
613
614       ENDDO
615
616!
617!--    WS3 as an intermediate step (top)
618       k         = nzt - 1
619       flux_d    = flux_t(k-1)
620       diss_d    = diss_t(k-1)
621       flux_t(k) = w(k,j,i) * (                                                &
622                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )            &
623                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )            &
624                              ) * adv_sca_3
625       diss_t(k) = -ABS( w(k,j,i) ) * (                                        &
626                                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )    &
627                                            - ( sk(k+2,j,i) - sk(k-1,j,i) )    &
628                                      ) * adv_sca_3
629
630       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
631                                   * ddzw(k)
632!                                 
633!--    2nd-order scheme (top)
634       k         = nzt 
635       flux_d    = flux_t(k-1)
636       diss_d    = diss_t(k-1)
637       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
638
639!       
640!--    sk(k+1) is referenced two times to avoid a segmentation fault at top
641       diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i), &
642                             sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) )
643
644       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
645                                   * ddzw(k)
646!       
647!--    Evaluation of statistics
648       SELECT CASE ( sk_char )
649
650          CASE ( 'pt' )
651
652             DO  k = nzb_s_inner(j,i), nzt
653                sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) +                      &
654                                       ( flux_t(k) + diss_t(k) )               &
655                                 * weight_substep(intermediate_timestep_count)
656             ENDDO
657             
658          CASE ( 'sa' )
659
660             DO  k = nzb_s_inner(j,i), nzt
661                sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +                      &
662                                       ( flux_t(k) + diss_t(k) )               &
663                                 * weight_substep(intermediate_timestep_count)
664             ENDDO
665             
666          CASE ( 'q' )
667
668             DO  k = nzb_s_inner(j,i), nzt
669                sums_wsqs_ws_l(k,tn)  = sums_wsqs_ws_l(k,tn) +                       &
670                                      ( flux_t(k) + diss_t(k) )                &
671                                 * weight_substep(intermediate_timestep_count)
672             ENDDO
673
674         END SELECT
675
676    END SUBROUTINE advec_s_ws_ij
677
678
679
680
681!------------------------------------------------------------------------------!
682! Advection of u-component - Call for grid point i,j
683!------------------------------------------------------------------------------!
684    SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
685
686       USE arrays_3d
687       USE constants
688       USE control_parameters
689       USE grid_variables
690       USE indices
691       USE statistics
692
693       IMPLICIT NONE
694
695       INTEGER ::  i, i_omp, j, k, tn
696       LOGICAL ::  degraded_l, degraded_s
697       REAL    ::  gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp
698       REAL, DIMENSION(nzb:nzt+1) ::  flux_t, diss_t, flux_r, diss_r,        &
699                                      flux_n, diss_n, u_comp
700                                       
701       degraded_l = .FALSE.
702       degraded_s = .FALSE.
703
704       gu = 2.0 * u_gtrans
705       gv = 2.0 * v_gtrans
706         
707       IF ( boundary_flags(j,i) /= 0 )  THEN
708!       
709!--      Degrade the order for Dirichlet bc. at the outflow boundary
710         SELECT CASE ( boundary_flags(j,i) )
711         
712            CASE ( 1 )
713               DO  k = nzb_u_inner(j,i)+1, nzt
714                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
715                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
716                              7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
717                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
718                  diss_r(k) = - ABS( u_comp(k) - gu ) * (                     &
719                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
720                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
721              ENDDO
722             
723            CASE ( 2 )
724               DO  k = nzb_u_inner(j,i)+1, nzt
725                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
726                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
727                                u(k,j,i+1) + u(k,j,i) ) * 0.25 
728                  diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i),     &
729                                        u(k,j,i-1), u(k,j,i-2), u_comp(k),    &
730                                        0.25, ddx )
731               ENDDO
732               
733            CASE ( 3 )
734               DO  k = nzb_u_inner(j,i)+1, nzt
735                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
736                  flux_n(k) = v_comp * (                                      &
737                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
738                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
739                  diss_n(k) = - ABS( v_comp ) * (                             &
740                              3.0 * ( u(k,j+1,i) - u(k,j,i)     )             &
741                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
742               ENDDO
743               
744            CASE ( 4 )
745               DO  k = nzb_u_inner(j,i)+1, nzt
746                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
747                  flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
748                  diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i), u(k,j,i),     &
749                                        u(k,j-1,i), u(k,j-2,i), v_comp,       &
750                                        0.25, ddy )
751               ENDDO
752               
753            CASE ( 5 )
754               DO  k = nzb_u_inner(j,i)+1, nzt
755!               
756!--               Compute leftside fluxes for the left boundary of PE domain
757                  u_comp(k)     = u(k,j,i) + u(k,j,i-1) - gu
758                  flux_l_u(k,j,tn) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
759                  diss_l_u(k,j,tn) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
760                                            u(k,j,i-1), u(k,j,i-1), u_comp(k),&
761                                            0.25, ddx )
762                 
763                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
764                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
765                              7.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
766                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
767                  diss_r(k) = - ABS( u_comp(k) -gu ) * (                      &
768                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
769                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
770               ENDDO
771               degraded_l = .TRUE.
772               
773            CASE ( 7 )
774               DO  k = nzb_u_inner(j,i)+1, nzt                           
775                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
776                  flux_n(k) = v_comp * (                                      &
777                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
778                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
779                  diss_n(k) = - ABS( v_comp ) * (                             &
780                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
781                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
782               ENDDO 
783               
784            CASE ( 8 )
785               DO  k = nzb_u_inner(j,i)+1, nzt
786!               
787!--              Compute southside fluxes for the south boundary of PE domain
788                  v_comp      = v(k,j,i) + v(k,j,i-1) - gv
789                  flux_s_u(k,tn) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25 
790                  diss_s_u(k,tn) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i),   &
791                                          u(k,j-1,i), u(k,j-1,i), v_comp,     &
792                                          0.25, ddy )
793                                 
794                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
795                  flux_n(k) = v_comp * (                                      &
796                              7.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
797                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
798                  diss_n(k) = - ABS( v_comp ) * (                             &
799                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
800                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
801               ENDDO 
802               degraded_s = .TRUE.
803               
804            CASE DEFAULT
805           
806         END SELECT
807!         
808!--      Compute the crosswise 5th order fluxes at the outflow
809         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
810              boundary_flags(j,i) == 5 )  THEN
811         
812             DO  k = nzb_u_inner(j,i)+1, nzt
813               v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
814               flux_n(k) = v_comp * (                                         &
815                           37.0 * ( u(k,j+1,i) + u(k,j,i)   )                 &
816                         -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                 &
817                         +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
818               diss_n(k) = - ABS ( v_comp ) * (                               &
819                           10.0 * ( u(k,j+1,i) - u(k,j,i)   )                 &
820                         -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                 &
821                         +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
822             ENDDO
823             
824         ELSE
825         
826            DO  k = nzb_u_inner(j,i)+1, nzt
827               u_comp(k) = u(k,j,i+1) + u(k,j,i)
828               flux_r(k) = ( u_comp(k) - gu ) * (                             &
829                           37.0 * ( u(k,j,i+1) + u(k,j,i)   )                 &
830                         -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                 &
831                         +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
832               diss_r(k) = - ABS( u_comp(k) - gu ) * (                        &
833                           10.0 * ( u(k,j,i+1) - u(k,j,i) )                   &
834                         -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                 &
835                         +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
836            ENDDO
837           
838         ENDIF
839                 
840       ELSE
841!       
842!--       Compute the fifth order fluxes for the interior of PE domain.
843          DO  k = nzb_u_inner(j,i)+1, nzt
844             u_comp(k) = u(k,j,i+1) + u(k,j,i)
845             flux_r(k) = ( u_comp(k) - gu ) * (                               &
846                         37.0 * ( u(k,j,i+1) + u(k,j,i)   )                   &
847                       -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )                   &
848                       +        ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
849             diss_r(k) = - ABS( u_comp(k) - gu ) * (                          &
850                         10.0 * ( u(k,j,i+1) - u(k,j,i)   )                   &
851                       -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )                   &
852                       +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
853
854             v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
855             flux_n(k) = v_comp * (                                           &
856                         37.0 * ( u(k,j+1,i) + u(k,j,i)   )                   &
857                       -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )                   &
858                       +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
859             diss_n(k) = - ABS( v_comp ) * (                                  &
860                         10.0 * ( u(k,j+1,i) - u(k,j,i)   )                   &
861                       -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )                   &
862                       +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
863          ENDDO 
864         
865       ENDIF
866!       
867!--    Compute left- and southside fluxes for the respective boundary of PE
868       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
869       
870          DO  k = nzb_u_inner(j,i)+1, nzt
871             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
872             flux_s_u(k,tn) = v_comp * (                                         &
873                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
874                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
875                         +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
876             diss_s_u(k,tn) = - ABS(v_comp) * (                                  &
877                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
878                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
879                         +        ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
880          ENDDO
881         
882       ENDIF
883       
884       IF ( i == i_omp .AND. ( .NOT. degraded_l ) )  THEN
885       
886          DO  k = nzb_u_inner(j,i)+1, nzt
887             u_comp_l      = u(k,j,i)+u(k,j,i-1)-gu
888             flux_l_u(k,j,tn) = u_comp_l * (                                     &
889                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )               &
890                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )               &
891                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
892             diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                              &
893                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )               &
894                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )               &
895                           +        ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
896          ENDDO
897         
898       ENDIF
899!       
900!--    Now compute the tendency terms for the horizontal parts.
901       DO  k = nzb_u_inner(j,i)+1, nzt                   
902          tend(k,j,i) = tend(k,j,i) - (                                       &
903                            ( flux_r(k) + diss_r(k)                           &
904                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx           &
905                          + ( flux_n(k) + diss_n(k)                           &
906                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy  )
907
908           flux_l_u(k,j,tn) = flux_r(k)
909           diss_l_u(k,j,tn) = diss_r(k)
910           flux_s_u(k,tn)   = flux_n(k)
911           diss_s_u(k,tn)   = diss_n(k)
912!
913!--        Statistical Evaluation of u'u'. The factor has to be applied for
914!--        right evaluation when gallilei_trans = .T. .
915           sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn)                                &
916             + ( flux_r(k) *                                                  &
917               ( u_comp(k) - 2.0 * hom(k,1,1,0) )                             &
918             / ( u_comp(k) - gu + 1.0E-20      )                              &
919             +   diss_r(k) *                                                  &
920                 ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )                        &
921             / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )                          &
922             *   weight_substep(intermediate_timestep_count)
923       ENDDO 
924       sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
925                                           
926
927!
928!--    Vertical advection, degradation of order near surface and top.
929!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
930!--    statistical evaluation the top flux at the surface should be 0
931       flux_t(nzb_u_inner(j,i)) = 0. !statistical reasons
932       diss_t(nzb_u_inner(j,i)) = 0.
933!
934!--    2nd order scheme (bottom)
935       k         = nzb_u_inner(j,i)+1
936       flux_d    = flux_t(k-1)
937       diss_d    = diss_t(k-1)
938       w_comp    = w(k,j,i) + w(k,j,i-1)
939       flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) *0.25
940       diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), 0., 0.,        &
941                             w_comp, 0.25, ddzw(k) )
942             
943       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
944                                 -   flux_d    - diss_d      ) * ddzw(k)
945!
946!--    WS3 as an intermediate step (bottom)
947       k         = nzb_u_inner(j,i)+2
948       flux_d    = flux_t(k-1)
949       diss_d    = diss_t(k-1)
950       w_comp    = w(k,j,i) + w(k,j,i-1)
951       flux_t(k) = w_comp * (                                                 &
952                   7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
953                 -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
954       diss_t(k) = -ABS( w_comp) * (                                          &
955                   3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          & 
956                 -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
957
958       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
959                                 -   flux_d    - diss_d      ) * ddzw(k)
960!
961!--    WS5
962       DO  k = nzb_u_inner(j,i)+3, nzt-2
963          flux_d    = flux_t(k-1)
964          diss_d    = diss_t(k-1)
965          w_comp    = w(k,j,i) + w(k,j,i-1)
966          flux_t(k) = w_comp * (                                              &
967                      37.0 * ( u(k+1,j,i) + u(k,j,i)   )                      &
968                    -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                      &
969                    +        ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
970          diss_t(k) = - ABS( w_comp ) * (                                     &
971                      10.0 *  ( u(k+1,j,i) - u(k,j,i)   )                     & 
972                    -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                      &
973                    +        ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
974
975          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
976                                    -   flux_d    - diss_d      ) * ddzw(k)
977       ENDDO
978!
979!--    WS3 as an intermediate step (top)
980       k         = nzt - 1
981       flux_d    = flux_t(k-1)
982       diss_d    = diss_t(k-1)
983       w_comp    = w(k,j,i) + w(k,j,i-1)
984       flux_t(k) = w_comp * (                                                 &
985                   7.0 * ( u(k+1,j,i) + u(k,j,i)   )                          &
986                 -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
987       diss_t(k) = - ABS( w_comp ) * (                                        &
988                   3.0 * ( u(k+1,j,i) - u(k,j,i)   )                          & 
989                 -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
990                   
991       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
992                                 -   flux_d    - diss_d      ) * ddzw(k)
993       
994!
995!--    2nd order scheme (top)
996       k         = nzt
997       flux_d    = flux_t(k-1)
998       diss_d    = diss_t(k-1)
999       w_comp    = w(k,j,i)+w(k,j,i-1)
1000       flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
1001       diss_t(k) = diss_2nd( u(k+1,j,i), u(k+1,j,i), u(k,j,i), u(k-1,j,i),    &
1002                             u(k-2,j,i), w_comp, 0.25, ddzw(k) )
1003
1004       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1005                                 -   flux_d    - diss_d      ) * ddzw(k)
1006
1007!
1008!--    sum up the vertical momentum fluxes
1009       DO  k = nzb_u_inner(j,i), nzt
1010          sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                               &
1011              + ( flux_t(k) + diss_t(k) )                                     &
1012              *   weight_substep(intermediate_timestep_count)
1013       ENDDO
1014
1015
1016    END SUBROUTINE advec_u_ws_ij
1017
1018
1019
1020
1021!-----------------------------------------------------------------------------!
1022! Advection of v-component - Call for grid point i,j
1023!-----------------------------------------------------------------------------!
1024   SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
1025
1026       USE arrays_3d
1027       USE constants
1028       USE control_parameters
1029       USE grid_variables
1030       USE indices
1031       USE statistics
1032
1033       IMPLICIT NONE
1034
1035       INTEGER  ::  i, i_omp, j, k, tn
1036       LOGICAL  ::  degraded_l, degraded_s
1037       REAL     ::  gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp
1038       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_n,                &
1039                                       diss_n, flux_r, diss_r, v_comp
1040                                       
1041       degraded_l = .FALSE.
1042       degraded_s = .FALSE.
1043     
1044       gu = 2.0 * u_gtrans
1045       gv = 2.0 * v_gtrans
1046       
1047       IF ( boundary_flags(j,i) /= 0 )  THEN
1048!       
1049!--       Degrade the order for Dirichlet bc. at the outflow boundary
1050          SELECT CASE ( boundary_flags(j,i) )
1051         
1052             CASE ( 1 )
1053                DO  k = nzb_v_inner(j,i)+1, nzt
1054                   u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1055                   flux_r(k) = u_comp * (                                     &
1056                               7.0 * (v(k,j,i+1) + v(k,j,i)    )              &
1057                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
1058                   diss_r(k) = - ABS( u_comp ) * (                            &
1059                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
1060                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
1061                ENDDO
1062               
1063             CASE ( 2 )
1064                DO  k = nzb_v_inner(j,i)+1, nzt
1065                   u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1066                   flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 
1067                   diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1), v(k,j,i),    &
1068                                         v(k,j,i-1), v(k,j,i-2), u_comp,      &
1069                                         0.25, ddx )
1070                ENDDO
1071               
1072             CASE ( 3 )
1073                DO  k = nzb_v_inner(j,i)+1, nzt
1074                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
1075                   flux_n(k) = ( v_comp(k)- gv ) * (                          &
1076                               7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
1077                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
1078                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
1079                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
1080                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
1081                ENDDO
1082               
1083             CASE ( 4 )
1084                DO  k = nzb_v_inner(j,i)+1, nzt
1085                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
1086                   flux_n(k) = ( v_comp(k) - gv ) *                           &
1087                               ( v(k,j+1,i) + v(k,j,i) ) * 0.25 
1088                   diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i), v(k,j,i),    &
1089                                         v(k,j-1,i), v(k,j-2,i), v_comp(k),   &
1090                                         0.25, ddy )
1091                ENDDO
1092               
1093             CASE ( 5 )
1094                DO  k = nzb_w_inner(j,i)+1, nzt
1095                   u_comp    = u(k,j-1,i) + u(k,j,i) - gu 
1096                   flux_r(k) = u_comp * (                                     &
1097                               7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
1098                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
1099                   diss_r(k) = - ABS( u_comp ) * (                            &
1100                               3.0 * ( w(k,j,i+1) - w(k,j,i)   )              &
1101                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
1102                ENDDO
1103               
1104             CASE ( 6 )
1105                DO  k = nzb_v_inner(j,i)+1, nzt
1106                   u_comp        = u(k,j-1,i) + u(k,j,i) - gu
1107                   flux_l_v(k,j,tn) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
1108                   diss_l_v(k,j,tn) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
1109                                             v(k,j,i-1), v(k,j,i-1), u_comp,  &
1110                                             0.25, ddx )
1111                                 
1112                   u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1113                   flux_r(k) = u_comp * (                                     &
1114                               7.0 * ( v(k,j,i+1) + v(k,j,i)   )              &
1115                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
1116                   diss_r(k) = - ABS( u_comp ) * (                            &
1117                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
1118                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 
1119                ENDDO
1120                degraded_l = .TRUE.
1121               
1122             CASE ( 7 )
1123                DO  k = nzb_v_inner(j,i)+1, nzt
1124                   v_comp(k)   = v(k,j,i) + v(k,j-1,i) - gv
1125                   flux_s_v(k,tn) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
1126                   diss_s_v(k,tn) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i),  &
1127                                           v(k,j-1,i), v(k,j-1,i), v_comp(k), &
1128                                           0.25, ddy )
1129                                 
1130                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
1131                   flux_n(k) = ( v_comp(k) - gv ) * (                         &
1132                               7.0 * ( v(k,j+1,i) + v(k,j,i)   )              &
1133                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
1134                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
1135                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
1136                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
1137                ENDDO
1138                degraded_s = .TRUE.
1139               
1140              CASE DEFAULT
1141             
1142          END SELECT
1143!         
1144!--       Compute the crosswise 5th order fluxes at the outflow
1145          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.   &
1146               boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
1147         
1148             DO  k = nzb_v_inner(j,i)+1, nzt
1149                v_comp(k) = v(k,j+1,i) + v(k,j,i)
1150                flux_n(k) = ( v_comp(k) - gv ) * (                            &
1151                            37.0 * ( v(k,j+1,i) + v(k,j,i)   )                &
1152                          -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                &
1153                          +        ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1154                diss_n(k) = - ABS( v_comp(k) - gv ) * (                       &
1155                             10.0 * ( v(k,j+1,i) - v(k,j,i)   )               &
1156                           -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )               &
1157                           +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
1158              ENDDO
1159             
1160          ELSE
1161         
1162             DO  k = nzb_v_inner(j,i)+1, nzt
1163                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1164                flux_r(k) = u_comp * (                                        &
1165                            37.0 * ( v(k,j,i+1) + v(k,j,i)   )                &
1166                          -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                &
1167                          +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1168                diss_r(k) = - ABS( u_comp ) * (                               &
1169                            10.0 * ( v(k,j,i+1) - v(k,j,i)   )                &
1170                          -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                &
1171                          +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1172             ENDDO
1173             
1174          ENDIF
1175         
1176       ELSE
1177!       
1178!--       Compute the fifth order fluxes for the interior of PE domain.                 
1179          DO  k = nzb_v_inner(j,i)+1, nzt
1180             u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1181             flux_r(k) = u_comp * (                                           &
1182                         37.0 * ( v(k,j,i+1) + v(k,j,i)   )                   &
1183                       -  8.0 * ( v(k,j,i+2) + v(k,j,i-1) )                   &
1184                       +        ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1185             diss_r(k) = - ABS( u_comp ) * (                                  &
1186                         10.0 * ( v(k,j,i+1) - v(k,j,i) )                     &
1187                       -  5.0 * ( v(k,j,i+2) - v(k,j,i-1) )                   &
1188                       +        ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1189
1190             v_comp(k) = v(k,j+1,i) + v(k,j,i)
1191             flux_n(k) = ( v_comp(k) - gv ) * (                               &
1192                         37.0 * ( v(k,j+1,i) + v(k,j,i)   )                   &
1193                       -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )                   &
1194                         +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1195             diss_n(k) = - ABS( v_comp(k) - gv ) * (                          &
1196                         10.0 * ( v(k,j+1,i) - v(k,j,i)   )                   &
1197                       -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )                   &
1198                       +        ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
1199          ENDDO 
1200         
1201       ENDIF
1202!       
1203!--    Compute left- and southside fluxes for the respective boundary       
1204       IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
1205
1206          DO  k = nzb_v_inner(j,i)+1, nzt
1207             u_comp        = u(k,j-1,i) + u(k,j,i) - gu
1208             flux_l_v(k,j,tn) = u_comp * (                                       &
1209                             37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
1210                           -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
1211                           +        ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
1212             diss_l_v(k,j,tn) = - ABS( u_comp ) * (                              &
1213                             10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
1214                           -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
1215                           +        ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
1216          ENDDO
1217         
1218       ENDIF
1219       
1220       IF ( j == nysv  .AND.  .NOT. degraded_s )  THEN
1221       
1222          DO  k = nzb_v_inner(j,i)+1, nzt
1223             v_comp_l    = v(k,j,i) + v(k,j-1,i) - gv
1224             flux_s_v(k,tn) = v_comp_l * (                                       &
1225                           37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
1226                         -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
1227                         +        ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
1228             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                              &
1229                           10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
1230                         -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
1231                         +        ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
1232          ENDDO
1233         
1234       ENDIF
1235!       
1236!--    Now compute the tendency terms for the horizontal parts.         
1237       DO  k = nzb_v_inner(j,i)+1, nzt                 
1238          tend(k,j,i) = tend(k,j,i) - (                                       &
1239                         ( flux_r(k) + diss_r(k)                              &
1240                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx            &
1241                       + ( flux_n(k) + diss_n(k)                              &
1242                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy     )
1243       
1244           flux_l_v(k,j,tn) = flux_r(k)
1245           diss_l_v(k,j,tn) = diss_r(k)
1246           flux_s_v(k,tn)   = flux_n(k)
1247           diss_s_v(k,tn)   = diss_n(k)
1248
1249!
1250!--        Statistical Evaluation of v'v'. The factor has to be applied for
1251!--        right evaluation when gallilei_trans = .T. .
1252
1253           sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn)                                &
1254             + ( flux_n(k)                                                    &
1255             * ( v_comp(k) - 2.0 * hom(k,1,2,0) )                             &
1256             / ( v_comp(k) - gv + 1.0E-20 )                                   &
1257             +   diss_n(k)                                                    &
1258             *   ABS( v_comp(k) - 2.0 * hom(k,1,2,0) )                        &
1259             / ( ABS( v_comp(k) - gv ) +1.0E-20 ) )                           &
1260             *   weight_substep(intermediate_timestep_count)
1261
1262       ENDDO
1263       sums_vs2_ws_l(nzb_v_inner(j,i),tn) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,tn) 
1264                                           
1265!
1266!--    Vertical advection, degradation of order near surface and top.
1267!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1268!--    statistical evaluation the top flux at the surface should be 0
1269       flux_t(nzb_v_inner(j,i)) = 0.0  !statistical reasons
1270       diss_t(nzb_v_inner(j,i)) = 0.0
1271!
1272!--    2nd order scheme (bottom)     
1273       k         = nzb_v_inner(j,i)+1
1274       flux_d    = flux_t(k-1)
1275       diss_d    = diss_t(k-1) 
1276       w_comp    = w(k,j-1,i) + w(k,j,i)
1277       flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
1278       diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0.0, 0.0,      &
1279                             w_comp, 0.25, ddzw(k) ) 
1280
1281       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1282                                 -   flux_d    - diss_d       ) * ddzw(k)
1283       
1284!
1285!--    WS3 as an intermediate step (bottom)
1286       k         = nzb_v_inner(j,i)+2
1287       flux_d    = flux_t(k-1)
1288       diss_d    = diss_t(k-1)
1289       w_comp    = w(k,j-1,i) + w(k,j,i)
1290       flux_t(k) = w_comp * (                                                 &
1291                   7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &
1292                 -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
1293       diss_t(k) = - ABS( w_comp ) * (                                        &
1294                   3.0 * ( v(k+1,j,i) - v(k,j,i)   )                          & 
1295                 -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
1296
1297       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1298                                 -   flux_d    - diss_d       ) * ddzw(k)
1299!                                 
1300!--    WS5
1301       DO  k = nzb_v_inner(j,i)+3, nzt-2
1302          flux_d    = flux_t(k-1)
1303          diss_d    = diss_t(k-1)
1304          w_comp    = w(k,j-1,i) + w(k,j,i)
1305          flux_t(k) = w_comp * (                                              &
1306                      37.0 * ( v(k+1,j,i) + v(k,j,i)   )                      &           
1307                    -  8.0 * ( v(k+2,j,i) + v(k-1,j,i) )                      &
1308                    +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
1309          diss_t(k) = - ABS( w_comp ) * (                                     &
1310                      10.0 * ( v(k+1,j,i) - v(k,j,i)   )                      & 
1311                    -  5.0 * ( v(k+2,j,i) - v(k-1,j,i) )                      &
1312                    +        ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
1313
1314          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
1315                                 -   flux_d    - diss_d       ) * ddzw(k)
1316       ENDDO
1317!
1318!--    WS3 as an intermediate step (top)
1319       k         = nzt - 1
1320       flux_d    = flux_t(k-1)
1321       diss_d    = diss_t(k-1)
1322       w_comp    = w(k,j-1,i) + w(k,j,i)
1323       flux_t(k) = w_comp * (                                                 &
1324                   7.0 * ( v(k+1,j,i) + v(k,j,i)   )                          &   
1325                 -       ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
1326       diss_t(k) = - ABS( w_comp ) * (                                        &
1327                   3.0 * ( v(k+1,j,i) - v(k,j,i) )                            & 
1328                 -       ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
1329       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1330                                 -   flux_d    - diss_d       ) * ddzw(k)
1331!
1332!--    2nd order scheme (top)
1333       k         = nzt
1334       flux_d    = flux_t(k-1)
1335       diss_d    = diss_t(k-1)
1336       w_comp    = w(k,j-1,i)+w(k,j,i)
1337       flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
1338       diss_t(k) = diss_2nd( v(k+1,j,i), v(k+1,j,i), v(k,j,i), v(k-1,j,i),    &
1339                             v(k-2,j,i), w_comp, 0.25, ddzw(k) )
1340 
1341       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1342                                 -   flux_d    - diss_d       ) * ddzw(k)
1343             
1344       DO  k = nzb_v_inner(j,i), nzt
1345          sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn)                               &
1346                 + ( flux_t(k) + diss_t(k) )                                  &
1347                 *   weight_substep(intermediate_timestep_count) 
1348       ENDDO
1349
1350    END SUBROUTINE advec_v_ws_ij
1351
1352
1353
1354!------------------------------------------------------------------------------!
1355! Advection of w-component - Call for grid point i,j
1356!------------------------------------------------------------------------------!
1357    SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
1358
1359       USE arrays_3d
1360       USE constants
1361       USE control_parameters
1362       USE grid_variables
1363       USE indices
1364       USE statistics
1365
1366       IMPLICIT NONE
1367
1368       INTEGER ::  i, i_omp, j, k, tn
1369       LOGICAL ::  degraded_l, degraded_s
1370       REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
1371       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_r, diss_r, flux_n, &
1372                                       diss_n
1373                                       
1374       degraded_l = .FALSE. 
1375       degraded_s = .FALSE.
1376
1377       gu = 2.0 * u_gtrans
1378       gv = 2.0 * v_gtrans
1379
1380             
1381       IF ( boundary_flags(j,i) /= 0 )  THEN
1382!       
1383!--      Degrade the order for Dirichlet bc. at the outflow boundary
1384         SELECT CASE ( boundary_flags(j,i) )
1385         
1386            CASE ( 1 )
1387               DO  k = nzb_w_inner(j,i)+1, nzt
1388                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1389                  flux_r(k) = u_comp * (                                      &
1390                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
1391                            -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
1392                  diss_r(k) = -ABS( u_comp ) * (                              &
1393                              3.0 * ( w(k,j,i+1) - w(k,j,i)   )               & 
1394                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 
1395               ENDDO
1396               
1397            CASE ( 2 )
1398               DO  k = nzb_w_inner(j,i)+1, nzt
1399                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1400                  flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 
1401                  diss_r(k) = diss_2nd( w(k,j,i+1), w(k,j,i+1), w(k,j,i),     &
1402                                        w(k,j,i-1), w(k,j,i-2), u_comp,       &
1403                                        0.25, ddx )
1404               ENDDO
1405               
1406            CASE ( 3 )
1407               DO  k = nzb_w_inner(j,i)+1, nzt
1408                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1409                  flux_n(k) = v_comp * (                                      &
1410                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
1411                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1412                  diss_n(k) = -ABS( v_comp ) * (                              &
1413                              3.0 * ( w(k,j+1,i) - w(k,j,i)   )               & 
1414                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1415               ENDDO
1416               
1417            CASE ( 4 )
1418               DO  k = nzb_w_inner(j,i)+1, nzt
1419                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1420                  flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 
1421                  diss_n(k) = diss_2nd( w(k,j+1,i), w(k,j+1,i), w(k,j,i),     &
1422                                        w(k,j-1,i), w(k,j-2,i), v_comp,       &
1423                                        0.25, ddy )
1424               ENDDO
1425               
1426            CASE ( 5 )
1427               DO  k = nzb_w_inner(j,i)+1, nzt
1428                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu 
1429                  flux_r(k) = u_comp * (                                      &
1430                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
1431                            -       ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
1432                  diss_r(k) = - ABS( u_comp ) * (                             &
1433                              3.0 * ( w(k,j,i+1) - w(k,j,i) )                 & 
1434                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
1435               ENDDO
1436               
1437            CASE ( 6 )
1438               DO  k = nzb_w_inner(j,i)+1, nzt
1439!               
1440!--               Compute leftside fluxes for the left boundary of PE domain
1441                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1442                  flux_r(k) = u_comp *(                                       &
1443                              7.0 * ( w(k,j,i+1) + w(k,j,i)   )               &
1444                            -       ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
1445                  diss_r(k) = - ABS( u_comp ) * (                             &
1446                              3.0 * ( w(k,j,i+1) - w(k,j,i) )                 & 
1447                            -       ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
1448                 
1449                  u_comp        = u(k+1,j,i) + u(k,j,i) - gu
1450                  flux_l_w(k,j,tn) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
1451                  diss_l_w(k,j,tn) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
1452                                            w(k,j,i-1), w(k,j,i-1), u_comp,   &
1453                                            0.25, ddx )
1454               ENDDO
1455               degraded_l = .TRUE.
1456               
1457            CASE ( 7 )
1458               DO  k = nzb_w_inner(j,i)+1, nzt
1459                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1460                  flux_n(k) = v_comp *(                                       &
1461                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
1462                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1463                  diss_n(k) = - ABS( v_comp ) * (                             &
1464                              3.0 * ( w(k,j+1,i) - w(k,j,i)   )               & 
1465                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1466                ENDDO
1467               
1468            CASE ( 8 )
1469               DO  k = nzb_w_inner(j,i)+1, nzt
1470                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1471                  flux_n(k) = v_comp * (                                      &
1472                              7.0 * ( w(k,j+1,i) + w(k,j,i)   )               &
1473                            -       ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1474                  diss_n(k) = - ABS( v_comp ) * (                             &
1475                              3.0 * ( w(k,j+1,i) - w(k,j,i) )                 & 
1476                            -       ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1477!                   
1478!--              Compute southside fluxes for the south boundary of PE domain           
1479                  v_comp      = v(k+1,j,i) + v(k,j,i) - gv
1480                  flux_s_w(k,tn) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
1481                  diss_s_w(k,tn) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i),   &
1482                                          w(k,j-1,i), w(k,j-1,i), v_comp,     &
1483                                          0.25, ddy )                 
1484               ENDDO 
1485               degraded_s = .TRUE.
1486               
1487            CASE DEFAULT
1488           
1489         END SELECT
1490!         
1491!--      Compute the crosswise 5th order fluxes at the outflow
1492         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2  .OR.    &
1493              boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
1494         
1495            DO  k = nzb_w_inner(j,i)+1, nzt
1496               v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1497               flux_n(k) = v_comp * (                                         &
1498                           37.0 * ( w(k,j+1,i) + w(k,j,i)   )                 &
1499                         -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                 &
1500                         +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1501               diss_n(k) = - ABS( v_comp ) * (                                &
1502                           10.0 * ( w(k,j+1,i) - w(k,j,i)   )                 &
1503                         -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                 &
1504                         +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1505            ENDDO
1506           
1507         ELSE
1508         
1509            DO  k = nzb_w_inner(j,i)+1, nzt         
1510               u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1511               flux_r(k) = u_comp * (                                         &
1512                           37.0 * ( w(k,j,i+1) + w(k,j,i) )                   &
1513                         -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                 &
1514                         +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1515               diss_r(k) = - ABS( u_comp ) * (                                &
1516                           10.0 * ( w(k,j,i+1) - w(k,j,i)   )                 &
1517                         -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                 &
1518                         +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1519            ENDDO
1520           
1521         ENDIF
1522         
1523       ELSE
1524!       
1525!--      Compute the fifth order fluxes for the interior of PE domain.               
1526         DO  k = nzb_w_inner(j,i)+1, nzt
1527            u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1528            flux_r(k) = u_comp * (                                            &
1529                        37.0 * ( w(k,j,i+1) + w(k,j,i)   )                    &
1530                      -  8.0 * ( w(k,j,i+2) + w(k,j,i-1) )                    &
1531                      +        ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1532            diss_r(k) = - ABS( u_comp ) * (                                   &
1533                        10.0 * ( w(k,j,i+1) - w(k,j,i)   )                    &
1534                      -  5.0 * ( w(k,j,i+2) - w(k,j,i-1) )                    &
1535                      +        ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1536                 
1537            v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1538            flux_n(k) = v_comp * (                                            &
1539                        37.0 * ( w(k,j+1,i) + w(k,j,i)   )                    &
1540                      -  8.0 * ( w(k,j+2,i) + w(k,j-1,i) )                    &
1541                      +        ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1542            diss_n(k) = - ABS( v_comp ) * (                                   &
1543                        10.0 * ( w(k,j+1,i) - w(k,j,i)   )                    &
1544                      -  5.0 * ( w(k,j+2,i) - w(k,j-1,i) )                    &
1545                      +        ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1546         ENDDO 
1547         
1548       ENDIF
1549!       
1550!--    Compute left- and southside fluxes for the respective boundary     
1551       IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
1552       
1553          DO  k = nzb_w_inner(j,i)+1, nzt
1554             v_comp      = v(k+1,j,i) + v(k,j,i) - gv
1555             flux_s_w(k,tn) = v_comp * (                                         &
1556                           37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
1557                         -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
1558                         +        ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
1559             diss_s_w(k,tn) = - ABS( v_comp ) * (                                &
1560                           10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
1561                         -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
1562                         +        ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
1563          ENDDO
1564         
1565       ENDIF
1566       
1567       IF ( i == i_omp  .AND.   .NOT. degraded_l ) THEN
1568         
1569          DO  k = nzb_w_inner(j,i)+1, nzt
1570            u_comp        = u(k+1,j,i) + u(k,j,i) - gu
1571            flux_l_w(k,j,tn) = u_comp * (                                        &
1572                            37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
1573                          -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                &
1574                          +        ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
1575            diss_l_w(k,j,tn) = - ABS( u_comp ) * (                               &
1576                            10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
1577                          -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
1578                          +        ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 
1579          ENDDO
1580         
1581       ENDIF
1582
1583!       
1584!--    Now compute the tendency terms for the horizontal parts.         
1585       DO  k = nzb_w_inner(j,i)+1, nzt
1586          tend(k,j,i) = tend(k,j,i) - (                                       &
1587                      ( flux_r(k) + diss_r(k)                                 &
1588                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx               &
1589                    + ( flux_n(k) + diss_n(k)                                 &
1590                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy  )
1591
1592          flux_l_w(k,j,tn) = flux_r(k)
1593          diss_l_w(k,j,tn) = diss_r(k)
1594          flux_s_w(k,tn) = flux_n(k)
1595          diss_s_w(k,tn) = diss_n(k)
1596       ENDDO
1597
1598!
1599!--    Vertical advection, degradation of order near surface and top.
1600!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1601!--    statistical evaluation the top flux at the surface should be 0
1602       flux_t(nzb_w_inner(j,i)) = 0.0  !statistical reasons
1603       diss_t(nzb_w_inner(j,i)) = 0.0
1604!
1605!--    2nd order scheme (bottom)       
1606       k         = nzb_w_inner(j,i)+1
1607       flux_d    = flux_t(k-1)
1608       diss_d    = diss_t(k-1)
1609       w_comp    = w(k+1,j,i) + w(k,j,i)
1610       flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
1611       diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0,        &
1612                             w_comp, 0.25, ddzu(k+1) )
1613                     
1614       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1615                                 -   flux_d    - diss_d      ) * ddzu(k+1)               
1616!
1617!--    WS3 as an intermediate step (bottom)
1618       k         = nzb_w_inner(j,i)+2
1619       flux_d    = flux_t(k-1)
1620       diss_d    = diss_t(k-1)
1621       w_comp    = w(k+1,j,i) + w(k,j,i)
1622       flux_t(k) = w_comp * (                                                 &
1623                   7.0 * ( w(k+1,j,i) + w(k,j,i) )                            &
1624                 -       ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
1625       diss_t(k) = - ABS( w_comp ) * (                                        &
1626                   3.0 * ( w(k+1,j,i) - w(k,j,i) )                            & 
1627                 -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
1628
1629       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1630                                 -   flux_d    - diss_d      ) * ddzu(k+1)
1631!                                 
1632!--    WS5
1633       DO  k = nzb_w_inner(j,i)+3, nzt-2 
1634          flux_d    = flux_t(k-1)
1635          diss_d    = diss_t(k-1)
1636          w_comp    = w(k+1,j,i) + w(k,j,i)
1637          flux_t(k) = w_comp * (                                              &
1638                      37.0 * ( w(k+1,j,i) + w(k,j,i)   )                      &
1639                    -  8.0 * ( w(k+2,j,i) + w(k-1,j,i) )                      &
1640                    +        ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
1641          diss_t(k) = - ABS( w_comp ) * (                                     &
1642                      10.0 * ( w(k+1,j,i) - w(k,j,i)   )                      & 
1643                    -  5.0 * ( w(k+2,j,i) - w(k-1,j,i) )                      & 
1644                    +        ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
1645
1646          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
1647                                    -   flux_d    - diss_d      ) * ddzu(k+1)
1648       ENDDO
1649!--    WS3 as an intermediate step (top)
1650       k         = nzt - 1
1651       flux_d    = flux_t(k-1)
1652       diss_d    = diss_t(k-1)
1653       w_comp    = w(k+1,j,i) + w(k,j,i)
1654       flux_t(k) = w_comp * (                                                 &
1655                   7.0 * ( w(k+1,j,i) + w(k,j,i)   )                          &
1656                 -       ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
1657       diss_t(k) = - ABS( w_comp ) * (                                        &
1658                   3.0 * ( w(k+1,j,i) - w(k,j,i)     )                        & 
1659                 -       ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
1660                   
1661       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1662                                 -   flux_d    - diss_d      ) * ddzu(k+1)
1663!
1664!--    2nd order scheme (top)
1665       k         = nzt
1666       flux_d    = flux_t(k-1)
1667       diss_d    = diss_t(k-1)
1668       w_comp    = w(k+1,j,i) + w(k,j,i)
1669       flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
1670       diss_t(k) = diss_2nd( w(k+1,j,i), w(k+1,j,i), w(k,j,i), w(k-1,j,i),    &
1671                             w(k-2,j,i), w_comp, 0.25, ddzu(k+1) )
1672
1673       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1674                                 -   flux_d    - diss_d      ) * ddzu(k+1)
1675       
1676       DO  k = nzb_w_inner(j,i), nzt
1677           sums_ws2_ws_l(k,tn)  = sums_ws2_ws_l(k,tn)                               &
1678                 + ( flux_t(k) + diss_t(k) )                                  &
1679                 *   weight_substep(intermediate_timestep_count)
1680       ENDDO
1681
1682    END SUBROUTINE advec_w_ws_ij
1683   
1684
1685!------------------------------------------------------------------------------!
1686! Scalar advection - Call for all grid points
1687!------------------------------------------------------------------------------!
1688    SUBROUTINE advec_s_ws( sk, sk_char )
1689
1690       USE arrays_3d
1691       USE constants
1692       USE control_parameters
1693       USE grid_variables
1694       USE indices
1695       USE statistics
1696
1697       IMPLICIT NONE
1698
1699       INTEGER ::  i, j, k, tn = 0
1700       REAL, DIMENSION(:,:,:), POINTER ::  sk
1701       REAL    :: flux_d, diss_d, u_comp, v_comp
1702       REAL, DIMENSION(nzb:nzt+1)  ::  flux_r, diss_r, flux_n, diss_n
1703       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, swap_diss_y_local,    &
1704                                     flux_t, diss_t
1705       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
1706                                             swap_diss_x_local
1707       CHARACTER (LEN = *), INTENT(IN) :: sk_char
1708
1709!
1710!--    Compute the fluxes for the whole left boundary of the processor domain.
1711       i = nxl
1712       DO  j = nys, nyn
1713          IF ( boundary_flags(j,i) == 6 )  THEN
1714         
1715             DO  k = nzb_s_inner(j,i)+1, nzt
1716                u_comp                 = u(k,j,i) - u_gtrans
1717                swap_flux_x_local(k,j) = u_comp * (                           &
1718                                         sk(k,j,i) + sk(k,j,i-1)) * 0.5
1719                swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2), sk(k,j,i+1),  &
1720                                                   sk(k,j,i), sk(k,j,i-1),    &
1721                                                   sk(k,j,i-1), u_comp,       &
1722                                                   0.5, ddx ) 
1723             ENDDO
1724             
1725          ELSE
1726         
1727             DO  k = nzb_s_inner(j,i)+1, nzt
1728                u_comp                 = u(k,j,i) - u_gtrans
1729                swap_flux_x_local(k,j) = u_comp*(                              &
1730                                         37.0 * ( sk(k,j,i)+sk(k,j,i-1)     )  &
1731                                       -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
1732                                       +        ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
1733                                       * adv_sca_5
1734                swap_diss_x_local(k,j) = - ABS( u_comp ) * (                   &
1735                                         10.0 * (sk(k,j,i) - sk(k,j,i-1)    )  &
1736                                       -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
1737                                       +        ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
1738                                       * adv_sca_5 
1739             ENDDO
1740          ENDIF
1741       ENDDO
1742!
1743!--    The following loop computes the horizontal fluxes for the interior of the
1744!--    processor domain plus south boundary points. Furthermore tendency terms
1745!--    are computed.
1746       DO  i = nxl, nxr
1747          j = nys
1748          IF ( boundary_flags(j,i) == 8 )  THEN
1749         
1750             DO  k = nzb_s_inner(j,i)+1, nzt
1751                v_comp               = v(k,j,i) - v_gtrans
1752                swap_flux_y_local(k) = v_comp *                                &
1753                                       ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 
1754                swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),     &
1755                                                 sk(k,j,i), sk(k,j-1,i),       &
1756                                                 sk(k,j-1,i), v_comp, 0.5, ddy )
1757             ENDDO
1758             
1759          ELSE
1760         
1761             DO  k = nzb_s_inner(j,i)+1, nzt
1762                v_comp               = v(k,j,i) - v_gtrans
1763                swap_flux_y_local(k) = v_comp * (                              &
1764                                       37.0 * ( sk(k,j,i) + sk(k,j-1,i)   )    &
1765                                     -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )    &
1766                                     +        ( sk(k,j+2,i) + sk(k,j-3,i) ) )  & 
1767                                     * adv_sca_5
1768                swap_diss_y_local(k)= - ABS( v_comp ) * (                      &
1769                                      10.0 * ( sk(k,j,i) - sk(k,j-1,i)   )     &
1770                                    -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
1771                                    +        ( sk(k,j+2,i)-sk(k,j-3,i) ) )     &
1772                                    * adv_sca_5
1773             ENDDO
1774             
1775          ENDIF
1776
1777          DO  j = nys, nyn
1778            IF ( boundary_flags(j,i) /= 0 )  THEN 
1779!
1780!--            Degrade the order for Dirichlet bc. at the outflow boundary
1781               SELECT CASE ( boundary_flags(j,i) )
1782               
1783                  CASE ( 1 )
1784                     DO  k = nzb_s_inner(j,i)+1, nzt
1785                        u_comp    = u(k,j,i+1) - u_gtrans
1786                        flux_r(k) = u_comp * (                                 &
1787                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
1788                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
1789                                  * adv_sca_3
1790                        diss_r(k) = - ABS( u_comp ) * (                        &
1791                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        & 
1792                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
1793                                  * adv_sca_3
1794                     ENDDO
1795                     
1796                  CASE ( 2 )
1797                     DO  k = nzb_s_inner(j,i)+1, nzt
1798                        u_comp    = u(k,j,i+1) - u_gtrans
1799                        flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
1800                        diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1),        &
1801                                              sk(k,j,i), sk(k,j,i-1),          &
1802                                              sk(k,j,i-2), u_comp, 0.5, ddx )
1803                     ENDDO
1804                     
1805                  CASE ( 3 )
1806                     DO  k = nzb_s_inner(j,i)+1, nzt
1807                        v_comp    = v(k,j+1,i) - v_gtrans
1808                        flux_n(k) = v_comp * (                                 &
1809                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
1810                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
1811                                  * adv_sca_3
1812                        diss_n(k) = - ABS( v_comp ) * (                        &
1813                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        & 
1814                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
1815                                  * adv_sca_3
1816                     ENDDO
1817                     
1818                  CASE ( 4 )
1819                     DO  k = nzb_s_inner(j,i)+1, nzt
1820                        v_comp    = v(k,j+1,i) - v_gtrans
1821                        flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
1822                        diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i),        &
1823                                              sk(k,j,i), sk(k,j-1,i),          &
1824                                              sk(k,j-2,i), v_comp, 0.5, ddy )     
1825                     ENDDO
1826                     
1827                  CASE ( 5 )
1828                     DO  k = nzb_w_inner(j,i)+1, nzt
1829                        u_comp    = u(k,j,i+1) - u_gtrans 
1830                        flux_r(k) = u_comp * (                                 &
1831                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
1832                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
1833                                  * adv_sca_3
1834                        diss_r(k) = - ABS( u_comp ) * (                        &
1835                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        & 
1836                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
1837                                  * adv_sca_3
1838                     ENDDO 
1839                     
1840                  CASE ( 6 )
1841                     DO  k = nzb_s_inner(j,i)+1, nzt
1842                        u_comp    = u(k,j,i+1) - u_gtrans
1843                        flux_r(k) = u_comp * (                                 &
1844                                    7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )        &
1845                                  -       ( sk(k,j,i+2) + sk(k,j,i-1) ) )      &
1846                                  * adv_sca_3
1847                        diss_r(k) = - ABS( u_comp ) * (                        &
1848                                    3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )        & 
1849                                  -       ( sk(k,j,i+2) - sk(k,j,i-1) ) )      &
1850                                  * adv_sca_3
1851                     ENDDO
1852                     
1853                  CASE ( 7 )
1854                     DO  k = nzb_s_inner(j,i)+1, nzt
1855                        v_comp    = v(k,j+1,i) - v_gtrans
1856                        flux_n(k) = v_comp * (                                 &
1857                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
1858                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
1859                                  * adv_sca_3
1860                        diss_n(k) = - ABS( v_comp ) * (                        &
1861                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        & 
1862                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
1863                                  * adv_sca_3
1864                     ENDDO
1865                     
1866                  CASE ( 8 )
1867                     DO  k = nzb_s_inner(j,i)+1, nzt
1868                        v_comp    = v(k,j+1,i) - v_gtrans
1869                        flux_n(k) = v_comp * (                                 &
1870                                    7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )        &
1871                                  -       ( sk(k,j+2,i) + sk(k,j-1,i) ) )      &
1872                                  * adv_sca_3
1873                        diss_n(k) = -  ABS( v_comp ) * (                       &
1874                                    3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )        & 
1875                                  -       ( sk(k,j+2,i) - sk(k,j-1,i) ) )      &
1876                                  * adv_sca_3
1877                     ENDDO 
1878                     
1879                  CASE DEFAULT
1880                 
1881               END SELECT
1882               
1883               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
1884                   boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )    &
1885               THEN
1886             
1887                  DO  k = nzb_s_inner(j,i)+1, nzt
1888                     v_comp    = v(k,j+1,i) - v_gtrans
1889                     flux_n(k) = v_comp * (                                    &
1890                                 37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )          &
1891                               -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )          &
1892                               +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )        &
1893                               * adv_sca_5
1894                     diss_n(k) = - ABS( v_comp ) * (                           &
1895                                 10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )          &
1896                               -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )          &
1897                               +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )        &
1898                               * adv_sca_5
1899                  ENDDO
1900                 
1901               ELSE
1902               
1903                  DO  k = nzb_s_inner(j,i)+1, nzt
1904                     u_comp    = u(k,j,i+1) - u_gtrans 
1905                     flux_r(k) = u_comp * (                                    & 
1906                                 37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )          &
1907                               -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )          &
1908                               +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )        &
1909                               * adv_sca_5
1910                     diss_r(k) = - ABS( u_comp ) * (                           &
1911                                 10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )          &
1912                               -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )          &
1913                               +        ( sk(k,j,i+3) - sk(k,j,i-2) ) )        &
1914                               * adv_sca_5
1915                  ENDDO
1916                 
1917               ENDIF     
1918           
1919            ELSE
1920           
1921               DO  k = nzb_s_inner(j,i)+1, nzt
1922                  u_comp    = u(k,j,i+1) - u_gtrans
1923                  flux_r(k) = u_comp * (                                       &
1924                              37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
1925                            -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )             &
1926                            +        ( sk(k,j,i+3) + sk(k,j,i-2) ) )           &
1927                            * adv_sca_5
1928                  diss_r(k) = - ABS( u_comp ) * (                              &
1929                              10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
1930                            -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )             &
1931                            +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )             &
1932                            * adv_sca_5
1933                 
1934                  v_comp    = v(k,j+1,i) - v_gtrans
1935                  flux_n(k) = v_comp * (                                       &
1936                              37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
1937                            -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )             &
1938                            +        ( sk(k,j+3,i) + sk(k,j-2,i) ) )           &
1939                            * adv_sca_5
1940                  diss_n(k) = - ABS( v_comp ) * (                              &
1941                              10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
1942                            -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )             &
1943                            +        ( sk(k,j+3,i) - sk(k,j-2,i) ) )           &
1944                            * adv_sca_5                 
1945               ENDDO
1946               
1947            ENDIF
1948           
1949            DO  k = nzb_s_inner(j,i)+1, nzt
1950               tend(k,j,i) = tend(k,j,i) - (                                  &
1951                  ( flux_r(k) + diss_r(k)                                     &
1952                -   swap_flux_x_local(k,j) - swap_diss_x_local(k,j)   ) * ddx &
1953                + ( flux_n(k) + diss_n(k)                                     &
1954                -   swap_flux_y_local(k) - swap_diss_y_local(k)       ) * ddy)
1955                   
1956                swap_flux_x_local(k,j) = flux_r(k)
1957                swap_diss_x_local(k,j) = diss_r(k)
1958                swap_flux_y_local(k)   = flux_n(k)
1959                swap_diss_y_local(k)   = diss_n(k)
1960            ENDDO
1961         ENDDO
1962      ENDDO
1963
1964!
1965!--   Vertical advection, degradation of order near surface and top.
1966!--   The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1967!--   statistical evaluation the top flux at the surface should be 0
1968      DO  i = nxl, nxr
1969         DO  j = nys, nyn
1970!
1971!--         2nd order scheme (bottom)
1972            k=nzb_s_inner(j,i)+1
1973!           
1974!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
1975!--         reasons the top flux at the surface should be 0.
1976            flux_t(nzb_s_inner(j,i)) = 0.0
1977            diss_t(nzb_s_inner(j,i)) = 0.0
1978            flux_d = flux_t(k-1)
1979            diss_d = diss_t(k-1)
1980            flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
1981            diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i),        &
1982                                  sk(k,j,i), sk(k,j,i), w(k,j,i),             &
1983                                  0.5, ddzw(k) )
1984            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
1985                                      -   flux_d    - diss_d       ) * ddzw(k)   
1986!
1987!--         WS3 as an intermediate step (bottom)
1988            k         = nzb_s_inner(j,i)+2
1989            flux_d    = flux_t(k-1)
1990            diss_d    = diss_t(k-1)
1991            flux_t(k) = w(k,j,i) * (                                          &
1992                        7.0 * ( sk(k+1,j,i) + sk(k,j,i) )                     &
1993                       - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 
1994            diss_t(k) = - ABS( w(k,j,i) ) * (                                 &
1995                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   & 
1996                      -       ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
1997            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
1998                                      -   flux_d    - diss_d       ) * ddzw(k)
1999!
2000!--         WS5
2001            DO  k = nzb_s_inner(j,i)+3, nzt-2
2002               flux_d    = flux_t(k-1)
2003               diss_d    = diss_t(k-1)
2004               flux_t(k) = w(k,j,i) * (                                       &
2005                           37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )               &
2006                         -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )               &
2007                         +        ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
2008               diss_t(k) = - ABS(w(k,j,i)) * (                                &
2009                           10.0 * ( sk(k+1,j,i) -sk(k,j,i)    )               &
2010                         -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )               &
2011                         +        ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
2012
2013               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
2014                                         -   flux_d    - diss_d     ) * ddzw(k)
2015            ENDDO
2016!
2017!--         WS3 as an intermediate step (top)
2018            k         = nzt - 1
2019            flux_d    = flux_t(k-1)
2020            diss_d    = diss_t(k-1)
2021            flux_t(k) = w(k,j,i) * (                                          &
2022                        7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )                   &
2023                      -       ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
2024            diss_t(k) = - ABS(w(k,j,i)) * (                                   &
2025                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )                   & 
2026                      -       ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
2027                       
2028            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2029                                      -   flux_d    - diss_d       ) * ddzw(k)
2030!
2031!--         2nd order scheme (top)
2032            k         = nzt
2033            flux_d    = flux_t(k-1)
2034            diss_d    = diss_t(k-1)
2035            flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
2036            diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i),        &
2037                                  sk(k-1,j,i), sk(k-2,j,i), w(k,j,i),         &
2038                                  0.5, ddzw(k) )
2039
2040            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2041                                      -   flux_d    - diss_d       ) * ddzw(k)
2042!
2043!--         evaluation of statistics
2044            SELECT CASE ( sk_char )
2045
2046               CASE ( 'pt' )
2047                 DO  k = nzb_s_inner(j,i), nzt
2048                   sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn)                    &
2049                      + ( flux_t(k) + diss_t(k) )                             &
2050                      *   weight_substep(intermediate_timestep_count)
2051                 ENDDO
2052               CASE ( 'sa' )
2053                 DO  k = nzb_s_inner(j,i), nzt
2054                   sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn)                    &
2055                      + ( flux_t(k) + diss_t(k) )                             &
2056                      *   weight_substep(intermediate_timestep_count)
2057                 ENDDO
2058               CASE ( 'q' )
2059                 DO  k = nzb_s_inner(j,i), nzt
2060                   sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn)                      &
2061                      + ( flux_t(k) + diss_t(k) )                             &
2062                      *   weight_substep(intermediate_timestep_count)
2063                 ENDDO
2064
2065            END SELECT
2066         ENDDO
2067      ENDDO
2068
2069
2070    END SUBROUTINE advec_s_ws
2071
2072
2073!------------------------------------------------------------------------------!
2074! Advection of u - Call for all grid points
2075!------------------------------------------------------------------------------!
2076    SUBROUTINE advec_u_ws
2077
2078       USE arrays_3d
2079       USE constants
2080       USE control_parameters
2081       USE grid_variables
2082       USE indices
2083       USE statistics
2084
2085       IMPLICIT NONE
2086
2087       INTEGER ::  i, j, k, tn = 0
2088       REAL    ::  gu, gv, flux_d, diss_d, v_comp, w_comp
2089       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u, swap_diss_y_local_u
2090       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u,             &
2091                                             swap_diss_x_local_u
2092       REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n,  &
2093                                     diss_n, u_comp
2094 
2095       gu = 2.0 * u_gtrans
2096       gv = 2.0 * v_gtrans
2097
2098!
2099!--    Compute the fluxes for the whole left boundary of the processor domain.
2100       i = nxlu
2101       DO  j = nys, nyn
2102          IF( boundary_flags(j,i) == 5 )  THEN
2103         
2104             DO  k = nzb_u_inner(j,i)+1, nzt
2105                u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
2106                swap_flux_x_local_u(k,j) = u_comp(k) *                         &
2107                                           ( u(k,j,i) + u(k,j,i-1) ) * 0.25
2108                swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1),   &
2109                                                     u(k,j,i), u(k,j,i-1),     &
2110                                                     u(k,j,i-1), u_comp(k),    &
2111                                                     0.25, ddx )
2112             ENDDO
2113             
2114          ELSE
2115         
2116             DO  k = nzb_u_inner(j,i)+1, nzt
2117                u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
2118                swap_flux_x_local_u(k,j) = u_comp(k) * (                       &
2119                                           37.0 * ( u(k,j,i) + u(k,j,i-1)   )  &
2120                                         -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )  &
2121                                         +        (u(k,j,i+2)+u(k,j,i-3) )  )  &
2122                                         * adv_mom_5
2123                swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * (                &
2124                                           10.0 * ( u(k,j,i) - u(k,j,i-1)   )  &
2125                                         -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )  &
2126                                         +        ( u(k,j,i+2) - u(k,j,i-3) ) )&
2127                                         * adv_mom_5
2128             ENDDO
2129             
2130          ENDIF
2131         
2132       ENDDO
2133
2134       DO i = nxlu, nxr
2135!       
2136!--      The following loop computes the fluxes for the south boundary points
2137         j = nys
2138         IF ( boundary_flags(j,i) == 8 )  THEN
2139!         
2140!--         Compute southside fluxes for the south boundary of PE domain
2141            DO  k = nzb_u_inner(j,i)+1, nzt
2142               v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2143               swap_flux_y_local_u(k) = v_comp *                              &
2144                                        ( u(k,j,i) + u(k,j-1,i) ) * 0.25 
2145               swap_diss_y_local_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i),     &
2146                                                  u(k,j,i), u(k,j-1,i),       &
2147                                                  u(k,j-1,i), v_comp,         &
2148                                                  0.25, ddy )
2149            ENDDO
2150           
2151         ELSE
2152         
2153            DO  k = nzb_u_inner(j,i)+1, nzt
2154               v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2155               swap_flux_y_local_u(k) = v_comp * (                             &
2156                                        37.0 * ( u(k,j,i) + u(k,j-1,i)   )     &
2157                                      -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )     &
2158                                      +        ( u(k,j+2,i) + u(k,j-3,i) ) )   &
2159                                      * adv_mom_5
2160               swap_diss_y_local_u(k) = - ABS( v_comp ) * (                    &
2161                                        10.0 * ( u(k,j,i) - u(k,j-1,i)   )     &
2162                                      -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )     &
2163                                      +        ( u(k,j+2,i) - u(k,j-3,i) ) )   &
2164                                      * adv_mom_5
2165            ENDDO
2166           
2167         ENDIF
2168!         
2169!--      Computation of interior fluxes and tendency terms
2170         DO  j = nys, nyn
2171            IF ( boundary_flags(j,i) /= 0 )  THEN
2172!           
2173!--            Degrade the order for Dirichlet bc. at the outflow boundary
2174               SELECT CASE ( boundary_flags(j,i) )
2175               
2176                  CASE ( 1 )
2177                     DO  k = nzb_u_inner(j,i)+1, nzt
2178                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
2179                        flux_r(k) = ( u_comp(k) - gu ) * (                     &
2180                                    7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
2181                                  -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
2182                                  * adv_mom_3
2183                        diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
2184                                    3.0 * ( u(k,j,i+1) - u(k,j,i) )            & 
2185                                  -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
2186                                  * adv_mom_3
2187                     ENDDO
2188                     
2189                  CASE ( 2 )
2190                     DO  k = nzb_u_inner(j,i)+1, nzt
2191                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
2192                        flux_r(k) = ( u_comp(k) - gu ) *                      &
2193                                    ( u(k,j,i+1) + u(k,j,i) ) * 0.25 
2194                        diss_r(k) = diss_2nd( u(k,j,i+1), u(k,j,i+1),         &
2195                                              u(k,j,i), u(k,j,i-1),           &
2196                                              u(k,j,i-2), u_comp(k) ,0.25 ,ddx)
2197                     ENDDO
2198                     
2199                  CASE ( 3 )
2200                     DO  k = nzb_u_inner(j,i)+1, nzt
2201                        v_comp   = v(k,j+1,i) + v(k,j+1,i-1) - gv
2202                        flux_n(k) = v_comp * (                                 &
2203                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
2204                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
2205                                  * adv_mom_3
2206                        diss_n(k) = - ABS( v_comp ) * (                        &
2207                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          & 
2208                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
2209                                  * adv_mom_3
2210                     ENDDO
2211                     
2212                  CASE ( 4 )
2213                     DO  k = nzb_u_inner(j,i)+1, nzt
2214                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2215                        flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
2216                        diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i),         &
2217                                              u(k,j,i), u(k,j-1,i),           &
2218                                              u(k,j-2,i), v_comp, 0.25, ddy )
2219                     ENDDO
2220                     
2221                  CASE ( 5 )
2222                     DO  k = nzb_u_inner(j,i)+1, nzt       
2223                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
2224                        flux_r(k) = ( u_comp(k) - gu ) * (                     &
2225                                    7.0 * ( u(k,j,i+1) + u(k,j,i)   )          &
2226                                  -       ( u(k,j,i+2) + u(k,j,i-1) ) )        &
2227                                  * adv_mom_3
2228                        diss_r(k) = - ABS( u_comp(k) - gu ) * (                &
2229                                    3.0 * ( u(k,j,i+1) - u(k,j,i)   )          & 
2230                                  -       ( u(k,j,i+2) - u(k,j,i-1) ) )        &
2231                                  * adv_mom_3
2232                     ENDDO
2233                     
2234                  CASE ( 7 )
2235                     DO  k = nzb_u_inner(j,i)+1, nzt                           
2236                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2237                        flux_n(k) = v_comp * (                                 &
2238                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
2239                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
2240                                  * adv_mom_3
2241                        diss_n(k) = - ABS( v_comp ) * (                        &
2242                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          & 
2243                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
2244                                  * adv_mom_3
2245                     ENDDO 
2246                     
2247                  CASE ( 8 )
2248                     DO  k = nzb_u_inner(j,i)+1, nzt
2249                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2250                        flux_n(k) = v_comp * (                                 &
2251                                    7.0 * ( u(k,j+1,i) + u(k,j,i)   )          &
2252                                  -       ( u(k,j+2,i) + u(k,j-1,i) ) )        &
2253                                  * adv_mom_3
2254                        diss_n(k) = - ABS( v_comp ) * (                        &
2255                                    3.0 * ( u(k,j+1,i) - u(k,j,i)   )          & 
2256                                  -       ( u(k,j+2,i) - u(k,j-1,i) ) )        &
2257                                  * adv_mom_3
2258                     ENDDO 
2259                     
2260                  CASE DEFAULT
2261                 
2262               END SELECT
2263!               
2264!--            Compute the crosswise 5th order fluxes at the outflow
2265               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
2266                    boundary_flags(j,i) == 5 )  THEN
2267             
2268                  DO  k = nzb_u_inner(j,i)+1, nzt
2269                     v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2270                     flux_n(k) = v_comp * (                                    &
2271                                 37.0 * ( u(k,j+1,i) + u(k,j,i)   )            &
2272                               -  8.0 * ( u(k,j+2,i) +u(k,j-1,i)  )            &
2273                               +        ( u(k,j+3,i) + u(k,j-2,i) ) )          &
2274                               * adv_mom_5
2275                     diss_n(k) = - ABS( v_comp ) * (                           &
2276                                 10.0 * ( u(k,j+1,i) - u(k,j,i)   )            &
2277                               -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )            &
2278                               +        ( u(k,j+3,i) - u(k,j-2,i) ) )          &
2279                               * adv_mom_5
2280                  ENDDO
2281                 
2282               ELSE
2283               
2284                  DO  k = nzb_u_inner(j,i)+1, nzt
2285                     u_comp(k) = u(k,j,i+1) + u(k,j,i)
2286                     flux_r(k) = ( u_comp(k) - gu ) * (                        &
2287                                 37.0 * ( u(k,j,i+1) + u(k,j,i)   )            &
2288                               -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )            &
2289                               +        ( u(k,j,i+3) + u(k,j,i-2) ) )          &
2290                               * adv_mom_5
2291                     diss_r(k) = - ABS( u_comp(k) - gu ) * (                   &
2292                                 10.0 * ( u(k,j,i+1) - u(k,j,i)   )            &
2293                               -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )            &
2294                               +        ( u(k,j,i+3) - u(k,j,i-2) ) )          &
2295                               * adv_mom_5
2296                  ENDDO
2297                 
2298               ENDIF
2299               
2300            ELSE
2301           
2302               DO  k = nzb_u_inner(j,i)+1, nzt
2303                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
2304                  flux_r(k) = ( u_comp(k) - gu ) * (                           &
2305                              37.0 * ( u(k,j,i+1) + u(k,j,i)   )               &
2306                            -  8.0 * ( u(k,j,i+2) + u(k,j,i-1) )               &
2307                            +        ( u(k,j,i+3) + u(k,j,i-2) ) )             &
2308                            * adv_mom_5
2309                  diss_r(k) = - ABS( u_comp(k) - gu ) * (                      &
2310                              10.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
2311                            -  5.0 * ( u(k,j,i+2) - u(k,j,i-1) )               &
2312                            +        ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2313
2314                  v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2315                  flux_n(k) = v_comp * (                                       &
2316                              37.0 * ( u(k,j+1,i) + u(k,j,i)   )               &
2317                            -  8.0 * ( u(k,j+2,i) + u(k,j-1,i) )               &
2318                            +        ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2319                  diss_n(k) = - ABS( v_comp ) * (                              &
2320                              10.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
2321                            -  5.0 * ( u(k,j+2,i) - u(k,j-1,i) )               &
2322                            +        ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2323                               
2324               ENDDO
2325               
2326            ENDIF
2327           
2328            DO  k = nzb_u_inner(j,i)+1, nzt
2329           
2330               tend(k,j,i) = tend(k,j,i) - (                                   &
2331              ( flux_r(k) + diss_r(k)                                          &
2332            -   swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j)   ) * ddx  &
2333            + ( flux_n(k) + diss_n(k)                                          &
2334            -   swap_flux_y_local_u(k) - swap_diss_y_local_u(k)       ) * ddy )
2335                 
2336               swap_flux_x_local_u(k,j) = flux_r(k)   
2337               swap_diss_x_local_u(k,j) = diss_r(k)
2338               swap_flux_y_local_u(k)   = flux_n(k)     
2339               swap_diss_y_local_u(k)   = diss_n(k)
2340                     
2341               sums_us2_ws_l(k,tn)  = sums_us2_ws_l(k,tn)                            &
2342                 + ( flux_r(k)                                                 &
2343                 * ( u_comp(k) - 2.0 * hom(k,1,1,0) )                          &
2344                 / ( u_comp(k) - gu + 1.0E-20 )                                & 
2345                 +   diss_r(k)                                                 &
2346                 *   ABS( u_comp(k) - 2.0 * hom(k,1,1,0) )                     &
2347                 / ( ABS( u_comp(k) - gu) + 1.0E-20) )                         &
2348                 *   weight_substep(intermediate_timestep_count)
2349            ENDDO
2350            sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
2351         ENDDO
2352       ENDDO
2353
2354!
2355!--   Vertical advection, degradation of order near surface and top.
2356!--   The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
2357!--   statistical evaluation the top flux at the surface should be 0
2358       DO i = nxlu, nxr
2359          DO  j = nys, nyn
2360             k = nzb_u_inner(j,i)+1
2361!             
2362!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
2363!--         reasons the top flux at the surface should be 0.
2364            flux_t(nzb_u_inner(j,i)) = 0.0
2365            diss_t(nzb_u_inner(j,i)) = 0.0
2366            flux_d = flux_t(k-1)
2367            diss_d = diss_t(k-1)             
2368!
2369!--         2nd order scheme (bottom)
2370             w_comp    = w(k,j,i) + w(k,j,i-1)
2371             flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
2372             diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i),           &
2373                                   0.0, 0.0, w_comp, 0.25, ddzw(k) )
2374                                   
2375             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2376                                       -   flux_d    -  diss_d       ) * ddzw(k)
2377!
2378!--         WS3 as an intermediate step (bottom)
2379            k         = nzb_u_inner(j,i)+2
2380            flux_d    = flux_t(k-1)
2381            diss_d    = diss_t(k-1)
2382            w_comp    = w(k,j,i) + w(k,j,i-1)
2383            flux_t(k) = w_comp * (                                            &
2384                        7.0 * ( u(k+1,j,i) + u(k,j,i)   )                     &
2385                      -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
2386            diss_t(k) = - ABS( w_comp ) * (                                   &
2387                        3.0 * ( u(k+1,j,i) - u(k,j,i)   )                     & 
2388                      -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
2389
2390            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2391                                      -   flux_d    -  diss_d      ) * ddzw(k)
2392!
2393!WS5
2394             DO  k = nzb_u_inner(j,i)+3, nzt-2
2395
2396                flux_d    = flux_t(k-1)
2397                diss_d    = diss_t(k-1)
2398                w_comp    = w(k,j,i) + w(k,j,i-1)
2399                flux_t(k) = w_comp * (                                        &
2400                            37.0 * ( u(k+1,j,i) + u(k,j,i)   )                &
2401                          -  8.0 * ( u(k+2,j,i) + u(k-1,j,i) )                &
2402                          +        ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
2403                diss_t(k) = - ABS( w_comp ) * (                               &
2404                            10.0 * ( u(k+1,j,i) - u(k,j,i)   )                & 
2405                          -  5.0 * ( u(k+2,j,i) - u(k-1,j,i) )                & 
2406                          +        ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
2407
2408                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
2409                                          -   flux_d    -  diss_d   ) * ddzw(k)
2410
2411             ENDDO
2412!
2413!--          WS3 as an intermediate step (top)
2414             k         = nzt-1
2415             flux_d    = flux_t(k-1)
2416             diss_d    = diss_t(k-1)
2417             w_comp    = w(k,j,i) + w(k,j,i-1)
2418             flux_t(k) = w_comp * (                                           &
2419                         7.0 * ( u(k+1,j,i) + u(k,j,i) )                      &
2420                       -       ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
2421             diss_t(k) = - ABS( w_comp ) * (                                  &
2422                         3.0 * ( u(k+1,j,i) - u(k,j,i) )                      & 
2423                       -       ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
2424                           
2425             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
2426                                       -   flux_d    -  diss_d      ) * ddzw(k)
2427!
2428!--         2nd order scheme (top)
2429             k         = nzt
2430             flux_d    = flux_t(k-1)
2431             diss_d    = diss_t(k-1)
2432             w_comp    = w(k,j,i) + w(k,j,i-1)
2433             flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
2434             diss_t(k) = diss_2nd( u(nzt+1,j,i), u(nzt+1,j,i), u(k,j,i),      &
2435                                   u(k-1,j,i), u(k-2,j,i), w_comp,            &
2436                                   0.25, ddzw(k)) 
2437
2438             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
2439                                       -   flux_d    -  diss_d      ) * ddzw(k)
2440!
2441!-- at last vertical momentum flux is accumulated
2442            DO  k = nzb_u_inner(j,i), nzt
2443               sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn)                           &
2444                              + ( flux_t(k) + diss_t(k) )                      &
2445                              *   weight_substep(intermediate_timestep_count)
2446            ENDDO
2447          ENDDO
2448       ENDDO
2449
2450
2451    END SUBROUTINE advec_u_ws
2452   
2453   
2454!------------------------------------------------------------------------------!
2455! Advection of v - Call for all grid points
2456!------------------------------------------------------------------------------!
2457    SUBROUTINE advec_v_ws
2458
2459       USE arrays_3d
2460       USE constants
2461       USE control_parameters
2462       USE grid_variables
2463       USE indices
2464       USE statistics
2465
2466       IMPLICIT NONE
2467
2468
2469       INTEGER ::  i, j, k, tn = 0
2470       REAL    ::  gu, gv, flux_l, flux_s, flux_d, diss_l, diss_s, diss_d,    &
2471                   u_comp, w_comp
2472       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v, swap_diss_y_local_v
2473       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v,             &
2474                                             swap_diss_x_local_v
2475       REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_n, diss_n, flux_r,  &
2476                                     diss_r, v_comp
2477
2478       gu = 2.0 * u_gtrans
2479       gv = 2.0 * v_gtrans
2480!
2481!-- First compute the whole left boundary of the processor domain
2482       i = nxl
2483       DO  j = nysv, nyn
2484       
2485          IF ( boundary_flags(j,i) == 6 )  THEN
2486             DO  k = nzb_v_inner(j,i)+1, nzt
2487                u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
2488                swap_flux_x_local_v(k,j) = u_comp *                           &
2489                                           ( v(k,j,i) + v(k,j,i-1)) * 0.25
2490                swap_diss_x_local_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1),  &
2491                                                     v(k,j,i), v(k,j,i-1),    &
2492                                                     v(k,j,i-1), u_comp,      &
2493                                                     0.25, ddx )
2494             ENDDO
2495             
2496          ELSE
2497         
2498             DO  k = nzb_v_inner(j,i)+1, nzt
2499                u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
2500                swap_flux_x_local_v(k,j) = u_comp * (                          &
2501                                           37.0 * ( v(k,j,i) + v(k,j,i-1)   )  &
2502                                         -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )  &
2503                                         +        ( v(k,j,i+2) + v(k,j,i-3) ) )&
2504                                         * adv_mom_5
2505                swap_diss_x_local_v(k,j) = - ABS( u_comp ) * (                 &
2506                                           10.0 * ( v(k,j,i) - v(k,j,i-1)   )  &
2507                                         -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )  &
2508                                         +        ( v(k,j,i+2) - v(k,j,i-3) ) )&
2509                                         * adv_mom_5
2510             ENDDO
2511             
2512          ENDIF
2513         
2514       ENDDO
2515
2516       DO i = nxl, nxr
2517         j = nysv
2518         IF ( boundary_flags(j,i) == 7 )  THEN
2519         
2520            DO  k = nzb_v_inner(j,i)+1, nzt
2521               v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
2522               swap_flux_y_local_v(k) = v_comp(k) *                           &
2523                                        ( v(k,j,i) + v(k,j-1,i) ) * 0.25
2524               swap_diss_y_local_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i),     &
2525                                                  v(k,j,i), v(k,j-1,i),       &
2526                                                  v(k,j-1,i), v_comp(k),      &
2527                                                  0.25, ddy )                               
2528            ENDDO
2529           
2530         ELSE
2531         
2532            DO  k = nzb_v_inner(j,i)+1, nzt
2533               v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
2534               swap_flux_y_local_v(k) = v_comp(k) * (                          &
2535                                        37.0 * ( v(k,j,i) + v(k,j-1,i)   )     &
2536                                      -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )     &
2537                                      +        ( v(k,j+2,i) + v(k,j-3,i) ) )   &
2538                                      * adv_mom_5
2539               swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * (                 &
2540                                        10.0 * ( v(k,j,i) - v(k,j-1,i)   )     &
2541                                      -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )     &
2542                                      +        ( v(k,j+2,i) - v(k,j-3,i) ) )   &
2543                                      * adv_mom_5
2544            ENDDO
2545           
2546         ENDIF
2547         
2548         DO  j = nysv, nyn
2549            IF ( boundary_flags(j,i) /= 0 )  THEN
2550!       
2551!--       Degrade the order for Dirichlet bc. at the outflow boundary
2552               SELECT CASE ( boundary_flags(j,i) )
2553         
2554                  CASE ( 1 )
2555                     DO  k = nzb_v_inner(j,i)+1, nzt
2556                        u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2557                        flux_r(k) = u_comp * (                                 &
2558                                    7.0 * (v(k,j,i+1) + v(k,j,i)    )          &
2559                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )        &
2560                                  * adv_mom_3
2561                        diss_r(k) = - ABS( u_comp ) * (                        &
2562                                    3.0 * ( v(k,j,i+1) - v(k,j,i)   )          & 
2563                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )        &
2564                                  * adv_mom_3
2565                     ENDDO
2566               
2567                  CASE ( 2 )
2568                     DO  k = nzb_v_inner(j,i)+1, nzt
2569                        u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2570                        flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 
2571                        diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1),         &
2572                                              v(k,j,i), v(k,j,i-1),           &
2573                                              v(k,j,i-2), u_comp, 0.25, ddx )
2574                     ENDDO
2575               
2576                  CASE ( 3 )
2577                     DO  k = nzb_v_inner(j,i)+1, nzt
2578                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
2579                        flux_n(k) = ( v_comp(k)- gv ) * (                     &
2580                                    7.0 * ( v(k,j+1,i) + v(k,j,i)   )         &
2581                                  -       ( v(k,j+2,i) + v(k,j-1,i) ) )       &
2582                                  * adv_mom_3
2583                        diss_n(k) = - ABS(v_comp(k) - gv) * (                 &
2584                                    3.0 * ( v(k,j+1,i) - v(k,j,i)   )         & 
2585                                  -       ( v(k,j+2,i) - v(k,j-1,i) ) )       &
2586                                  * adv_mom_3
2587                     ENDDO
2588               
2589                  CASE ( 4 )
2590                     DO  k = nzb_v_inner(j,i)+1, nzt
2591                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
2592                        flux_n(k) = ( v_comp(k) - gv ) *                      &
2593                                   ( v(k,j+1,i) + v(k,j,i) ) * 0.25 
2594                        diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i),         &
2595                                              v(k,j,i), v(k,j-1,i),           &
2596                                              v(k,j-2,i), v_comp(k), 0.25, ddy)
2597                     ENDDO
2598               
2599                  CASE ( 5 )
2600                     DO  k = nzb_w_inner(j,i)+1, nzt
2601                        u_comp    = u(k,j-1,i) + u(k,j,i) - gu 
2602                        flux_r(k) = u_comp * (                                &
2603                                    7.0 * ( v(k,j,i+1) + v(k,j,i)   )         &
2604                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )       &
2605                                  * adv_mom_3
2606                        diss_r(k) = - ABS (u_comp ) * (                       &
2607                                    3.0 * ( w(k,j,i+1) - w(k,j,i)   )         & 
2608                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )       &
2609                                  * adv_mom_3
2610                     ENDDO
2611                 
2612                  CASE ( 6 )
2613                     DO  k = nzb_v_inner(j,i)+1, nzt
2614                                 
2615                        u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2616                        flux_r(k) = u_comp * (                                &
2617                                    7.0 * ( v(k,j,i+1) + v(k,j,i)   )         &
2618                                  -       ( v(k,j,i+2) + v(k,j,i-1) ) )       &
2619                                  * adv_mom_3
2620                        diss_r(k) = - ABS( u_comp ) * (                       &
2621                                    3.0 * ( v(k,j,i+1) - v(k,j,i)   )         & 
2622                                  -       ( v(k,j,i+2) - v(k,j,i-1) ) )       &
2623                                  * adv_mom_3 
2624                     ENDDO
2625               
2626                  CASE ( 7 )
2627                     DO  k = nzb_v_inner(j,i)+1, nzt                                 
2628                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
2629                        flux_n(k) = ( v_comp(k) - gv ) * (                    &
2630                                    7.0 * ( v(k,j+1,i) + v(k,j,i)   )         &
2631                                  -       ( v(k,j+2,i) + v(k,j-1,i) ) )       &
2632                                  * adv_mom_3
2633                        diss_n(k) = - ABS( v_comp(k) - gv ) * (               &
2634                                    3.0 * ( v(k,j+1,i) - v(k,j,i)   )         & 
2635                                  -       ( v(k,j+2,i) - v(k,j-1,i) ) )       &
2636                                  * adv_mom_3
2637                     ENDDO
2638               
2639                  CASE DEFAULT
2640             
2641               END SELECT
2642!         
2643!--            Compute the crosswise 5th order fluxes at the outflow
2644               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
2645                    boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )   &
2646               THEN
2647         
2648                  DO  k = nzb_v_inner(j,i)+1, nzt
2649                     v_comp(k) = v(k,j+1,i) + v(k,j,i)
2650                     flux_n(k) = ( v_comp(k) - gv ) * (                       &
2651                                 37.0 * ( v(k,j+1,i) + v(k,j,i)   )           &
2652                               -  8.0 * ( v(k,j+2,i) + v(k,j-1,i) )           &
2653                               +        ( v(k,j+3,i) + v(k,j-2,i) ) )         &
2654                               * adv_mom_5
2655                     diss_n(k) = - ABS( v_comp(k) - gv ) * (                  &
2656                                 10.0 * ( v(k,j+1,i) - v(k,j,i)   )           &
2657                               -  5.0 * ( v(k,j+2,i) - v(k,j-1,i) )           &
2658                               +        ( v(k,j+3,i) - v(k,j-2,i) ) )         &
2659                               * adv_mom_5
2660                   ENDDO
2661             
2662               ELSE
2663           
2664                  DO  k = nzb_v_inner(j,i)+1, nzt
2665                     u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2666                     flux_r(k) = u_comp * (                                   &
2667                                 37.0 * ( v(k,j,i+1) + v(k,j,i)   )           &
2668                              -   8.0 * ( v(k,j,i+2) + v(k,j,i-1) )           &
2669                              +         ( v(k,j,i+3) + v(k,j,i-2) ) )         &
2670                              * adv_mom_5
2671                     diss_r(k) = - ABS( u_comp ) * (                          &
2672