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

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

reformatted advec_ws.f90, reformulate constants as broken numbers, bugfix in vertical advection of vertical velocity (concerning vector optimized routine)

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