Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 1616)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 1628)
@@ -493,7 +493,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
-             ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
-             ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
+             ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
+             ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
+             ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
 
              v_comp                  = v(k,j,i) - v_gtrans
@@ -554,7 +554,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
-             ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
-             ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
+             ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
+             ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
+             ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
 
              u_comp                     = u(k,j,i) - u_gtrans
@@ -874,7 +874,23 @@
 !--       correction is needed to overcome numerical instabilities caused
 !--       by a not sufficient reduction of divergences near topography. 
-          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
-                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
-                        + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
+          div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
+                          - u(k,j,i)   * ( IBITS(wall_flags_0(k,j,i-1),0,1)    &
+                                         + IBITS(wall_flags_0(k,j,i-1),1,1)    &
+                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
+                                         )                                     &
+                          ) * ddx                                              &
+                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
+                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
+                                         + IBITS(wall_flags_0(k,j-1,i),4,1)    &
+                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
+                                         )                                     &
+                          ) * ddy                                              &
+                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
+                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
+                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
+                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
+                                         )                                     &      
+                          ) * ddzw(k)
+
 
           tend(k,j,i) = tend(k,j,i) - (                                       &
@@ -1247,7 +1263,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
-             ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
-             ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
+             ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
+             ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
+             ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
 
              v_comp      = v(k,j,i) + v(k,j,i-1) - gv
@@ -1305,7 +1321,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
-             ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
-             ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
+             ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
+             ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
+             ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
 
              u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
@@ -1481,8 +1497,27 @@
 !--       correction is needed to overcome numerical instabilities introduced
 !--       by a not sufficient reduction of divergences near topography. 
-          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
-               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
-               +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)   &
+          div = ( ( u_comp(k)       * ( ibit9 + ibit10 + ibit11 )             &
+                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
+                                    * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
+                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
+                                      )                                       &   
+                  ) * ddx                                                     &
+               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
+                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
+                                    * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),13,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
+                                      )                                       &
+                  ) * ddy                                                     &
+               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
+                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
+                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
+                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
+                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
+                                      )                                       &  
+                  ) * ddzw(k)   &
                 ) * 0.5_wp
+
 
           tend(k,j,i) = tend(k,j,i) - (                                       &
@@ -1702,7 +1737,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
-             ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
-             ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
+             ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
+             ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
+             ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
 
              u_comp           = u(k,j-1,i) + u(k,j,i) - gu
@@ -1760,7 +1795,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
-             ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
-             ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
+             ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
+             ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
+             ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
 
              v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
@@ -1937,8 +1972,30 @@
 !--       correction is needed to overcome numerical instabilities introduced
 !--       by a not sufficient reduction of divergences near topography. 
-          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
-               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
-               +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)   &
+          div = ( ( ( u_comp     + gu )                                       &
+                                       * ( ibit18 + ibit19 + ibit20 )         &
+                  - ( u(k,j-1,i) + u(k,j,i) )                                 &
+                                       * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
+                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
+                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
+                                         )                                    &   
+                  ) * ddx                                                     &
+               +  ( v_comp(k)                                                 &
+                                       * ( ibit21 + ibit22 + ibit23 )         &
+                - ( v(k,j,i)     + v(k,j-1,i) )                               &
+                                       * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
+                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
+                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
+                                         )                                    &   
+                  ) * ddy                                                     &
+               +  ( w_comp                                                    &
+                                       * ( ibit24 + ibit25 + ibit26 )         &
+                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
+                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
+                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
+                                         + IBITS(wall_flags_0(k-1,j,i),26,1)  &
+                                         )                                    &
+                   ) * ddzw(k)   &
                 ) * 0.5_wp
+
 
           tend(k,j,i) = tend(k,j,i) - (                                       &
@@ -2162,7 +2219,7 @@
 
           DO  k = nzb+1, nzb_max
-             ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
-             ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
-             ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
+             ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1)
+             ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
+             ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
 
              v_comp         = v(k+1,j,i) + v(k,j,i) - gv
