Index: /palm/trunk/SOURCE/poismg_fast_mod.f90
===================================================================
--- /palm/trunk/SOURCE/poismg_fast_mod.f90	(revision 1897)
+++ /palm/trunk/SOURCE/poismg_fast_mod.f90	(revision 1898)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: bottom and top boundary condition in resid_fast
+! Bugfix: restriction at nzb+1
+! formatting adjustments, variable descriptions added in some declaration blocks
 ! 
 ! Former revisions:
@@ -247,7 +249,7 @@
        USE control_parameters,                                                 &
            ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
-                  inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l,  &
-                  nest_bound_n, nest_bound_r, nest_bound_s, outflow_l,         &
-                  outflow_n, outflow_r, outflow_s, topography
+                  inflow_n, inflow_r, inflow_s, nest_bound_l, nest_bound_n,    &
+                  nest_bound_r, nest_bound_s, outflow_l, outflow_n, outflow_r, &
+                  outflow_s
 
        USE grid_variables,                                                     &
@@ -255,27 +257,24 @@
 
        USE indices,                                                            &
-           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
-                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
-                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
-                  nys_mg, nyn_mg, nzb, nzt_mg
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i
-       INTEGER(iwp) ::  j
-       INTEGER(iwp) ::  k
-       INTEGER(iwp) ::  l
-       INTEGER(iwp) :: km1      !<
-       INTEGER(iwp) :: kp1      !<
+       INTEGER(iwp) ::  i        !< index variable along x
+       INTEGER(iwp) ::  j        !< index variable along y
+       INTEGER(iwp) ::  k        !< index variable along z
+       INTEGER(iwp) ::  l        !< index indicating grid level
+       INTEGER(iwp) ::  km1      !< index variable along z dimension (k-1)
+       INTEGER(iwp) ::  kp1      !< index variable along z dimension (k+1)
 
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
-                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !< velocity divergence
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
-                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !< perturbation pressure
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
-                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !<
+                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !< residuum of perturbation pressure
 
 !
@@ -284,116 +283,37 @@
 
        CALL cpu_log( log_point_s(53), 'resid', 'start' )
-       IF ( topography == 'flat'  .OR.  masking_method )  THEN
-          !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
-          !$OMP DO
-          DO  i = nxl_mg(l), nxr_mg(l)
-             DO  j = nys_mg(l), nyn_mg(l)
+       !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
+       !$OMP DO
+       DO  i = nxl_mg(l), nxr_mg(l)
+          DO  j = nys_mg(l), nyn_mg(l)
                 !DIR$ IVDEP
-                DO k = ind_even_odd+1, nzt_mg(l)
-                   km1 = k-ind_even_odd-1
-                   kp1 = k-ind_even_odd
-                   r(k,j,i) = f_mg(k,j,i)                                 &
-                      - ddx2_mg(l) *                                      &
-                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                 &
-                      - ddy2_mg(l) *                                      &
-                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                  &
-                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                      &
-                      - f3_mg_b(k,l) * p_mg(km1,j,i)                      &
+             DO k = ind_even_odd+1, nzt_mg(l)
+                km1 = k-ind_even_odd-1
+                kp1 = k-ind_even_odd
+                r(k,j,i) = f_mg(k,j,i)                                         &
+                      - ddx2_mg(l) *                                           &
+                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
+                      - ddy2_mg(l) *                                           &
+                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
+                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
+                      - f3_mg_b(k,l) * p_mg(km1,j,i)                           &
                       + f1_mg_b(k,l) * p_mg(k,j,i)
-                ENDDO
-                !DIR$ IVDEP
-                DO k = nzb+1, ind_even_odd
-                   km1 = k+ind_even_odd
-                   kp1 = k+ind_even_odd+1
-                   r(k,j,i) = f_mg(k,j,i)                                 &
-                      - ddx2_mg(l) *                                      &
-                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                 &
-                      - ddy2_mg(l) *                                      &
-                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                  &
-                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                      &
-                      - f3_mg_b(k,l) * p_mg(km1,j,i)                      &
+             ENDDO
+             !DIR$ IVDEP
+             DO k = nzb+1, ind_even_odd
+                km1 = k+ind_even_odd
+                kp1 = k+ind_even_odd+1
+                r(k,j,i) = f_mg(k,j,i)                                         &
+                      - ddx2_mg(l) *                                           &
+                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                      &
+                      - ddy2_mg(l) *                                           &
+                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                       &
+                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                           &
+                      - f3_mg_b(k,l) * p_mg(km1,j,i)                           &
                       + f1_mg_b(k,l) * p_mg(k,j,i)
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-       ELSE
-!
-!--       Choose flag array of this level
-          SELECT CASE ( l )
-             CASE ( 1 )
-                flags => wall_flags_1
-             CASE ( 2 )
-                flags => wall_flags_2
-             CASE ( 3 )
-                flags => wall_flags_3
-             CASE ( 4 )
-                flags => wall_flags_4
-             CASE ( 5 )
-                flags => wall_flags_5
-             CASE ( 6 )
-                flags => wall_flags_6
-             CASE ( 7 )
-                flags => wall_flags_7
-             CASE ( 8 )
-                flags => wall_flags_8
-             CASE ( 9 )
-                flags => wall_flags_9
-             CASE ( 10 )
-                flags => wall_flags_10
-          END SELECT
-
-          !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
-          !$OMP DO
-          DO  i = nxl_mg(l), nxr_mg(l)
-             DO  j = nys_mg(l), nyn_mg(l)
-
-                !DIR$ IVDEP
-                DO  k = ind_even_odd+1, nzt_mg(l)
-                   km1 = k-ind_even_odd-1
-                   kp1 = k-ind_even_odd
-
-                   r(k,j,i) = f_mg(k,j,i)                                         &
-                    - ddx2_mg(l) *                                                &
-                    ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 ))   &
-                    + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
-                    - ddy2_mg(l)                                                  &
-                    * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
-                    + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
-                    - f2_mg(k,l) * p_mg(kp1,j,i)                                  &
-                    - f3_mg(k,l) *                                                &
-                      MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 ))   &
-                    + f1_mg(k,l) * p_mg(k,j,i)
-                ENDDO
-
-                !DIR$ IVDEP
-                DO  k = nzb+1, ind_even_odd
-                   km1 = k+ind_even_odd
-                   kp1 = k+ind_even_odd+1
-
-                   r(k,j,i) = f_mg(k,j,i)                                         &
-                    - ddx2_mg(l) *                                                &
-                    ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 ))   &
-                    + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
-                    - ddy2_mg(l)                                                  &
-                    * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
-                    + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
-                    - f2_mg(k,l) * p_mg(kp1,j,i)                                  &
-                    - f3_mg(k,l) *                                                &
-                      MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 ))   &
-                    + f1_mg(k,l) * p_mg(k,j,i)
-                ENDDO
-!
-!--             The residual within topography should be zero
-                WHERE( BTEST(flags(nzb+1:nzt_mg(l),j,i), 6) )
-                   r(nzb+1:nzt_mg(l),j,i) = 0.0
-                END WHERE
-
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-
-       ENDIF
-
+             ENDDO
+          ENDDO
+       ENDDO
+       !$OMP END PARALLEL
 !
 !--    Horizontal boundary conditions
