Index: palm/trunk/SOURCE/poismg.f90
===================================================================
--- palm/trunk/SOURCE/poismg.f90	(revision 77)
+++ palm/trunk/SOURCE/poismg.f90	(revision 114)
@@ -8,5 +8,8 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Boundary conditions at walls are implicitly set using flag arrays. Only
+! Neumann BC is allowed. Upper walls are still not realized.
+! Bottom and top BCs for array f_mg in restrict removed because boundary
+! values are not needed (right hand side of SOR iteration)
 !
 ! Former revisions:
@@ -150,4 +153,29 @@
     l = grid_level
 
+!
+!-- 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)
 !$OMP DO
@@ -155,10 +183,23 @@
        DO  j = nys_mg(l), nyn_mg(l) 
           DO  k = nzb+1, nzt_mg(l)
-             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(k,l) * p_mg(k+1,j,i)                     &
-                        - f3_mg(k,l) * p_mg(k-1,j,i)                     &
+             r(k,j,i) = f_mg(k,j,i)                                         &
+                        - 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(k+1,j,i)                        &
+                        - f3_mg(k,l) *                                      &
+                            ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                          ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
                         + f1_mg(k,l) * p_mg(k,j,i)
+!
+!--          Residual within topography should be zero
+             r(k,j,i) = r(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) )
           ENDDO
        ENDDO
@@ -181,11 +222,6 @@
 
 !
-!-- Bottom and top boundary conditions
-    IF ( ibc_p_b == 1 )  THEN
-       r(nzb,:,: ) = r(nzb+1,:,:)
-    ELSE
-       r(nzb,:,: ) = 0.0
-    ENDIF
-
+!-- Top boundary condition
+!-- A Neumann boundary condition for r is implicitly set in routine restrict
     IF ( ibc_p_t == 1 )  THEN
        r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
@@ -217,4 +253,8 @@
     INTEGER ::  i, ic, j, jc, k, kc, l
 
+    REAL ::  rkjim, rkjip, rkjmi, rkjmim, rkjmip, rkjpi, rkjpim, rkjpip,       &
+             rkmji, rkmjim, rkmjip, rkmjmi, rkmjmim, rkmjmip, rkmjpi, rkmjpim, &
+             rkmjpip
+
     REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
                     nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
@@ -229,4 +269,29 @@
     l = grid_level
 
+!
+!-- Choose flag array of the upper 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,ic,jc,kc)
 !$OMP DO
@@ -237,15 +302,53 @@
           DO  kc = nzb+1, nzt_mg(l)
              k = 2*kc-1
+!
+!--          Use implicit Neumann BCs if the respective gridpoint is inside
+!--          the building
+             rkjim   = r(k,j,i-1)     + IBITS( flags(k,j,i-1), 6, 1 ) *     &
+                                        ( r(k,j,i) - r(k,j,i-1) )
+             rkjip   = r(k,j,i+1)     + IBITS( flags(k,j,i+1), 6, 1 ) *     &
+                                        ( r(k,j,i) - r(k,j,i+1) )
+             rkjpi   = r(k,j+1,i)     + IBITS( flags(k,j+1,i), 6, 1 ) *     &
+                                        ( r(k,j,i) - r(k,j+1,i) )
+             rkjmi   = r(k,j-1,i)     + IBITS( flags(k,j-1,i), 6, 1 ) *     &
+                                        ( r(k,j,i) - r(k,j-1,i) )
+             rkjmim  = r(k,j-1,i-1)   + IBITS( flags(k,j-1,i-1), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k,j-1,i-1) )
+             rkjpim  = r(k,j+1,i-1)   + IBITS( flags(k,j+1,i-1), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k,j+1,i-1) )
+             rkjmip  = r(k,j-1,i+1)   + IBITS( flags(k,j-1,i+1), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k,j-1,i+1) )
+             rkjpip  = r(k,j+1,i+1)   + IBITS( flags(k,j+1,i+1), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k,j+1,i+1) )
+             rkmji   = r(k-1,j,i)     + IBITS( flags(k-1,j,i), 6, 1 ) *     &
+                                        ( r(k,j,i) - r(k-1,j,i) )
+             rkmjim  = r(k-1,j,i-1)   + IBITS( flags(k-1,j,i-1), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k-1,j,i-1) )
+             rkmjip  = r(k-1,j,i+1)   + IBITS( flags(k-1,j,i+1), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k-1,j,i+1) )
+             rkmjpi  = r(k-1,j+1,i)   + IBITS( flags(k-1,j+1,i), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k-1,j+1,i) )
+             rkmjmi  = r(k-1,j-1,i)   + IBITS( flags(k-1,j-1,i), 6, 1 ) *   &
+                                        ( r(k,j,i) - r(k-1,j-1,i) )
+             rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * &
+                                        ( r(k,j,i) - r(k-1,j-1,i-1) )
+             rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * &
+                                        ( r(k,j,i) - r(k-1,j+1,i-1) )
+             rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * &
+                                        ( r(k,j,i) - r(k-1,j-1,i+1) )
+             rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * &
+                                        ( r(k,j,i) - r(k-1,j+1,i+1) )
+
              f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
                               8.0 * r(k,j,i)                            &