@@ -2220,7 +2277,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
-             ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
-             ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
+             ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
+             ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
+             ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
 
              u_comp           = u(k+1,j,i) + u(k,j,i) - gu
@@ -2402,8 +2459,27 @@
 !--       correction is needed to overcome numerical instabilities introduced
 !--       by a not sufficient reduction of divergences near topography. 
-          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
-              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
-              +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
+          div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )            &
+                  - ( u(k+1,j,i) + u(k,j,i)   )                               & 
+                                    * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),28,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
+                                      )                                       &
+                  ) * ddx                                                     &
+              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
+                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
+                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),31,1)     &
+                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
+                                      )                                       &
+                  ) * ddy                                                     &
+              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
+                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
+                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
+                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
+                                      + IBITS(wall_flags_00(k-1,j,i),3,1)     &
+                                      )                                       & 
+                  ) * ddzu(k+1)   &
                 ) * 0.5_wp
+
 
           tend(k,j,i) = tend(k,j,i) - (                                       &
@@ -2631,7 +2707,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
-             ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
-             ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
+             ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
+             ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
+             ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
 
              u_comp                 = u(k,j,i) - u_gtrans
@@ -2692,7 +2768,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
-             ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
-             ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
+             ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
+             ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
+             ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
 
              v_comp               = v(k,j,i) - v_gtrans
@@ -3006,7 +3082,23 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
-                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
-                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
+                div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
+                          - u(k,j,i)   * ( IBITS(wall_flags_0(k,j,i-1),0,1)    &
+                                         + IBITS(wall_flags_0(k,j,i-1),1,1)    &
+                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
+                                         )                                     &
+                          ) * ddx                                              &
+                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
+                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
+                                         + IBITS(wall_flags_0(k,j-1,i),4,1)    &
+                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
+                                         )                                     &
+                          ) * ddy                                              &
+                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
+                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
+                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
+                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
+                                         )                                     &      
+                          ) * ddzw(k)
+
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
@@ -3388,7 +3480,7 @@
              DO  k = nzb+1, nzt
 
-                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
-                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
-                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
+                ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
+                ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
+                ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
 
                 u_comp              = u(k,j,i) - u_gtrans
@@ -3423,4 +3515,8 @@
                                                         )
 
+                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
+                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
+                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
+
                 u_comp    = u(k,j,i+1) - u_gtrans
                 flux_r    = u_comp * (                                        &
@@ -3454,7 +3550,7 @@
                                              )
 
-                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
-                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
-                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
+                ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
+                ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
+                ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
 
                 v_comp    = v(k,j,i) - v_gtrans
@@ -3489,4 +3585,7 @@
                                              )
 
+                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
+                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
+                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
 
                 v_comp    = v(k,j+1,i) - v_gtrans
@@ -3732,7 +3831,23 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
-                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
-                              + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
+                div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
+                          - u(k,j,i)   * ( IBITS(wall_flags_0(k,j,i-1),0,1)    &
+                                         + IBITS(wall_flags_0(k,j,i-1),1,1)    &
+                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
+                                         )                                     &
+                          ) * ddx                                              &
+                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
+                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
+                                         + IBITS(wall_flags_0(k,j-1,i),4,1)    &
+                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
+                                         )                                     &
+                          ) * ddy                                              &
+                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
+                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
+                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
+                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
+                                         )                                     &      
+                          ) * ddzw(k)
+
 
                 tend(k,j,i) = - (                                             &
@@ -3852,7 +3967,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
-             ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
-             ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
+             ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
+             ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
+             ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
 
              u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
@@ -3910,7 +4025,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
-             ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
-             ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
+             ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
+             ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
+             ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
 
              v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
@@ -4085,9 +4200,28 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
-                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
-                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
-                                                                    * ddzw(k)  &
-                      ) * 0.5_wp
+                div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 )             &
+                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
+                                    * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
+                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
+                                      )                                       &   
+                  ) * ddx                                                     &
+               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
+                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
+                                    * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),13,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
+                                      )                                       &
+                  ) * ddy                                                     &
+               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
+                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
+                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
+                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
+                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
+                                      )                                       &  
+                  ) * ddzw(k)   &
+                ) * 0.5_wp
+
+
 
                 tend(k,j,i) = tend(k,j,i) - (                                  &
@@ -4316,7 +4450,7 @@
              DO  k = nzb+1, nzt
 
-                ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
-                ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
-                ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
+                ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
+                ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
+                ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
 
                 u_comp_l           = u(k,j,i) + u(k,j,i-1) - gu
@@ -4351,4 +4485,8 @@
                                                          )
 
