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

Last change on this file since 857 was 857, checked in by suehring, 9 years ago

last commit documented

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