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

Last change on this file since 709 was 709, checked in by raasch, 10 years ago

formatting adjustments

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