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

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

Calculation of turbulent fluxes with WS-scheme only for the whole model domain, not for user-defined subregions

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