@@ -419,9 +339,10 @@
 
 !
-!--    Boundary conditions at bottom and top of the domain.outflow_l, outflow_n,
-!--    These points are not handled by the above loop. Points may be within
+!--    Boundary conditions at bottom and top of the domain. Points may be within
 !--    buildings, but that doesn't matter.
        IF ( ibc_p_b == 1 )  THEN
-          r(nzb,:,: ) = r(nzb+1,:,:)
+!
+!--       equivalent to r(nzb,:,: ) = r(nzb+1,:,:)
+          r(nzb,:,: ) = r(ind_even_odd+1,:,:)
        ELSE
           r(nzb,:,: ) = 0.0_wp
@@ -429,5 +350,7 @@
 
        IF ( ibc_p_t == 1 )  THEN
-          r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
+!
+!--       equivalent to r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
+          r(nzt_mg(l)+1,:,: ) = r(ind_even_odd,:,:)
        ELSE
           r(nzt_mg(l)+1,:,: ) = 0.0_wp
@@ -450,51 +373,33 @@
        USE control_parameters,                                                 &
            ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
-                  inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l,  &
-                  nest_bound_n, nest_bound_r, nest_bound_s, outflow_l,         &
-                  outflow_n, outflow_r, outflow_s, topography
+                  inflow_n, inflow_r, inflow_s, nest_bound_l, nest_bound_n,    &
+                  nest_bound_r, nest_bound_s, outflow_l, outflow_n, outflow_r, &
+                  outflow_s
 
        USE indices,                                                            &
-           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
-                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
-                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
-                  nys_mg, nyn_mg, nzb, nzt_mg
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i    !<
-       INTEGER(iwp) ::  ic   !<
-       INTEGER(iwp) ::  j    !<
-       INTEGER(iwp) ::  jc   !<
-       INTEGER(iwp) ::  k    !<
-       INTEGER(iwp) ::  kc   !<
-       INTEGER(iwp) ::  l    !<
-       INTEGER(iwp) ::  km1  !<
-       INTEGER(iwp) ::  kp1  !<
-
-       REAL(wp) ::  rkjim    !<
-       REAL(wp) ::  rkjip    !<
-       REAL(wp) ::  rkjmi    !<
-       REAL(wp) ::  rkjmim   !<
-       REAL(wp) ::  rkjmip   !<
-       REAL(wp) ::  rkjpi    !<
-       REAL(wp) ::  rkjpim   !<
-       REAL(wp) ::  rkjpip   !<
-       REAL(wp) ::  rkmji    !<
-       REAL(wp) ::  rkmjim   !<
-       REAL(wp) ::  rkmjip   !<
-       REAL(wp) ::  rkmjmi   !<
-       REAL(wp) ::  rkmjmim  !<
-       REAL(wp) ::  rkmjmip  !<
-       REAL(wp) ::  rkmjpi   !<
-       REAL(wp) ::  rkmjpim  !<
-       REAL(wp) ::  rkmjpip  !<
+       INTEGER(iwp) ::  i    !< index variable along x on finer grid
+       INTEGER(iwp) ::  ic   !< index variable along x on coarser grid
+       INTEGER(iwp) ::  j    !< index variable along y on finer grid
+       INTEGER(iwp) ::  jc   !< index variable along y on coarser grid
+       INTEGER(iwp) ::  k    !< index variable along z on finer grid
+       INTEGER(iwp) ::  kc   !< index variable along z on coarser grid
+       INTEGER(iwp) ::  l    !< index indicating finer grid level
+       INTEGER(iwp) ::  km1  !< index variable along z dimension (k-1 on finer level) 
+       INTEGER(iwp) ::  kp1  !< index variable along z dimension (k+1 on finer level)
+
 
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                            nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
-                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
+                                         f_mg  !< Residual on coarser grid level
 
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,                         &
                            nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,      &
-                           nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r !<
+                           nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::   &
+                                         r !< Residual on finer grid level
 
 !
@@ -503,152 +408,50 @@
 
        CALL cpu_log( log_point_s(54), 'restrict', 'start' )
