Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 110)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 114)
@@ -1,4 +1,11 @@
 New:
 ---
+
+Pressure boundary conditions for vertical walls added to the multigrid solver.
+They are applied using new wall flag arrays (wall_flags_..) which are defined
+for each grid level. New argument gls added to routine user_init_grid
+(user_interface).
+
+init_grid, init_pegrid, modules, user_interface
 
 
@@ -9,2 +16,9 @@
 Errors:
 ------
+
+Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
+numbers along y (advec_particles)
+
+Bugfix: model_string needed a default value (combine_plot_fields)
+
+advec_particles, combine_plot_fields
Index: palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- palm/trunk/SOURCE/advec_particles.f90	(revision 110)
+++ palm/trunk/SOURCE/advec_particles.f90	(revision 114)
@@ -4,4 +4,6 @@
 ! Actual revisions:
 ! -----------------
+! Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
+! numbers along y
 ! TEST: PRINT statements on unit 9 (commented out)
 !
@@ -2773,6 +2775,6 @@
           IF ( use_particle_tails )  THEN
 
-             CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, pleft,  0, &
-                                trnpt_count_recv, 1, MPI_INTEGER, pright, 0, &
+             CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
+                                trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
                                 comm2d, status, ierr )
 
@@ -2850,6 +2852,6 @@
           IF ( use_particle_tails )  THEN
 
-             CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pright, 0, &
-                                trspt_count_recv, 1, MPI_INTEGER, pleft,  0, &
+             CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
+                                trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
                                 comm2d, status, ierr )
 
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 110)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 114)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Multigrid solver allows topography
 !
 ! Former revisions:
@@ -286,5 +286,5 @@
           WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
        ENDIF
-       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  THEN
+       IF ( psolver == 'sor' )  THEN
           WRITE( action, '(A,A)' )  'psolver = ', psolver
        ENDIF
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 110)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 114)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation of wall flag arrays
 !
 ! Former revisions:
@@ -44,6 +44,6 @@
     IMPLICIT NONE
 
-    INTEGER ::  bh, blx, bly, bxl, bxr, byn, bys, i, i_center, j, j_center, k, &
-                l, nzb_si, nzt_l, vi
+    INTEGER ::  bh, blx, bly, bxl, bxr, byn, bys, gls, i, inc, i_center, j, &
+                j_center, k, l, nxl_l, nxr_l, nyn_l, nys_l, nzb_si, nzt_l, vi
 
     INTEGER, DIMENSION(:), ALLOCATABLE   ::  vertical_influence
@@ -217,9 +217,12 @@
 !
 !-- Allocate outer and inner index arrays for topography and set
-!-- defaults
-    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), &
-              corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), &
-              nzb_local(-1:ny+1,-1:nx+1), nzb_tmp(-1:ny+1,-1:nx+1),   &
-              wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),       &
+!-- defaults.
+!-- nzb_local has to contain additional layers of ghost points for calculating
+!-- the flag arrays needed for the multigrid method
+    gls = 2**( maximum_grid_level )
+    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr),       &
+              corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr),       &
+              nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &
+              wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),             &
               wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
     ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), &
@@ -388,8 +391,11 @@
              ENDDO
           ENDDO
-          nzb_local(-1,0:nx)   = nzb_local(ny,0:nx)
-          nzb_local(ny+1,0:nx) = nzb_local(0,0:nx)
-          nzb_local(:,-1)      = nzb_local(:,nx)
-          nzb_local(:,nx+1)    = nzb_local(:,0)
+!
+!--       Add cyclic boundaries (additional layers are for calculating flag
+!--       arrays needed for the multigrid sover)
+          nzb_local(-gls:-1,0:nx)     = nzb_local(ny-gls+1:ny,0:nx)
+          nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
+          nzb_local(:,-gls:-1)        = nzb_local(:,nx-gls+1:nx)
+          nzb_local(:,nx+1:nx+gls)    = nzb_local(:,0:gls-1)
       
           GOTO 12
