Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 1932)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 1942)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Topography filter implemented to fill holes resolved by only one grid point.
+! Initialization of flags for ws-scheme moved to advec_ws.  
 ! 
 ! Former revisions:
@@ -171,4 +172,6 @@
  SUBROUTINE init_grid
  
+    USE advec_ws,                                                              &
+        ONLY:  ws_init_flags
 
     USE arrays_3d,                                                             &
@@ -185,6 +188,5 @@
                io_group, inflow_l, inflow_n, inflow_r, inflow_s,               &
                masking_method, maximum_grid_level, message_string,             &
-               momentum_advec, nest_domain, nest_bound_l, nest_bound_n,        &
-               nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n,        &
+               momentum_advec, nest_domain, ocean, outflow_l, outflow_n,       &
                outflow_r, outflow_s, psolver, scalar_advec, topography,        &
                topography_grid_convention, use_surface_fluxes, use_top_fluxes, &
@@ -197,6 +199,6 @@
         
     USE indices,                                                               &
-        ONLY:  flags, nbgp, nx, nxl, nxlg, nxlu, nxl_mg, nxr, nxrg, nxr_mg,    &
-               ny, nyn, nyng, nyn_mg, nys, nysv, nys_mg, nysg, nz, nzb,        &
+        ONLY:  flags, nbgp, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg,          &
+               ny, nyn, nyng, nyn_mg, nys, nys_mg, nysg, nz, nzb,              &
                nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u,       &
                nzb_diff_v, nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,     &
