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

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

Bugfix concerning OpenMP parallelization. Calculation of turbulent fluxes in advec_ws is now thread-safe.

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