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

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

Bugfix in ws-scheme concerning OpenMP paralellization.

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