-
-       IF ( topography == 'flat'  .OR.  masking_method )  THEN
-!
-!--       No wall treatment
-          !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1)
-          DO  ic = nxl_mg(l), nxr_mg(l)
-             i = 2*ic
-             !$OMP DO SCHEDULE( STATIC )
-             DO  jc = nys_mg(l), nyn_mg(l)
-!
-!--             Calculation for the first point along k
-                j  = 2*jc
-                k   = ind_even_odd+1
+!
+!--    No wall treatment
+       !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1)
+       DO  ic = nxl_mg(l), nxr_mg(l)
+          i = 2*ic
+          !$OMP DO SCHEDULE( STATIC )
+          DO  jc = nys_mg(l), nyn_mg(l)
+!
+!--          Calculation for the first point along k
+             j  = 2*jc
+!
+!--          Calculation for the other points along k
+             !DIR$ IVDEP
+             DO k = ind_even_odd+1, nzt_mg(l+1)    ! Fine grid at this point
+                km1 = k-ind_even_odd-1
                 kp1 = k-ind_even_odd
-                kc  = k-ind_even_odd
-                f_mg(kc,jc,ic) =  1.0_wp / 64.0_wp * (                         &
-                                  8.0_wp * r(k,j,i)                            &
-                               +  4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
-                                             r(k,j+1,i)     + r(k,j-1,i)     ) &
-                               +  2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
-                                             r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
-                               + 16.0_wp * r(k,j,i)                            &
-                               +  4.0_wp * r(kp1,j,i)                          &
-                               +  2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
-                                             r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
-                               +           ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
-                                             r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
-                                                      )
-!
-!--             Calculation for the other points along k
-                !DIR$ IVDEP
-                DO k = ind_even_odd+2, nzt_mg(l+1)    ! Fine grid at this point
-                   km1 = k-ind_even_odd-1
-                   kp1 = k-ind_even_odd
-                   kc  = k-ind_even_odd               ! Coarse grid index
-
-                   f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
-                                  8.0_wp * r(k,j,i)                            &
-                                + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
-                                             r(k,j+1,i)     + r(k,j-1,i)     ) &
-                                + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
-                                             r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
-                                + 4.0_wp * r(km1,j,i)                          &
-                                + 2.0_wp * ( r(km1,j,i-1)   + r(km1,j,i+1)   + &
-                                             r(km1,j+1,i)   + r(km1,j-1,i)   ) &
-                                +          ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + &
-                                             r(km1,j-1,i+1) + r(km1,j+1,i+1) ) &
-                                + 4.0_wp * r(kp1,j,i)                          &
-                                + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
-                                             r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
-                                +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
-                                             r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
-                                           )
-                ENDDO
-             ENDDO
-             !$OMP ENDDO nowait
-          ENDDO
-          !$OMP END PARALLEL
-
-       ELSE
-!
-!--       Choose flag array of the upper level
-          SELECT CASE ( l+1 )
-             CASE ( 1 )
-                flags => wall_flags_1
-             CASE ( 2 )
-                flags => wall_flags_2
-             CASE ( 3 )
-                flags => wall_flags_3
-             CASE ( 4 )
-                flags => wall_flags_4
-             CASE ( 5 )
-                flags => wall_flags_5
-             CASE ( 6 )
-                flags => wall_flags_6
-             CASE ( 7 )
-                flags => wall_flags_7
-             CASE ( 8 )
-                flags => wall_flags_8
-             CASE ( 9 )
-                flags => wall_flags_9
-             CASE ( 10 )
-                flags => wall_flags_10
-          END SELECT
-    
-          !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1,rkjim,rkjip,rkjpi,    &
-          !$OMP                   rkjmim,rkjpim,rkjmip,rkjpip,rkmji,rkmjim,    &
-          !$OMP                   rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim,        &
-          !$OMP                   rkmjmip,rkmjpip)
-          !$OMP DO
-          DO  ic = nxl_mg(l), nxr_mg(l)
-             i = 2*ic
-             DO  jc = nys_mg(l), nyn_mg(l)
-                j = 2*jc
-                !DIR$ IVDEP
-                DO k = ind_even_odd+1, nzt_mg(l+1)      !fine grid at this point
-
-                    km1 = k-ind_even_odd-1
-                    kp1 = k-ind_even_odd
-                    kc  = k-ind_even_odd                !Coarse grid index
-!
-!--                 Use implicit Neumann BCs if the respective gridpoint is
-!--                 inside the building
-                    rkjim   = MERGE(r(k,j,i),r(k,j,i-1),BTEST( flags(k,j,i-1), 6))
-                    rkjip   = MERGE(r(k,j,i),r(k,j,i+1),BTEST( flags(k,j,i+1), 6))
-                    rkjpi   = MERGE(r(k,j,i),r(k,j+1,i),BTEST( flags(k,j+1,i), 6))
-                    rkjmi   = MERGE(r(k,j,i),r(k,j-1,i),BTEST( flags(k,j-1,i), 6))
-
-                    rkjmim  = MERGE(r(k,j,i),r(k,j-1,i-1),BTEST( flags(k,j-1,i-1), 6))
-                    rkjpim  = MERGE(r(k,j,i),r(k,j+1,i-1),BTEST( flags(k,j+1,i-1), 6))
-                    rkjmip  = MERGE(r(k,j,i),r(k,j-1,i+1),BTEST( flags(k,j-1,i+1), 6))
-                    rkjpip  = MERGE(r(k,j,i),r(k,j+1,i+1),BTEST( flags(k,j+1,i+1), 6))
-
-                    rkmji   = MERGE(r(k,j,i),r(km1,j,i)  ,BTEST( flags(km1,j,i)  , 6))
-                    rkmjim  = MERGE(r(k,j,i),r(km1,j,i-1),BTEST( flags(km1,j,i-1), 6))
-                    rkmjip  = MERGE(r(k,j,i),r(km1,j,i+1),BTEST( flags(km1,j,i+1), 6))
-                    rkmjpi  = MERGE(r(k,j,i),r(km1,j+1,i),BTEST( flags(km1,j+1,i), 6))
-                    rkmjmi  = MERGE(r(k,j,i),r(km1,j-1,i),BTEST( flags(km1,j-1,i), 6))
-
-                    rkmjmim = MERGE(r(k,j,i),r(km1,j-1,i-1),BTEST( flags(km1,j-1,i-1), 6))
-                    rkmjpim = MERGE(r(k,j,i),r(km1,j+1,i-1),BTEST( flags(km1,j+1,i-1), 6))
-                    rkmjmip = MERGE(r(k,j,i),r(km1,j-1,i+1),BTEST( flags(km1,j-1,i+1), 6))
-                    rkmjpip = MERGE(r(k,j,i),r(km1,j+1,i+1),BTEST( flags(km1,j+1,i+1), 6))
-
-                    f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
-                           8.0_wp * r(k,j,i)                                   &
-                         + 4.0_wp * ( rkjim   + rkjip   +  rkjpi   + rkjmi   ) &
-                         + 2.0_wp * ( rkjmim  + rkjpim  +  rkjmip  + rkjpip  ) &
-                         + 4.0_wp * rkmji                                      &
-                         + 2.0_wp * ( rkmjim  + rkmjip  +  rkmjpi  + rkmjmi  ) &
-                         +          ( rkmjmim + rkmjpim +  rkmjmip + rkmjpip ) &
-                         + 4.0_wp * r(kp1,j,i)                                 &
-                         + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   +        &
-                                      r(kp1,j+1,i)   + r(kp1,j-1,i)   )        &
-                         +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) +        &
-                                      r(kp1,j-1,i+1) + r(kp1,j+1,i+1) )        &
-                                                        )
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-
-       ENDIF
-
+                kc  = k-ind_even_odd               ! Coarse grid index
+
+                f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
+                               8.0_wp * r(k,j,i)                            &
+                             + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
+                                          r(k,j+1,i)     + r(k,j-1,i)     ) &
+                             + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
+                                          r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
+                             + 4.0_wp * r(km1,j,i)                          &
+                             + 2.0_wp * ( r(km1,j,i-1)   + r(km1,j,i+1)   + &
+                                          r(km1,j+1,i)   + r(km1,j-1,i)   ) &
+                             +          ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + &
+                                          r(km1,j-1,i+1) + r(km1,j+1,i+1) ) &
+                             + 4.0_wp * r(kp1,j,i)                          &
+                             + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
+                                          r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
+                             +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
+                                          r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
+                                        )
+             ENDDO
+          ENDDO
+          !$OMP ENDDO nowait
+       ENDDO
+       !$OMP END PARALLEL
+
+!
+!--    Ghost point exchange
+       CALL exchange_horiz( f_mg, 1)
 !
 !--    Horizontal boundary conditions
