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

Last change on this file since 744 was 744, checked in by suehring, 13 years ago

last commit documented

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