+                ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
+                ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
+                ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
+
                 u_comp    = u(k,j,i+1) + u(k,j,i)
                 flux_r    = ( u_comp   - gu ) * (                           &
@@ -4382,7 +4520,7 @@
                                                      )
 
-                ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
-                ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
-                ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
+                ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
+                ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
+                ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
 
                 v_comp_s                 = v(k,j,i) + v(k,j,i-1) - gv
@@ -4418,4 +4556,8 @@
 
 
+                ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
+                ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
+                ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
+
                 v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
                 flux_n    = v_comp * (                                       &
@@ -4532,9 +4674,27 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
-                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
-                     +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
-                                                                    * ddzw(k)  &
-                      ) * 0.5_wp
+                div = ( ( u_comp    * ( ibit9 + ibit10 + ibit11 )             &
+                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
+                                    * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
+                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
+                                      )                                       &   
+                  ) * ddx                                                     &
+               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
+                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
+                                    * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),13,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
+                                      )                                       &
+                  ) * ddy                                                     &
+               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
+                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
+                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
+                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
+                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
+                                      )                                       &  
+                  ) * ddzw(k)   &
+                ) * 0.5_wp
+
 
                 tend(k,j,i) = - (                                              &
@@ -4646,7 +4806,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
-             ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
-             ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
+             ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
+             ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
+             ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
 
              u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
@@ -4704,7 +4864,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
-             ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
-             ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
+             ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
+             ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
+             ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
 
              v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
@@ -4878,9 +5038,30 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
-                     +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
-                     +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) )         &
-                                                                  ) * ddzw(k) &
-                      ) * 0.5_wp
+                div = ( ( ( u_comp     + gu )                                 &
+                                       * ( ibit18 + ibit19 + ibit20 )         &
+                - ( u(k,j-1,i)   + u(k,j,i) )                                 &
+                                       * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
+                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
+                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
+                                         )                                    &   
+                  ) * ddx                                                     &
+               +  ( v_comp(k)                                                 &
+                                       * ( ibit21 + ibit22 + ibit23 )         &
+                - ( v(k,j,i)     + v(k,j-1,i) )                               &
+                                       * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
+                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
+                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
+                                         )                                    &   
+                  ) * ddy                                                     &
+               +  ( w_comp                                                    &
+                                       * ( ibit24 + ibit25 + ibit26 )         &
+                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
+                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
+                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
+                                         + IBITS(wall_flags_0(k-1,j,i),26,1)  &
+                                         )                                    &
+                   ) * ddzw(k)   &
+                ) * 0.5_wp
+
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
@@ -5120,7 +5301,7 @@
              DO  k = nzb+1, nzt
 
-                ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
-                ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
-                ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
+                ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
+                ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
+                ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
 
                 u_comp_l                 = u(k,j-1,i) + u(k,j,i) - gu
@@ -5155,4 +5336,8 @@
                                                            )
 
+                ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
+                ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
+                ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
+
                 u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
                 flux_r    = u_comp * (                                       &
@@ -5186,7 +5371,7 @@
                                               )
 
-                ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
-                ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
-                ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
+                ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
+                ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
+                ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
 
 
@@ -5222,4 +5407,8 @@
                                                           )
 