-       CALL exchange_horiz( f_mg, 1)
-
        IF ( .NOT. bc_lr_cyc )  THEN
           IF ( inflow_l .OR. outflow_l .OR. nest_bound_l )  THEN
@@ -672,5 +475,7 @@
 !--    Boundary conditions at bottom and top of the domain.
 !--    These points are not handled by the above loop. Points may be within
-!--    buildings, but that doesn't matter.
+!--    buildings, but that doesn't matter. Remark: f_mg is ordered sequentielly
+!--    after interpolation on coarse grid (is ordered in odd-even blocks further
+!--    below). 
        IF ( ibc_p_b == 1 )  THEN
           f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
@@ -687,5 +492,6 @@
        CALL cpu_log( log_point_s(54), 'restrict', 'stop' )
 !
-!--    Why do we need a sorting here?
+!--    Since residual is in sequential order after interpolation, an additional 
+!--    sorting in odd-even blocks along z dimension is required at this point.
        CALL sort_k_to_even_odd_blocks( f_mg , l)
 
@@ -713,19 +519,21 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i   !<
-       INTEGER(iwp) ::  j   !<
-       INTEGER(iwp) ::  k   !<
-       INTEGER(iwp) ::  l   !<
-       INTEGER(iwp) ::  kp1 !<
-       INTEGER(iwp) ::  ke  !< Index for prolog even
-       INTEGER(iwp) ::  ko  !< Index for prolog odd
+       INTEGER(iwp) ::  i   !< index variable along x on coarser grid level
+       INTEGER(iwp) ::  j   !< index variable along y on coarser grid level
+       INTEGER(iwp) ::  k   !< index variable along z on coarser grid level
+       INTEGER(iwp) ::  l   !< index indicating finer grid level
+       INTEGER(iwp) ::  kp1 !< index variable along z
+       INTEGER(iwp) ::  ke  !< index for prolog even
+       INTEGER(iwp) ::  ko  !< index for prolog odd
 
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
                            nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
-                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p  !<
+                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  &
+                               p     !< perturbation pressure on coarser grid level
 
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                            nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
-                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  temp  !<
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
+                               temp  !< perturbation pressure on finer grid level
 
 
@@ -772,4 +580,5 @@
                                                  p(kp1,j,i)   + p(kp1,j,i+1) + &
                                                  p(kp1,j+1,i) + p(kp1,j+1,i+1) )
+
              ENDDO
 
@@ -804,4 +613,5 @@
                                                  p(kp1,j,i)   + p(kp1,j,i+1) + &
                                                  p(kp1,j+1,i) + p(kp1,j+1,i+1) )
+
              ENDDO
 
@@ -836,4 +646,6 @@
 !--    Bottom and top boundary conditions
        IF ( ibc_p_b == 1 )  THEN
+!
+!--       equivalent to temp(nzb,:,: ) = temp(nzb+1,:,:)
           temp(nzb,:,: ) = temp(ind_even_odd+1,:,:)
        ELSE
@@ -842,4 +654,6 @@
 
        IF ( ibc_p_t == 1 )  THEN
+!
+!--       equivalent to temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
           temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:)
        ELSE