@@ -214,32 +216,34 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  bh      !< temporary vertical index of building height
-    INTEGER(iwp) ::  blx     !< grid point number of building size along x
-    INTEGER(iwp) ::  bly     !< grid point number of building size along y
-    INTEGER(iwp) ::  bxl     !< index for left building wall
-    INTEGER(iwp) ::  bxr     !< index for right building wall
-    INTEGER(iwp) ::  byn     !< index for north building wall
-    INTEGER(iwp) ::  bys     !< index for south building wall
-    INTEGER(iwp) ::  ch      !< temporary vertical index for canyon height
-    INTEGER(iwp) ::  cwx     !< grid point number of canyon size along x
-    INTEGER(iwp) ::  cwy     !< grid point number of canyon size along y
-    INTEGER(iwp) ::  cxl     !< index for left canyon wall
-    INTEGER(iwp) ::  cxr     !< index for right canyon wall
-    INTEGER(iwp) ::  cyn     !< index for north canyon wall
-    INTEGER(iwp) ::  cys     !< index for south canyon wall
-    INTEGER(iwp) ::  gls     !< number of lateral ghost points at total model domain boundaries required for multigrid solver
-    INTEGER(iwp) ::  i       !< index variable along x
-    INTEGER(iwp) ::  ii      !< loop variable for reading topography file
-    INTEGER(iwp) ::  inc     !< incremental parameter for coarsening grid level
-    INTEGER(iwp) ::  j       !< index variable along y
-    INTEGER(iwp) ::  k       !< index variable along z
-    INTEGER(iwp) ::  l       !< loop variable 
-    INTEGER(iwp) ::  nxl_l   !< index of left PE boundary for multigrid level
-    INTEGER(iwp) ::  nxr_l   !< index of right PE boundary for multigrid level
-    INTEGER(iwp) ::  nyn_l   !< index of north PE boundary for multigrid level
-    INTEGER(iwp) ::  nys_l   !< index of south PE boundary for multigrid level
-    INTEGER(iwp) ::  nzb_si  !< dummy index for local nzb_s_inner
-    INTEGER(iwp) ::  nzt_l   !< index of top PE boundary for multigrid level
-    INTEGER(iwp) ::  vi      !< dummy for vertical influence
+    INTEGER(iwp) ::  bh       !< temporary vertical index of building height
+    INTEGER(iwp) ::  blx      !< grid point number of building size along x
+    INTEGER(iwp) ::  bly      !< grid point number of building size along y
+    INTEGER(iwp) ::  bxl      !< index for left building wall
+    INTEGER(iwp) ::  bxr      !< index for right building wall
+    INTEGER(iwp) ::  byn      !< index for north building wall
+    INTEGER(iwp) ::  bys      !< index for south building wall
+    INTEGER(iwp) ::  ch       !< temporary vertical index for canyon height
+    INTEGER(iwp) ::  cwx      !< grid point number of canyon size along x
+    INTEGER(iwp) ::  cwy      !< grid point number of canyon size along y
+    INTEGER(iwp) ::  cxl      !< index for left canyon wall
+    INTEGER(iwp) ::  cxr      !< index for right canyon wall
+    INTEGER(iwp) ::  cyn      !< index for north canyon wall
+    INTEGER(iwp) ::  cys      !< index for south canyon wall
+    INTEGER(iwp) ::  gls      !< number of lateral ghost points at total model domain boundaries required for multigrid solver
+    INTEGER(iwp) ::  i        !< index variable along x
+    INTEGER(iwp) ::  ii       !< loop variable for reading topography file
+    INTEGER(iwp) ::  inc      !< incremental parameter for coarsening grid level
+    INTEGER(iwp) ::  j        !< index variable along y
+    INTEGER(iwp) ::  k        !< index variable along z
+    INTEGER(iwp) ::  l        !< loop variable 
+    INTEGER(iwp) ::  nxl_l    !< index of left PE boundary for multigrid level
+    INTEGER(iwp) ::  nxr_l    !< index of right PE boundary for multigrid level
+    INTEGER(iwp) ::  nyn_l    !< index of north PE boundary for multigrid level
+    INTEGER(iwp) ::  nys_l    !< index of south PE boundary for multigrid level
+    INTEGER(iwp) ::  nzb_si   !< dummy index for local nzb_s_inner
+    INTEGER(iwp) ::  nzt_l    !< index of top PE boundary for multigrid level
+    INTEGER(iwp) ::  num_hole !< number of holes (in topography) resolved by only one grid point 
+    INTEGER(iwp) ::  num_wall !< number of surrounding vertical walls for a single grid point
+    INTEGER(iwp) ::  vi       !< dummy for vertical influence
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::                               &
@@ -259,5 +263,5 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp    !< dummy to calculate topography indices on u- and v-grid
 
-    LOGICAL  :: flag_set = .FALSE.  !< steering variable for advection flags
+    LOGICAL  ::  hole = .FALSE.  !< flag to check if any holes resolved by only 1 grid point were filled
 
     REAL(wp) ::  dx_l          !< grid spacing along x on different multigrid level
@@ -744,5 +748,58 @@
 
           DEALLOCATE ( topo_height )
