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

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

last commit documented

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