@@ -866,7 +680,7 @@
        USE control_parameters,                                                 &
            ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
-                  inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l,  &
-                  nest_bound_n, nest_bound_r, nest_bound_s, ngsrb,             &
-                  outflow_l, outflow_n, outflow_r, outflow_s, topography
+                  inflow_n, inflow_r, inflow_s, nest_bound_l, nest_bound_n,    &
+                  nest_bound_r, nest_bound_s, ngsrb, outflow_l, outflow_n,     &
+                  outflow_r, outflow_s
 
        USE grid_variables,                                                     &
@@ -874,71 +688,34 @@
 
        USE indices,                                                            &
-           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
-                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
-                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
-                  nys_mg, nyn_mg, nzb, nzt_mg
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
 
        IMPLICIT NONE
 
-       INTEGER(iwp) :: color    !<
-       INTEGER(iwp) :: i        !<
-       INTEGER(iwp) :: ic       !<
-       INTEGER(iwp) :: j        !<
-       INTEGER(iwp) :: jc       !<
-       INTEGER(iwp) :: jj       !<
-       INTEGER(iwp) :: k        !<
-       INTEGER(iwp) :: l        !<
-       INTEGER(iwp) :: n        !<
-       INTEGER(iwp) :: km1      !<
-       INTEGER(iwp) :: kp1      !<
-
-       LOGICAL :: unroll        !<
-
-       REAL(wp) ::  wall_left   !<
-       REAL(wp) ::  wall_north  !<
-       REAL(wp) ::  wall_right  !<
-       REAL(wp) ::  wall_south  !<
-       REAL(wp) ::  wall_total  !<
-       REAL(wp) ::  wall_top    !<
+       INTEGER(iwp) :: color    !< grid point color, either red or black
+       INTEGER(iwp) :: i        !< index variable along x
+       INTEGER(iwp) :: ic       !< index variable along x
+       INTEGER(iwp) :: j        !< index variable along y
+       INTEGER(iwp) :: jc       !< index variable along y
+       INTEGER(iwp) :: jj       !< index variable along y
+       INTEGER(iwp) :: k        !< index variable along z
+       INTEGER(iwp) :: l        !< grid level
+       INTEGER(iwp) :: n        !< loop variable GauÃ-Seidel iterations
+       INTEGER(iwp) :: km1      !< index variable (k-1)
+       INTEGER(iwp) :: kp1      !< index variable (k+1)
+
+       LOGICAL      :: unroll   !< flag indicating whether loop unrolling is possible
 
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                            nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
-                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
+                                      f_mg  !< residual of perturbation pressure 
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                            nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
-                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
+                                      p_mg  !< perturbation pressure
 
        l = grid_level
 
-!
-!--    Choose flag array of this level
-       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
-
-          SELECT CASE ( l )
-             CASE ( 1 )
-                flags => wall_flags_1
-             CASE ( 2 )
-                flags => wall_flags_2
-             CASE ( 3 )
-                flags => wall_flags_3
-             CASE ( 4 )
-                flags => wall_flags_4
-             CASE ( 5 )
-                flags => wall_flags_5
-             CASE ( 6 )
-                flags => wall_flags_6
-             CASE ( 7 )
-                flags => wall_flags_7
-             CASE ( 8 )
-                flags => wall_flags_8
-             CASE ( 9 )
-                flags => wall_flags_9
-             CASE ( 10 )
-                flags => wall_flags_10
-          END SELECT
-
-       ENDIF
-
-       unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
+       unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND.                  &
                   MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
 
@@ -946,61 +723,20 @@
        
           DO  color = 1, 2
-!
-!--          No wall treatment in case of flat topography
-             IF ( topography == 'flat'  .OR.  masking_method )  THEN
-
-                IF ( .NOT. unroll )  THEN
-
-                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' )
-!
-!--                Without unrolling of loops, no cache optimization
-                   !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
-                   !$OMP DO
-                   DO  i = nxl_mg(l), nxr_mg(l), 2
-                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
-                         !DIR$ IVDEP
-                         DO  k = ind_even_odd+1, nzt_mg(l)
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
-                               + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
-                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
-                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-                      ENDDO
-                   ENDDO
-   
-                   !$OMP DO
-                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
-                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
-                          !DIR$ IVDEP
-                          DO  k = ind_even_odd+1, nzt_mg(l)
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
-                               + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
-                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
-                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-                      ENDDO
-                   ENDDO
-  
-                   !$OMP DO
-                   DO  i = nxl_mg(l), nxr_mg(l), 2
-                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
+
+             IF ( .NOT. unroll )  THEN
+
+                CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' )
+!
+!--             Without unrolling of loops, no cache optimization
+                !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
+                !$OMP DO
+                DO  i = nxl_mg(l), nxr_mg(l), 2
+                   DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+                      !DIR$ IVDEP
+                      DO  k = ind_even_odd+1, nzt_mg(l)
+                         km1 = k-ind_even_odd-1
+                         kp1 = k-ind_even_odd
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
@@ -1008,18 +744,18 @@
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
                                + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
+                               - f_mg(k,j,i)                   )
                       ENDDO
                    ENDDO
-
-                   !$OMP DO
-                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
-                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
+                ENDDO
+   
+                !$OMP DO
+                DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+                   DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+                       !DIR$ IVDEP
+                       DO  k = ind_even_odd+1, nzt_mg(l)
+                         km1 = k-ind_even_odd-1
+                         kp1 = k-ind_even_odd
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
@@ -1027,30 +763,18 @@
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
                                + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
+                               - f_mg(k,j,i)                   )
                       ENDDO
                    ENDDO