+                ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
+                ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
+                ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
+
                 v_comp = v(k,j+1,i) + v(k,j,i)
                 flux_n = ( v_comp - gv ) * (                                 &
@@ -5336,9 +5525,30 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
-                     +  ( v_comp      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
-                     +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) )         &
-                                                                  ) * ddzw(k) &
-                      ) * 0.5_wp
+                div = ( ( ( u_comp     + gu )                                 &
+                                       * ( ibit18 + ibit19 + ibit20 )         &
+                - ( u(k,j-1,i)   + u(k,j,i) )                                 &
+                                       * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
+                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
+                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
+                                         )                                    &   
+                  ) * ddx                                                     &
+               +  ( v_comp                                                    &
+                                       * ( ibit21 + ibit22 + ibit23 )         &
+                - ( v(k,j,i)     + v(k,j-1,i) )                               &
+                                       * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
+                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
+                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
+                                         )                                    &   
+                  ) * ddy                                                     &
+               +  ( w_comp                                                    &
+                                       * ( ibit24 + ibit25 + ibit26 )         &
+                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
+                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
+                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
+                                         + IBITS(wall_flags_0(k-1,j,i),26,1)  &
+                                         )                                    &
+                   ) * ddzw(k)   &
+                ) * 0.5_wp
+
 
                 tend(k,j,i) = - (                                              &
@@ -5452,7 +5662,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
-             ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
-             ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
+             ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
+             ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
+             ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
 
              u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
@@ -5510,7 +5720,7 @@
           DO  k = nzb+1, nzb_max
 
-             ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
-             ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
-             ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
+             ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1)
+             ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
+             ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
 
              v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
@@ -5690,9 +5900,28 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
-                    +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
-                    +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
-                                                                 * ddzu(k+1) &
-                      ) * 0.5_wp
+                div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
+                  - ( u(k+1,j,i) + u(k,j,i)   )                               & 
+                                    * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),28,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
+                                      )                                       &
+                  ) * ddx                                                     &
+              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
+                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
+                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),31,1)     &
+                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
+                                      )                                       &
+                  ) * ddy                                                     &
+              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
+                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
+                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
+                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
+                                      + IBITS(wall_flags_00(k-1,j,i),3,1)     &
+                                      )                                       & 
+                  ) * ddzu(k+1)   &
+                ) * 0.5_wp
+
+
 
                 tend(k,j,i) = tend(k,j,i) - (                                 &
@@ -5903,7 +6132,7 @@
              DO  k = nzb+1, nzt
 
-                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
-                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
-                ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
+                ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
+                ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
+                ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
 
                 u_comp_l                 = u(k+1,j,i) + u(k,j,i) - gu
@@ -5938,4 +6167,8 @@
                                                             )
 
+                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
+                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
+                ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
+
                 u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
                 flux_r    = u_comp * (                                       &
@@ -5968,7 +6201,7 @@
                                  ( w(k,j,i+3) - w(k,j,i-2) )                 &
                                               )
-                ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
-                ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
-                ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
+                ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1)
+                ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
+                ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
 
                 v_comp_s               = v(k+1,j,i) + v(k,j,i) - gv
@@ -6003,4 +6236,8 @@
                                                         )
 
+                ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
+                ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
+                ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
+
                 v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
                 flux_n    = v_comp * (                                       &
@@ -6118,9 +6355,27 @@
 !--             correction is needed to overcome numerical instabilities caused
 !--             by a not sufficient reduction of divergences near topography.
-                div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
-                    +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
-                    +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
-                                                                 * ddzu(k+1) &
-                      ) * 0.5_wp 
+                div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
+                  - ( u(k+1,j,i) + u(k,j,i)   )                               & 
+                                    * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),28,1)     &
+                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
+                                      )                                       &
+                  ) * ddx                                                     &
+              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            & 
+                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
+                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
+                                      + IBITS(wall_flags_0(k,j-1,i),31,1)     &
+                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
+                                      )                                       &
+                  ) * ddy                                                     &
+              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
+                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
+                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
+                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
+                                      + IBITS(wall_flags_00(k-1,j,i),3,1)     &
+                                      )                                       & 
+                  ) * ddzu(k+1)   &
+                ) * 0.5_wp
+
 
                 tend(k,j,i) = - (                                                &
