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

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

last commit documented

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