-
+!
+!--       Filter topography, i.e. fill holes resolved by only one grid point.  
+!--       Such holes are suspected to lead to velocity blow-ups as continuity 
+!--       equation on discrete grid cannot be fulfilled in such case.
+!--       For now, check only for holes and fill them to the lowest height level
+!--       of the directly adjoining grid points along x- and y- direction. 
+!--       Before checking for holes, set lateral boundary conditions for 
+!--       topography. After hole-filling, boundary conditions must be set again!
+          IF ( bc_ns_cyc )  THEN
+             nzb_local(-1,:)   = nzb_local(ny,:)
+             nzb_local(ny+1,:) = nzb_local(0,:)
+          ELSE
+             nzb_local(-1,:)   = nzb_local(0,:)
+             nzb_local(ny+1,:) = nzb_local(ny,:)
+          ENDIF
+
+          IF ( bc_lr_cyc )  THEN
+             nzb_local(:,-1)   = nzb_local(:,nx)
+             nzb_local(:,nx+1) = nzb_local(:,0)
+          ELSE
+             nzb_local(:,-1)   = nzb_local(:,0)
+             nzb_local(:,nx+1) = nzb_local(:,nx)
+          ENDIF
+
+          num_hole = 0
+          DO i = 0, nx
+             DO j = 0, ny
+
+                num_wall = 0
+
+                IF ( nzb_local(j-1,i) > nzb_local(j,i) )                       &
+                   num_wall = num_wall + 1
+                IF ( nzb_local(j+1,i) > nzb_local(j,i) )                       &
+                   num_wall = num_wall + 1
+                IF ( nzb_local(j,i-1) > nzb_local(j,i) )                       &
+                   num_wall = num_wall + 1
+                IF ( nzb_local(j,i+1) > nzb_local(j,i) )                       &
+                   num_wall = num_wall + 1
+
+                IF ( num_wall == 4 )  THEN
+                   hole           = .TRUE.
+                   nzb_local(j,i) = MIN( nzb_local(j-1,i), nzb_local(j+1,i),   &
+                                         nzb_local(j,i-1), nzb_local(j,i+1) )
+                   num_hole       = num_hole + 1
+                ENDIF
+             ENDDO
+          ENDDO
+!
+!--       Create an informative message if any hole was removed.
+          IF ( hole )  THEN
+             WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '//&
+                                                  'one grid point were filled'
+             CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 )
+          ENDIF
 !
 !--       Add cyclic or Neumann boundary conditions (additional layers are for
@@ -1237,15 +1294,4 @@
           ENDIF
 
-!
-!--       Test output of flag arrays
-!          i = nxl_l
-!          WRITE (9,*)  ' '
-!          WRITE (9,*)  '*** mg level ', l, ' ***', mg_switch_to_pe0_level
-!          WRITE (9,*)  '    inc=', inc, '  i =', nxl_l
-!          WRITE (9,*)  '    nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
-!          DO  k = nzt_l+1, nzb, -1
-!             WRITE (9,'(194(1X,I2))')  ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
-!          ENDDO
-
           inc = inc * 2
 
@@ -1254,423 +1300,16 @@
     ENDIF
 !
-!-- Allocate flags needed for masking walls.
+!-- Allocate flags needed for masking walls. Even though these flags are only
+!-- required in the ws-scheme, the arrays need to be allocated as they are 
+!-- used in OpenACC directives. 
     ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
               wall_flags_00(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     wall_flags_0  = 0
     wall_flags_00 = 0
-
-    IF ( scalar_advec == 'ws-scheme' .OR.                                      &
-         scalar_advec == 'ws-scheme-mono' )  THEN
-!
-!--    Set flags to steer the degradation of the advection scheme in advec_ws
-!--    near topography, inflow- and outflow boundaries as well as bottom and
-!--    top of model domain. wall_flags_0 remains zero for all non-prognostic
-!--    grid points.
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb_s_inner(j,i)+1, nzt
-!
-!--             scalar - x-direction
-!--             WS1 (0), WS3 (1), WS5 (2)
-                IF ( ( k <= nzb_s_inner(j,i+1) .OR. k <= nzb_s_inner(j,i+2)    &   
-                  .OR. k <= nzb_s_inner(j,i-1) )                               &
-                    .OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i == nxl   )    .OR.                          &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr   ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
-                ELSEIF ( ( k <= nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i+1)&
-                                                   .AND. k > nzb_s_inner(j,i+2)&
-                                                   .AND. k > nzb_s_inner(j,i-1)&
-                         )                       .OR.                          &
-                         ( k <= nzb_s_inner(j,i-2) .AND. k > nzb_s_inner(j,i+1)&
-                                                   .AND. k > nzb_s_inner(j,i+2)&
-                                                   .AND. k > nzb_s_inner(j,i-1)&
-                         )                                                     &
-                                                 .OR.                          &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr-1 )    .OR.                          &
-                         ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i == nxlu  ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
-                ELSEIF ( k > nzb_s_inner(j,i+1) .AND. k > nzb_s_inner(j,i+2)   &
-                   .AND. k > nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i-1)   &
-                   .AND. k > nzb_s_inner(j,i-2) )                              &
-                THEN 
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
-                ENDIF
-!
-!--             scalar - y-direction
-!--             WS1 (3), WS3 (4), WS5 (5)
-                IF ( ( k <= nzb_s_inner(j+1,i)  .OR. k <= nzb_s_inner(j+2,i)   &   
-                                                .OR. k <= nzb_s_inner(j-1,i) ) &
-                                                 .OR.                          &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j == nys   )    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn   ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
-!
-!--             WS3
-                ELSEIF ( ( k <= nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j+1,i)&
-                                                   .AND. k > nzb_s_inner(j+2,i)&
-                                                   .AND. k > nzb_s_inner(j-1,i)&
-                         )                           .OR.                      &
-                         ( k <= nzb_s_inner(j-2,i) .AND. k > nzb_s_inner(j+1,i)&
-                                                   .AND. k > nzb_s_inner(j+2,i)&
-                                                   .AND. k > nzb_s_inner(j-1,i)&
-                         )                                                     &
-                                                     .OR.                      &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j == nysv  )    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn-1 ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
-!
-!--             WS5
-                ELSEIF ( k > nzb_s_inner(j+1,i) .AND. k > nzb_s_inner(j+2,i)   &
-                   .AND. k > nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j-1,i)   &
-                   .AND. k > nzb_s_inner(j-2,i) )                              &
-                THEN 
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
-                ENDIF
-!
-!--             scalar - z-direction
-!--             WS1 (6), WS3 (7), WS5 (8)
-                flag_set = .FALSE.
-                IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 )
-                   flag_set = .TRUE.
-                ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 )
-                   flag_set = .TRUE.
-                ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 )
-                ENDIF
-
-             ENDDO
-          ENDDO
-       ENDDO
-    ENDIF
-
-    IF ( momentum_advec == 'ws-scheme' )  THEN
-!
-!--    Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws
-!--    near topography, inflow- and outflow boundaries as well as bottom and
-!--    top of model domain. wall_flags_0 remains zero for all non-prognostic
-!--    grid points.
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-!
-!--             At first, set flags to WS1. 
-!--             Since fluxes are swapped in advec_ws.f90, this is necessary to 
-!--             in order to handle the left/south flux.
-!--             near vertical walls. 
-                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
-                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
-!
-!--             u component - x-direction
-!--             WS1 (9), WS3 (10), WS5 (11)
-                IF ( k <= nzb_u_inner(j,i+1)     .OR.                          &
-                         ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i <= nxlu  )    .OR.                          &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr   ) )                                &
-                THEN
-                    wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
-                ELSEIF ( ( k <= nzb_u_inner(j,i+2) .AND.                       &
-                           k >  nzb_u_inner(j,i+1) ) .OR.                      &
-                           k <= nzb_u_inner(j,i-1)                             &
-                                                     .OR.                      &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr-1 )    .OR.                          &
-                         ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i == nxlu+1) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 )
-                ELSEIF ( k > nzb_u_inner(j,i+1) .AND. k > nzb_u_inner(j,i+2)   &
-                                                .AND. k > nzb_u_inner(j,i-1) ) &
-                THEN    
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 )
-                ENDIF
-!
-!--             u component - y-direction
-!--             WS1 (12), WS3 (13), WS5 (14)
-                IF ( k <= nzb_u_inner(j+1,i) .OR.                              &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j == nys   )    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn   ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
-                ELSEIF ( ( k <= nzb_u_inner(j+2,i) .AND.                       &
-                           k >  nzb_u_inner(j+1,i) ) .OR.                      &
-                           k <= nzb_u_inner(j-1,i)                             &
-                                                     .OR.                      &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j == nysv  )    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn-1 ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 )
-                ELSEIF ( k > nzb_u_inner(j+1,i) .AND. k > nzb_u_inner(j+2,i)   &
-                                                .AND. k > nzb_u_inner(j-1,i) ) &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 )
-                ENDIF
-!
-!--             u component - z-direction
-!--             WS1 (15), WS3 (16), WS5 (17)
-                flag_set = .FALSE.
-                IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 )
-                   flag_set = .TRUE.
-                ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 )
-                   flag_set = .TRUE.
-                ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 )
-                ENDIF
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-!
-!--             At first, set flags to WS1. 
-!--             Since fluxes are swapped in advec_ws.f90, this is necessary to 
-!--             in order to handle the left/south flux.
-                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
-                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
-!
-!--             v component - x-direction
-!--             WS1 (18), WS3 (19), WS5 (20)
-                IF ( k <= nzb_v_inner(j,i+1) .OR.                              &
-                         ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i == nxl   )    .OR.                          &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr   ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
-!
-!--             WS3
-                ELSEIF ( ( k <= nzb_v_inner(j,i+2) .AND.                       &
-                           k >  nzb_v_inner(j,i+1) ) .OR.                      &
-                           k <= nzb_v_inner(j,i-1)                             &
-                                                 .OR.                          &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr-1 )    .OR.                          &
-                         ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i == nxlu  ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 )
-                ELSEIF ( k > nzb_v_inner(j,i+1) .AND. k > nzb_v_inner(j,i+2)   &
-                                                .AND. k > nzb_v_inner(j,i-1) ) &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 )
-                ENDIF
-!
-!--             v component - y-direction
-!--             WS1 (21), WS3 (22), WS5 (23)
-                IF ( k <= nzb_v_inner(j+1,i) .OR.                              &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j <= nysv  )    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn   ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
-                ELSEIF ( ( k <= nzb_v_inner(j+2,i) .AND.                       &
-                           k >  nzb_v_inner(j+1,i) ) .OR.                      &
-                           k <= nzb_v_inner(j-1,i)                             &
-                                                     .OR.                      &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j == nysv+1)    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn-1 ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 )
-                ELSEIF ( k > nzb_v_inner(j+1,i) .AND. k > nzb_v_inner(j+2,i)   &
-                                                .AND. k > nzb_v_inner(j-1,i) ) &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 )
-                ENDIF
-!
-!--             v component - z-direction
-!--             WS1 (24), WS3 (25), WS5 (26)
-                flag_set = .FALSE.
-                IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
-                   flag_set = .TRUE.
-                ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 )
-                   flag_set = .TRUE.
-                ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set )  THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 )
-                ENDIF
-
-             ENDDO
-          ENDDO
-       ENDDO
-       DO  i = nxl, nxr
-          DO  j = nys, nyn
-             DO  k = nzb+1, nzt
-!
-!--             At first, set flags to WS1. 
-!--             Since fluxes are swapped in advec_ws.f90, this is necessary to 
-!--             in order to handle the left/south flux.
-                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
-                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
-!
-!--             w component - x-direction
-!--             WS1 (27), WS3 (28), WS5 (29)
-                IF ( k <= nzb_w_inner(j,i+1) .OR.                              &
-                         ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i == nxl   )    .OR.                          &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr   ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
-                ELSEIF ( ( k <= nzb_w_inner(j,i+2) .AND.                       &
-                           k >  nzb_w_inner(j,i+1) ) .OR.                      &
-                           k <= nzb_w_inner(j,i-1)                             &
-                                                     .OR.                      &
-                         ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
-                           .AND. i == nxr-1 )    .OR.                          &
-                         ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
-                           .AND. i == nxlu  ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 )
-                ELSEIF ( k > nzb_w_inner(j,i+1) .AND. k > nzb_w_inner(j,i+2)   &
-                                                .AND. k > nzb_w_inner(j,i-1) ) &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 )
-                ENDIF
-!
-!--             w component - y-direction
-!--             WS1 (30), WS3 (31), WS5 (32)
-                IF ( k <= nzb_w_inner(j+1,i) .OR.                              &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j == nys   )    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn   ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
-                ELSEIF ( ( k <= nzb_w_inner(j+2,i) .AND.                       &
-                           k >  nzb_w_inner(j+1,i) ) .OR.                      &
-                           k <= nzb_w_inner(j-1,i)                             &
-                                                     .OR.                      &
-                         ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
-                           .AND. j == nysv  )    .OR.                          &
-                         ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
-                           .AND. j == nyn-1 ) )                                &
-                THEN
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 )
-                ELSEIF ( k > nzb_w_inner(j+1,i) .AND. k > nzb_w_inner(j+2,i)   &
-                                                .AND. k > nzb_w_inner(j-1,i) ) &
-                THEN
-                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 0 )
-!
-!--                Clear flag for WS1
-                   wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 )
-                ENDIF
-!
-!--             w component - z-direction
-!--             WS1 (33), WS3 (34), WS5 (35)
-                flag_set = .FALSE.
-                IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1      &
-                                           .OR. k == nzt )  THEN
-!
-!--                Please note, at k == nzb_w_inner(j,i) a flag is explictely
-!--                set, although this is not a prognostic level. However, 
-!--                contrary to the advection of u,v and s this is necessary
-!--                because flux_t(nzb_w_inner(j,i)) is used for the tendency
-!--                at k == nzb_w_inner(j,i)+1.
-                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 1 )
-                   flag_set = .TRUE.
-                ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
-                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 2 )
-                   flag_set = .TRUE.
-                ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set )  THEN
-                   wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 3 )
-                ENDIF
-
-             ENDDO
-          ENDDO
-       ENDDO
-
-    ENDIF
-
-!
-!-- Exchange 3D integer wall_flags.
-    IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme'     &
-    .OR. scalar_advec == 'ws-scheme-mono' )  THEN  
-!
-!--    Exchange ghost points for advection flags
-       CALL exchange_horiz_int( wall_flags_0,  nbgp )
-       CALL exchange_horiz_int( wall_flags_00, nbgp )
-!
-!--    Set boundary flags at inflow and outflow boundary in case of 
-!--    non-cyclic boundary conditions.
-       IF ( inflow_l .OR. outflow_l .OR. nest_bound_l )  THEN
-          wall_flags_0(:,:,nxl-1)  = wall_flags_0(:,:,nxl)
-          wall_flags_00(:,:,nxl-1) = wall_flags_00(:,:,nxl)
-       ENDIF
-
-       IF ( inflow_r .OR. outflow_r .OR. nest_bound_r )  THEN
-          wall_flags_0(:,:,nxr+1)  = wall_flags_0(:,:,nxr)
-          wall_flags_00(:,:,nxr+1) = wall_flags_00(:,:,nxr)
-       ENDIF
-
-       IF ( inflow_n .OR. outflow_n .OR. nest_bound_n )  THEN
-          wall_flags_0(:,nyn+1,:)  = wall_flags_0(:,nyn,:)
-          wall_flags_00(:,nyn+1,:) = wall_flags_00(:,nyn,:)
-       ENDIF
-
-       IF ( inflow_s .OR. outflow_s  .OR. nest_bound_s )  THEN
-          wall_flags_0(:,nys-1,:)  = wall_flags_0(:,nys,:)
-          wall_flags_00(:,nys-1,:) = wall_flags_00(:,nys,:)
-       ENDIF
-
+!
+!-- Init flags for ws-scheme to degrade order near walls
+    IF ( momentum_advec == 'ws-scheme'  .OR.  scalar_advec == 'ws-scheme'  .OR.&
+         scalar_advec   == 'ws-scheme-mono' )  THEN
+       CALL ws_init_flags
     ENDIF
 