-                   !$OMP END PARALLEL
-
-                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' )
-
-                ELSE
-!
-!--                Loop unrolling along y, only one i loop for better cache use
-                   CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' )
-
-                   !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
-                   !$OMP DO
-                   DO  ic = nxl_mg(l), nxr_mg(l), 2
-                      DO  jc = nys_mg(l), nyn_mg(l), 4
-                         i  = ic
-                         jj = jc+2-color
-                         !DIR$ IVDEP
-                         DO  k = ind_even_odd+1, nzt_mg(l)
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-                            j   = jj
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
+                ENDDO
+  
+                !$OMP DO
+                DO  i = nxl_mg(l), nxr_mg(l), 2
+                   DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
+                      !DIR$ IVDEP
+                      DO  k = nzb+1, ind_even_odd
+                         km1 = k+ind_even_odd
+                         kp1 = k+ind_even_odd+1
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
@@ -1058,8 +782,18 @@
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
                                + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
+                               - f_mg(k,j,i)                   )
+                      ENDDO
+                   ENDDO
+                ENDDO
+
+                !$OMP DO
+                DO  i = nxl_mg(l)+1, nxr_mg(l), 2
+                   DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
+                      !DIR$ IVDEP
+                      DO  k = nzb+1, ind_even_odd
+                         km1 = k+ind_even_odd
+                         kp1 = k+ind_even_odd+1
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
@@ -1067,51 +801,30 @@
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
                                + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-
-                         i  = ic+1
-                         jj = jc+color-1
-                         !DIR$ IVDEP
-                         DO  k = ind_even_odd+1, nzt_mg(l)
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-                            j   = jj
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
-                               + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
-                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
-                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
-                               + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
-                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
-                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-
-                         i  = ic
-                         jj = jc+color-1
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-                            j   = jj
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
-                               + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
-                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
-                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
+                               - f_mg(k,j,i)                   )
+                      ENDDO
+                   ENDDO
+                ENDDO
+                !$OMP END PARALLEL
+
+                CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' )
+
+             ELSE
+!
+!--              Loop unrolling along y, only one i loop for better cache use
+                CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' )
+
+                !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
+                !$OMP DO
+                DO  ic = nxl_mg(l), nxr_mg(l), 2
+                   DO  jc = nys_mg(l), nyn_mg(l), 4
+                      i  = ic
+                      jj = jc+2-color
+                      !DIR$ IVDEP
+                      DO  k = ind_even_odd+1, nzt_mg(l)
+                         km1 = k-ind_even_odd-1
+                         kp1 = k-ind_even_odd
+                         j   = jj
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
@@ -1119,324 +832,99 @@
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
                                + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-
-                         i  = ic+1
-                         jj = jc+2-color
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-                            j   = jj
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
+                               - f_mg(k,j,i)                   )
+                         j = jj+2
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
                                ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
                                + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
-                               ddx2_mg(l) *                                    &
+                               - f_mg(k,j,i)                   )
+                      ENDDO
+
+                      i  = ic+1
+                      jj = jc+color-1
+                      !DIR$ IVDEP
+                      DO  k = ind_even_odd+1, nzt_mg(l)
+                         km1 = k-ind_even_odd-1
+                         kp1 = k-ind_even_odd
+                         j   = jj
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
                                ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
                                + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
                                + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-
+                               - f_mg(k,j,i)                   )
+                         j = jj+2
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)                   )
                       ENDDO
-                   ENDDO
-                   !$OMP END PARALLEL
-
-                   CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' )
-
-                ENDIF
-
-             ELSE
-
-                IF ( .NOT. unroll )  THEN
-
-                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
-!
-!--                Without unrolling of loops, no cache optimization /
-!--                vectorization. Therefore, the non-vectorized IBITS function
-!--                is still used.
-                   !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
-                   !$OMP DO
-                   DO  i = nxl_mg(l), nxr_mg(l), 2
-                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
-                         !DIR$ IVDEP
-                         DO  k = ind_even_odd+1, nzt_mg(l)
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-
-                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
-                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
+
+                      i  = ic
+                      jj = jc+color-1
+                      !DIR$ IVDEP
+                      DO  k = nzb+1, ind_even_odd
+                         km1 = k+ind_even_odd
+                         kp1 = k+ind_even_odd+1
+                         j   = jj
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
-                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
-                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
-                               + f3_mg(k,l) *                                  &
-                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
-                               - f_mg(k,j,i)               )
-                         ENDDO
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)                   )
+                         j = jj+2
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)                   )
                       ENDDO
-                   ENDDO
-   
-                   !$OMP DO
-                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
-                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
-                          !DIR$ IVDEP
-                          DO  k = ind_even_odd+1, nzt_mg(l)
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-
-                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
-                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
+
+                      i  = ic+1
+                      jj = jc+2-color
+                      !DIR$ IVDEP
+                      DO  k = nzb+1, ind_even_odd
+                         km1 = k+ind_even_odd
+                         kp1 = k+ind_even_odd+1
+                         j   = jj
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
                                + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
-                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
-                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
-                               + f3_mg(k,l) *                                  &
-                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
-                               - f_mg(k,j,i)               )
-                         ENDDO
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)                   )
+                         j = jj+2
+                         p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (               &
+                                 ddx2_mg(l) *                                  &
+                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
+                               + ddy2_mg(l) *                                  &
+                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
+                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
+                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
+                               - f_mg(k,j,i)                   )
                       ENDDO