@@ -413,7 +419,14 @@
 !--       case in the user interface. There, the subroutine user_init_grid 
 !--       checks which of these two conditions applies.
-          CALL user_init_grid( nzb_local )
+          CALL user_init_grid( gls, nzb_local )
 
     END SELECT
+
+!
+!-- Test output of nzb_local -1:ny+1,-1:nx+1
+    WRITE (9,*)  '*** nzb_local ***'
+    DO  j = ny+1, -1, -1
+       WRITE (9,'(194(1X,I2))')  ( nzb_local(j,i), i = -1, nx+1 )
+    ENDDO
 
 !
@@ -478,5 +491,5 @@
 !--    nzb_s_outer: 
 !--    extend nzb_local east-/westwards first, then north-/southwards
-       nzb_tmp = nzb_local
+       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
        DO  j = -1, ny + 1
           DO  i = 0, nx
@@ -510,5 +523,5 @@
 !--    nzb_u_inner: 
 !--    extend nzb_local rightwards only
-       nzb_tmp = nzb_local
+       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
        DO  j = -1, ny + 1
           DO  i = 0, nx + 1
@@ -542,5 +555,5 @@
 !--    nzb_v_inner:
 !--    extend nzb_local northwards only
-       nzb_tmp = nzb_local
+       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
        DO  i = -1, nx + 1
           DO  j = 0, ny + 1
@@ -728,4 +741,118 @@
 
 !
+!-- Calculate wall flag arrays for the multigrid method
+    IF ( psolver == 'multigrid' )  THEN
+!
+!--    Gridpoint increment of the current level
+       inc = 1
+
+       DO  l = maximum_grid_level, 1 , -1
+
+          nxl_l = nxl_mg(l)
+          nxr_l = nxr_mg(l)
+          nys_l = nys_mg(l)
+          nyn_l = nyn_mg(l)
+          nzt_l = nzt_mg(l)
+
+!
+!--       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
+
+!
+!--       Depending on the grid level, set the respective bits in case of
+!--       neighbouring walls
+!--       Bit 0:  wall to the bottom
+!--       Bit 1:  wall to the top (not realized in remaining PALM code so far)
+!--       Bit 2:  wall to the south
+!--       Bit 3:  wall to the north
+!--       Bit 4:  wall to the left
+!--       Bit 5:  wall to the right
+!--       Bit 6:  inside / outside building (1/0)
+
+          flags = 0
+
+          DO  i = nxl_l-1, nxr_l+1
+             DO  j = nys_l-1, nyn_l+1
+                DO  k = nzb, nzt_l+1
+                         
+!
+!--                Inside/outside building (inside building does not need
+!--                further tests for walls)
+                   IF ( k*inc <= nzb_local(j*inc,i*inc) )  THEN
+
+                      flags(k,j,i) = IBSET( flags(k,j,i), 6 )
+
+                   ELSE
+!
+!--                   Bottom wall
+                      IF ( (k-1)*inc <= nzb_local(j*inc,i*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 0 )
+                      ENDIF
+!
+!--                   South wall
+                      IF ( k*inc <= nzb_local((j-1)*inc,i*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 2 )
+                      ENDIF
+!
+!--                   North wall
+                      IF ( k*inc <= nzb_local((j+1)*inc,i*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 3 )
+                      ENDIF
+!
+!--                   Left wall
+                      IF ( k*inc <= nzb_local(j*inc,(i-1)*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 4 )
+                      ENDIF
+!
+!--                   Right wall
+                      IF ( k*inc <= nzb_local(j*inc,(i+1)*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 5 )
+                      ENDIF
+
+                   ENDIF
+                            
+                ENDDO
+             ENDDO
+          ENDDO 
+
+!
+!--       Test output of flag arrays
+          i = nxl_l
+          WRITE (9,*)  ' '
+          WRITE (9,*)  '*** mg level ', l, ' ***', mg_switch_to_pe0_level
+          WRITE (9,*)  '    inc=', inc, '  i =', nxl_l
+          WRITE (9,*)  '    nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
+          DO  k = nzt_l+1, nzb, -1
+             WRITE (9,'(194(1X,I2))')  ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
+          ENDDO
+
+          inc = inc * 2
+
+       ENDDO
+
+    ENDIF
+
+!
 !-- In case of topography: limit near-wall mixing length l_wall further:
 !-- Go through all points of the subdomain one by one and look for the closest
Index: palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- palm/trunk/SOURCE/init_pegrid.f90	(revision 110)
+++ palm/trunk/SOURCE/init_pegrid.f90	(revision 114)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! Allocation of wall flag arrays for multigrid solver
 ! TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values
 !
@@ -915,3 +916,72 @@
     ENDIF
 
+!
+!-- Allocate wall flag arrays used in the multigrid solver
+    IF ( psolver == 'multigrid' )  THEN
+
+       DO  i = maximum_grid_level, 1, -1
+
+           SELECT CASE ( i )
+
+              CASE ( 1 )
+                 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 2 )
+                 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 3 )
+                 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 4 )
+                 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 5 )
+                 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 6 )
+                 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 7 )
+                 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 8 )
+                 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 9 )
+                 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1,         &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE ( 10 )
+                 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1,        &
+                                        nys_mg(i)-1:nyn_mg(i)+1, &
+                                        nxl_mg(i)-1:nxr_mg(i)+1) )
+
+              CASE DEFAULT
+                 IF ( myid == 0 )  PRINT*, '+++ init_pegrid: more than 10 ', &
+                                           ' multigrid levels'
+                 CALL local_stop
+
+          END SELECT
+
+       ENDDO
+
+    ENDIF
+
  END SUBROUTINE init_pegrid
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 110)
+++ palm/trunk/SOURCE/modules.f90	(revision 114)
@@ -5,5 +5,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +flags, wall_flags_1..10
 !
 ! Former revisions:
@@ -589,4 +589,11 @@
                 nzb_v_outer, nzb_w_inner, nzb_w_outer, nzb_2d
 
+    INTEGER, DIMENSION(:,:,:), POINTER ::  flags
+
+    INTEGER, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  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
+
     SAVE
 
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 110)
+++ palm/trunk/SOURCE/palm.f90	(revision 114)
@@ -66,5 +66,5 @@
     INTEGER           ::  i, run_description_header_i(80)
 
-    version = 'PALM 3.4'
+    version = 'PALM 3.4a'
 
 #if defined( __parallel )
Index: palm/trunk/SOURCE/poismg.f90
===================================================================
--- palm/trunk/SOURCE/poismg.f90	(revision 110)
+++ 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
Index: palm/trunk/SOURCE/user_interface.f90
===================================================================
--- palm/trunk/SOURCE/user_interface.f90	(revision 110)
+++ palm/trunk/SOURCE/user_interface.f90	(revision 114)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +argument gls in user_init_grid
 !
 ! Former revisions:
@@ -238,5 +238,5 @@
 
 
- SUBROUTINE user_init_grid( nzb_local )
+ SUBROUTINE user_init_grid( gls, nzb_local )
 
 !------------------------------------------------------------------------------!
@@ -245,4 +245,6 @@
 ! ------------
 ! Execution of user-defined grid initializing actions
+! First argument gls contains the number of ghost layers, which is > 1 if the
+! multigrid method for the pressure solver is used
 !------------------------------------------------------------------------------!
 
@@ -253,5 +255,7 @@
     IMPLICIT NONE
 
-    INTEGER, DIMENSION(-1:ny+1,-1:nx+1) ::  nzb_local
+    INTEGER ::  gls
+
+    INTEGER, DIMENSION(-gls:ny+gls,-gls:nx+gls) ::  nzb_local
 
 !