-                            + 4.0 * ( r(k,j,i-1)     + r(k,j,i+1)     + &
-                                      r(k,j+1,i)     + r(k,j-1,i)     ) &
-                            + 2.0 * ( 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 * r(k-1,j,i)                          &
-                            + 2.0 * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
-                                      r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
-                            +       ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
-                                      r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
+                            + 4.0 * ( rkjim   + rkjip   +               &
+                                      rkjpi   + rkjmi   )               &
+                            + 2.0 * ( rkjmim  + rkjpim  +               &
+                                      rkjmip  + rkjpip  )               &
+                            + 4.0 * rkmji                               &
+                            + 2.0 * ( rkmjim  + rkmjim  +               &
+                                      rkmjpi  + rkmjmi  )               &
+                            +       ( rkmjmim + rkmjpim +               &
+                                      rkmjmip + rkmjpip )               &
                             + 4.0 * r(k+1,j,i)                          &
                             + 2.0 * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
@@ -254,4 +357,22 @@
                                       r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
                                            ) 
+
+!             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
+!                              8.0 * r(k,j,i)                            &
+!                            + 4.0 * ( r(k,j,i-1)     + r(k,j,i+1)     + &
+!                                      r(k,j+1,i)     + r(k,j-1,i)     ) &
+!                            + 2.0 * ( 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 * r(k-1,j,i)                          &
+!                            + 2.0 * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
+!                                      r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
+!                            +       ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
+!                                      r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
+!                            + 4.0 * r(k+1,j,i)                          &
+!                            + 2.0 * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
+!                                      r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
+!                            +       ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
+!                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
+!                                           ) 
           ENDDO
        ENDDO
@@ -275,15 +396,15 @@
 !
 !-- Bottom and top boundary conditions
-    IF ( ibc_p_b == 1 )  THEN
-       f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
-    ELSE
-       f_mg(nzb,:,: ) = 0.0
-    ENDIF
-
-    IF ( ibc_p_t == 1 )  THEN
-       f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
-    ELSE
-       f_mg(nzt_mg(l)+1,:,: ) = 0.0
-    ENDIF
+!    IF ( ibc_p_b == 1 )  THEN
+!       f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
+!    ELSE
+!       f_mg(nzb,:,: ) = 0.0
+!    ENDIF
+!
+!    IF ( ibc_p_t == 1 )  THEN
+!       f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
+!    ELSE
+!       f_mg(nzt_mg(l)+1,:,: ) = 0.0
+!    ENDIF
 
 
@@ -412,4 +533,6 @@
     LOGICAL :: unroll
 
+    REAL ::  wall_left, wall_north, wall_right, wall_south, wall_total, wall_top
+
     REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                 &
                     nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                &
@@ -418,4 +541,29 @@
 
     l = grid_level
+
+!
+!-- 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
 
     unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
@@ -434,10 +582,27 @@
                 DO  j = nys_mg(l) + 2 - colour, nyn_mg(l), 2 
                    DO  k = nzb+1, nzt_mg(l), 2
+!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
+!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
+!                                                       )
+
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
                 ENDDO
@@ -448,9 +613,19 @@
                    DO  k = nzb+1, nzt_mg(l), 2 
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
                 ENDDO
@@ -461,9 +636,19 @@
                    DO  k = nzb+2, nzt_mg(l), 2
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
                 ENDDO
@@ -474,9 +659,19 @@
                    DO  k = nzb+2, nzt_mg(l), 2
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
                 ENDDO
@@ -496,30 +691,36 @@
                       j = jj
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                       j = jj+2
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
-!                      j = jj+4
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                    )
-!                      j = jj+6
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
    
@@ -529,30 +730,36 @@
                       j =jj
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                       j = jj+2
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
-!                      j = jj+4
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                       )
-!                      j = jj+6
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
 
@@ -562,30 +769,36 @@
                       j =jj
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                       j = jj+2
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
-!                      j = jj+4
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                       )
-!                      j = jj+6
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
 
@@ -595,30 +808,36 @@
                       j =jj
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                       j = jj+2
                       p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-                                                       )
-!                      j = jj+4
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                       )
-!                      j = jj+6
-!                      p_mg(k,j,i) = 1.0 / f1_mg(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(k,l) * p_mg(k+1,j,i)                     &
-!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
-!                                                       )
+                             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(k+1,j,i)                        &
+                           + f3_mg(k,l) *                                      &
+                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
+                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
+                           - f_mg(k,j,i)               )
                    ENDDO
 
@@ -669,4 +888,42 @@
     ENDDO
 
+!
+!-- Set pressure within topography and at the topography surfaces
+!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
+!$OMP DO
+    DO  i = nxl_mg(l), nxr_mg(l)
+       DO  j = nys_mg(l), nyn_mg(l) 
+          DO  k = nzb, nzt_mg(l)
+!
+!--          First, set pressure inside topography to zero
+             p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0 - 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(k+1,j,i), 0, 1 )
+             wall_total = wall_left + wall_right + wall_south + wall_north + &
+                          wall_top
+
+             IF ( wall_total > 0.0 )  THEN
+                p_mg(k,j,i) = 1.0 / 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(k+1,j,i) )
+             ENDIF
+          ENDDO
+       ENDDO
+    ENDDO
+!$OMP END PARALLEL 
+
+!
+!-- One more time horizontal boundary conditions
+    CALL exchange_horiz( p_mg )
 
  END SUBROUTINE redblack