-                   ENDDO
-
-                   !$OMP DO
-                   DO  i = nxl_mg(l), nxr_mg(l), 2
-                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
-                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
-                               + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
-                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
-                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
-                               + f3_mg(k,l) *                                  &
-                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) ) &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-                      ENDDO
-                   ENDDO
-
-                   !$OMP DO
-                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
-                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-
-                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
-                               ddx2_mg(l) *                                    &
-                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
-                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
-                               + ddy2_mg(l) *                                  &
-                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
-                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
-                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
-                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
-                               + f3_mg(k,l) *                                  &
-                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
-                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
-                               - f_mg(k,j,i)               )
-                         ENDDO
-                      ENDDO
-                   ENDDO
-                   !$OMP END PARALLEL
-
-                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
-
-                ELSE
-!
-!--                Loop unrolling along y, only one i loop for better cache use
-                   CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
-
-                   !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
-                   !$OMP DO
-                   DO  ic = nxl_mg(l), nxr_mg(l), 2
-                      DO  jc = nys_mg(l), nyn_mg(l), 4
-                         i  = ic
-                         jj = jc+2-color
-                         !DIR$ IVDEP
-                         DO  k = ind_even_odd+1, nzt_mg(l)
-
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-
-                            j = jj
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-
-                         ENDDO
-
-                         i  = ic+1
-                         jj = jc+color-1
-                         !DIR$ IVDEP
-                         DO  k = ind_even_odd+1, nzt_mg(l)
-
-                            km1 = k-ind_even_odd-1
-                            kp1 = k-ind_even_odd
-
-                            j =jj
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-
-                         ENDDO
-
-                         i  = ic
-                         jj = jc+color-1
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-
-                            j =jj
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-
-                         ENDDO
-
-                         i  = ic+1
-                         jj = jc+2-color
-                         !DIR$ IVDEP
-                         DO  k = nzb+1, ind_even_odd
-
-                            km1 = k+ind_even_odd
-                            kp1 = k+ind_even_odd+1
-
-                            j =jj
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-
-                            j = jj+2
-                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
-                                ddx2_mg(l) *                                                 &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
-                                + ddy2_mg(l) *                                               &
-                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
-                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
-                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
-                                + f3_mg(k,l) *                                               &
-                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
-                                - f_mg(k,j,i)               )
-
-                         ENDDO
-
-                      ENDDO
-                   ENDDO
-                   !$OMP END PARALLEL
-
-                   CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
-
-                ENDIF
+
+                   ENDDO
+                ENDDO
+                !$OMP END PARALLEL
+
+                CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' )
 
              ENDIF
@@ -1467,4 +955,6 @@
 !--          Bottom and top boundary conditions
              IF ( ibc_p_b == 1 )  THEN
+!
+!--             equivalent to p_mg(nzb,:,: ) = p_mg(nzb+1,:,:) 
                 p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:)
              ELSE
@@ -1473,4 +963,6 @@
 
              IF ( ibc_p_t == 1 )  THEN
+!
+!--             equivalent to p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
                 p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:)
              ELSE
@@ -1479,84 +971,6 @@
 
           ENDDO
+
        ENDDO
-!
-!--    Set pressure within topography and at the topography surfaces
-       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
-
-          !$OMP PARALLEL PRIVATE (i,j,k,kp1,wall_left,wall_north,wall_right,   &
-          !$OMP                   wall_south,wall_top,wall_total)
-          !$OMP DO
-          DO  i = nxl_mg(l), nxr_mg(l)
-             DO  j = nys_mg(l), nyn_mg(l)
-                !DIR$ IVDEP
-                DO  k = ind_even_odd+1, nzt_mg(l)
-                   km1 = k-ind_even_odd-1
-                   kp1 = k-ind_even_odd
-!
-!--                First, set pressure inside topography to zero
-                   p_mg(k,j,i) = p_mg(k,j,i) *                                 &
-                                 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
-!
-!--                Second, determine if the gridpoint inside topography is
-!--                adjacent to a wall and set its value to a value given by the
-!--                average of those values obtained from Neumann boundary
-!--                condition
-                   wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
-                   wall_right = IBITS( flags(k,j,i+1), 4, 1 )
-                   wall_south = IBITS( flags(k,j-1,i), 3, 1 )
-                   wall_north = IBITS( flags(k,j+1,i), 2, 1 )
-                   wall_top   = IBITS( flags(kp1,j,i), 0, 1 )
-                   wall_total = wall_left + wall_right + wall_south +          &
-                                wall_north + wall_top
-
-                   IF ( wall_total > 0.0_wp )  THEN
-                      p_mg(k,j,i) = 1.0_wp / wall_total *                      &
-                                    ( wall_left  * p_mg(k,j,i-1) +             &
-                                      wall_right * p_mg(k,j,i+1) +             &
-                                      wall_south * p_mg(k,j-1,i) +             &
-                                      wall_north * p_mg(k,j+1,i) +             &
-                                      wall_top   * p_mg(kp1,j,i) )
-                   ENDIF
-                ENDDO
-
-                !DIR$ IVDEP
-                DO  k = nzb+1, ind_even_odd
-                   km1 = k+ind_even_odd
-                   kp1 = k+ind_even_odd+1
-!
-!--                First, set pressure inside topography to zero
-                   p_mg(k,j,i) = p_mg(k,j,i) *                                 &
-                                 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
-!
-!--                Second, determine if the gridpoint inside topography is
-!--                adjacent to a wall and set its value to a value given by the
-!--                average of those values obtained from Neumann boundary
-!--                condition
-                   wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
-                   wall_right = IBITS( flags(k,j,i+1), 4, 1 )
-                   wall_south = IBITS( flags(k,j-1,i), 3, 1 )
-                   wall_north = IBITS( flags(k,j+1,i), 2, 1 )
-                   wall_top   = IBITS( flags(kp1,j,i), 0, 1 )
-                   wall_total = wall_left + wall_right + wall_south +          &
-                                wall_north + wall_top
-
-                   IF ( wall_total > 0.0_wp )  THEN
-                      p_mg(k,j,i) = 1.0_wp / wall_total *                      &
-                                    ( wall_left  * p_mg(k,j,i-1) +             &
-                                      wall_right * p_mg(k,j,i+1) +             &
-                                      wall_south * p_mg(k,j-1,i) +             &
-                                      wall_north * p_mg(k,j+1,i) +             &
-                                      wall_top   * p_mg(kp1,j,i) )
-                   ENDIF
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-
-!
-!--       One more time horizontal boundary conditions
-          CALL exchange_horiz( p_mg, 1)
-
-       ENDIF
 
     END SUBROUTINE redblack_fast
