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

Last change on this file since 673 was 673, checked in by suehring, 13 years ago

Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.

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