source: palm/tags/release-3.8/SOURCE/advec_ws.f90 @ 4343

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

last commit documented

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