@@ -1581,17 +995,18 @@
        IMPLICIT NONE
 
-       INTEGER(iwp), INTENT(IN) ::  glevel
+       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level 
 
        REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1,                               &
                            nys_mg(glevel)-1:nyn_mg(glevel)+1,                  &
-                           nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::  p_mg  !<
+                           nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::               &
+                                      p_mg  !< array to be sorted
 !
 !--    Local variables
-       INTEGER(iwp) :: i        !<
-       INTEGER(iwp) :: j        !<
-       INTEGER(iwp) :: k        !<
-       INTEGER(iwp) :: l        !<
-       INTEGER(iwp) :: ind      !<
-       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !<
+       INTEGER(iwp) :: i        !< index variable along x
+       INTEGER(iwp) :: j        !< index variable along y
+       INTEGER(iwp) :: k        !< index variable along z
+       INTEGER(iwp) :: l        !< grid level
+       INTEGER(iwp) :: ind      !< index variable along z
+       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< odd-even sorted temporary array
 
 
@@ -1646,13 +1061,13 @@
        IMPLICIT NONE
 
-       INTEGER(iwp), INTENT(IN) ::  glevel
-
-       REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) ::  f_mg
-       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  f_mg_b
+       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
+
+       REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) ::  f_mg    !< 1D input array
+       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  f_mg_b  !< 1D output array
 
 !
 !--    Local variables
-       INTEGER(iwp) :: ind   !<
-       INTEGER(iwp) :: k     !<
+       INTEGER(iwp) :: ind   !< index variable along z
+       INTEGER(iwp) :: k     !< index variable along z 
 
 
@@ -1696,17 +1111,18 @@
        IMPLICIT NONE
 
-       INTEGER(iwp), INTENT(IN) ::  glevel
+       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
 
        INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
                                nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
-                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::  i_mg  !<
+                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::           &
+                                    i_mg    !< array to be sorted 
 !
 !--    Local variables
-       INTEGER(iwp) :: i        !<
-       INTEGER(iwp) :: j        !<
-       INTEGER(iwp) :: k        !<
-       INTEGER(iwp) :: l        !<
-       INTEGER(iwp) :: ind      !<
-       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp
+       INTEGER(iwp) :: i        !< index variabel along x
+       INTEGER(iwp) :: j        !< index variable along y
+       INTEGER(iwp) :: k        !< index variable along z
+       INTEGER(iwp) :: l        !< grid level 
+       INTEGER(iwp) :: ind      !< index variable along z
+       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< temporary odd-even sorted array
 
 
@@ -1770,12 +1186,13 @@
        REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
                            nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
-                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
+                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::       &
+                                     p_mg  !< array to be sorted
 !
 !--    Local variables
-       INTEGER(iwp) :: i        !<
-       INTEGER(iwp) :: j        !<
-       INTEGER(iwp) :: k        !<
-       INTEGER(iwp) :: l        !<
-       INTEGER(iwp) :: ind      !<
+       INTEGER(iwp) :: i        !< index variable along x
+       INTEGER(iwp) :: j        !< index variable along y
+       INTEGER(iwp) :: k        !< index variable along z
+       INTEGER(iwp) :: l        !< grid level
+       INTEGER(iwp) :: ind      !< index variable along z
 
        REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) ::  tmp
@@ -1958,7 +1375,7 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i            !<
-       INTEGER(iwp) ::  j            !<
-       INTEGER(iwp) ::  k            !<
+       INTEGER(iwp) ::  i            !< index variable along x
+       INTEGER(iwp) ::  j            !< index variable along y
+       INTEGER(iwp) ::  k            !< index variable along z
        INTEGER(iwp) ::  nxl_mg_save  !<
        INTEGER(iwp) ::  nxr_mg_save  !<
@@ -2308,5 +1725,5 @@
 
        USE control_parameters,                                                 &
-           ONLY:  grid_level, masking_method, maximum_grid_level, topography
+           ONLY:  grid_level, maximum_grid_level
 
        USE indices,                                                            &
@@ -2314,13 +1731,10 @@
 
        USE indices,                                                            &
-           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
-                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
-                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
-                  nys_mg, nyn_mg, nzb, nzt_mg
+           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
 
        IMPLICIT NONE
 !
 !--    Local variables
-       INTEGER(iwp) ::  i     !<
+       INTEGER(iwp) ::  i     !<  
        INTEGER(iwp) ::  l     !<
 
@@ -2357,40 +1771,4 @@
        ENDDO
 
-!
-!--    Sort flags array for all levels to block (even/odd) structure
-       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
-
-          DO  l = maximum_grid_level, 1 , -1
-!
-!--          Assign the flag level to be calculated
-             SELECT CASE ( l )
-                CASE ( 1 )
-                   flags => wall_flags_1
-                CASE ( 2 )
-                   flags => wall_flags_2
-                CASE ( 3 )
-                   flags => wall_flags_3
-                CASE ( 4 )
-                   flags => wall_flags_4
-                CASE ( 5 )
-                   flags => wall_flags_5
-                CASE ( 6 )
-                   flags => wall_flags_6
-                CASE ( 7 )
-                   flags => wall_flags_7
-                CASE ( 8 )
-                   flags => wall_flags_8
-                CASE ( 9 )
-                   flags => wall_flags_9
-                CASE ( 10 )
-                   flags => wall_flags_10
-             END SELECT
-
-             CALL sort_k_to_even_odd_blocks( flags, l )
-
-          ENDDO
-
-       ENDIF
-
        lfirst = .FALSE.
 
