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

Last change on this file since 668 was 668, checked in by suehring, 11 years ago

last commit documented

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